Back to index

lightning-sunbird  0.9+nobinonly
Functions | Variables
e_sqrt.c File Reference
#include "fdlibm.h"

Go to the source code of this file.

Functions

double __ieee754_sqrt (double x)

Variables

static double one = 1.0
static double tiny = 1.0e-300

Function Documentation

double __ieee754_sqrt ( double  x)

Definition at line 137 of file e_sqrt.c.

{
        fd_twoints u;
       double z;
       int    sign = (int)0x80000000; 
       unsigned r,t1,s1,ix1,q1;
       int ix0,s0,q,m,t,i;

        u.d = x;
       ix0 = __HI(u);                     /* high word of x */
       ix1 = __LO(u);              /* low word of x */

    /* take care of Inf and NaN */
       if((ix0&0x7ff00000)==0x7ff00000) {               
           return x*x+x;           /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
                                      sqrt(-inf)=sNaN */
       } 
    /* take care of zero */
       if(ix0<=0) {
           if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
           else if(ix0<0)
              return (x-x)/(x-x);         /* sqrt(-ve) = sNaN */
       }
    /* normalize x */
       m = (ix0>>20);
       if(m==0) {                         /* subnormal x */
           while(ix0==0) {
              m -= 21;
              ix0 |= (ix1>>11); ix1 <<= 21;
           }
           for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
           m -= i-1;
           ix0 |= (ix1>>(32-i));
           ix1 <<= i;
       }
       m -= 1023;    /* unbias exponent */
       ix0 = (ix0&0x000fffff)|0x00100000;
       if(m&1){      /* odd m, double x to make it even */
           ix0 += ix0 + ((ix1&sign)>>31);
           ix1 += ix1;
       }
       m >>= 1;      /* m = [m/2] */

    /* generate sqrt(x) bit by bit */
       ix0 += ix0 + ((ix1&sign)>>31);
       ix1 += ix1;
       q = q1 = s0 = s1 = 0;       /* [q,q1] = sqrt(x) */
       r = 0x00200000;             /* r = moving bit from right to left */

       while(r!=0) {
           t = s0+r; 
           if(t<=ix0) { 
              s0   = t+r; 
              ix0 -= t; 
              q   += r; 
           } 
           ix0 += ix0 + ((ix1&sign)>>31);
           ix1 += ix1;
           r>>=1;
       }

       r = sign;
       while(r!=0) {
           t1 = s1+r; 
           t  = s0;
           if((t<ix0)||((t==ix0)&&(t1<=ix1))) { 
              s1  = t1+r;
              if(((int)(t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
              ix0 -= t;
              if (ix1 < t1) ix0 -= 1;
              ix1 -= t1;
              q1  += r;
           }
           ix0 += ix0 + ((ix1&sign)>>31);
           ix1 += ix1;
           r>>=1;
       }

    /* use floating add to find out rounding direction */
       if((ix0|ix1)!=0) {
           z = one-tiny; /* trigger inexact flag */
           if (z>=one) {
               z = one+tiny;
               if (q1==(unsigned)0xffffffff) { q1=0; q += 1;}
              else if (z>one) {
                  if (q1==(unsigned)0xfffffffe) q+=1;
                  q1+=2; 
              } else
                   q1 += (q1&1);
           }
       }
       ix0 = (q>>1)+0x3fe00000;
       ix1 =  q1>>1;
       if ((q&1)==1) ix1 |= sign;
       ix0 += (m <<20);
        u.d = z;
       __HI(u) = ix0;
       __LO(u) = ix1;
        z = u.d;
       return z;
}

Here is the caller graph for this function:


Variable Documentation

double one = 1.0 [static]

Definition at line 131 of file e_sqrt.c.

double tiny = 1.0e-300 [static]

Definition at line 131 of file e_sqrt.c.