Update fallback rsqrt implementation to use optimal constants.

Improves max relative error from 0.00175126 to 0.000650197.
Also add unit tests to check error bounds.

BUG=chromium:511458

Review URL: https://codereview.chromium.org/1251423002

2 files changed, 39 insertions, 2 deletions

diff --git a/include/core/SkFloatingPoint.h b/include/core/SkFloatingPoint.h
index 7c34706f7a..13e8963c14 100644
--- a/include/core/SkFloatingPoint.h
+++ b/include/core/SkFloatingPoint.h
@@ -151,6 +151,8 @@ static inline float sk_float_rsqrt(const float x) {
 //
 // We do one step of Newton's method to refine the estimates in the NEON and null paths.  No
 // refinement is faster, but very innacurate.  Two steps is more accurate, but slower than 1/sqrt.
+//
+// Optimized constants in the null path courtesy of http://rrrola.wz.cz/inv_sqrt.html
 #if SK_CPU_SSE_LEVEL >= SK_CPU_SSE_LEVEL_SSE1
     return _mm_cvtss_f32(_mm_rsqrt_ss(_mm_set_ss(x)));
 #elif defined(SK_ARM_HAS_NEON)
@@ -165,12 +167,12 @@ static inline float sk_float_rsqrt(const float x) {
 #else
     // Get initial estimate.
     int i = *SkTCast<int*>(&x);
-    i = 0x5f3759df - (i>>1);
+    i = 0x5F1FFFF9 - (i>>1);
     float estimate = *SkTCast<float*>(&i);
 
     // One step of Newton's method to refine.
     const float estimate_sq = estimate*estimate;
-    estimate *= (1.5f-0.5f*x*estimate_sq);
+    estimate *= 0.703952253f*(2.38924456f-x*estimate_sq);
     return estimate;
 #endif
 }
diff --git a/tests/MathTest.cpp b/tests/MathTest.cpp
index 772bade185..3c2bc64407 100644
--- a/tests/MathTest.cpp
+++ b/tests/MathTest.cpp
@@ -382,6 +382,40 @@ static void unittest_half(skiatest::Reporter* reporter) {
 
 }
 
+static void test_rsqrt(skiatest::Reporter* reporter) {
+    const float maxRelativeError = 6.50196699e-4f;
+
+    // test close to 0 up to 1
+    float input = 0.000001f;
+    for (int i = 0; i < 1000; ++i) {
+        float exact = 1.0f/sk_float_sqrt(input);
+        float estimate = sk_float_rsqrt(input);
+        float relativeError = sk_float_abs(exact - estimate)/exact;
+        REPORTER_ASSERT(reporter, relativeError <= maxRelativeError);
+        input += 0.001f;
+    }
+
+    // test 1 to ~100
+    input = 1.0f;
+    for (int i = 0; i < 1000; ++i) {
+        float exact = 1.0f/sk_float_sqrt(input);
+        float estimate = sk_float_rsqrt(input);
+        float relativeError = sk_float_abs(exact - estimate)/exact;
+        REPORTER_ASSERT(reporter, relativeError <= maxRelativeError);
+        input += 0.01f;
+    }
+
+    // test some big numbers
+    input = 1000000.0f;
+    for (int i = 0; i < 100; ++i) {
+        float exact = 1.0f/sk_float_sqrt(input);
+        float estimate = sk_float_rsqrt(input);
+        float relativeError = sk_float_abs(exact - estimate)/exact;
+        REPORTER_ASSERT(reporter, relativeError <= maxRelativeError);
+        input += 754326.f;
+    }
+}
+
 static void test_muldiv255(skiatest::Reporter* reporter) {
     for (int a = 0; a <= 255; a++) {
         for (int b = 0; b <= 255; b++) {
@@ -521,6 +555,7 @@ DEF_TEST(Math, reporter) {
     unittest_fastfloat(reporter);
     unittest_isfinite(reporter);
     unittest_half(reporter);
+    test_rsqrt(reporter);
 
     for (i = 0; i < 10000; i++) {
         SkFixed numer = rand.nextS();