ReactOS  0.4.15-dev-1150-g593bcce
sqrt.c
Go to the documentation of this file.
1 /*
2  * COPYRIGHT: BSD - See COPYING.ARM in the top level directory
3  * PROJECT: ReactOS CRT library
4  * PURPOSE: Portable implementation of sqrt
5  * PROGRAMMER: Timo Kreuzer (timo.kreuzer@reactos.org)
6  */
7 
8 #include <math.h>
9 #include <assert.h>
10 
11 double
12 __cdecl
14  double x)
15 {
16  const double threehalfs = 1.5;
17  const double x2 = x * 0.5;
18  long long bits;
19  double inv, y;
20 
21  /* Handle special cases */
22  if (x == 0.0)
23  {
24  return x;
25  }
26  else if (x < 0.0)
27  {
28  return -NAN;
29  }
30 
31  /* Convert into a 64 bit integer */
32  bits = *(long long *)&x;
33 
34  /* Check for !finite(x) */
35  if ((bits & 0x7ff7ffffffffffffLL) == 0x7ff0000000000000LL)
36  {
37  return x;
38  }
39 
40  /* Step 1: quick approximation of 1/sqrt(x) with bit magic
41  See http://en.wikipedia.org/wiki/Fast_inverse_square_root */
42  bits = 0x5fe6eb50c7b537a9ll - (bits >> 1);
43  inv = *(double*)&bits;
44 
45  /* Step 2: 3 Newton iterations to approximate 1 / sqrt(x) */
46  inv = inv * (threehalfs - (x2 * inv * inv));
47  inv = inv * (threehalfs - (x2 * inv * inv));
48  inv = inv * (threehalfs - (x2 * inv * inv));
49 
50  /* Step 3: 1 additional Heron iteration has shown to maximize the precision.
51  Normally the formula would be: y = (y + (x / y)) * 0.5;
52  Instead we use the inverse sqrt directly */
53  y = ((1 / inv) + (x * inv)) * 0.5;
54 
55  //assert(y == (double)((y + (x / y)) * 0.5));
56  /* GCC BUG: While the C-Standard requires that an explicit cast to
57  double converts the result of a computation to the appropriate
58  64 bit value, our GCC ignores this and uses an 80 bit FPU register
59  in an intermediate value, so we need to make sure it is stored in
60  a memory location before comparison */
61 //#if DBG
62 // {
63 // volatile double y1 = y, y2;
64 // y2 = (y + (x / y)) * 0.5;
65 // assert(y1 == y2);
66 // }
67 //#endif
68 
69  return y;
70 }
#define LL
Definition: tui.h:84
#define __cdecl
Definition: accygwin.h:79
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
GLenum GLint GLenum GLsizei GLsizei GLsizei GLint GLsizei const GLvoid * bits
Definition: glext.h:10929
w ll
Definition: byte_order.h:166
GLint GLint GLint GLint GLint GLint y
Definition: gl.h:1548
_In_ CLIPOBJ _In_ BRUSHOBJ _In_ LONG _In_ LONG _In_ LONG x2
Definition: winddi.h:3706
#define NAN
Definition: effect.c:41
double sqrt(double x)
Definition: sqrt.c:5