From 6d85192ee672a87bfda61fc369bdb5e33965a327 Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Wed, 28 Sep 2011 17:05:00 +0000
Subject: [PATCH] test: minor updates to the Pi program (now almost deprecated)
 and the Remez exchange program.

---
 test/math/pi.cpp    | 36 +++++++++++++++++-------------------
 test/math/remez.cpp | 36 ++++++++++++++++++++++++++++--------
 2 files changed, 45 insertions(+), 27 deletions(-)

diff --git a/test/math/pi.cpp b/test/math/pi.cpp
index 3c3ecf7e..953caf12 100644
--- a/test/math/pi.cpp
+++ b/test/math/pi.cpp
@@ -18,30 +18,28 @@ using namespace lol;
 
 int main(int argc, char **argv)
 {
-    /* Approximate Pi using Machin's formula: 16*atan(1/5)-4*atan(1/239) */
-    real sum = 0.0, x0 = 5.0, x1 = 239.0;
-    real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16.0, r4 = 4.0;
+    printf("Pi: "); real::R_PI.print();
+    printf("e: "); real::R_E.print();
+    printf("ln(2): "); real::R_LN2.print();
+    printf("sqrt(2): "); real::R_SQRT2.print();
+    printf("sqrt(1/2): "); real::R_SQRT1_2.print();
 
-    /* Degree 240 is required for 512-bit mantissa precision */
-    for (int i = 1; i < 240; i += 2)
-    {
-        sum += r16 / (x0 * (real)i) - r4 / (x1 * (real)i);
-        x0 *= m0;
-        x1 *= m1;
-    }
-
-    sum.print();
-
-    /* Bonus: compute e for fun. */
-    sum = 0.0;
-    x0 = 1.0;
+    /* Gauss-Legendre computation of Pi -- doesn't work well at all,
+     * probably because we use finite precision. */
+    real a = 1.0, b = sqrt((real)0.5), t = 0.25, p = 1.0;
 
-    for (int i = 1; i < 100; i++)
+    for (int i = 0; i < 3; i++)
     {
-        sum += fres(x0);
-        x0 *= (real)i;
+        real a2 = (a + b) * (real)0.5;
+        real b2 = sqrt(a * b);
+        real tmp = a - a2;
+        real t2 = t - p * tmp * tmp;
+        real p2 = p + p;
+        a = a2; b = b2; t = t2; p = p2;
     }
 
+    real sum = a + b;
+    sum = sum * sum / ((real)4 * t);
     sum.print();
 
     return EXIT_SUCCESS;
diff --git a/test/math/remez.cpp b/test/math/remez.cpp
index faaff77f..3fabf190 100644
--- a/test/math/remez.cpp
+++ b/test/math/remez.cpp
@@ -19,13 +19,33 @@
 using namespace lol;
 
 /* The order of the approximation we're looking for */
-static int const ORDER = 4;
+static int const ORDER = 8;
+
+static real compute_pi()
+{
+    /* Approximate Pi using Machin's formula: 16*atan(1/5)-4*atan(1/239) */
+    real sum = 0.0, x0 = 5.0, x1 = 239.0;
+    real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16.0, r4 = 4.0;
+
+    /* Degree 240 is required for 512-bit mantissa precision */
+    for (int i = 1; i < 240; i += 2)
+    {
+        sum += r16 / (x0 * (real)i) - r4 / (x1 * (real)i);
+        x0 *= m0;
+        x1 *= m1;
+    }
+    return sum;
+}
 
 /* The function we want to approximate */
-real myfun(real const &x)
+static real myfun(real const &x)
 {
+    static real const half_pi = compute_pi() * (real)0.5;
     static real const one = 1.0;
-    return cos(x) - one;
+    if (x == (real)0.0 || x == (real)-0.0)
+        return half_pi;
+    return sin(x * half_pi) / x;
+    //return cos(x) - one;
     //return exp(x);
 }
 
@@ -46,7 +66,7 @@ template<int N> struct Matrix
 
     Matrix<N> inv() const
     {
-        Matrix a = *this, b(real(1.0));
+        Matrix a = *this, b((real)1.0);
 
         /* Inversion method: iterate through all columns and make sure
          * all the terms are 1 on the diagonal and 0 everywhere else */
@@ -72,7 +92,7 @@ template<int N> struct Matrix
 
             /* Now we know the diagonal term is non-zero. Get its inverse
              * and use that to nullify all other terms in the column */
-            real x = fres(a.m[i][i]);
+            real x = (real)1.0 / a.m[i][i];
             for (int j = 0; j < N; j++)
             {
                 if (j == i)
@@ -162,7 +182,7 @@ static void remez_init(real *coeff, real *zeroes)
     real fxn[ORDER + 1];
     for (int i = 0; i < ORDER + 1; i++)
     {
-        zeroes[i] = real(2 * i - ORDER) / real(ORDER + 1);
+        zeroes[i] = (real)(2 * i - ORDER) / (real)(ORDER + 1);
         fxn[i] = myfun(zeroes[i]);
     }
 
@@ -183,7 +203,7 @@ static void remez_init(real *coeff, real *zeroes)
             real sum = 0.0;
             for (int k = 0; k < ORDER + 1; k++)
                 if (cheby[n][k])
-                    sum += real(cheby[n][k]) * powers[k];
+                    sum += (real)cheby[n][k] * powers[k];
             mat.m[i][n] = sum;
         }
     }
@@ -307,7 +327,7 @@ static void remez_step(real *coeff, real *control)
             real sum = 0.0;
             for (int k = 0; k < ORDER + 1; k++)
                 if (cheby[n][k])
-                    sum += real(cheby[n][k]) * powers[k];
+                    sum += (real)cheby[n][k] * powers[k];
             mat.m[i][n] = sum;
         }
         mat.m[i][ORDER + 1] = (real)(-1 + (i & 1) * 2);