ReactOS 0.4.15-dev-7887-g64a59a1
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
11double
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 __cdecl
Definition: accygwin.h:79
double sqrt(double x)
Definition: sqrt.c:5
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
GLint GLint GLint GLint GLint GLint y
Definition: gl.h:1548
GLenum GLint GLenum GLsizei GLsizei GLsizei GLint GLsizei const GLvoid * bits
Definition: glext.h:10929
#define bits
Definition: infblock.c:15
#define NAN
Definition: mesh.c:39
_In_ CLIPOBJ _In_ BRUSHOBJ _In_ LONG _In_ LONG _In_ LONG x2
Definition: winddi.h:3710