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 double | one = 1.0
static double | tiny = 1.0e-300

## Function Documentation

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

Definition at line 102 of file e_sqrt.c.

References EXTRACT_WORDS, i, INSERT_WORDS, one, and tiny.

{
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:2072
int32_t
signed int int32_t
Definition: SDL_config_windows.h:61
x
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1567
q
GLdouble GLdouble GLdouble GLdouble q
Definition: SDL_opengl.h:2080
m
const GLfloat * m
Definition: SDL_opengl_glext.h:6092
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:29
EXTRACT_WORDS
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: math_private.h:93
INSERT_WORDS
#define INSERT_WORDS(d, ix0, ix1)
Definition: math_private.h:121
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:42
z
GLdouble GLdouble z
Definition: SDL_opengl_glext.h:404
tiny
static double tiny
Definition: e_sqrt.c:94
t
GLdouble GLdouble t
Definition: SDL_opengl.h:2064
one
static double one
Definition: e_sqrt.c:94

## Variable Documentation

double one = 1.0

static

Definition at line 94 of file e_sqrt.c.

Referenced by __ieee754_sqrt().

double tiny = 1.0e-300

static

Definition at line 94 of file e_sqrt.c.

Referenced by __ieee754_sqrt().