Ignore:
Timestamp:
Oct 10, 2011, 3:33:15 AM (9 years ago)
Author:
sam
Message:

test: various improvements to the Remez exchange solver.

File:
1 edited

Legend:

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

    r1010 r1014  
    2121    }
    2222
    23     void Run(real a, real b, RealFunc *func, RealFunc *error, int steps)
     23    void Run(real a, real b, RealFunc *func, RealFunc *weight, int steps)
    2424    {
    2525        m_func = func;
    26         m_error = error;
     26        m_weight = weight;
    2727        m_k1 = (b + a) >> 1;
    2828        m_k2 = (b - a) >> 1;
     
    3636        for (int n = 0; n < steps; n++)
    3737        {
    38             FindError();
     38            FindExtrema();
    3939            Step();
    4040
     
    4444        }
    4545
    46         FindError();
     46        FindExtrema();
    4747        Step();
    4848
     
    111111    void FindZeroes()
    112112    {
    113         for (int i = 0; i < ORDER + 1; i++)
    114         {
    115             real a = control[i];
    116             real ea = ChebyEval(a) - Value(a);
    117             real b = control[i + 1];
    118             real eb = ChebyEval(b) - Value(b);
    119 
    120             while (fabs(a - b) > (real)1e-140)
     113        /* Find ORDER + 1 zeroes of the error function. No need to
     114         * compute the relative error: its zeroes are at the same
     115         * place as the absolute error! */
     116        for (int i = 0; i < ORDER + 1; i++)
     117        {
     118            struct { real value, error; } left, right, mid;
     119
     120            left.value = control[i];
     121            left.error = ChebyEval(left.value) - Value(left.value);
     122            right.value = control[i + 1];
     123            right.error = ChebyEval(right.value) - Value(right.value);
     124
     125            static real limit = real::R_1 >> 500;
     126            while (fabs(left.value - right.value) > limit)
    121127            {
    122                 real c = (a + b) * (real)0.5;
    123                 real ec = ChebyEval(c) - Value(c);
    124 
    125                 if ((ea < (real)0 && ec < (real)0)
    126                      || (ea > (real)0 && ec > (real)0))
    127                 {
    128                     a = c;
    129                     ea = ec;
    130                 }
     128                mid.value = (left.value + right.value) >> 1;
     129                mid.error = ChebyEval(mid.value) - Value(mid.value);
     130
     131                if ((left.error < real::R_0 && mid.error < real::R_0)
     132                     || (left.error > real::R_0 && mid.error > real::R_0))
     133                    left = mid;
    131134                else
    132                 {
    133                     b = c;
    134                     eb = ec;
    135                 }
     135                    right = mid;
    136136            }
    137137
    138             zeroes[i] = a;
    139         }
    140     }
    141 
    142     void FindError()
    143     {
     138            zeroes[i] = mid.value;
     139        }
     140    }
     141
     142    void FindExtrema()
     143    {
     144        /* Find ORDER + 2 extrema of the error function. We need to
     145         * compute the relative error, since its extrema are at slightly
     146         * different locations than the absolute error’s. */
    144147        real final = 0;
    145148
     
    155158            for (;;)
    156159            {
    157                 real c = a, delta = (b - a) / (real)10.0;
     160                real c = a, delta = (b - a) >> 3;
    158161                real maxerror = 0;
     162                real maxweight = 0;
    159163                int best = -1;
    160                 for (int k = 0; k <= 10; k++)
     164                for (int k = 1; k <= 7; k++)
    161165                {
    162                     real e = fabs(ChebyEval(c) - Value(c));
    163                     if (e > maxerror)
     166                    real error = ChebyEval(c) - Value(c);
     167                    real weight = Weight(c);
     168                    if (fabs(error * maxweight) >= fabs(maxerror * weight))
    164169                    {
    165                         maxerror = e;
     170                        maxerror = error;
     171                        maxweight = weight;
    166172                        best = k;
    167173                    }
     
    169175                }
    170176
    171                 if (best == 0)
    172                     best = 1;
    173                 if (best == 10)
    174                     best = 9;
    175 
    176177                b = a + (real)(best + 1) * delta;
    177178                a = a + (real)(best - 1) * delta;
    178179
    179                 if (b - a < (real)1e-15)
     180                if (b - a < (real)1e-18)
    180181                {
    181                     if (maxerror > final)
    182                         final = maxerror;
    183                     control[i] = (a + b) * (real)0.5;
    184                     printf("%g (at %g)\n", (double)maxerror, (double)control[i]);
     182                    real e = maxerror / maxweight;
     183                    if (e > final)
     184                        final = e;
     185                    control[i] = (a + b) >> 1;
     186                    printf("%g (at %g)\n", (double)e, (double)control[i]);
    185187                    break;
    186188                }
     
    219221            }
    220222            if (i & 1)
    221                 mat.m[i][ORDER + 1] = fabs(Error(control[i]));
     223                mat.m[i][ORDER + 1] = fabs(Weight(control[i]));
    222224            else
    223                 mat.m[i][ORDER + 1] = -fabs(Error(control[i]));
     225                mat.m[i][ORDER + 1] = -fabs(Weight(control[i]));
    224226        }
    225227
     
    289291        }
    290292
     293        printf("Final polynomial:\n");
    291294        for (int j = 0; j < ORDER + 1; j++)
    292             printf("%s%18.16gx^%i", j && (an[j] >= real::R_0) ? "+" : "", (double)an[j], j);
     295        {
     296            if (j)
     297                printf("+");
     298            printf("x^%i*", j);
     299            an[j].print(40);
     300        }
    293301        printf("\n");
    294302    }
     
    299307    }
    300308
    301     real Error(real const &x)
    302     {
    303         return m_error(x * m_k2 + m_k1);
     309    real Weight(real const &x)
     310    {
     311        return m_weight(x * m_k2 + m_k1);
    304312    }
    305313
     
    312320
    313321private:
    314     RealFunc *m_func, *m_error;
     322    RealFunc *m_func, *m_weight;
    315323    real m_k1, m_k2, m_invk1, m_invk2;
    316324};
Note: See TracChangeset for help on using the changeset viewer.