Changeset 1010


Ignore:
Timestamp:
Oct 5, 2011, 7:01:33 PM (9 years ago)
Author:
sam
Message:

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

Location:
trunk/test/math
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/math/remez-solver.h

    r1009 r1010  
    1919    RemezSolver()
    2020    {
    21         ChebyInit();
    22     }
    23 
    24     void Run(RealFunc *func, RealFunc *error, int steps)
     21    }
     22
     23    void Run(real a, real b, RealFunc *func, RealFunc *error, int steps)
    2524    {
    2625        m_func = func;
    2726        m_error = error;
     27        m_k1 = (b + a) >> 1;
     28        m_k2 = (b - a) >> 1;
     29        m_invk2 = re(m_k2);
     30        m_invk1 = -m_k1 * m_invk2;
    2831
    2932        Init();
    3033
    31         ChebyCoeff();
    32         for (int j = 0; j < ORDER + 1; j++)
    33             printf("%s%14.12gx^%i", j && (bn[j] >= real::R_0) ? "+" : "", (double)bn[j], j);
    34         printf("\n");
     34        PrintPoly();
    3535
    3636        for (int n = 0; n < steps; n++)
     
    3939            Step();
    4040
    41             ChebyCoeff();
    42             for (int j = 0; j < ORDER + 1; j++)
    43                 printf("%s%14.12gx^%i", j && (bn[j] >= real::R_0) ? "+" : "", (double)bn[j], j);
    44             printf("\n");
     41            PrintPoly();
    4542
    4643            FindZeroes();
     
    5047        Step();
    5148
    52         ChebyCoeff();
    53         for (int j = 0; j < ORDER + 1; j++)
    54             printf("%s%14.12gx^%i", j && (bn[j] >= real::R_0) ? "+" : "", (double)bn[j], j);
    55         printf("\n");
    56     }
    57 
    58     /* Fill the Chebyshev tables */
    59     void ChebyInit()
    60     {
    61         memset(cheby, 0, sizeof(cheby));
    62 
    63         cheby[0][0] = 1;
    64         cheby[1][1] = 1;
    65 
    66         for (int i = 2; i < ORDER + 1; i++)
    67         {
    68             cheby[i][0] = -cheby[i - 2][0];
    69             for (int j = 1; j < ORDER + 1; j++)
    70                 cheby[i][j] = 2 * cheby[i - 1][j - 1] - cheby[i - 2][j];
    71         }
    72     }
    73 
    74     void ChebyCoeff()
    75     {
    76         for (int i = 0; i < ORDER + 1; i++)
    77         {
    78             bn[i] = 0;
    79             for (int j = 0; j < ORDER + 1; j++)
    80             if (cheby[j][i])
    81                     bn[i] += coeff[j] * (real)cheby[j][i];
    82         }
     49        PrintPoly();
    8350    }
    8451
     
    9158            real mul = 0;
    9259            for (int j = 0; j < ORDER + 1; j++)
    93             if (cheby[j][i])
    94                     mul += coeff[j] * (real)cheby[j][i];
     60                mul += coeff[j] * (real)Cheby(j, i);
    9561            ret += mul * xn;
    9662            xn *= x;
     
    10773        {
    10874            zeroes[i] = (real)(2 * i - ORDER) / (real)(ORDER + 1);
    109             fxn[i] = m_func(zeroes[i]);
     75            fxn[i] = Value(zeroes[i]);
    11076        }
    11177
     
    12692                real sum = 0.0;
    12793                for (int k = 0; k < ORDER + 1; k++)
    128                     if (cheby[n][k])
    129                         sum += (real)cheby[n][k] * powers[k];
     94                    sum += (real)Cheby(n, k) * powers[k];
    13095                mat.m[i][n] = sum;
    13196            }
     
    149114        {
    150115            real a = control[i];
    151             real ea = ChebyEval(a) - m_func(a);
     116            real ea = ChebyEval(a) - Value(a);
    152117            real b = control[i + 1];
    153             real eb = ChebyEval(b) - m_func(b);
     118            real eb = ChebyEval(b) - Value(b);
    154119
    155120            while (fabs(a - b) > (real)1e-140)
    156121            {
    157122                real c = (a + b) * (real)0.5;
    158                 real ec = ChebyEval(c) - m_func(c);
     123                real ec = ChebyEval(c) - Value(c);
    159124
    160125                if ((ea < (real)0 && ec < (real)0)
     
    195160                for (int k = 0; k <= 10; k++)
    196161                {
    197                     real e = fabs(ChebyEval(c) - m_func(c));
     162                    real e = fabs(ChebyEval(c) - Value(c));
    198163                    if (e > maxerror)
    199164                    {
     
    223188        }
    224189
    225         printf("Final error: %g\n", (double)final);
     190        printf("Final error: ");
     191        final.print(40);
    226192    }
    227193
     
    231197        real fxn[ORDER + 2];
    232198        for (int i = 0; i < ORDER + 2; i++)
    233             fxn[i] = m_func(control[i]);
     199            fxn[i] = Value(control[i]);
    234200
    235201        /* We build a matrix of Chebishev evaluations: row i contains the
     
    249215                real sum = 0.0;
    250216                for (int k = 0; k < ORDER + 1; k++)
    251                     if (cheby[n][k])
    252                         sum += (real)cheby[n][k] * powers[k];
     217                    sum += (real)Cheby(n, k) * powers[k];
    253218                mat.m[i][n] = sum;
    254219            }
    255220            if (i & 1)
    256                 mat.m[i][ORDER + 1] = fabs(m_error(control[i]));
     221                mat.m[i][ORDER + 1] = fabs(Error(control[i]));
    257222            else
    258                 mat.m[i][ORDER + 1] = -fabs(m_error(control[i]));
     223                mat.m[i][ORDER + 1] = -fabs(Error(control[i]));
    259224        }
    260225
     
    276241    }
    277242
    278     int cheby[ORDER + 1][ORDER + 1];
    279 
    280     /* ORDER + 1 chebyshev coefficients and 1 error value */
     243    int Cheby(int n, int k)
     244    {
     245        if (k > n || k < 0)
     246            return 0;
     247        if (n <= 1)
     248            return (n ^ k ^ 1) & 1;
     249        return 2 * Cheby(n - 1, k - 1) - Cheby(n - 2, k);
     250    }
     251
     252    int Comb(int n, int k)
     253    {
     254        if (k == 0 || k == n)
     255            return 1;
     256        return Comb(n - 1, k - 1) + Comb(n - 1, k);
     257    }
     258
     259    void PrintPoly()
     260    {
     261        /* Transform Chebyshev polynomial weights into powers of X^i
     262         * in the [-1..1] range. */
     263        real bn[ORDER + 1];
     264
     265        for (int i = 0; i < ORDER + 1; i++)
     266        {
     267            bn[i] = 0;
     268            for (int j = 0; j < ORDER + 1; j++)
     269                bn[i] += coeff[j] * (real)Cheby(j, i);
     270        }
     271
     272        /* Transform a polynomial in the [-1..1] range into a polynomial
     273         * in the [a..b] range. */
     274        real k1p[ORDER + 1], k2p[ORDER + 1];
     275        real an[ORDER + 1];
     276
     277        for (int i = 0; i < ORDER + 1; i++)
     278        {
     279            k1p[i] = i ? k1p[i - 1] * m_invk1 : real::R_1;
     280            k2p[i] = i ? k2p[i - 1] * m_invk2 : real::R_1;
     281        }
     282
     283        for (int i = 0; i < ORDER + 1; i++)
     284        {
     285            an[i] = 0;
     286            for (int j = i; j < ORDER + 1; j++)
     287                an[i] += (real)Comb(j, i) * k1p[j - i] * bn[j];
     288            an[i] *= k2p[i];
     289        }
     290
     291        for (int j = 0; j < ORDER + 1; j++)
     292            printf("%s%18.16gx^%i", j && (an[j] >= real::R_0) ? "+" : "", (double)an[j], j);
     293        printf("\n");
     294    }
     295
     296    real Value(real const &x)
     297    {
     298        return m_func(x * m_k2 + m_k1);
     299    }
     300
     301    real Error(real const &x)
     302    {
     303        return m_error(x * m_k2 + m_k1);
     304    }
     305
     306    /* ORDER + 1 Chebyshev coefficients and 1 error value */
    281307    real coeff[ORDER + 2];
    282308    /* ORDER + 1 zeroes of the error function */
     
    285311    real control[ORDER + 2];
    286312
    287     real bn[ORDER + 1];
    288 
    289313private:
    290     RealFunc *m_func;
    291     RealFunc *m_error;
     314    RealFunc *m_func, *m_error;
     315    real m_k1, m_k2, m_invk1, m_invk2;
    292316};
    293317
  • trunk/test/math/remez.cpp

    r1009 r1010  
    2727static real myfun(real const &x)
    2828{
    29     return exp(x);
     29    real y = sqrt(x);
     30    if (!y)
     31        return real::R_PI_2;
     32    return sin(real::R_PI_2 * y) / y;
    3033}
    3134
    32 static real myerror(real const &x)
     35static real myerr(real const &x)
    3336{
    34     return myfun(x);
     37    return real::R_1;
    3538}
    3639
     
    3841{
    3942    RemezSolver<4> solver;
    40 
    41     solver.Run(myfun, myerror, 10);
     43    solver.Run(0, 1, myfun, myfun, 15);
     44    //solver.Run(-1, 1, myfun, myfun, 15);
     45    //solver.Run(0, real::R_PI * real::R_PI >> 4, myfun, myfun, 15);
    4246
    4347    return EXIT_SUCCESS;
Note: See TracChangeset for help on using the changeset viewer.