Changeset 1009


Ignore:
Timestamp:
Oct 5, 2011, 9:01:44 AM (11 years ago)
Author:
sam
Message:

test: some refactoring in the Remez solver to prepare multiple function
solving.

Location:
trunk/test/math
Files:
2 added
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/test/math/Makefile.am

    r985 r1009  
    1717pi_DEPENDENCIES = $(top_builddir)/src/liblol.a
    1818
    19 remez_SOURCES = remez.cpp
     19remez_SOURCES = remez.cpp remez-matrix.h remez-solver.h
    2020remez_CPPFLAGS = @LOL_CFLAGS@ @PIPI_CFLAGS@
    2121remez_LDFLAGS = $(top_builddir)/src/liblol.a @LOL_LIBS@ @PIPI_LIBS@

  • trunk/test/math/remez.cpp

    r1007 r1009  
    2121using namespace std;
    2222
    23 /* The order of the approximation we're looking for */
    24 static int const ORDER = 4;
     23#include "remez-matrix.h"
     24#include "remez-solver.h"
    2525
    2626/* The function we want to approximate */
    … …  
    2828{
    2929    return exp(x);
    30     //if (!x)
    31     //    return real::R_PI_2;
    32     //return sin(x * real::R_PI_2) / x;
    3330}
    3431
    … …  
    3633{
    3734    return myfun(x);
    38     //return real::R_1;
    39 }
    40 
    41 /* Naive matrix inversion */
    42 template<int N> struct Matrix
    43 {
    44     inline Matrix() {}
    45 
    46     Matrix(real x)
    47     {
    48         for (int j = 0; j < N; j++)
    49             for (int i = 0; i < N; i++)
    50                 if (i == j)
    51                     m[i][j] = x;
    52                 else
    53                     m[i][j] = 0;
    54     }
    55 
    56     Matrix<N> inv() const
    57     {
    58         Matrix a = *this, b((real)1.0);
    59 
    60         /* Inversion method: iterate through all columns and make sure
    61          * all the terms are 1 on the diagonal and 0 everywhere else */
    62         for (int i = 0; i < N; i++)
    63         {
    64             /* If the expected coefficient is zero, add one of
    65              * the other lines. The first we meet will do. */
    66             if ((double)a.m[i][i] == 0.0)
    67             {
    68                 for (int j = i + 1; j < N; j++)
    69                 {
    70                     if ((double)a.m[i][j] == 0.0)
    71                         continue;
    72                     /* Add row j to row i */
    73                     for (int n = 0; n < N; n++)
    74                     {
    75                         a.m[n][i] += a.m[n][j];
    76                         b.m[n][i] += b.m[n][j];
    77                     }
    78                     break;
    79                 }
    80             }
    81 
    82             /* Now we know the diagonal term is non-zero. Get its inverse
    83              * and use that to nullify all other terms in the column */
    84             real x = (real)1.0 / a.m[i][i];
    85             for (int j = 0; j < N; j++)
    86             {
    87                 if (j == i)
    88                     continue;
    89                 real mul = x * a.m[i][j];
    90                 for (int n = 0; n < N; n++)
    91                 {
    92                     a.m[n][j] -= mul * a.m[n][i];
    93                     b.m[n][j] -= mul * b.m[n][i];
    94                 }
    95             }
    96 
    97             /* Finally, ensure the diagonal term is 1 */
    98             for (int n = 0; n < N; n++)
    99             {
    100                 a.m[n][i] *= x;
    101                 b.m[n][i] *= x;
    102             }
    103         }
    104 
    105         return b;
    106     }
    107 
    108     void print() const
    109     {
    110         for (int j = 0; j < N; j++)
    111         {
    112             for (int i = 0; i < N; i++)
    113                 printf("%9.5f ", (double)m[j][i]);
    114             printf("\n");
    115         }
    116     }
    117 
    118     real m[N][N];
    119 };
    120 
    121 
    122 static int cheby[ORDER + 1][ORDER + 1];
    123 
    124 /* Fill the Chebyshev tables */
    125 static void cheby_init()
    126 {
    127     memset(cheby, 0, sizeof(cheby));
    128 
    129     cheby[0][0] = 1;
    130     cheby[1][1] = 1;
    131 
    132     for (int i = 2; i < ORDER + 1; i++)
    133     {
    134         cheby[i][0] = -cheby[i - 2][0];
    135         for (int j = 1; j < ORDER + 1; j++)
    136             cheby[i][j] = 2 * cheby[i - 1][j - 1] - cheby[i - 2][j];
    137     }
    138 }
    139 
    140 static void cheby_coeff(real *coeff, real *bn)
    141 {
    142     for (int i = 0; i < ORDER + 1; i++)
    143     {
    144         bn[i] = 0;
    145         for (int j = 0; j < ORDER + 1; j++)
    146             if (cheby[j][i])
    147                 bn[i] += coeff[j] * (real)cheby[j][i];
    148     }
    149 }
    150 
    151 static real cheby_eval(real *coeff, real const &x)
    152 {
    153     real ret = 0.0, xn = 1.0;
    154 
    155     for (int i = 0; i < ORDER + 1; i++)
    156     {
    157         real mul = 0;
    158         for (int j = 0; j < ORDER + 1; j++)
    159             if (cheby[j][i])
    160                 mul += coeff[j] * (real)cheby[j][i];
    161         ret += mul * xn;
    162         xn *= x;
    163     }
    164 
    165     return ret;
    166 }
    167 
    168 static void remez_init(real *coeff, real *zeroes)
    169 {
    170     /* Pick up x_i where error will be 0 and compute f(x_i) */
    171     real fxn[ORDER + 1];
    172     for (int i = 0; i < ORDER + 1; i++)
    173     {
    174         zeroes[i] = (real)(2 * i - ORDER) / (real)(ORDER + 1);
    175         fxn[i] = myfun(zeroes[i]);
    176     }
    177 
    178     /* We build a matrix of Chebishev evaluations: row i contains the
    179      * evaluations of x_i for polynomial order n = 0, 1, ... */
    180     Matrix<ORDER + 1> mat;
    181     for (int i = 0; i < ORDER + 1; i++)
    182     {
    183         /* Compute the powers of x_i */
    184         real powers[ORDER + 1];
    185         powers[0] = 1.0;
    186         for (int n = 1; n < ORDER + 1; n++)
    187              powers[n] = powers[n - 1] * zeroes[i];
    188 
    189         /* Compute the Chebishev evaluations at x_i */
    190         for (int n = 0; n < ORDER + 1; n++)
    191         {
    192             real sum = 0.0;
    193             for (int k = 0; k < ORDER + 1; k++)
    194                 if (cheby[n][k])
    195                     sum += (real)cheby[n][k] * powers[k];
    196             mat.m[i][n] = sum;
    197         }
    198     }
    199 
    200     /* Solve the system */
    201     mat = mat.inv();
    202 
    203     /* Compute interpolation coefficients */
    204     for (int j = 0; j < ORDER + 1; j++)
    205     {
    206         coeff[j] = 0;
    207         for (int i = 0; i < ORDER + 1; i++)
    208             coeff[j] += mat.m[j][i] * fxn[i];
    209     }
    210 }
    211 
    212 static void remez_findzeroes(real *coeff, real *zeroes, real *control)
    213 {
    214     for (int i = 0; i < ORDER + 1; i++)
    215     {
    216         real a = control[i];
    217         real ea = cheby_eval(coeff, a) - myfun(a);
    218         real b = control[i + 1];
    219         real eb = cheby_eval(coeff, b) - myfun(b);
    220 
    221         while (fabs(a - b) > (real)1e-140)
    222         {
    223             real c = (a + b) * (real)0.5;
    224             real ec = cheby_eval(coeff, c) - myfun(c);
    225 
    226             if ((ea < (real)0 && ec < (real)0)
    227                  || (ea > (real)0 && ec > (real)0))
    228             {
    229                 a = c;
    230                 ea = ec;
    231             }
    232             else
    233             {
    234                 b = c;
    235                 eb = ec;
    236             }
    237         }
    238 
    239         zeroes[i] = a;
    240     }
    241 }
    242 
    243 static void remez_finderror(real *coeff, real *zeroes, real *control)
    244 {
    245     real final = 0;
    246 
    247     for (int i = 0; i < ORDER + 2; i++)
    248     {
    249         real a = -1, b = 1;
    250         if (i > 0)
    251             a = zeroes[i - 1];
    252         if (i < ORDER + 1)
    253             b = zeroes[i];
    254 
    255         printf("Error for [%g..%g]: ", (double)a, (double)b);
    256         for (;;)
    257         {
    258             real c = a, delta = (b - a) / (real)10.0;
    259             real maxerror = 0;
    260             int best = -1;
    261             for (int k = 0; k <= 10; k++)
    262             {
    263                 real e = fabs(cheby_eval(coeff, c) - myfun(c));
    264                 if (e > maxerror)
    265                 {
    266                     maxerror = e;
    267                     best = k;
    268                 }
    269                 c += delta;
    270             }
    271 
    272             if (best == 0)
    273                 best = 1;
    274             if (best == 10)
    275                 best = 9;
    276 
    277             b = a + (real)(best + 1) * delta;
    278             a = a + (real)(best - 1) * delta;
    279 
    280             if (b - a < (real)1e-15)
    281             {
    282                 if (maxerror > final)
    283                     final = maxerror;
    284                 control[i] = (a + b) * (real)0.5;
    285                 printf("%g (at %g)\n", (double)maxerror, (double)control[i]);
    286                 break;
    287             }
    288         }
    289     }
    290 
    291     printf("Final error: %g\n", (double)final);
    292 }
    293 
    294 static void remez_step(real *coeff, real *control)
    295 {
    296     /* Pick up x_i where error will be 0 and compute f(x_i) */
    297     real fxn[ORDER + 2];
    298     for (int i = 0; i < ORDER + 2; i++)
    299         fxn[i] = myfun(control[i]);
    300 
    301     /* We build a matrix of Chebishev evaluations: row i contains the
    302      * evaluations of x_i for polynomial order n = 0, 1, ... */
    303     Matrix<ORDER + 2> mat;
    304     for (int i = 0; i < ORDER + 2; i++)
    305     {
    306         /* Compute the powers of x_i */
    307         real powers[ORDER + 1];
    308         powers[0] = 1.0;
    309         for (int n = 1; n < ORDER + 1; n++)
    310              powers[n] = powers[n - 1] * control[i];
    311 
    312         /* Compute the Chebishev evaluations at x_i */
    313         for (int n = 0; n < ORDER + 1; n++)
    314         {
    315             real sum = 0.0;
    316             for (int k = 0; k < ORDER + 1; k++)
    317                 if (cheby[n][k])
    318                     sum += (real)cheby[n][k] * powers[k];
    319             mat.m[i][n] = sum;
    320         }
    321         if (i & 1)
    322             mat.m[i][ORDER + 1] = fabs(myerror(control[i]));
    323         else
    324             mat.m[i][ORDER + 1] = -fabs(myerror(control[i]));
    325     }
    326 
    327     /* Solve the system */
    328     mat = mat.inv();
    329 
    330     /* Compute interpolation coefficients */
    331     for (int j = 0; j < ORDER + 1; j++)
    332     {
    333         coeff[j] = 0;
    334         for (int i = 0; i < ORDER + 2; i++)
    335             coeff[j] += mat.m[j][i] * fxn[i];
    336     }
    337 
    338     /* Compute the error */
    339     real error = 0;
    340     for (int i = 0; i < ORDER + 2; i++)
    341         error += mat.m[ORDER + 1][i] * fxn[i];
    34235}
    34336
    34437int main(int argc, char **argv)
    34538{
    346     cheby_init();
     39    RemezSolver<4> solver;
    34740
    348     /* ORDER + 1 chebyshev coefficients and 1 error value */
    349     real coeff[ORDER + 2];
    350     /* ORDER + 1 zeroes of the error function */
    351     real zeroes[ORDER + 1];
    352     /* ORDER + 2 control points */
    353     real control[ORDER + 2];
    354 
    355     real bn[ORDER + 1];
    356 
    357     remez_init(coeff, zeroes);
    358 
    359     cheby_coeff(coeff, bn);
    360     for (int j = 0; j < ORDER + 1; j++)
    361         printf("%s%12.10gx^%i", j ? "+" : "", (double)bn[j], j);
    362     printf("\n");
    363 
    364     for (int n = 0; n < 200; n++)
    365     {
    366         remez_finderror(coeff, zeroes, control);
    367         remez_step(coeff, control);
    368 
    369         cheby_coeff(coeff, bn);
    370         for (int j = 0; j < ORDER + 1; j++)
    371             printf("%s%12.10gx^%i", j ? "+" : "", (double)bn[j], j);
    372         printf("\n");
    373 
    374         remez_findzeroes(coeff, zeroes, control);
    375     }
    376 
    377     remez_finderror(coeff, zeroes, control);
    378     remez_step(coeff, control);
    379 
    380     cheby_coeff(coeff, bn);
    381     for (int j = 0; j < ORDER + 1; j++)
    382         printf("%s%12.10gx^%i", j ? "+" : "", (double)bn[j], j);
    383     printf("\n");
     41    solver.Run(myfun, myerror, 10);
    38442
    38543    return EXIT_SUCCESS;
Note: See TracChangeset for help on using the changeset viewer.