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 13 sqrt( 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 } 71