From 5fff7db5143ef9cba4ba983a624ea8c10c7396b9 Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Wed, 22 Mar 2006 14:56:50 +0000
Subject: [PATCH]   * Use ln(x) = 2 * (t + t^3/3 + t^5/5 + ...) with t =
 (x-1)/(x+1).

---
 cucul/bitmap.c | 25 +++++++++++++++----------
 1 file changed, 15 insertions(+), 10 deletions(-)

diff --git a/cucul/bitmap.c b/cucul/bitmap.c
index 548239e..0abdb95 100644
--- a/cucul/bitmap.c
+++ b/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++)
     {