Changeset 3766


Ignore:
Timestamp:
Jan 12, 2015, 7:38:52 PM (7 years ago)
Author:
sam
Message:

lolremez: greatly improve root search times by using simple regula falsi.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/tools/lolremez/solver.cpp

    r3765 r3766  
    180180    for (int i = 0; i < m_order + 1; i++)
    181181    {
    182         struct { real value, error; } left, right, mid;
    183 
    184         left.value = m_control[i];
    185         left.error = EvalEstimate(left.value) - EvalFunc(left.value);
    186         right.value = m_control[i + 1];
    187         right.error = EvalEstimate(right.value) - EvalFunc(right.value);
     182        struct { real value, error; } a, b, c;
     183
     184        a.value = m_control[i];
     185        a.error = EvalEstimate(a.value) - EvalFunc(a.value);
     186        b.value = m_control[i + 1];
     187        b.error = EvalEstimate(b.value) - EvalFunc(b.value);
    188188
    189189        static real limit = ldexp((real)1, -500);
    190190        static real zero = (real)0;
    191         while (fabs(left.value - right.value) > limit)
     191        while (fabs(a.value - b.value) > limit)
    192192        {
    193             mid.value = (left.value + right.value) / 2;
    194             mid.error = EvalEstimate(mid.value) - EvalFunc(mid.value);
    195 
    196             if ((left.error <= zero && mid.error <= zero)
    197                  || (left.error >= zero && mid.error >= zero))
    198                 left = mid;
     193            /* Interpolate linearly instead of taking the midpoint, this
     194             * leads to far better convergence (6:1 speedups). */
     195            real t = abs(b.error) / (abs(a.error) + abs(b.error));
     196            real newc = b.value + t * (a.value - b.value);
     197
     198            /* If the third point didn't change since last iteration,
     199             * we may be at an inflection point. Use the midpoint to get
     200             * out of this situation. */
     201            c.value = newc == c.value ? (a.value + b.value) / 2 : newc;
     202            c.error = EvalEstimate(c.value) - EvalFunc(c.value);
     203
     204            if (c.error == zero)
     205                break;
     206
     207            if ((a.error < zero && c.error < zero)
     208                 || (a.error > zero && c.error > zero))
     209                a = c;
    199210            else
    200                 right = mid;
     211                b = c;
    201212        }
    202213
    203         m_zeroes[i] = mid.value;
     214        m_zeroes[i] = c.value;
    204215    }
    205216
Note: See TracChangeset for help on using the changeset viewer.