Back to index

glibc  2.9
e_hypotf.c
Go to the documentation of this file.
00001 /* e_hypotf.c -- float version of e_hypot.c.
00002  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
00003  */
00004 
00005 /*
00006  * ====================================================
00007  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00008  *
00009  * Developed at SunPro, a Sun Microsystems, Inc. business.
00010  * Permission to use, copy, modify, and distribute this
00011  * software is freely granted, provided that this notice
00012  * is preserved.
00013  * ====================================================
00014  */
00015 
00016 #if defined(LIBM_SCCS) && !defined(lint)
00017 static char rcsid[] = "$NetBSD: e_hypotf.c,v 1.5 1995/05/12 04:57:30 jtc Exp $";
00018 #endif
00019 
00020 #include "math.h"
00021 #include "math_private.h"
00022 
00023 #ifdef __STDC__
00024        float __ieee754_hypotf(float x, float y)
00025 #else
00026        float __ieee754_hypotf(x,y)
00027        float x, y;
00028 #endif
00029 {
00030        float a,b,t1,t2,y1,y2,w;
00031        int32_t j,k,ha,hb;
00032 
00033        GET_FLOAT_WORD(ha,x);
00034        ha &= 0x7fffffff;
00035        GET_FLOAT_WORD(hb,y);
00036        hb &= 0x7fffffff;
00037        if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
00038        SET_FLOAT_WORD(a,ha);       /* a <- |a| */
00039        SET_FLOAT_WORD(b,hb);       /* b <- |b| */
00040        if((ha-hb)>0xf000000) {return a+b;} /* x/y > 2**30 */
00041        k=0;
00042        if(ha > 0x58800000) {       /* a>2**50 */
00043           if(ha >= 0x7f800000) {   /* Inf or NaN */
00044               w = a+b;                    /* for sNaN */
00045               if(ha == 0x7f800000) w = a;
00046               if(hb == 0x7f800000) w = b;
00047               return w;
00048           }
00049           /* scale a and b by 2**-60 */
00050           ha -= 0x1e000000; hb -= 0x1e000000;    k += 60;
00051           SET_FLOAT_WORD(a,ha);
00052           SET_FLOAT_WORD(b,hb);
00053        }
00054        if(hb < 0x26800000) {       /* b < 2**-50 */
00055            if(hb <= 0x007fffff) {  /* subnormal b or 0 */
00056                if(hb==0) return a;
00057               SET_FLOAT_WORD(t1,0x7e800000);     /* t1=2^126 */
00058               b *= t1;
00059               a *= t1;
00060               k -= 126;
00061            } else {         /* scale a and b by 2^60 */
00062                ha += 0x1e000000;   /* a *= 2^60 */
00063               hb += 0x1e000000;    /* b *= 2^60 */
00064               k -= 60;
00065               SET_FLOAT_WORD(a,ha);
00066               SET_FLOAT_WORD(b,hb);
00067            }
00068        }
00069     /* medium size a and b */
00070        w = a-b;
00071        if (w>b) {
00072            SET_FLOAT_WORD(t1,ha&0xfffff000);
00073            t2 = a-t1;
00074            w  = __ieee754_sqrtf(t1*t1-(b*(-b)-t2*(a+t1)));
00075        } else {
00076            a  = a+a;
00077            SET_FLOAT_WORD(y1,hb&0xfffff000);
00078            y2 = b - y1;
00079            SET_FLOAT_WORD(t1,ha+0x00800000);
00080            t2 = a - t1;
00081            w  = __ieee754_sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
00082        }
00083        if(k!=0) {
00084            SET_FLOAT_WORD(t1,0x3f800000+(k<<23));
00085            return t1*w;
00086        } else return w;
00087 }