Experiments on calculating reciprocal of square root

include/core/SkMath.h

Index: include/core/SkMath.h
diff --git a/include/core/SkMath.h b/include/core/SkMath.h
index 078c8fccdc2206975c07ed830222fe82f7e710af..66df00f40813514fe347a41df418392276431986 100644
--- a/include/core/SkMath.h
+++ b/include/core/SkMath.h
@@ -173,4 +173,13 @@ static inline U8CPU SkMulDiv255Round(U16CPU a, U16CPU b) {
     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