SDL  2.0
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.

105 {
106  double z;
107  int32_t sign = (int) 0x80000000;
108  int32_t ix0, s0, q, m, t, i;
109  u_int32_t r, t1, s1, ix1, q1;
110 
111  EXTRACT_WORDS(ix0, ix1, x);
112 
113  /* take care of Inf and NaN */
114  if ((ix0 & 0x7ff00000) == 0x7ff00000) {
115  return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
116  sqrt(-inf)=sNaN */
117  }
118  /* take care of zero */
119  if (ix0 <= 0) {
120  if (((ix0 & (~sign)) | ix1) == 0)
121  return x; /* sqrt(+-0) = +-0 */
122  else if (ix0 < 0)
123  return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
124  }
125  /* normalize x */
126  m = (ix0 >> 20);
127  if (m == 0) { /* subnormal x */
128  while (ix0 == 0) {
129  m -= 21;
130  ix0 |= (ix1 >> 11);
131  ix1 <<= 21;
132  }
133  for (i = 0; (ix0 & 0x00100000) == 0; i++)
134  ix0 <<= 1;
135  m -= i - 1;
136  ix0 |= (ix1 >> (32 - i));
137  ix1 <<= i;
138  }
139  m -= 1023; /* unbias exponent */
140  ix0 = (ix0 & 0x000fffff) | 0x00100000;
141  if (m & 1) { /* odd m, double x to make it even */
142  ix0 += ix0 + ((ix1 & sign) >> 31);
143  ix1 += ix1;
144  }
145  m >>= 1; /* m = [m/2] */
146 
147  /* generate sqrt(x) bit by bit */
148  ix0 += ix0 + ((ix1 & sign) >> 31);
149  ix1 += ix1;
150  q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
151  r = 0x00200000; /* r = moving bit from right to left */
152 
153  while (r != 0) {
154  t = s0 + r;
155  if (t <= ix0) {
156  s0 = t + r;
157  ix0 -= t;
158  q += r;
159  }
160  ix0 += ix0 + ((ix1 & sign) >> 31);
161  ix1 += ix1;
162  r >>= 1;
163  }
164 
165  r = sign;
166  while (r != 0) {
167  t1 = s1 + r;
168  t = s0;
169  if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
170  s1 = t1 + r;
171  if (((t1 & sign) == sign) && (s1 & sign) == 0)
172  s0 += 1;
173  ix0 -= t;
174  if (ix1 < t1)
175  ix0 -= 1;
176  ix1 -= t1;
177  q1 += r;
178  }
179  ix0 += ix0 + ((ix1 & sign) >> 31);
180  ix1 += ix1;
181  r >>= 1;
182  }
183 
184  /* use floating add to find out rounding direction */
185  if ((ix0 | ix1) != 0) {
186  z = one - tiny; /* trigger inexact flag */
187  if (z >= one) {
188  z = one + tiny;
189  if (q1 == (u_int32_t) 0xffffffff) {
190  q1 = 0;
191  q += 1;
192  } else if (z > one) {
193  if (q1 == (u_int32_t) 0xfffffffe)
194  q += 1;
195  q1 += 2;
196  } else
197  q1 += (q1 & 1);
198  }
199  }
200  ix0 = (q >> 1) + 0x3fe00000;
201  ix1 = q1 >> 1;
202  if ((q & 1) == 1)
203  ix1 |= sign;
204  ix0 += (m << 20);
205  INSERT_WORDS(z, ix0, ix1);
206  return z;
207 }
GLdouble GLdouble GLdouble r
Definition: SDL_opengl.h:2079
GLdouble GLdouble z
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
signed int int32_t
GLdouble GLdouble GLdouble GLdouble q
Definition: SDL_opengl.h:2087
const GLfloat * m
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
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
unsigned int u_int32_t
Definition: math_private.h:29
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: math_private.h:93
#define INSERT_WORDS(d, ix0, ix1)
Definition: math_private.h:121
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s1
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s0
static double tiny
Definition: e_sqrt.c:94
GLdouble GLdouble t
Definition: SDL_opengl.h:2071
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().