Changeset 989 for trunk/test/math


Ignore:
Timestamp:
Sep 27, 2011, 7:09:35 PM (9 years ago)
Author:
sam
Message:

test: more work on the Remez exchange algorithm.

Location:
trunk/test/math
Files:
2 edited

Legend:

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

    r984 r989  
    3333    sum.print();
    3434
     35    /* Bonus: compute e for fun. */
     36    sum = 0.0;
     37    x0 = 1.0;
     38
     39    for (int i = 1; i < 100; i++)
     40    {
     41        sum += fres(x0);
     42        x0 *= (real)i;
     43    }
     44
     45    sum.print();
     46
    3547    return EXIT_SUCCESS;
    3648}
  • trunk/test/math/remez.cpp

    r985 r989  
    1919using namespace lol;
    2020
    21 static int const ORDER = 12;
     21/* The order of the approximation we're looking for */
     22static int const ORDER = 3;
    2223
    23 static int cheby[ORDER][ORDER];
     24/* The function we want to approximate */
     25double myfun(double x)
     26{
     27    return exp(x);
     28}
    2429
     30real myfun(real const &x)
     31{
     32    return exp(x);
     33}
     34
     35/* Naive matrix inversion */
     36template<int N> class Matrix
     37{
     38public:
     39    inline Matrix() {}
     40
     41    Matrix(real x)
     42    {
     43        for (int j = 0; j < N; j++)
     44            for (int i = 0; i < N; i++)
     45                if (i == j)
     46                    m[i][j] = x;
     47                else
     48                    m[i][j] = 0;
     49    }
     50
     51    Matrix<N> inv() const
     52    {
     53        Matrix a = *this, b(real(1.0));
     54
     55        /* Inversion method: iterate through all columns and make sure
     56         * all the terms are zero except on the diagonal where it is one */
     57        for (int i = 0; i < N; i++)
     58        {
     59            /* If the expected coefficient is zero, add one of
     60             * the other lines. The first we meet will do. */
     61            if ((double)a.m[i][i] == 0.0)
     62            {
     63                for (int j = i + 1; j < N; j++)
     64                {
     65                    if ((double)a.m[i][j] == 0.0)
     66                        continue;
     67                    /* Add row j to row i */
     68                    for (int n = 0; n < N; n++)
     69                    {
     70                        a.m[n][i] += a.m[n][j];
     71                        b.m[n][i] += b.m[n][j];
     72                    }
     73                    break;
     74                }
     75            }
     76
     77            /* Now we know the diagonal term is non-zero. Get its inverse
     78             * and use that to nullify all other terms in the column */
     79            real x = fres(a.m[i][i]);
     80            for (int j = 0; j < N; j++)
     81            {
     82                if (j == i)
     83                    continue;
     84                real mul = x * a.m[i][j];
     85                for (int n = 0; n < N; n++)
     86                {
     87                    a.m[n][j] -= mul * a.m[n][i];
     88                    b.m[n][j] -= mul * b.m[n][i];
     89                }
     90            }
     91
     92            /* Finally, ensure the diagonal term is 1 */
     93            for (int n = 0; n < N; n++)
     94            {
     95                a.m[n][i] *= x;
     96                b.m[n][i] *= x;
     97            }
     98        }
     99
     100        return b;
     101    }
     102
     103    void print() const
     104    {
     105        for (int j = 0; j < N; j++)
     106        {
     107            for (int i = 0; i < N; i++)
     108                printf("%9.5f ", (double)m[j][i]);
     109            printf("\n");
     110        }
     111    }
     112
     113    real m[N][N];
     114};
     115
     116
     117static int cheby[ORDER + 1][ORDER + 1];
     118
     119/* Fill the Chebyshev tables */
    25120static void make_table()
    26121{
    27     memset(cheby[0], 0, ORDER * sizeof(int));
     122    memset(cheby, 0, sizeof(cheby));
     123
    28124    cheby[0][0] = 1;
    29     memset(cheby[1], 0, ORDER * sizeof(int));
    30125    cheby[1][1] = 1;
    31126
    32     for (int i = 2; i < ORDER; i++)
     127    for (int i = 2; i < ORDER + 1; i++)
    33128    {
    34129        cheby[i][0] = -cheby[i - 2][0];
    35         for (int j = 1; j < ORDER; j++)
     130        for (int j = 1; j < ORDER + 1; j++)
    36131            cheby[i][j] = 2 * cheby[i - 1][j - 1] - cheby[i - 2][j];
    37132    }
     
    42137    make_table();
    43138
    44 for (int i = 0; i < ORDER; i++)
    45 {
    46     for (int j = 0; j < ORDER; j++)
    47         printf("% 5i ", cheby[i][j]);
    48     printf("\n");
    49 }
     139    /* We start with ORDER+1 points and their images through myfun() */
     140    real xn[ORDER + 1];
     141    real fxn[ORDER + 1];
     142    for (int i = 0; i < ORDER + 1; i++)
     143    {
     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. */
     152    Matrix<ORDER + 1> mat;
     153    for (int j = 0; j < ORDER + 1; j++)
     154    {
     155        /* Compute the powers of x_j */
     156        real powers[ORDER + 1];
     157        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++)
     163        {
     164            real sum = 0.0;
     165            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 */
     173    mat = mat.inv();
     174    real an[ORDER + 1];
     175    for (int j = 0; j < ORDER + 1; j++)
     176    {
     177        an[j] = 0;
     178        for (int i = 0; i < ORDER + 1; i++)
     179            an[j] += mat.m[j][i] * fxn[i];
     180        an[j].print(10);
     181    }
    50182
    51183    return EXIT_SUCCESS;
Note: See TracChangeset for help on using the changeset viewer.