From afd2cfd69d3d7581effadaa6d60d4beb1344f318 Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Tue, 27 Sep 2011 17:09:35 +0000
Subject: [PATCH] test: more work on the Remez exchange algorithm.
---
test/math/pi.cpp | 12 ++++
test/math/remez.cpp | 156 ++++++++++++++++++++++++++++++++++++++++----
2 files changed, 156 insertions(+), 12 deletions(-)
diff --git a/test/math/pi.cpp b/test/math/pi.cpp
index e00379c1..3c3ecf7e 100644
--- a/test/math/pi.cpp
+++ b/test/math/pi.cpp
@@ -32,6 +32,18 @@ int main(int argc, char **argv)
sum.print();
+ /* Bonus: compute e for fun. */
+ sum = 0.0;
+ x0 = 1.0;
+
+ for (int i = 1; i < 100; i++)
+ {
+ sum += fres(x0);
+ x0 *= (real)i;
+ }
+
+ sum.print();
+
return EXIT_SUCCESS;
}
diff --git a/test/math/remez.cpp b/test/math/remez.cpp
index 44d8a1cc..2babb661 100644
--- a/test/math/remez.cpp
+++ b/test/math/remez.cpp
@@ -18,21 +18,116 @@
using namespace lol;
-static int const ORDER = 12;
+/* The order of the approximation we're looking for */
+static int const ORDER = 3;
-static int cheby[ORDER][ORDER];
+/* The function we want to approximate */
+double myfun(double x)
+{
+ return exp(x);
+}
+
+real myfun(real const &x)
+{
+ return exp(x);
+}
+
+/* Naive matrix inversion */
+template<int N> class Matrix
+{
+public:
+ inline Matrix() {}
+
+ Matrix(real x)
+ {
+ for (int j = 0; j < N; j++)
+ for (int i = 0; i < N; i++)
+ if (i == j)
+ m[i][j] = x;
+ else
+ m[i][j] = 0;
+ }
+
+ Matrix<N> inv() const
+ {
+ Matrix a = *this, b(real(1.0));
+
+ /* Inversion method: iterate through all columns and make sure
+ * all the terms are zero except on the diagonal where it is one */
+ for (int i = 0; i < N; i++)
+ {
+ /* If the expected coefficient is zero, add one of
+ * the other lines. The first we meet will do. */
+ if ((double)a.m[i][i] == 0.0)
+ {
+ for (int j = i + 1; j < N; j++)
+ {
+ if ((double)a.m[i][j] == 0.0)
+ continue;
+ /* Add row j to row i */
+ for (int n = 0; n < N; n++)
+ {
+ a.m[n][i] += a.m[n][j];
+ b.m[n][i] += b.m[n][j];
+ }
+ break;
+ }
+ }
+
+ /* 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]);
+ for (int j = 0; j < N; j++)
+ {
+ if (j == i)
+ continue;
+ real mul = x * a.m[i][j];
+ for (int n = 0; n < N; n++)
+ {
+ a.m[n][j] -= mul * a.m[n][i];
+ b.m[n][j] -= mul * b.m[n][i];
+ }
+ }
+ /* Finally, ensure the diagonal term is 1 */
+ for (int n = 0; n < N; n++)
+ {
+ a.m[n][i] *= x;
+ b.m[n][i] *= x;
+ }
+ }
+
+ return b;
+ }
+
+ void print() const
+ {
+ for (int j = 0; j < N; j++)
+ {
+ for (int i = 0; i < N; i++)
+ printf("%9.5f ", (double)m[j][i]);
+ printf("\n");
+ }
+ }
+
+ real m[N][N];
+};
+
+
+static int cheby[ORDER + 1][ORDER + 1];
+
+/* Fill the Chebyshev tables */
static void make_table()
{
- memset(cheby[0], 0, ORDER * sizeof(int));
+ memset(cheby, 0, sizeof(cheby));
+
cheby[0][0] = 1;
- memset(cheby[1], 0, ORDER * sizeof(int));
cheby[1][1] = 1;
- for (int i = 2; i < ORDER; i++)
+ for (int i = 2; i < ORDER + 1; i++)
{
cheby[i][0] = -cheby[i - 2][0];
- for (int j = 1; j < ORDER; j++)
+ for (int j = 1; j < ORDER + 1; j++)
cheby[i][j] = 2 * cheby[i - 1][j - 1] - cheby[i - 2][j];
}
}
@@ -41,12 +136,49 @@ int main(int argc, char **argv)
{
make_table();
-for (int i = 0; i < ORDER; i++)
-{
- for (int j = 0; j < ORDER; j++)
- printf("% 5i ", cheby[i][j]);
- printf("\n");
-}
+ /* We start with ORDER+1 points and their images through myfun() */
+ real xn[ORDER + 1];
+ real fxn[ORDER + 1];
+ for (int i = 0; i < ORDER + 1; i++)
+ {
+ //xn[i] = real(2 * i - ORDER) / real(ORDER + 1);
+ xn[i] = real(2 * i - ORDER + 1) / real(ORDER - 1);
+ fxn[i] = myfun(xn[i]);
+ }
+
+ /* We build a matrix of Chebishev evaluations: one row per point
+ * in our array, and column i is the evaluation of the ith order
+ * polynomial. */
+ Matrix<ORDER + 1> mat;
+ for (int j = 0; j < ORDER + 1; j++)
+ {
+ /* Compute the powers of x_j */
+ real powers[ORDER + 1];
+ powers[0] = 1.0;
+ for (int i = 1; i < ORDER + 1; i++)
+ powers[i] = powers[i - 1] * xn[j];
+
+ /* Compute the Chebishev evaluations at x_j */
+ for (int i = 0; i < ORDER + 1; i++)
+ {
+ real sum = 0.0;
+ for (int k = 0; k < ORDER + 1; k++)
+ if (cheby[i][k])
+ sum += real(cheby[i][k]) * powers[k];
+ mat.m[j][i] = sum;
+ }
+ }
+
+ /* Invert the matrix and build interpolation coefficients */
+ mat = mat.inv();
+ real an[ORDER + 1];
+ for (int j = 0; j < ORDER + 1; j++)
+ {
+ an[j] = 0;
+ for (int i = 0; i < ORDER + 1; i++)
+ an[j] += mat.m[j][i] * fxn[i];
+ an[j].print(10);
+ }
return EXIT_SUCCESS;
}