* Use ln(x) = 2 * (t + t^3/3 + t^5/5 + ...) with t = (x-1)/(x+1).

cucul/bitmap.c

@@ -183,7 +183,7 @@ static float gammapow(float x, float y)
     register double logx;
     register long double v, e;
 #else
-    register float tmp, t, r;
+    register float tmp, t, t2, r;
     int i;
 #endif
 
@@ -208,20 +208,25 @@ static float gammapow(float x, float y)
                  : "=t" (v) : "0" (v), "u" (e));
     return v;
 #else
-    /* Compute ln(x) as ln(1+t) where t = x - 1, x ∈ ]0,1[
-     *   ln(1+t) = t - t^2/2 + t^3/3 - t^4/4 + t^5/5 ...
-     * The convergence is quite slow, especially when x is near 0. */
-    tmp = t = r = x - 1.0;
-    for(i = 2; i < 30; i++)
+    /* Compute ln(x) for x ∈ ]0,1]
+     *   ln(x) = 2 * (t + t^3/3 + t^5/5 + ...) with t = (x-1)/(x+1)
+     * The convergence is a bit slow, especially when x is near 0. */
+    t = (x - 1.0) / (x + 1.0);
+    t2 = t * t;
+    tmp = r = t;
+    for(i = 3; i < 20; i += 2)
     {
-        r *= -t;
+        r *= t2;
         tmp += r / i;
     }
 
-    /* Compute x^-y as e^t where t = y*ln(x):
+    /* Compute -y*ln(x) */
+    tmp = - y * 2.0 * tmp;
+
+    /* Compute x^-y as e^t where t = -y*ln(x):
      *   e^t = 1 + t/1! + t^2/2! + t^3/3! + t^4/4! + t^5/5! ...
-     * The convergence is a lot faster here, thanks to the factorial. */
-    r = t = - y * tmp;
+     * The convergence is quite faster here, thanks to the factorial. */
+    r = t = tmp;
     tmp = 1.0 + t;
     for(i = 2; i < 16; i++)
     {