Changeset 1010Tweet

Timestamp:

Oct 5, 2011, 7:01:33 PM (11 years ago)

Author:

sam

Message:

test: allow to perform Remez solving on an arbitrary range.

Location:

trunk/test/math

Files:

• remez-solver.h (modified) (13 diffs)
• remez.cpp (modified) (2 diffs)

trunk/test/math/remez-solver.h

@@ -19 +19 @@
     RemezSolver()
     {
-        ChebyInit();
-    }
-
-    void Run(RealFunc *func, RealFunc *error, int steps)
+    }
+
+    void Run(real a, real b, RealFunc *func, RealFunc *error, int steps)
     {
         m_func = func;
         m_error = error;
+        m_k1 = (b + a) >> 1;
+        m_k2 = (b - a) >> 1;
+        m_invk2 = re(m_k2);
+        m_invk1 = -m_k1 * m_invk2;
 
         Init();
 
-        ChebyCoeff();
-        for (int j = 0; j < ORDER + 1; j++)
-            printf("%s%14.12gx^%i", j && (bn[j] >= real::R_0) ? "+" : "", (double)bn[j], j);
-        printf("\n");
+        PrintPoly();
 
         for (int n = 0; n < steps; n++)
@@ -39 +39 @@
             Step();
 
-            ChebyCoeff();
-            for (int j = 0; j < ORDER + 1; j++)
-                printf("%s%14.12gx^%i", j && (bn[j] >= real::R_0) ? "+" : "", (double)bn[j], j);
-            printf("\n");
+            PrintPoly();
 
             FindZeroes();
@@ -50 +47 @@
         Step();
 
-        ChebyCoeff();
-        for (int j = 0; j < ORDER + 1; j++)
-            printf("%s%14.12gx^%i", j && (bn[j] >= real::R_0) ? "+" : "", (double)bn[j], j);
-        printf("\n");
-    }
-
-    /* Fill the Chebyshev tables */
-    void ChebyInit()
-    {
-        memset(cheby, 0, sizeof(cheby));
-
-        cheby[0][0] = 1;
-        cheby[1][1] = 1;
-
-        for (int i = 2; i < ORDER + 1; i++)
-        {
-            cheby[i][0] = -cheby[i - 2][0];
-            for (int j = 1; j < ORDER + 1; j++)
-                cheby[i][j] = 2 * cheby[i - 1][j - 1] - cheby[i - 2][j];
-        }
-    }
-
-    void ChebyCoeff()
-    {
-        for (int i = 0; i < ORDER + 1; i++)
-        {
-            bn[i] = 0;
-            for (int j = 0; j < ORDER + 1; j++)
-            if (cheby[j][i])
-                    bn[i] += coeff[j] * (real)cheby[j][i];
-        }
+        PrintPoly();
     }
 
@@ -91 +58 @@
             real mul = 0;
             for (int j = 0; j < ORDER + 1; j++)
-            if (cheby[j][i])
-                    mul += coeff[j] * (real)cheby[j][i];
+                mul += coeff[j] * (real)Cheby(j, i);
             ret += mul * xn;
             xn *= x;
@@ -107 +73 @@
         {
             zeroes[i] = (real)(2 * i - ORDER) / (real)(ORDER + 1);
-            fxn[i] = m_func(zeroes[i]);
+            fxn[i] = Value(zeroes[i]);
         }
 
@@ -126 +92 @@
                 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;
             }
@@ -149 +114 @@
         {
             real a = control[i];
-            real ea = ChebyEval(a) - m_func(a);
+            real ea = ChebyEval(a) - Value(a);
             real b = control[i + 1];
-            real eb = ChebyEval(b) - m_func(b);
+            real eb = ChebyEval(b) - Value(b);
 
             while (fabs(a - b) > (real)1e-140)
             {
                 real c = (a + b) * (real)0.5;
-                real ec = ChebyEval(c) - m_func(c);
+                real ec = ChebyEval(c) - Value(c);
 
                 if ((ea < (real)0 && ec < (real)0)
@@ -195 +160 @@
                 for (int k = 0; k <= 10; k++)
                 {
-                    real e = fabs(ChebyEval(c) - m_func(c));
+                    real e = fabs(ChebyEval(c) - Value(c));
                     if (e > maxerror)
                     {
@@ -223 +188 @@
         }
 
-        printf("Final error: %g\n", (double)final);
+        printf("Final error: ");
+        final.print(40);
     }
 
@@ -231 +197 @@
         real fxn[ORDER + 2];
         for (int i = 0; i < ORDER + 2; i++)
-            fxn[i] = m_func(control[i]);
+            fxn[i] = Value(control[i]);
 
         /* We build a matrix of Chebishev evaluations: row i contains the
@@ -249 +215 @@
                 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;
             }
             if (i & 1)
-                mat.m[i][ORDER + 1] = fabs(m_error(control[i]));
+                mat.m[i][ORDER + 1] = fabs(Error(control[i]));
             else
-                mat.m[i][ORDER + 1] = -fabs(m_error(control[i]));
+                mat.m[i][ORDER + 1] = -fabs(Error(control[i]));
         }
 
@@ -276 +241 @@
     }
 
-    int cheby[ORDER + 1][ORDER + 1];
-
-    /* ORDER + 1 chebyshev coefficients and 1 error value */
+    int Cheby(int n, int k)
+    {
+        if (k > n || k < 0)
+            return 0;
+        if (n <= 1)
+            return (n ^ k ^ 1) & 1;
+        return 2 * Cheby(n - 1, k - 1) - Cheby(n - 2, k);
+    }
+
+    int Comb(int n, int k)
+    {
+        if (k == 0 || k == n)
+            return 1;
+        return Comb(n - 1, k - 1) + Comb(n - 1, k);
+    }
+
+    void PrintPoly()
+    {
+        /* Transform Chebyshev polynomial weights into powers of X^i
+         * in the [-1..1] range. */
+        real bn[ORDER + 1];
+
+        for (int i = 0; i < ORDER + 1; i++)
+        {
+            bn[i] = 0;
+            for (int j = 0; j < ORDER + 1; j++)
+                bn[i] += coeff[j] * (real)Cheby(j, i);
+        }
+
+        /* Transform a polynomial in the [-1..1] range into a polynomial
+         * in the [a..b] range. */
+        real k1p[ORDER + 1], k2p[ORDER + 1];
+        real an[ORDER + 1];
+
+        for (int i = 0; i < ORDER + 1; i++)
+        {
+            k1p[i] = i ? k1p[i - 1] * m_invk1 : real::R_1;
+            k2p[i] = i ? k2p[i - 1] * m_invk2 : real::R_1;
+        }
+
+        for (int i = 0; i < ORDER + 1; i++)
+        {
+            an[i] = 0;
+            for (int j = i; j < ORDER + 1; j++)
+                an[i] += (real)Comb(j, i) * k1p[j - i] * bn[j];
+            an[i] *= k2p[i];
+        }
+
+        for (int j = 0; j < ORDER + 1; j++)
+            printf("%s%18.16gx^%i", j && (an[j] >= real::R_0) ? "+" : "", (double)an[j], j);
+        printf("\n");
+    }
+
+    real Value(real const &x)
+    {
+        return m_func(x * m_k2 + m_k1);
+    }
+
+    real Error(real const &x)
+    {
+        return m_error(x * m_k2 + m_k1);
+    }
+
+    /* ORDER + 1 Chebyshev coefficients and 1 error value */
     real coeff[ORDER + 2];
     /* ORDER + 1 zeroes of the error function */
@@ -285 +311 @@
     real control[ORDER + 2];
 
-    real bn[ORDER + 1];
-
 private:
-    RealFunc *m_func;
-    RealFunc *m_error;
+    RealFunc *m_func, *m_error;
+    real m_k1, m_k2, m_invk1, m_invk2;
 };
 

trunk/test/math/remez.cpp

@@ -27 +27 @@
 static real myfun(real const &x)
 {
-    return exp(x);
+    real y = sqrt(x);
+    if (!y)
+        return real::R_PI_2;
+    return sin(real::R_PI_2 * y) / y;
 }
 
-static real myerror(real const &x)
+static real myerr(real const &x)
 {
-    return myfun(x);
+    return real::R_1;
 }
 
@@ -38 +41 @@
 {
     RemezSolver<4> solver;
-
-    solver.Run(myfun, myerror, 10);
+    solver.Run(0, 1, myfun, myfun, 15);
+    //solver.Run(-1, 1, myfun, myfun, 15);
+    //solver.Run(0, real::R_PI * real::R_PI >> 4, myfun, myfun, 15);
 
     return EXIT_SUCCESS;