SDL  2.0
e_sqrt.c
Go to the documentation of this file.
1 /*
2  * ====================================================
3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4  *
5  * Developed at SunPro, a Sun Microsystems, Inc. business.
6  * Permission to use, copy, modify, and distribute this
7  * software is freely granted, provided that this notice
8  * is preserved.
9  * ====================================================
10  */
11 
12 /* __ieee754_sqrt(x)
13  * Return correctly rounded sqrt.
14  * ------------------------------------------
15  * | Use the hardware sqrt if you have one |
16  * ------------------------------------------
17  * Method:
18  * Bit by bit method using integer arithmetic. (Slow, but portable)
19  * 1. Normalization
20  * Scale x to y in [1,4) with even powers of 2:
21  * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
22  * sqrt(x) = 2^k * sqrt(y)
23  * 2. Bit by bit computation
24  * Let q = sqrt(y) truncated to i bit after binary point (q = 1),
25  * i 0
26  * i+1 2
27  * s = 2*q , and y = 2 * ( y - q ). (1)
28  * i i i i
29  *
30  * To compute q from q , one checks whether
31  * i+1 i
32  *
33  * -(i+1) 2
34  * (q + 2 ) <= y. (2)
35  * i
36  * -(i+1)
37  * If (2) is false, then q = q ; otherwise q = q + 2 .
38  * i+1 i i+1 i
39  *
40  * With some algebric manipulation, it is not difficult to see
41  * that (2) is equivalent to
42  * -(i+1)
43  * s + 2 <= y (3)
44  * i i
45  *
46  * The advantage of (3) is that s and y can be computed by
47  * i i
48  * the following recurrence formula:
49  * if (3) is false
50  *
51  * s = s , y = y ; (4)
52  * i+1 i i+1 i
53  *
54  * otherwise,
55  * -i -(i+1)
56  * s = s + 2 , y = y - s - 2 (5)
57  * i+1 i i+1 i i
58  *
59  * One may easily use induction to prove (4) and (5).
60  * Note. Since the left hand side of (3) contain only i+2 bits,
61  * it does not necessary to do a full (53-bit) comparison
62  * in (3).
63  * 3. Final rounding
64  * After generating the 53 bits result, we compute one more bit.
65  * Together with the remainder, we can decide whether the
66  * result is exact, bigger than 1/2ulp, or less than 1/2ulp
67  * (it will never equal to 1/2ulp).
68  * The rounding mode can be detected by checking whether
69  * huge + tiny is equal to huge, and whether huge - tiny is
70  * equal to huge for some floating point number "huge" and "tiny".
71  *
72  * Special cases:
73  * sqrt(+-0) = +-0 ... exact
74  * sqrt(inf) = inf
75  * sqrt(-ve) = NaN ... with invalid signal
76  * sqrt(NaN) = NaN ... with invalid signal for signaling NaN
77  *
78  * Other methods : see the appended file at the end of the program below.
79  *---------------
80  */
81 
82 #include "math_libm.h"
83 #include "math_private.h"
84 
85 static const double one = 1.0, tiny = 1.0e-300;
86 
88 {
89  double z;
90  int32_t sign = (int)0x80000000;
91  int32_t ix0,s0,q,m,t,i;
92  u_int32_t r,t1,s1,ix1,q1;
93 
94  EXTRACT_WORDS(ix0,ix1,x);
95 
96  /* take care of Inf and NaN */
97  if((ix0&0x7ff00000)==0x7ff00000) {
98  return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
99  sqrt(-inf)=sNaN */
100  }
101  /* take care of zero */
102  if(ix0<=0) {
103  if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
104  else if(ix0<0)
105  return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
106  }
107  /* normalize x */
108  m = (ix0>>20);
109  if(m==0) { /* subnormal x */
110  while(ix0==0) {
111  m -= 21;
112  ix0 |= (ix1>>11); ix1 <<= 21;
113  }
114  for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
115  m -= i-1;
116  ix0 |= (ix1>>(32-i));
117  ix1 <<= i;
118  }
119  m -= 1023; /* unbias exponent */
120  ix0 = (ix0&0x000fffff)|0x00100000;
121  if(m&1){ /* odd m, double x to make it even */
122  ix0 += ix0 + ((ix1&sign)>>31);
123  ix1 += ix1;
124  }
125  m >>= 1; /* m = [m/2] */
126 
127  /* generate sqrt(x) bit by bit */
128  ix0 += ix0 + ((ix1&sign)>>31);
129  ix1 += ix1;
130  q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
131  r = 0x00200000; /* r = moving bit from right to left */
132 
133  while(r!=0) {
134  t = s0+r;
135  if(t<=ix0) {
136  s0 = t+r;
137  ix0 -= t;
138  q += r;
139  }
140  ix0 += ix0 + ((ix1&sign)>>31);
141  ix1 += ix1;
142  r>>=1;
143  }
144 
145  r = sign;
146  while(r!=0) {
147  t1 = s1+r;
148  t = s0;
149  if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
150  s1 = t1+r;
151  if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
152  ix0 -= t;
153  if (ix1 < t1) ix0 -= 1;
154  ix1 -= t1;
155  q1 += r;
156  }
157  ix0 += ix0 + ((ix1&sign)>>31);
158  ix1 += ix1;
159  r>>=1;
160  }
161 
162  /* use floating add to find out rounding direction */
163  if((ix0|ix1)!=0) {
164  z = one-tiny; /* trigger inexact flag */
165  if (z>=one) {
166  z = one+tiny;
167  if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
168  else if (z>one) {
169  if (q1==(u_int32_t)0xfffffffe) q+=1;
170  q1+=2;
171  } else
172  q1 += (q1&1);
173  }
174  }
175  ix0 = (q>>1)+0x3fe00000;
176  ix1 = q1>>1;
177  if ((q&1)==1) ix1 |= sign;
178  ix0 += (m <<20);
179  INSERT_WORDS(z,ix0,ix1);
180  return z;
181 }
182 
183 /*
184  * wrapper sqrt(x)
185  */
186 #ifndef _IEEE_LIBM
187 double sqrt(double x)
188 {
189  double z = __ieee754_sqrt(x);
190  if (_LIB_VERSION == _IEEE_ || isnan(x))
191  return z;
192  if (x < 0.0)
193  return __kernel_standard(x, x, 26); /* sqrt(negative) */
194  return z;
195 }
196 #else
198 #endif
199 libm_hidden_def(sqrt)
200 
201 
202 /*
203 Other methods (use floating-point arithmetic)
204 -------------
205 (This is a copy of a drafted paper by Prof W. Kahan
206 and K.C. Ng, written in May, 1986)
207 
208  Two algorithms are given here to implement sqrt(x)
209  (IEEE double precision arithmetic) in software.
210  Both supply sqrt(x) correctly rounded. The first algorithm (in
211  Section A) uses newton iterations and involves four divisions.
212  The second one uses reciproot iterations to avoid division, but
213  requires more multiplications. Both algorithms need the ability
214  to chop results of arithmetic operations instead of round them,
215  and the INEXACT flag to indicate when an arithmetic operation
216  is executed exactly with no roundoff error, all part of the
217  standard (IEEE 754-1985). The ability to perform shift, add,
218  subtract and logical AND operations upon 32-bit words is needed
219  too, though not part of the standard.
220 
221 A. sqrt(x) by Newton Iteration
222 
223  (1) Initial approximation
224 
225  Let x0 and x1 be the leading and the trailing 32-bit words of
226  a floating point number x (in IEEE double format) respectively
227 
228  1 11 52 ...widths
229  ------------------------------------------------------
230  x: |s| e | f |
231  ------------------------------------------------------
232  msb lsb msb lsb ...order
233 
234 
235  ------------------------ ------------------------
236  x0: |s| e | f1 | x1: | f2 |
237  ------------------------ ------------------------
238 
239  By performing shifts and subtracts on x0 and x1 (both regarded
240  as integers), we obtain an 8-bit approximation of sqrt(x) as
241  follows.
242 
243  k := (x0>>1) + 0x1ff80000;
244  y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits
245  Here k is a 32-bit integer and T1[] is an integer array containing
246  correction terms. Now magically the floating value of y (y's
247  leading 32-bit word is y0, the value of its trailing word is 0)
248  approximates sqrt(x) to almost 8-bit.
249 
250  Value of T1:
251  static int T1[32]= {
252  0, 1024, 3062, 5746, 9193, 13348, 18162, 23592,
253  29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215,
254  83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581,
255  16499, 12183, 8588, 5674, 3403, 1742, 661, 130,};
256 
257  (2) Iterative refinement
258 
259  Apply Heron's rule three times to y, we have y approximates
260  sqrt(x) to within 1 ulp (Unit in the Last Place):
261 
262  y := (y+x/y)/2 ... almost 17 sig. bits
263  y := (y+x/y)/2 ... almost 35 sig. bits
264  y := y-(y-x/y)/2 ... within 1 ulp
265 
266 
267  Remark 1.
268  Another way to improve y to within 1 ulp is:
269 
270  y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x)
271  y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x)
272 
273  2
274  (x-y )*y
275  y := y + 2* ---------- ...within 1 ulp
276  2
277  3y + x
278 
279 
280  This formula has one division fewer than the one above; however,
281  it requires more multiplications and additions. Also x must be
282  scaled in advance to avoid spurious overflow in evaluating the
283  expression 3y*y+x. Hence it is not recommended uless division
284  is slow. If division is very slow, then one should use the
285  reciproot algorithm given in section B.
286 
287  (3) Final adjustment
288 
289  By twiddling y's last bit it is possible to force y to be
290  correctly rounded according to the prevailing rounding mode
291  as follows. Let r and i be copies of the rounding mode and
292  inexact flag before entering the square root program. Also we
293  use the expression y+-ulp for the next representable floating
294  numbers (up and down) of y. Note that y+-ulp = either fixed
295  point y+-1, or multiply y by nextafter(1,+-inf) in chopped
296  mode.
297 
298  I := FALSE; ... reset INEXACT flag I
299  R := RZ; ... set rounding mode to round-toward-zero
300  z := x/y; ... chopped quotient, possibly inexact
301  If(not I) then { ... if the quotient is exact
302  if(z=y) {
303  I := i; ... restore inexact flag
304  R := r; ... restore rounded mode
305  return sqrt(x):=y.
306  } else {
307  z := z - ulp; ... special rounding
308  }
309  }
310  i := TRUE; ... sqrt(x) is inexact
311  If (r=RN) then z=z+ulp ... rounded-to-nearest
312  If (r=RP) then { ... round-toward-+inf
313  y = y+ulp; z=z+ulp;
314  }
315  y := y+z; ... chopped sum
316  y0:=y0-0x00100000; ... y := y/2 is correctly rounded.
317  I := i; ... restore inexact flag
318  R := r; ... restore rounded mode
319  return sqrt(x):=y.
320 
321  (4) Special cases
322 
323  Square root of +inf, +-0, or NaN is itself;
324  Square root of a negative number is NaN with invalid signal.
325 
326 
327 B. sqrt(x) by Reciproot Iteration
328 
329  (1) Initial approximation
330 
331  Let x0 and x1 be the leading and the trailing 32-bit words of
332  a floating point number x (in IEEE double format) respectively
333  (see section A). By performing shifs and subtracts on x0 and y0,
334  we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
335 
336  k := 0x5fe80000 - (x0>>1);
337  y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits
338 
339  Here k is a 32-bit integer and T2[] is an integer array
340  containing correction terms. Now magically the floating
341  value of y (y's leading 32-bit word is y0, the value of
342  its trailing word y1 is set to zero) approximates 1/sqrt(x)
343  to almost 7.8-bit.
344 
345  Value of T2:
346  static int T2[64]= {
347  0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
348  0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
349  0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
350  0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
351  0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
352  0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
353  0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
354  0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};
355 
356  (2) Iterative refinement
357 
358  Apply Reciproot iteration three times to y and multiply the
359  result by x to get an approximation z that matches sqrt(x)
360  to about 1 ulp. To be exact, we will have
361  -1ulp < sqrt(x)-z<1.0625ulp.
362 
363  ... set rounding mode to Round-to-nearest
364  y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x)
365  y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
366  ... special arrangement for better accuracy
367  z := x*y ... 29 bits to sqrt(x), with z*y<1
368  z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x)
369 
370  Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
371  (a) the term z*y in the final iteration is always less than 1;
372  (b) the error in the final result is biased upward so that
373  -1 ulp < sqrt(x) - z < 1.0625 ulp
374  instead of |sqrt(x)-z|<1.03125ulp.
375 
376  (3) Final adjustment
377 
378  By twiddling y's last bit it is possible to force y to be
379  correctly rounded according to the prevailing rounding mode
380  as follows. Let r and i be copies of the rounding mode and
381  inexact flag before entering the square root program. Also we
382  use the expression y+-ulp for the next representable floating
383  numbers (up and down) of y. Note that y+-ulp = either fixed
384  point y+-1, or multiply y by nextafter(1,+-inf) in chopped
385  mode.
386 
387  R := RZ; ... set rounding mode to round-toward-zero
388  switch(r) {
389  case RN: ... round-to-nearest
390  if(x<= z*(z-ulp)...chopped) z = z - ulp; else
391  if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
392  break;
393  case RZ:case RM: ... round-to-zero or round-to--inf
394  R:=RP; ... reset rounding mod to round-to-+inf
395  if(x<z*z ... rounded up) z = z - ulp; else
396  if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
397  break;
398  case RP: ... round-to-+inf
399  if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
400  if(x>z*z ...chopped) z = z+ulp;
401  break;
402  }
403 
404  Remark 3. The above comparisons can be done in fixed point. For
405  example, to compare x and w=z*z chopped, it suffices to compare
406  x1 and w1 (the trailing parts of x and w), regarding them as
407  two's complement integers.
408 
409  ...Is z an exact square root?
410  To determine whether z is an exact square root of x, let z1 be the
411  trailing part of z, and also let x0 and x1 be the leading and
412  trailing parts of x.
413 
414  If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0
415  I := 1; ... Raise Inexact flag: z is not exact
416  else {
417  j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2
418  k := z1 >> 26; ... get z's 25-th and 26-th
419  fraction bits
420  I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
421  }
422  R:= r ... restore rounded mode
423  return sqrt(x):=z.
424 
425  If multiplication is cheaper then the foregoing red tape, the
426  Inexact flag can be evaluated by
427 
428  I := i;
429  I := (z*z!=x) or I.
430 
431  Note that z*z can overwrite I; this value must be sensed if it is
432  True.
433 
434  Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
435  zero.
436 
437  --------------------
438  z1: | f2 |
439  --------------------
440  bit 31 bit 0
441 
442  Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
443  or even of logb(x) have the following relations:
444 
445  -------------------------------------------------
446  bit 27,26 of z1 bit 1,0 of x1 logb(x)
447  -------------------------------------------------
448  00 00 odd and even
449  01 01 even
450  10 10 odd
451  10 00 even
452  11 01 even
453  -------------------------------------------------
454 
455  (4) Special cases (see (4) of Section A).
456 
457  */
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
static const double tiny
Definition: e_sqrt.c:85
const GLfloat * m
#define strong_alias(x, y)
Definition: math_private.h:28
static const double one
Definition: e_sqrt.c:85
#define attribute_hidden
Definition: math_private.h:25
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
unsigned int u_int32_t
Definition: math_private.h:31
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: math_private.h:98
#define INSERT_WORDS(d, ix0, ix1)
Definition: math_private.h:126
double attribute_hidden __ieee754_sqrt(double x)
Definition: e_sqrt.c:87
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:50
libm_hidden_def(scalbln)
Definition: s_scalbn.c:62
GLdouble GLdouble t
Definition: SDL_opengl.h:2071