Changeset 1126


Ignore:
Timestamp:
Jan 9, 2012, 2:30:13 PM (8 years ago)
Author:
sam
Message:

math: write a faster factorial method for use in exp() and sin(). These
functions are now about 20% faster.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r1123 r1126  
    694694}
    695695
     696static real fast_fact(int x)
     697{
     698    real ret = real::R_1;
     699    int i = 1, multiplier = 1, exponent = 0;
     700
     701    for (;;)
     702    {
     703        if (i++ >= x)
     704            /* Multiplication is a no-op if multiplier == 1 */
     705            return ldexp(ret * multiplier, exponent);
     706
     707        int tmp = i;
     708        while ((tmp & 1) == 0)
     709        {
     710            tmp >>= 1;
     711            exponent++;
     712        }
     713        if (multiplier * tmp / tmp != multiplier)
     714        {
     715            ret *= multiplier;
     716            multiplier = 1;
     717        }
     718        multiplier *= tmp;
     719    }
     720}
     721
    696722real gamma(real const &x)
    697723{
     
    803829     * domain of validity. The argument y is used for cases where we
    804830     * don't want the leading 1 in the Taylor series. */
    805     real ret = real::R_1 - y, fact = real::R_1, xn = x;
    806 
    807     for (int i = 1; ; i++)
     831    real ret = real::R_1 - y, xn = x;
     832    int i = 1;
     833
     834    for (;;)
    808835    {
    809836        real newret = ret + xn;
    810837        if (newret == ret)
    811838            break;
    812         ret = newret;
    813         real mul = (i + 1);
    814         fact *= mul;
    815         ret *= mul;
     839        ret = newret * ++i;
    816840        xn *= x;
    817841    }
    818     ret /= fact;
    819 
    820     return ret;
     842
     843    return ret / fast_fact(i);
    821844}
    822845
     
    10051028
    10061029    real ret = real::R_0, fact = real::R_1, xn = absx, mx2 = -absx * absx;
    1007     for (int i = 1; ; i += 2)
     1030    int i = 1;
     1031    for (;;)
    10081032    {
    10091033        real newret = ret + xn;
    10101034        if (newret == ret)
    10111035            break;
    1012         ret = newret;
    1013         real mul = (i + 1) * (i + 2);
    1014         fact *= mul;
    1015         ret *= mul;
     1036        ret = newret * ((i + 1) * (i + 2));
    10161037        xn *= mx2;
    1017     }
    1018     ret /= fact;
     1038        i += 2;
     1039    }
     1040    ret /= fast_fact(i);
    10191041
    10201042    /* Propagate sign */
     
    12871309{
    12881310    /* Approximate Pi using Machin's formula: 16*atan(1/5)-4*atan(1/239) */
    1289     real ret = 0.0, x0 = 5.0, x1 = 239.0;
    1290     real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16.0, r4 = 4.0;
     1311    real ret = 0, x0 = 5, x1 = 239;
     1312    real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16, r4 = 4;
    12911313
    12921314    for (int i = 1; ; i += 2)
     
    13091331real const real::R_10       = (real)10.0;
    13101332
     1333/*
     1334 * Initialisation order is important here:
     1335 *  - fast_log() requires R_1
     1336 *  - log() requires R_LN2
     1337 *  - re() require R_2
     1338 *  - exp() requires R_0, R_1, R_LN2
     1339 *  - sqrt() requires R_3
     1340 */
    13111341real const real::R_LN2      = fast_log(R_2);
    13121342real const real::R_LN10     = log(R_10);
  • trunk/test/math/remez.cpp

    r1119 r1126  
    2525/* See the tutorial at http://lol.zoy.org/wiki/doc/maths/remez */
    2626real f(real const &x) { return exp(x); }
     27real g(real const &x) { return exp(x); }
    2728
    2829int main(int argc, char **argv)
    2930{
    3031    RemezSolver<4, real> solver;
    31     solver.Run(-1, 1, f, 30);
     32    solver.Run(-1, 1, f, g, 30);
    3233    return 0;
    3334}
Note: See TracChangeset for help on using the changeset viewer.