Ignore:
Timestamp:
Sep 9, 2012, 4:55:26 PM (8 years ago)
Author:
sam
Message:

math: refactor real number constant declarations so that they are only
computed on demand with static initialisation.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/math/real.cpp

    r1667 r1893  
    3737namespace lol
    3838{
     39
     40/*
     41 * First handle explicit specialisation of our templates.
     42 *
     43 * Initialisation order is not important because everything is
     44 * done on demand, but here is the dependency list anyway:
     45 *  - fast_log() requires R_1
     46 *  - log() requires R_LN2
     47 *  - re() require R_2
     48 *  - exp() requires R_0, R_1, R_LN2
     49 *  - sqrt() requires R_3
     50 */
     51
     52static real fast_log(real const &x);
     53static real fast_pi();
     54
     55#define LOL_CONSTANT_GETTER(name, value) \
     56    template<> real const& real::name() \
     57    { \
     58        static real const ret = value; \
     59        return ret; \
     60    }
     61
     62LOL_CONSTANT_GETTER(R_0,        (real)0.0);
     63LOL_CONSTANT_GETTER(R_1,        (real)1.0);
     64LOL_CONSTANT_GETTER(R_2,        (real)2.0);
     65LOL_CONSTANT_GETTER(R_3,        (real)3.0);
     66LOL_CONSTANT_GETTER(R_10,       (real)10.0);
     67
     68LOL_CONSTANT_GETTER(R_LN2,      fast_log(R_2()));
     69LOL_CONSTANT_GETTER(R_LN10,     log(R_10()));
     70LOL_CONSTANT_GETTER(R_LOG2E,    re(R_LN2()));
     71LOL_CONSTANT_GETTER(R_LOG10E,   re(R_LN10()));
     72LOL_CONSTANT_GETTER(R_E,        exp(R_1()));
     73LOL_CONSTANT_GETTER(R_PI,       fast_pi());
     74LOL_CONSTANT_GETTER(R_PI_2,     R_PI() / 2);
     75LOL_CONSTANT_GETTER(R_PI_3,     R_PI() / R_3());
     76LOL_CONSTANT_GETTER(R_PI_4,     R_PI() / 4);
     77LOL_CONSTANT_GETTER(R_1_PI,     re(R_PI()));
     78LOL_CONSTANT_GETTER(R_2_PI,     R_1_PI() * 2);
     79LOL_CONSTANT_GETTER(R_2_SQRTPI, re(sqrt(R_PI())) * 2);
     80LOL_CONSTANT_GETTER(R_SQRT2,    sqrt(R_2()));
     81LOL_CONSTANT_GETTER(R_SQRT3,    sqrt(R_3()));
     82LOL_CONSTANT_GETTER(R_SQRT1_2,  R_SQRT2() / 2);
     83
     84#undef LOL_CONSTANT_GETTER
     85
     86/*
     87 * Now carry on with the rest of the Real class.
     88 */
    3989
    4090template<> real::Real()
     
    190240
    191241    if (exponent)
    192         ret *= pow(R_10, (real)exponent);
     242        ret *= pow(R_10(), (real)exponent);
    193243
    194244    if (negative)
     
    602652     * convergence, but this hasn't been checked seriously. */
    603653    for (int i = 1; i <= real::BIGITS; i *= 2)
    604         ret = ret * (real::R_2 - ret * x);
     654        ret = ret * (real::R_2() - ret * x);
    605655
    606656    return ret;
     
    647697    for (int i = 1; i <= real::BIGITS; i *= 2)
    648698    {
    649         ret = ret * (real::R_3 - ret * ret * x);
     699        ret = ret * (real::R_3() - ret * ret * x);
    650700        ret.m_signexp--;
    651701    }
     
    684734    for (int i = 1; i <= real::BIGITS; i *= 2)
    685735    {
    686         static real third = re(real::R_3);
     736        static real third = re(real::R_3());
    687737        ret = third * (x / (ret * ret) + (ret / 2));
    688738    }
     
    694744{
    695745    if (!y)
    696         return real::R_1;
     746        return real::R_1();
    697747    if (!x)
    698         return real::R_0;
    699     if (x > real::R_0)
     748        return real::R_0();
     749    if (x > real::R_0())
    700750        return exp(y * log(x));
    701751    else /* x < 0 */
     
    710760
    711761        /* FIXME: negative nth root */
    712         return real::R_0;
     762        return real::R_0();
    713763    }
    714764}
     
    716766static real fast_fact(int x)
    717767{
    718     real ret = real::R_1;
     768    real ret = real::R_1();
    719769    int i = 1, multiplier = 1, exponent = 0;
    720770
     
    749799    int a = ceilf(logf(2) / logf(2 * M_PI) * real::BIGITS * real::BIGIT_BITS);
    750800
    751     real ret = sqrt(real::R_PI * 2);
    752     real fact_k_1 = real::R_1;
     801    real ret = sqrt(real::R_PI() * 2);
     802    real fact_k_1 = real::R_1();
    753803
    754804    for (int k = 1; k < a; k++)
     
    761811    }
    762812
    763     ret *= pow(x + (real)(a - 1), x - (real::R_1 / 2));
     813    ret *= pow(x + (real)(a - 1), x - (real::R_1() / 2));
    764814    ret *= exp(-x - (real)(a - 1));
    765815
     
    792842     * sqrt() call. */
    793843    real y = sqrt(x);
    794     real z = (y - real::R_1) / (y + real::R_1), z2 = z * z, zn = z2;
    795     real sum = real::R_1;
     844    real z = (y - real::R_1()) / (y + real::R_1()), z2 = z * z, zn = z2;
     845    real sum = real::R_1();
    796846
    797847    for (int i = 3; ; i += 2)
     
    819869    }
    820870    tmp.m_signexp = (1 << 30) - 1;
    821     return (real)(int)(x.m_signexp - (1 << 30) + 1) * real::R_LN2
     871    return (real)(int)(x.m_signexp - (1 << 30) + 1) * real::R_LN2()
    822872           + fast_log(tmp);
    823873}
     
    835885    tmp.m_signexp = (1 << 30) - 1;
    836886    return (real)(int)(x.m_signexp - (1 << 30) + 1)
    837            + fast_log(tmp) * real::R_LOG2E;
     887           + fast_log(tmp) * real::R_LOG2E();
    838888}
    839889
    840890template<> real log10(real const &x)
    841891{
    842     return log(x) * real::R_LOG10E;
     892    return log(x) * real::R_LOG10E();
    843893}
    844894
     
    849899     * domain of validity. The argument y is used for cases where we
    850900     * don't want the leading 1 in the Taylor series. */
    851     real ret = real::R_1 - y, xn = x;
     901    real ret = real::R_1() - y, xn = x;
    852902    int i = 1;
    853903
     
    882932     *  return x1 * 2^E0
    883933     */
    884     int e0 = x / real::R_LN2;
    885     real x0 = x - (real)e0 * real::R_LN2;
    886     real x1 = fast_exp_sub(x0, real::R_0);
     934    int e0 = x / real::R_LN2();
     935    real x0 = x - (real)e0 * real::R_LN2();
     936    real x1 = fast_exp_sub(x0, real::R_0());
    887937    x1.m_signexp += e0;
    888938    return x1;
     
    894944    int e0 = x;
    895945    real x0 = x - (real)e0;
    896     real x1 = fast_exp_sub(x0 * real::R_LN2, real::R_0);
     946    real x1 = fast_exp_sub(x0 * real::R_LN2(), real::R_0());
    897947    x1.m_signexp += e0;
    898948    return x1;
     
    904954     * accuracy near zero. We only use this identity for |x|>0.5. If
    905955     * |x|<=0.5, we compute exp(x)-1 and exp(-x)-1 instead. */
    906     bool near_zero = (fabs(x) < real::R_1 / 2);
    907     real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x);
    908     real x2 = near_zero ? fast_exp_sub(-x, real::R_1) : exp(-x);
     956    bool near_zero = (fabs(x) < real::R_1() / 2);
     957    real x1 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x);
     958    real x2 = near_zero ? fast_exp_sub(-x, real::R_1()) : exp(-x);
    909959    return (x1 - x2) / 2;
    910960}
     
    913963{
    914964    /* See sinh() for the strategy here */
    915     bool near_zero = (fabs(x) < real::R_1 / 2);
    916     real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x);
    917     real x2 = near_zero ? fast_exp_sub(-x, real::R_1) : exp(-x);
    918     real x3 = near_zero ? x1 + x2 + real::R_2 : x1 + x2;
     965    bool near_zero = (fabs(x) < real::R_1() / 2);
     966    real x1 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x);
     967    real x2 = near_zero ? fast_exp_sub(-x, real::R_1()) : exp(-x);
     968    real x3 = near_zero ? x1 + x2 + real::R_2() : x1 + x2;
    919969    return (x1 - x2) / x3;
    920970}
     
    9611011template<> real ulp(real const &x)
    9621012{
    963     real ret = real::R_1;
     1013    real ret = real::R_1();
    9641014    if (x)
    9651015        ret.m_signexp = x.m_signexp + 1 - real::BIGITS * real::BIGIT_BITS;
     
    9951045     *  - otherwise, if e is the exponent, clear all bits except the
    9961046     *    first e. */
    997     if (x < -real::R_0)
     1047    if (x < -real::R_0())
    9981048        return -ceil(-x);
    9991049    if (!x)
    10001050        return x;
    1001     if (x < real::R_1)
    1002         return real::R_0;
     1051    if (x < real::R_1())
     1052        return real::R_0();
    10031053
    10041054    real ret = x;
     
    10241074     *  - if x == floor(x), return x
    10251075     *  - otherwise, return floor(x) + 1 */
    1026     if (x < -real::R_0)
     1076    if (x < -real::R_0())
    10271077        return -floor(-x);
    10281078    real ret = floor(x);
     
    10301080        return ret;
    10311081    else
    1032         return ret + real::R_1;
     1082        return ret + real::R_1();
    10331083}
    10341084
    10351085template<> real round(real const &x)
    10361086{
    1037     if (x < real::R_0)
     1087    if (x < real::R_0())
    10381088        return -round(-x);
    10391089
    1040     return floor(x + (real::R_1 / 2));
     1090    return floor(x + (real::R_1() / 2));
    10411091}
    10421092
     
    10441094{
    10451095    if (!y)
    1046         return real::R_0; /* FIXME: return NaN */
     1096        return real::R_0(); /* FIXME: return NaN */
    10471097
    10481098    if (!x)
     
    10571107    int switch_sign = x.m_signexp & 0x80000000u;
    10581108
    1059     real absx = fmod(fabs(x), real::R_PI * 2);
    1060     if (absx > real::R_PI)
    1061     {
    1062         absx -= real::R_PI;
     1109    real absx = fmod(fabs(x), real::R_PI() * 2);
     1110    if (absx > real::R_PI())
     1111    {
     1112        absx -= real::R_PI();
    10631113        switch_sign = !switch_sign;
    10641114    }
    10651115
    1066     if (absx > real::R_PI_2)
    1067         absx = real::R_PI - absx;
    1068 
    1069     real ret = real::R_0, fact = real::R_1, xn = absx, mx2 = -absx * absx;
     1116    if (absx > real::R_PI_2())
     1117        absx = real::R_PI() - absx;
     1118
     1119    real ret = real::R_0(), fact = real::R_1(), xn = absx, mx2 = -absx * absx;
    10701120    int i = 1;
    10711121    for (;;)
     
    10881138template<> real cos(real const &x)
    10891139{
    1090     return sin(real::R_PI_2 - x);
     1140    return sin(real::R_PI_2() - x);
    10911141}
    10921142
     
    10941144{
    10951145    /* Constrain input to [-π,π] */
    1096     real y = fmod(x, real::R_PI);
     1146    real y = fmod(x, real::R_PI());
    10971147
    10981148    /* Constrain input to [-π/2,π/2] */
    1099     if (y < -real::R_PI_2)
    1100         y += real::R_PI;
    1101     else if (y > real::R_PI_2)
    1102         y -= real::R_PI;
     1149    if (y < -real::R_PI_2())
     1150        y += real::R_PI();
     1151    else if (y > real::R_PI_2())
     1152        y -= real::R_PI();
    11031153
    11041154    /* In [-π/4,π/4] return sin/cos */
    1105     if (fabs(y) <= real::R_PI_4)
     1155    if (fabs(y) <= real::R_PI_4())
    11061156        return sin(y) / cos(y);
    11071157
    11081158    /* Otherwise, return cos/sin */
    1109     if (y > real::R_0)
    1110         y = real::R_PI_2 - y;
     1159    if (y > real::R_0())
     1160        y = real::R_PI_2() - y;
    11111161    else
    1112         y = -real::R_PI_2 - y;
     1162        y = -real::R_PI_2() - y;
    11131163
    11141164    return cos(y) / sin(y);
     
    11231173     * lose the precision around x=1. */
    11241174    real absx = fabs(x);
    1125     int around_zero = (absx < (real::R_1 / 2));
     1175    int around_zero = (absx < (real::R_1() / 2));
    11261176
    11271177    if (!around_zero)
    1128         absx = sqrt((real::R_1 - absx) / 2);
     1178        absx = sqrt((real::R_1() - absx) / 2);
    11291179
    11301180    real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1;
     
    11451195
    11461196    if (around_zero)
    1147         ret = is_asin ? ret : real::R_PI_2 - ret;
     1197        ret = is_asin ? ret : real::R_PI_2() - ret;
    11481198    else
    11491199    {
    1150         real adjust = is_negative ? real::R_PI : real::R_0;
     1200        real adjust = is_negative ? real::R_PI() : real::R_0();
    11511201        if (is_asin)
    1152             ret = real::R_PI_2 - adjust - ret * 2;
     1202            ret = real::R_PI_2() - adjust - ret * 2;
    11531203        else
    11541204            ret = adjust + ret * 2;
     
    11941244    real absx = fabs(x);
    11951245
    1196     if (absx < (real::R_1 / 2))
     1246    if (absx < (real::R_1() / 2))
    11971247    {
    11981248        real ret = x, xn = x, mx2 = -x * x;
     
    12101260    real ret = 0;
    12111261
    1212     if (absx < (real::R_3 / 2))
    1213     {
    1214         real y = real::R_1 - absx;
     1262    if (absx < (real::R_3() / 2))
     1263    {
     1264        real y = real::R_1() - absx;
    12151265        real yn = y, my2 = -y * y;
    12161266        for (int i = 0; ; i += 2)
     
    12261276            yn *= my2;
    12271277        }
    1228         ret = real::R_PI_4 - ret;
    1229     }
    1230     else if (absx < real::R_2)
    1231     {
    1232         real y = (absx - real::R_SQRT3) / 2;
     1278        ret = real::R_PI_4() - ret;
     1279    }
     1280    else if (absx < real::R_2())
     1281    {
     1282        real y = (absx - real::R_SQRT3()) / 2;
    12331283        real yn = y, my2 = -y * y;
    12341284        for (int i = 1; ; i += 6)
     
    12361286            real newret = ret + ((yn / (real)i) / 2);
    12371287            yn *= y;
    1238             newret -= (real::R_SQRT3 / 2) * yn / (real)(i + 1);
     1288            newret -= (real::R_SQRT3() / 2) * yn / (real)(i + 1);
    12391289            yn *= y;
    12401290            newret += yn / (real)(i + 2);
    12411291            yn *= y;
    1242             newret -= (real::R_SQRT3 / 2) * yn / (real)(i + 3);
     1292            newret -= (real::R_SQRT3() / 2) * yn / (real)(i + 3);
    12431293            yn *= y;
    12441294            newret += (yn / (real)(i + 4)) / 2;
     
    12481298            yn *= my2;
    12491299        }
    1250         ret = real::R_PI_3 + ret;
     1300        ret = real::R_PI_3() + ret;
    12511301    }
    12521302    else
     
    12631313            ret = newret;
    12641314        }
    1265         ret = real::R_PI_2 - ret;
     1315        ret = real::R_PI_2() - ret;
    12661316    }
    12671317
     
    12781328            return y;
    12791329        if (y.m_signexp >> 31)
    1280             return -real::R_PI;
    1281         return real::R_PI;
     1330            return -real::R_PI();
     1331        return real::R_PI();
    12821332    }
    12831333
     
    12851335    {
    12861336        if (y.m_signexp >> 31)
    1287             return -real::R_PI;
    1288         return real::R_PI;
     1337            return -real::R_PI();
     1338        return real::R_PI();
    12891339    }
    12901340
     
    12921342    real z = y / x;
    12931343    real ret = atan(z);
    1294     if (x < real::R_0)
    1295         ret += (y > real::R_0) ? real::R_PI : -real::R_PI;
     1344    if (x < real::R_0())
     1345        ret += (y > real::R_0()) ? real::R_PI() : -real::R_PI();
    12961346    return ret;
    12971347}
     
    13341384    /* FIXME: better use int64_t when the cast is implemented */
    13351385    int exponent = ceil(log10(x));
    1336     x /= pow(R_10, (real)exponent);
    1337 
    1338     if (x < R_1)
    1339     {
    1340         x *= R_10;
     1386    x /= pow(R_10(), (real)exponent);
     1387
     1388    if (x < R_1())
     1389    {
     1390        x *= R_10();
    13411391        exponent--;
    13421392    }
     
    13501400            *str++ = '.';
    13511401        x -= real(digit);
    1352         x *= R_10;
     1402        x *= R_10();
    13531403    }
    13541404
     
    13791429}
    13801430
    1381 template<> real const real::R_0        = (real)0.0;
    1382 template<> real const real::R_1        = (real)1.0;
    1383 template<> real const real::R_2        = (real)2.0;
    1384 template<> real const real::R_3        = (real)3.0;
    1385 template<> real const real::R_10       = (real)10.0;
    1386 
    1387 /*
    1388  * Initialisation order is important here:
    1389  *  - fast_log() requires R_1
    1390  *  - log() requires R_LN2
    1391  *  - re() require R_2
    1392  *  - exp() requires R_0, R_1, R_LN2
    1393  *  - sqrt() requires R_3
    1394  */
    1395 template<> real const real::R_LN2      = fast_log(R_2);
    1396 template<> real const real::R_LN10     = log(R_10);
    1397 template<> real const real::R_LOG2E    = re(R_LN2);
    1398 template<> real const real::R_LOG10E   = re(R_LN10);
    1399 template<> real const real::R_E        = exp(R_1);
    1400 template<> real const real::R_PI       = fast_pi();
    1401 template<> real const real::R_PI_2     = R_PI / 2;
    1402 template<> real const real::R_PI_3     = R_PI / R_3;
    1403 template<> real const real::R_PI_4     = R_PI / 4;
    1404 template<> real const real::R_1_PI     = re(R_PI);
    1405 template<> real const real::R_2_PI     = R_1_PI * 2;
    1406 template<> real const real::R_2_SQRTPI = re(sqrt(R_PI)) * 2;
    1407 template<> real const real::R_SQRT2    = sqrt(R_2);
    1408 template<> real const real::R_SQRT3    = sqrt(R_3);
    1409 template<> real const real::R_SQRT1_2  = R_SQRT2 / 2;
    1410 
    14111431} /* namespace lol */
    14121432
Note: See TracChangeset for help on using the changeset viewer.