Changeset 991 for trunk/test/math


Ignore:
Timestamp:
Sep 28, 2011, 2:59:42 AM (11 years ago)
Author:
sam
Message:

test: the Remez algorithm is now almost functional.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/math/remez.cpp

    r989 r991  
    2020
    2121/* The order of the approximation we're looking for */
    22 static int const ORDER = 3;
     22static int const ORDER = 4;
    2323
    2424/* The function we want to approximate */
    25 double myfun(double x)
    26 {
    27     return exp(x);
    28 }
    29 
    3025real myfun(real const &x)
    3126{
    32     return exp(x);
     27    static real const one = 1.0;
     28    return cos(x) - one;
     29    //return exp(x);
    3330}
    3431
    3532/* Naive matrix inversion */
    36 template<int N> class Matrix
    37 {
    38 public:
     33template<int N> struct Matrix
     34{
    3935    inline Matrix() {}
    4036
     
    5450
    5551        /* Inversion method: iterate through all columns and make sure
    56          * all the terms are zero except on the diagonal where it is one */
     52         * all the terms are 1 on the diagonal and 0 everywhere else */
    5753        for (int i = 0; i < N; i++)
    5854        {
     
    118114
    119115/* Fill the Chebyshev tables */
    120 static void make_table()
     116static void cheby_init()
    121117{
    122118    memset(cheby, 0, sizeof(cheby));
     
    133129}
    134130
    135 int main(int argc, char **argv)
    136 {
    137     make_table();
    138 
    139     /* We start with ORDER+1 points and their images through myfun() */
    140     real xn[ORDER + 1];
     131static void cheby_coeff(real *coeff, real *bn)
     132{
     133    for (int i = 0; i < ORDER + 1; i++)
     134    {
     135        bn[i] = 0;
     136        for (int j = 0; j < ORDER + 1; j++)
     137            if (cheby[j][i])
     138                bn[i] += coeff[j] * (real)cheby[j][i];
     139    }
     140}
     141
     142static real cheby_eval(real *coeff, real const &x)
     143{
     144    real ret = 0.0, xn = 1.0;
     145
     146    for (int i = 0; i < ORDER + 1; i++)
     147    {
     148        real mul = 0;
     149        for (int j = 0; j < ORDER + 1; j++)
     150            if (cheby[j][i])
     151                mul += coeff[j] * (real)cheby[j][i];
     152        ret += mul * xn;
     153        xn *= x;
     154    }
     155
     156    return ret;
     157}
     158
     159static void remez_init(real *coeff, real *zeroes)
     160{
     161    /* Pick up x_i where error will be 0 and compute f(x_i) */
    141162    real fxn[ORDER + 1];
    142163    for (int i = 0; i < ORDER + 1; i++)
    143164    {
    144         //xn[i] = real(2 * i - ORDER) / real(ORDER + 1);
    145         xn[i] = real(2 * i - ORDER + 1) / real(ORDER - 1);
    146         fxn[i] = myfun(xn[i]);
    147     }
    148 
    149     /* We build a matrix of Chebishev evaluations: one row per point
    150      * in our array, and column i is the evaluation of the ith order
    151      * polynomial. */
     165        zeroes[i] = real(2 * i - ORDER) / real(ORDER + 1);
     166        fxn[i] = myfun(zeroes[i]);
     167    }
     168
     169    /* We build a matrix of Chebishev evaluations: row i contains the
     170     * evaluations of x_i for polynomial order n = 0, 1, ... */
    152171    Matrix<ORDER + 1> mat;
    153     for (int j = 0; j < ORDER + 1; j++)
    154     {
    155         /* Compute the powers of x_j */
     172    for (int i = 0; i < ORDER + 1; i++)
     173    {
     174        /* Compute the powers of x_i */
    156175        real powers[ORDER + 1];
    157176        powers[0] = 1.0;
    158         for (int i = 1; i < ORDER + 1; i++)
    159              powers[i] = powers[i - 1] * xn[j];
    160 
    161         /* Compute the Chebishev evaluations at x_j */
    162         for (int i = 0; i < ORDER + 1; i++)
     177        for (int n = 1; n < ORDER + 1; n++)
     178             powers[n] = powers[n - 1] * zeroes[i];
     179
     180        /* Compute the Chebishev evaluations at x_i */
     181        for (int n = 0; n < ORDER + 1; n++)
    163182        {
    164183            real sum = 0.0;
    165184            for (int k = 0; k < ORDER + 1; k++)
    166                 if (cheby[i][k])
    167                     sum += real(cheby[i][k]) * powers[k];
    168             mat.m[j][i] = sum;
    169         }
    170     }
    171 
    172     /* Invert the matrix and build interpolation coefficients */
     185                if (cheby[n][k])
     186                    sum += real(cheby[n][k]) * powers[k];
     187            mat.m[i][n] = sum;
     188        }
     189    }
     190
     191    /* Solve the system */
    173192    mat = mat.inv();
    174     real an[ORDER + 1];
     193
     194    /* Compute interpolation coefficients */
    175195    for (int j = 0; j < ORDER + 1; j++)
    176196    {
    177         an[j] = 0;
     197        coeff[j] = 0;
    178198        for (int i = 0; i < ORDER + 1; i++)
    179             an[j] += mat.m[j][i] * fxn[i];
    180         an[j].print(10);
    181     }
     199            coeff[j] += mat.m[j][i] * fxn[i];
     200    }
     201}
     202
     203static void remez_findzeroes(real *coeff, real *zeroes, real *control)
     204{
     205    /* FIXME: this is fake for now */
     206    for (int i = 0; i < ORDER + 1; i++)
     207    {
     208        real a = control[i];
     209        real ea = cheby_eval(coeff, a) - myfun(a);
     210        real b = control[i + 1];
     211        real eb = cheby_eval(coeff, b) - myfun(b);
     212
     213        while (fabs(a - b) > (real)1e-140)
     214        {
     215            real c = (a + b) * (real)0.5;
     216            real ec = cheby_eval(coeff, c) - myfun(c);
     217
     218            if ((ea < (real)0 && ec < (real)0)
     219                 || (ea > (real)0 && ec > (real)0))
     220            {
     221                a = c;
     222                ea = ec;
     223            }
     224            else
     225            {
     226                b = c;
     227                eb = ec;
     228            }
     229        }
     230
     231        zeroes[i] = a;
     232    }
     233}
     234
     235static void remez_finderror(real *coeff, real *zeroes, real *control)
     236{
     237    real final = 0;
     238
     239    for (int i = 0; i < ORDER + 2; i++)
     240    {
     241        real a = -1, b = 1;
     242        if (i > 0)
     243            a = zeroes[i - 1];
     244        if (i < ORDER + 1)
     245            b = zeroes[i];
     246
     247        printf("Error for [%g..%g]: ", (double)a, (double)b);
     248        for (;;)
     249        {
     250            real c = a, delta = (b - a) / (real)10.0;
     251            real maxerror = 0;
     252            int best = -1;
     253            for (int k = 0; k <= 10; k++)
     254            {
     255                real e = fabs(cheby_eval(coeff, c) - myfun(c));
     256                if (e > maxerror)
     257                {
     258                    maxerror = e;
     259                    best = k;
     260                }
     261                c += delta;
     262            }
     263
     264            if (best == 0)
     265                best = 1;
     266            if (best == 10)
     267                best = 9;
     268
     269            b = a + (real)(best + 1) * delta;
     270            a = a + (real)(best - 1) * delta;
     271
     272            if (b - a < (real)1e-15)
     273            {
     274                if (maxerror > final)
     275                    final = maxerror;
     276                control[i] = (a + b) * (real)0.5;
     277                printf("%g (in %g)\n", (double)maxerror, (double)control[i]);
     278                break;
     279            }
     280        }
     281    }
     282
     283    printf("Final error: %g\n", (double)final);
     284}
     285
     286static void remez_step(real *coeff, real *control)
     287{
     288    /* Pick up x_i where error will be 0 and compute f(x_i) */
     289    real fxn[ORDER + 2];
     290    for (int i = 0; i < ORDER + 2; i++)
     291        fxn[i] = myfun(control[i]);
     292
     293    /* We build a matrix of Chebishev evaluations: row i contains the
     294     * evaluations of x_i for polynomial order n = 0, 1, ... */
     295    Matrix<ORDER + 2> mat;
     296    for (int i = 0; i < ORDER + 2; i++)
     297    {
     298        /* Compute the powers of x_i */
     299        real powers[ORDER + 1];
     300        powers[0] = 1.0;
     301        for (int n = 1; n < ORDER + 1; n++)
     302             powers[n] = powers[n - 1] * control[i];
     303
     304        /* Compute the Chebishev evaluations at x_i */
     305        for (int n = 0; n < ORDER + 1; n++)
     306        {
     307            real sum = 0.0;
     308            for (int k = 0; k < ORDER + 1; k++)
     309                if (cheby[n][k])
     310                    sum += real(cheby[n][k]) * powers[k];
     311            mat.m[i][n] = sum;
     312        }
     313        mat.m[i][ORDER + 1] = (real)(-1 + (i & 1) * 2);
     314    }
     315
     316    /* Solve the system */
     317    mat = mat.inv();
     318
     319    /* Compute interpolation coefficients */
     320    for (int j = 0; j < ORDER + 1; j++)
     321    {
     322        coeff[j] = 0;
     323        for (int i = 0; i < ORDER + 2; i++)
     324            coeff[j] += mat.m[j][i] * fxn[i];
     325    }
     326
     327    /* Compute the error */
     328    real error = 0;
     329    for (int i = 0; i < ORDER + 2; i++)
     330        error += mat.m[ORDER + 1][i] * fxn[i];
     331}
     332
     333int main(int argc, char **argv)
     334{
     335    cheby_init();
     336
     337    /* ORDER + 1 chebyshev coefficients and 1 error value */
     338    real coeff[ORDER + 2];
     339    /* ORDER + 1 zeroes of the error function */
     340    real zeroes[ORDER + 1];
     341    /* ORDER + 2 control points */
     342    real control[ORDER + 2];
     343
     344    real bn[ORDER + 1];
     345
     346    remez_init(coeff, zeroes);
     347
     348    cheby_coeff(coeff, bn);
     349    for (int j = 0; j < ORDER + 1; j++)
     350        printf("%s%12.10gx^%i", j ? "+" : "", (double)bn[j], j);
     351    printf("\n");
     352
     353    for (int n = 0; n < 200; n++)
     354    {
     355        remez_finderror(coeff, zeroes, control);
     356        remez_step(coeff, control);
     357
     358        cheby_coeff(coeff, bn);
     359        for (int j = 0; j < ORDER + 1; j++)
     360            printf("%s%12.10gx^%i", j ? "+" : "", (double)bn[j], j);
     361        printf("\n");
     362
     363        remez_findzeroes(coeff, zeroes, control);
     364    }
     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");
    182373
    183374    return EXIT_SUCCESS;
Note: See TracChangeset for help on using the changeset viewer.