Changeset 1013


Ignore:
Timestamp:
Oct 10, 2011, 3:33:09 AM (9 years ago)
Author:
sam
Message:

core: remove most dependencies on real number size in the various math
functions.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r1006 r1013  
    488488    ret.m_signexp |= (exponent - 1) & 0x7fffffffu;
    489489
    490     /* Five steps of Newton-Raphson seems enough for 32-bigit reals. */
    491     real two = 2;
    492     ret = ret * (two - ret * x);
    493     ret = ret * (two - ret * x);
    494     ret = ret * (two - ret * x);
    495     ret = ret * (two - ret * x);
    496     ret = ret * (two - ret * x);
     490    /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for
     491     * convergence, but this hasn't been checked seriously. */
     492    for (int i = 1; i < real::BIGITS; i *= 2)
     493        ret = ret * (real::R_2 - ret * x);
    497494
    498495    return ret;
     
    533530
    534531    uint32_t exponent = (x.m_signexp & 0x7fffffffu);
    535     exponent = ((1 << 30) + (1 << 29) -1) - (exponent + 1) / 2;
     532    exponent = ((1 << 30) + (1 << 29) - 1) - (exponent + 1) / 2;
    536533    exponent = exponent + (v.x >> 23) - (u.x >> 23);
    537534    ret.m_signexp |= exponent & 0x7fffffffu;
    538535
    539     /* Five steps of Newton-Raphson seems enough for 32-bigit reals. */
    540     real three = 3;
    541     ret = ret * (three - ret * ret * x);
    542     ret.m_signexp--;
    543     ret = ret * (three - ret * ret * x);
    544     ret.m_signexp--;
    545     ret = ret * (three - ret * ret * x);
    546     ret.m_signexp--;
    547     ret = ret * (three - ret * ret * x);
    548     ret.m_signexp--;
    549     ret = ret * (three - ret * ret * x);
    550     ret.m_signexp--;
     536    /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for
     537     * convergence, but this hasn't been checked seriously. */
     538    for (int i = 1; i < real::BIGITS; i *= 2)
     539    {
     540        ret = ret * (real::R_3 - ret * ret * x);
     541        ret.m_signexp--;
     542    }
    551543
    552544    return ret * x;
     
    560552}
    561553
    562 static real fastlog(real const &x)
     554static real fast_log(real const &x)
    563555{
    564556    /* This fast log method is tuned to work on the [1..2] range and
     
    578570     * sqrt() call. */
    579571    real y = sqrt(x);
    580     real z = (y - (real)1) / (y + (real)1), z2 = z * z, zn = z2;
    581     real sum = 1.0;
    582 
    583     for (int i = 3; i < 200; i += 2)
    584     {
    585         sum += zn / (real)i;
     572    real z = (y - real::R_1) / (y + real::R_1), z2 = z * z, zn = z2;
     573    real sum = real::R_1;
     574
     575    for (int i = 3; ; i += 2)
     576    {
     577        real newsum = sum + zn / (real)i;
     578        if (newsum == sum)
     579            break;
     580        sum = newsum;
    586581        zn *= z2;
    587582    }
     
    590585}
    591586
    592 static real LOG_2 = fastlog((real)2);
    593 
    594587real log(real const &x)
    595588{
    596589    /* Strategy for log(x): if x = 2^E*M then log(x) = E log(2) + log(M),
    597      * with the property that M is in [1..2[, so fastlog() applies here. */
     590     * with the property that M is in [1..2[, so fast_log() applies here. */
    598591    real tmp = x;
    599592    if (x.m_signexp >> 31 || x.m_signexp == 0)
     
    604597    }
    605598    tmp.m_signexp = (1 << 30) - 1;
    606     return (real)(x.m_signexp - (1 << 30) + 1) * LOG_2 + fastlog(tmp);
     599    return (real)(x.m_signexp - (1 << 30) + 1) * real::R_LN2 + fast_log(tmp);
    607600}
    608601
     
    625618     *  return x1 * 2^E0
    626619     */
    627     int e0 = x / LOG_2;
    628     real x0 = x - (real)e0 * LOG_2;
    629     real x1 = 1.0, fact = 1.0, xn = x0;
    630 
    631     for (int i = 1; i < 100; i++)
     620    int e0 = x / real::R_LN2;
     621    real x0 = x - (real)e0 * real::R_LN2;
     622    real x1 = real::R_1, fact = real::R_1, xn = x0;
     623
     624    for (int i = 1; ; i++)
    632625    {
    633626        fact *= (real)i;
    634         x1 += xn / fact;
     627        real newx1 = x1 + xn / fact;
     628        if (newx1 == x1)
     629            break;
     630        x1 = newx1;
    635631        xn *= x0;
    636632    }
     
    720716        absx = real::R_PI - absx;
    721717
    722     real ret = 0.0, fact = 1.0, xn = absx, x2 = absx * absx;
     718    real ret = real::R_0, fact = real::R_1, xn = absx, x2 = absx * absx;
    723719    for (int i = 1; ; i += 2)
    724720    {
    725721        real newret = ret + xn / fact;
    726         if (ret == newret)
     722        if (newret == ret)
    727723            break;
    728724        ret = newret;
     
    756752
    757753    real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1;
    758     for (int i = 1; i < 280; i++)
     754    for (int i = 1; ; i++)
    759755    {
    760756        xn *= x2;
    761         ret += (fact1 * xn / ((real)(2 * i + 1) * fact2)) >> (i * 2);
     757        real mul = (real)(2 * i + 1);
     758        real newret = ret + ((fact1 * xn / (mul * fact2)) >> (i * 2));
     759        if (newret == ret)
     760            break;
     761        ret = newret;
    762762        fact1 *= (real)((2 * i + 1) * (2 * i + 2));
    763763        fact2 *= (real)((i + 1) * (i + 1));
     
    820820    {
    821821        real ret = x, xn = x, mx2 = -x * x;
    822         for (int i = 3; i < 100; i += 2)
     822        for (int i = 3; ; i += 2)
    823823        {
    824824            xn *= mx2;
    825             ret += xn / (real)i;
     825            real newret = ret + xn / (real)i;
     826            if (newret == ret)
     827                break;
     828            ret = newret;
    826829        }
    827830        return ret;
     
    834837        real y = real::R_1 - absx;
    835838        real yn = y, my2 = -y * y;
    836         for (int i = 0; i < 200; i += 2)
     839        for (int i = 0; ; i += 2)
    837840        {
    838             ret += (yn / (real)(2 * i + 1)) >> (i + 1);
     841            real newret = ret + ((yn / (real)(2 * i + 1)) >> (i + 1));
    839842            yn *= y;
    840             ret += (yn / (real)(2 * i + 2)) >> (i + 1);
     843            newret += (yn / (real)(2 * i + 2)) >> (i + 1);
    841844            yn *= y;
    842             ret += (yn / (real)(2 * i + 3)) >> (i + 2);
     845            newret += (yn / (real)(2 * i + 3)) >> (i + 2);
     846            if (newret == ret)
     847                break;
     848            ret = newret;
    843849            yn *= my2;
    844850        }
     
    849855        real y = (absx - real::R_SQRT3) >> 1;
    850856        real yn = y, my2 = -y * y;
    851         for (int i = 1; i < 200; i += 6)
     857        for (int i = 1; ; i += 6)
    852858        {
    853             ret += (yn / (real)i) >> 1;
     859            real newret = ret + ((yn / (real)i) >> 1);
    854860            yn *= y;
    855             ret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 1);
     861            newret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 1);
    856862            yn *= y;
    857             ret += yn / (real)(i + 2);
     863            newret += yn / (real)(i + 2);
    858864            yn *= y;
    859             ret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 3);
     865            newret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 3);
    860866            yn *= y;
    861             ret += (yn / (real)(i + 4)) >> 1;
     867            newret += (yn / (real)(i + 4)) >> 1;
     868            if (newret == ret)
     869                break;
     870            ret = newret;
    862871            yn *= my2;
    863872        }
     
    869878        real yn = y, my2 = -y * y;
    870879        ret = y;
    871         for (int i = 3; i < 120; i += 2)
     880        for (int i = 3; ; i += 2)
    872881        {
    873882            yn *= my2;
    874             ret += yn / (real)i;
     883            real newret = ret + yn / (real)i;
     884            if (newret == ret)
     885                break;
     886            ret = newret;
    875887        }
    876888        ret = real::R_PI_2 - ret;
     
    945957    real const m0 = -x0 * x0, m1 = -x1 * x1, r16 = 16.0, r4 = 4.0;
    946958
    947     /* Degree 240 is required for 512-bit mantissa precision */
    948     for (int i = 1; i < 240; i += 2)
    949     {
    950         ret += r16 / (x0 * (real)i) - r4 / (x1 * (real)i);
     959    for (int i = 1; ; i += 2)
     960    {
     961        real newret = ret + r16 / (x0 * (real)i) - r4 / (x1 * (real)i);
     962        if (newret == ret)
     963            break;
     964        ret = newret;
    951965        x0 *= m0;
    952966        x1 *= m1;
     
    962976real const real::R_10       = (real)10.0;
    963977
    964 real const real::R_E        = exp(R_1);
    965 real const real::R_LN2      = log(R_2);
     978real const real::R_LN2      = fast_log(R_2);
    966979real const real::R_LN10     = log(R_10);
    967980real const real::R_LOG2E    = re(R_LN2);
    968981real const real::R_LOG10E   = re(R_LN10);
     982real const real::R_E        = exp(R_1);
    969983real const real::R_PI       = fast_pi();
    970984real const real::R_PI_2     = R_PI >> 1;
Note: See TracChangeset for help on using the changeset viewer.