Changeset 1127


Ignore:
Timestamp:
Jan 9, 2012, 11:49:28 PM (7 years ago)
Author:
sam
Message:

math: significant performance improvements in the Remez solver.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/lol/math/remez.h

    r1122 r1127  
    3333    }
    3434
    35     void Run(T a, T b, RealFunc *func, RealFunc *weight, int steps)
     35    void Run(T a, T b, RealFunc *func, RealFunc *weight, int decimals)
    3636    {
    3737        m_func = func;
     
    4141        m_invk2 = re(m_k2);
    4242        m_invk1 = -m_k1 * m_invk2;
     43        m_decimals = decimals;
     44        m_epsilon = pow((T)10, (T)-(decimals + 2));
    4345
    4446        Init();
     
    4648        PrintPoly();
    4749
    48         for (int n = 0; n < steps; n++)
    49         {
    50             FindExtrema();
     50        T error = -1;
     51
     52        for (int n = 0; ; n++)
     53        {
     54            T newerror = FindExtrema();
     55            printf("Step %i error: ", n);
     56            newerror.print(m_decimals);
     57
    5158            Step();
    5259
     60            if (error >= (T)0 && fabs(newerror - error) < error * m_epsilon)
     61                break;
     62            error = newerror;
     63
    5364            PrintPoly();
    5465
     
    5667        }
    5768
    58         FindExtrema();
    59         Step();
    60 
    6169        PrintPoly();
    6270    }
    6371
    64     inline void Run(T a, T b, RealFunc *func, int steps)
    65     {
    66         Run(a, b, func, NULL, steps);
     72    inline void Run(T a, T b, RealFunc *func, int decimals)
     73    {
     74        Run(a, b, func, NULL, decimals);
    6775    }
    6876
     
    158166    }
    159167
    160     void FindExtrema()
     168    real FindExtrema()
    161169    {
    162170        using std::printf;
     
    175183                b = zeroes[i];
    176184
    177             for (;;)
     185            T maxerror = 0, maxweight = 0;
     186            int best = -1;
     187
     188            for (int round = 0; ; round++)
    178189            {
    179                 T c = a, delta = (b - a) / 8;
    180                 T maxerror = 0;
    181                 T maxweight = 0;
    182                 int best = -1;
    183                 for (int k = 1; k <= 7; k++)
     190                T c = a, delta = (b - a) / 4;
     191                for (int k = 0; k <= 4; k++)
    184192                {
    185                     T error = ChebyEval(c) - Value(c);
    186                     T weight = Weight(c);
    187                     if (fabs(error * maxweight) >= fabs(maxerror * weight))
     193                    if (round == 0 || (k & 1))
    188194                    {
    189                         maxerror = error;
    190                         maxweight = weight;
    191                         best = k;
     195                        T error = ChebyEval(c) - Value(c);
     196                        T weight = Weight(c);
     197                        if (fabs(error * maxweight) >= fabs(maxerror * weight))
     198                        {
     199                            maxerror = error;
     200                            maxweight = weight;
     201                            best = k;
     202                        }
    192203                    }
    193204                    c += delta;
    194205                }
    195206
    196                 b = a + (T)(best + 1) * delta;
    197                 a = a + (T)(best - 1) * delta;
    198 
    199                 if (b - a < (T)1e-18)
     207                switch (best)
     208                {
     209                case 0:
     210                    b = a + delta * 2;
     211                    break;
     212                case 4:
     213                    a = b - delta * 2;
     214                    break;
     215                default:
     216                    b = a + delta * (best + 1);
     217                    a = a + delta * (best - 1);
     218                    best = 2;
     219                    break;
     220                }
     221
     222                if (delta < m_epsilon)
    200223                {
    201224                    T e = maxerror / maxweight;
     
    208231        }
    209232
    210         printf("Final error: ");
    211         final.print(40);
     233        return final;
    212234    }
    213235
     
    317339                printf("+");
    318340            printf("x**%i*", j);
    319             an[j].print(40);
     341            an[j].print(m_decimals);
    320342        }
    321343        printf("\n");
     
    343365private:
    344366    RealFunc *m_func, *m_weight;
    345     T m_k1, m_k2, m_invk1, m_invk2;
     367    T m_k1, m_k2, m_invk1, m_invk2, m_epsilon;
     368    int m_decimals;
    346369};
    347370
  • trunk/test/math/remez.cpp

    r1126 r1127  
    3030{
    3131    RemezSolver<4, real> solver;
    32     solver.Run(-1, 1, f, g, 30);
     32    solver.Run(-1, 1, f, g, 40);
    3333    return 0;
    3434}
Note: See TracChangeset for help on using the changeset viewer.