[tkreuzer] 67956: [CRT] Implement a portable version of sqrt. It's using an optimized combination of different algorithms, resulting in code that calculates the result in the maximum possible (when...

Author: tkreuzer
Date: Fri May 29 20:10:48 2015
New Revision: 67956
URL: http://svn.reactos.org/svn/reactos?rev=67956&view=rev
Log:
[CRT]
Implement a portable version of sqrt. It's using an optimized combination of different
algorithms, resulting in code that calculates the result in the maximum possible (when
internally using double) precision with only 11 multiplications and a single division.
Added:
    trunk/reactos/lib/sdk/crt/math/sqrt.c   (with props)
Removed:
    trunk/reactos/lib/sdk/crt/math/arm/sqrt.s
Modified:
    trunk/reactos/lib/sdk/crt/crt.cmake
    trunk/reactos/lib/sdk/crt/libcntpr.cmake
Modified: trunk/reactos/lib/sdk/crt/crt.cmake
URL:
http://svn.reactos.org/svn/reactos/trunk/reactos/lib/sdk/crt/crt.cmake?rev=…
==============================================================================
--- trunk/reactos/lib/sdk/crt/crt.cmake [iso-8859-1] (original)
+++ trunk/reactos/lib/sdk/crt/crt.cmake [iso-8859-1] Fri May 29 20:10:48 2015
@@ -477,6 +477,7 @@
     list(APPEND CRT_SOURCE
         except/arm/ehandler.c
         math/fabsf.c
+        math/sqrt.c
         math/arm/__rt_sdiv.c
         math/arm/__rt_sdiv64_worker.c
         math/arm/__rt_udiv.c
Modified: trunk/reactos/lib/sdk/crt/libcntpr.cmake
URL:
http://svn.reactos.org/svn/reactos/trunk/reactos/lib/sdk/crt/libcntpr.cmake…
==============================================================================
--- trunk/reactos/lib/sdk/crt/libcntpr.cmake    [iso-8859-1] (original)
+++ trunk/reactos/lib/sdk/crt/libcntpr.cmake    [iso-8859-1] Fri May 29 20:10:48 2015
@@ -165,7 +165,6 @@
         math/arm/log.s
         math/arm/log10.s
         math/arm/pow.s
-        math/arm/sqrt.s
         math/arm/tan.s
         math/arm/__dtoi64.s
         math/arm/__dtou64.s
Removed: trunk/reactos/lib/sdk/crt/math/arm/sqrt.s
URL:
http://svn.reactos.org/svn/reactos/trunk/reactos/lib/sdk/crt/math/arm/sqrt.…
==============================================================================
--- trunk/reactos/lib/sdk/crt/math/arm/sqrt.s   [iso-8859-1] (original)
+++ trunk/reactos/lib/sdk/crt/math/arm/sqrt.s   (removed)
@@ -1,24 +0,0 @@
-/*
- * COPYRIGHT:         BSD - See COPYING.ARM in the top level directory
- * PROJECT:           ReactOS CRT library
- * PURPOSE:           Implementation of sqrt
- * PROGRAMMER:        Timo Kreuzer (timo.kreuzer(a)reactos.org)
- */
-
-/* INCLUDES ******************************************************************/
-
-#include <kxarm.h>
-
-/* CODE **********************************************************************/
-
-    TEXTAREA
-
-    LEAF_ENTRY sqrt
-
-       __assertfail
-       bx      lr
-
-    LEAF_END sqrt
-
-    END
-/* EOF */
Added: trunk/reactos/lib/sdk/crt/math/sqrt.c
URL: http://svn.reactos.org/svn/reactos/trunk/reactos/lib/sdk/crt/math/sqrt.c?re…
==============================================================================
--- trunk/reactos/lib/sdk/crt/math/sqrt.c       (added)
+++ trunk/reactos/lib/sdk/crt/math/sqrt.c       [iso-8859-1] Fri May 29 20:10:48 2015
@@ -0,0 +1,70 @@
+/*
+ * COPYRIGHT:       BSD - See COPYING.ARM in the top level directory
+ * PROJECT:         ReactOS CRT library
+ * PURPOSE:         Portable implementation of sqrt
+ * PROGRAMMER:      Timo Kreuzer (timo.kreuzer(a)reactos.org)
+ */
+
+#include <math.h>
+#include <assert.h>
+
+double
+__cdecl
+sqrt(
+    double x)
+{
+    const double threehalfs = 1.5;
+    const double x2 = x * 0.5;
+    long long bits;
+    double inv, y;
+
+    /* Handle special cases */
+    if (x == 0.0)
+    {
+        return x;
+    }
+    else if (x < 0.0)
+    {
+        return -NAN;
+    }
+
+    /* Convert into a 64  bit integer */
+    bits = *(long long *)&x;
+
+    /* Check for !finite(x) */
+    if ((bits & 0x7ff7ffffffffffffLL) == 0x7ff0000000000000LL)
+    {
+        return x;
+    }
+
+    /* Step 1: quick approximation of 1/sqrt(x) with bit magic
+       See http://en.wikipedia.org/wiki/Fast_inverse_square_root */
+    bits = 0x5fe6eb50c7b537a9ll - (bits >> 1);
+    inv = *(double*)&bits;
+
+    /* Step 2: 3 Newton iterations to approximate 1 / sqrt(x) */
+    inv = inv * (threehalfs - (x2 * inv * inv));
+    inv = inv * (threehalfs - (x2 * inv * inv));
+    inv = inv * (threehalfs - (x2 * inv * inv));
+
+    /* Step 3: 1 additional Heron iteration has shown to maximize the precision.
+       Normally the formula would be: y = (y + (x / y)) * 0.5;
+       Instead we use the inverse sqrt directly */
+    y = ((1 / inv) + (x * inv)) * 0.5;
+
+    //assert(y == (double)((y + (x / y)) * 0.5));
+    /* GCC BUG: While the C-Standard requires that an explicit cast to
+       double converts the result of a computation to the appropriate
+       64 bit value, our GCC ignores this and uses an 80 bit FPU register
+       in an intermediate value, so we need to make sure it is stored in
+       a memory location before comparison */
+#if DBG
+    {
+        volatile double y1 = y, y2;
+        y2 = (y + (x / y)) * 0.5;
+        assert(y1 == y2);
+    }
+#endif
+
+    return y;
+}
Propchange: trunk/reactos/lib/sdk/crt/math/sqrt.c
------------------------------------------------------------------------------
    svn:eol-style = native