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:2072
signed int int32_t
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1567
GLdouble GLdouble GLdouble GLdouble q
Definition: SDL_opengl.h:2080
const GLfloat * m
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
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
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
GLdouble GLdouble z
static double tiny
Definition: e_sqrt.c:94
GLdouble GLdouble t
Definition: SDL_opengl.h:2064
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().