e_sqrt.c File Reference
`#include "math_libm.h"`
`#include "math_private.h"`
Include dependency graph for e_sqrt.c:

Go to the source code of this file.

## Functions

double attribute_hidden | __ieee754_sqrt (double x)

## Variables

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

## Function Documentation

## __ieee754_sqrt()

double attribute_hidden __ieee754_sqrt | ( | double | x | ) |

Definition at line 87 of file e_sqrt.c.

References __ieee754_sqrt(), EXTRACT_WORDS, i, INSERT_WORDS, libm_hidden_def(), one, strong_alias, and tiny.

Referenced by __ieee754_sqrt().

{
double z;
int32_t sign = (int)0x80000000;
int32_t ix0,s0,q,m,t,i;
u_int32_t r,t1,s1,ix1,q1;
EXTRACT_WORDS(ix0,ix1,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(((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==(u_int32_t)0xffffffff) { q1=0; q += 1;}
else if (z>one) {
if (q1==(u_int32_t)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);
INSERT_WORDS(z,ix0,ix1);
return z;
}
r
GLdouble GLdouble GLdouble r
Definition: SDL_opengl.h:2079
z
GLdouble GLdouble z
Definition: SDL_opengl_glext.h:404
x
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
int32_t
signed int int32_t
Definition: SDL_config_windows.h:62
q
GLdouble GLdouble GLdouble GLdouble q
Definition: SDL_opengl.h:2087
tiny
static const double tiny
Definition: e_sqrt.c:85
m
const GLfloat * m
Definition: SDL_opengl_glext.h:6092
one
static const double one
Definition: e_sqrt.c:85
t1
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
Definition: SDL_opengl_glext.h:8583
u_int32_t
unsigned int u_int32_t
Definition: math_private.h:31
EXTRACT_WORDS
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: math_private.h:98
INSERT_WORDS
#define INSERT_WORDS(d, ix0, ix1)
Definition: math_private.h:126
s1
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s1
Definition: SDL_opengl_glext.h:8583
s0
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s0
Definition: SDL_opengl_glext.h:8583
i
return Display return Display Bool Bool int int int return Display XEvent Bool(*) XPointer return Display return Display Drawable _Xconst char unsigned int unsigned int return Display Pixmap Pixmap XColor XColor unsigned int unsigned int return Display _Xconst char char int char return Display Visual unsigned int int int char unsigned int unsigned int in i)
Definition: SDL_x11sym.h:50
t
GLdouble GLdouble t
Definition: SDL_opengl.h:2071

## Variable Documentation

## one

const double one = 1.0

static

Definition at line 85 of file e_sqrt.c.

Referenced by __ieee754_sqrt().

## tiny

const double tiny = 1.0e-300

static

Definition at line 85 of file e_sqrt.c.

Referenced by __ieee754_sqrt().