Changeset 2394


Ignore:
Timestamp:
Feb 11, 2013, 4:27:49 PM (7 years ago)
Author:
sam
Message:

math: minor improvements to the Remez exchange algorithm.

Location:
trunk/src
Files:
2 edited

Legend:

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

    r2183 r2394  
    5757            printf("Step %i error: ", n);
    5858            newerror.print(m_decimals);
     59            printf("\n");
    5960
    6061            Step();
     
    7778    }
    7879
    79     T ChebyEval(T const &x)
     80    T EvalCheby(T const &x)
    8081    {
    8182        T ret = 0.0, xn = 1.0;
     
    100101        {
    101102            zeroes[i] = (T)(2 * i - ORDER) / (T)(ORDER + 1);
    102             fxn[i] = Value(zeroes[i]);
     103            fxn[i] = EvalFunc(zeroes[i]);
    103104        }
    104105
     
    146147
    147148            left.value = control[i];
    148             left.error = ChebyEval(left.value) - Value(left.value);
     149            left.error = EvalCheby(left.value) - EvalFunc(left.value);
    149150            right.value = control[i + 1];
    150             right.error = ChebyEval(right.value) - Value(right.value);
     151            right.error = EvalCheby(right.value) - EvalFunc(right.value);
    151152
    152153            static T limit = ldexp((T)1, -500);
     
    155156            {
    156157                mid.value = (left.value + right.value) / 2;
    157                 mid.error = ChebyEval(mid.value) - Value(mid.value);
    158 
    159                 if ((left.error < zero && mid.error < zero)
    160                      || (left.error > zero && mid.error > zero))
     158                mid.error = EvalCheby(mid.value) - EvalFunc(mid.value);
     159
     160                if ((left.error <= zero && mid.error <= zero)
     161                     || (left.error >= zero && mid.error >= zero))
    161162                    left = mid;
    162163                else
     
    185186                b = zeroes[i];
    186187
    187             T maxerror = 0, maxweight = 0;
    188             int best = -1;
    189 
    190188            for (int round = 0; ; round++)
    191189            {
     190                T maxerror = 0, maxweight = 0;
     191                int best = -1;
     192
    192193                T c = a, delta = (b - a) / 4;
    193194                for (int k = 0; k <= 4; k++)
     
    195196                    if (round == 0 || (k & 1))
    196197                    {
    197                         T error = ChebyEval(c) - Value(c);
    198                         T weight = Weight(c);
    199                         if (fabs(error * maxweight) >= fabs(maxerror * weight))
     198                        T error = fabs(EvalCheby(c) - EvalFunc(c));
     199                        T weight = fabs(Weight(c));
     200                        /* if error/weight >= maxerror/maxweight */
     201                        if (error * maxweight >= maxerror * weight)
    200202                        {
    201203                            maxerror = error;
     
    218220                    b = a + delta * (best + 1);
    219221                    a = a + delta * (best - 1);
    220                     best = 2;
    221222                    break;
    222223                }
     
    224225                if (delta < m_epsilon)
    225226                {
    226                     T e = maxerror / maxweight;
     227                    T e = fabs(maxerror / maxweight);
    227228                    if (e > final)
    228229                        final = e;
     
    241242        T fxn[ORDER + 2];
    242243        for (int i = 0; i < ORDER + 2; i++)
    243             fxn[i] = Value(control[i]);
     244            fxn[i] = EvalFunc(control[i]);
    244245
    245246        /* We build a matrix of Chebishev evaluations: row i contains the
     
    335336        }
    336337
    337         printf("Polynomial estimate:\n");
     338        printf("Polynomial estimate: ");
    338339        for (int j = 0; j < ORDER + 1; j++)
    339340        {
    340341            if (j)
    341                 printf("+");
    342             printf("x**%i*", j);
     342                printf(" + x**%i * ", j);
    343343            an[j].print(m_decimals);
    344344        }
    345         printf("\n");
    346     }
    347 
    348     T Value(T const &x)
     345        printf("\n\n");
     346    }
     347
     348    T EvalFunc(T const &x)
    349349    {
    350350        return m_func(x * m_k2 + m_k1);
  • trunk/src/math/real.cpp

    r2322 r2394  
    13621362    for (int i = 0; i < BIGITS; i++)
    13631363        std::printf(" %08x", m_mantissa[i]);
    1364     std::printf("\n");
    13651364}
    13661365
     
    13711370    char *buf = new char[ndigits + 32 + 10];
    13721371    real::sprintf(buf, ndigits);
    1373     std::printf("%s\n", buf);
     1372    std::printf("%s", buf);
    13741373    delete[] buf;
    13751374}
     
    13871386    if (!x)
    13881387    {
    1389         std::strcpy(str, "0.0\n");
     1388        std::strcpy(str, "0.0");
    13901389        return;
    13911390    }
     
    14151414    /* Print exponent information */
    14161415    if (exponent)
    1417         str += std::sprintf(str, "e%c%i", exponent > 0 ? '+' : '-', -exponent);
     1416        str += std::sprintf(str, "e%c%i",
     1417                            exponent >= 0 ? '+' : '-', lol::abs(exponent));
    14181418
    14191419    *str++ = '\0';
Note: See TracChangeset for help on using the changeset viewer.