Experiments on calculating reciprocal of square root

include/core/SkMath.h

 
 /*
  * Copyright 2006 The Android Open Source Project
  *
  * Use of this source code is governed by a BSD-style license that can be
  * found in the LICENSE file.
  */
 
 
 #ifndef SkMath_DEFINED
@@ skipping 155 matching lines @@
  *  Return a*b/255, rounding any fractional bits.
  *  Only valid if a and b are unsigned and <= 32767.
  */
 static inline U8CPU SkMulDiv255Round(U16CPU a, U16CPU b) {
     SkASSERT(a <= 32767);
     SkASSERT(b <= 32767);
     unsigned prod = SkMulS16(a, b) + 128;
     return (prod + (prod >> 8)) >> 8;
 }
 
+static inline float SkFloatInvSqrt(float x) {

[reed1, 2013/08/12 14:05:54]

Need dox for a public API

1. Should probably be renamed to somehow reflect that this is an approximation
of the inverse.
2. Why is this public right now? What is the driving case for it?
3. Possibly document why this works at all

[Yang Gu, 2013/08/13 03:33:31]

We may add "fast" into name (e.g., SkFloatFastInvSqrt) as this algorithm is
often referred as fast
invsqrt(http://en.wikipedia.org/wiki/Fast_inverse_square_root). I want to keep
"float" in name as we may need a version of double (which uses Newton's method
twice). However, I think current name is OK (SkFloatFastInvSqrt or so is a bit
too long), and many libraries just use invsqrt or rsqrt as name. For example, in
code of Quake 3 (origin of the algorithm though it has deeper root), the snippet
is
float InvSqrt (float x){
    float xhalf = 0.5f*x;
    int i = *(int*)&x;
    i = 0x5f3759df - (i>>1);
    x = *(float*)&i;
    x = x*(1.5f - xhalf*x*x);
    return x;
}

And in hackersdelight, the code is
float rsqrt(float x0) {
   union {int ix; float x;};
// union {int ihalf; float xhalf;}; // For alternative halving step.

   x = x0;                      // x can be viewed as int.
// ihalf = ix - 0x00800000;     // Alternative to line below, for x not a
denorm.
   float xhalf = 0.5f*x;
// ix = 0x5f3759df - (ix >> 1); // Initial guess (traditional),
//                                 but slightly better:
   ix = 0x5f375a82 - (ix >> 1); // Initial guess.
   x = x*(1.5f - xhalf*x*x);    // Newton step.
// x = x*(1.5008908 - xhalf*x*x);  // Newton step for a balanced error.
   return x;
}

What's more, we may also evolve this API to be a bit more complicated. According
to my test result, the performance of intrinsic solution for IA is on par with
this fast solution, but the intrinsic solution may have better precision. 
We may have this API be used everywhere in Skia to calculate the reciprocal
square root, without caring if it will lose precision. I'm working on this right
now. Is it safe to do this everywhere, that is, is there any place that requires
high precision that this algorithm can't meet?
The driven reason is performance. Actually we already had such benchmark in
MathBench, and @tomhudson created an issue at
https://code.google.com/p/skia/issues/detail?id=885. After some investigation, I
think there are some issues with original benchmark. After the correction of
benchmark, I saw an obvious speedup on all platforms (IA desktop, IA phone and
ARM phone (Nexus4)). Please refer my previous post for detailed data.
How about adding following comment on this API? (I referred a paper
http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf)
/**
 *  Return the approximate reciprocal square root of an IEEE
 *  float (single precision).
 *  The algorithm has obvious speedup than regular floating
 *  point division, with a small sacrifice of precision.
 *  According to some paper, its maximal relative error is
 *  0.175228. 
 */

Should I add a patch set for this?

+    float xhalf = 0.5f * x;
+    int i = *SkTCast<int*>(&x);
+    i = 0x5f3759df - (i >> 1);
+    x = *SkTCast<float*>(&i);
+    x = x * (1.5f - xhalf * x * x);
+    return x;
+}
+
 #endif