Changeset 1019


Ignore:
Timestamp:
Oct 12, 2011, 12:01:51 AM (9 years ago)
Author:
sam
Message:

core: implement real methods cbrt(), log2(), exp2(), and copysign().

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r1018 r1019  
    552552
    553553    return ret * x;
     554}
     555
     556real cbrt(real const &x)
     557{
     558    /* if zero, return x */
     559    if (!(x.m_signexp << 1))
     560        return x;
     561
     562    /* Use the system's float inversion to approximate cbrt(x). First
     563     * we construct a float in the [1..8[ range that has roughly the same
     564     * mantissa as our real. Its exponent is 0, 1 or 2, depending on the
     565     * value of x. The final exponent is 0 or 1 (special case). We use
     566     * the final exponent and final mantissa to pre-fill the result. */
     567    union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f };
     568    v.x += ((x.m_signexp % 3) << 23);
     569    v.x |= x.m_mantissa[0] >> 9;
     570    v.f = powf(v.f, 0.33333333333333333f);
     571
     572    real ret;
     573    ret.m_mantissa[0] = v.x << 9;
     574
     575    uint32_t sign = x.m_signexp & 0x80000000u;
     576    ret.m_signexp = sign;
     577
     578    int exponent = (x.m_signexp & 0x7fffffffu) - (1 << 30) + 1;
     579    exponent = exponent / 3 + (v.x >> 23) - (u.x >> 23);
     580    ret.m_signexp |= (exponent + (1 << 30) - 1) & 0x7fffffffu;
     581
     582    /* FIXME: 1+log2(BIGITS) steps of Newton-Raphson seems to be enough for
     583     * convergence, but this hasn't been checked seriously. */
     584    for (int i = 1; i <= real::BIGITS; i *= 2)
     585    {
     586        static real third = re(real::R_3);
     587        ret = third * (x / (ret * ret) + (ret << 1));
     588    }
     589
     590    return ret;
    554591}
    555592
     
    609646}
    610647
     648real log2(real const &x)
     649{
     650    /* Strategy for log2(x): see log(x). */
     651    real tmp = x;
     652    if (x.m_signexp >> 31 || x.m_signexp == 0)
     653    {
     654        tmp.m_signexp = 0xffffffffu;
     655        tmp.m_mantissa[0] = 0xffffffffu;
     656        return tmp;
     657    }
     658    tmp.m_signexp = (1 << 30) - 1;
     659    return (real)(x.m_signexp - (1 << 30) + 1) + fast_log(tmp) * real::R_LOG2E;
     660}
     661
     662static real fast_exp(real const &x)
     663{
     664    /* This fast exp method is tuned to work on the [0..1] range and
     665     * no effort whatsoever was made to improve convergence outside this
     666     * domain of validity. */
     667    real ret = real::R_1, fact = real::R_1, xn = x;
     668
     669    for (int i = 1; ; i++)
     670    {
     671        fact *= (real)i;
     672        real newret = ret + xn / fact;
     673        if (newret == ret)
     674            break;
     675        ret = newret;
     676        xn *= x;
     677    }
     678
     679    return ret;
     680}
     681
    611682real exp(real const &x)
    612683{
     
    629700    int e0 = x / real::R_LN2;
    630701    real x0 = x - (real)e0 * real::R_LN2;
    631     real x1 = real::R_1, fact = real::R_1, xn = x0;
    632 
    633     for (int i = 1; ; i++)
    634     {
    635         fact *= (real)i;
    636         real newx1 = x1 + xn / fact;
    637         if (newx1 == x1)
    638             break;
    639         x1 = newx1;
    640         xn *= x0;
    641     }
    642 
     702    real x1 = fast_exp(x0);
    643703    x1.m_signexp += e0;
    644704    return x1;
     705}
     706
     707real exp2(real const &x)
     708{
     709    /* Strategy for exp2(x): see strategy in exp(). */
     710    int e0 = x;
     711    real x0 = x - (real)e0;
     712    real x1 = fast_exp(x0 * real::R_LN2);
     713    x1.m_signexp += e0;
     714    return x1;
     715}
     716
     717real copysign(real const &x, real const &y)
     718{
     719    real ret = x;
     720    ret.m_signexp &= 0x7fffffffu;
     721    ret.m_signexp |= y.m_signexp & 0x80000000u;
     722    return ret;
    645723}
    646724
     
    9291007void real::print(int ndigits) const
    9301008{
    931     real const r1 = 1, r10 = 10;
    9321009    real x = *this;
    9331010
     
    9421019    if (x.m_signexp)
    9431020    {
    944         for (real div = r1, newdiv; true; div = newdiv)
     1021        for (real div = R_1, newdiv; true; div = newdiv)
    9451022        {
    946             newdiv = div * r10;
     1023            newdiv = div * R_10;
    9471024            if (x < newdiv)
    9481025            {
     
    9521029            exponent++;
    9531030        }
    954         for (real mul = 1, newx; true; mul *= r10)
     1031        for (real mul = 1, newx; true; mul *= R_10)
    9551032        {
    9561033            newx = x * mul;
    957             if (newx >= r1)
     1034            if (newx >= R_1)
    9581035            {
    9591036                x = newx;
     
    9721049            printf(".");
    9731050        x -= real(digit);
    974         x *= r10;
     1051        x *= R_10;
    9751052    }
    9761053
  • trunk/src/real.h

    r1018 r1019  
    6363    operator bool() const;
    6464
    65     friend real fabs(real const &x);
    66 
    67     friend real re(real const &x);
    68     friend real sqrt(real const &x);
    69     friend real log(real const &x);
    70     friend real exp(real const &x);
    71 
    72     friend real floor(real const &x);
    73     friend real ceil(real const &x);
    74     friend real round(real const &x);
    75     friend real fmod(real const &x, real const &y);
    76 
     65    /* Trigonometric functions */
    7766    friend real sin(real const &x);
    7867    friend real cos(real const &x);
     
    8170    friend real acos(real const &x);
    8271    friend real atan(real const &x);
     72    /* FIXME: atan2 */
     73
     74    /* Hyperbolic functions */
     75    /* FIXME: sinh cosh tanh */
     76
     77    /* Exponential and logarithmic functions */
     78    friend real exp(real const &x);
     79    friend real exp2(real const &x);
     80    friend real log(real const &x);
     81    friend real log2(real const &x);
     82    /* FIXME: frexp ldexp log10 modf */
     83
     84    /* Power functions */
     85    friend real re(real const &x);
     86    friend real sqrt(real const &x);
     87    friend real cbrt(real const &x);
     88    /* FIXME: pow */
     89
     90    /* Rounding, absolute value, remainder etc. */
     91    friend real ceil(real const &x);
     92    friend real copysign(real const &x, real const &y);
     93    friend real floor(real const &x);
     94    friend real fabs(real const &x);
     95    friend real round(real const &x);
     96    friend real fmod(real const &x, real const &y);
    8397
    8498    void print(int ndigits = 150) const;
  • trunk/test/unit/real.cpp

    r1017 r1019  
    3535        LOLUNIT_ASSERT_EQUAL(a10, 10.0);
    3636
    37         double b0 = log(real::R_E);
    38         LOLUNIT_ASSERT_EQUAL(b0, 1.0);
    39 
    40         double b1 = exp(re(real::R_LOG2E));
    41         LOLUNIT_ASSERT_EQUAL(b1, 2.0);
    42 
    43         double b2 = exp(re(real::R_LOG10E));
    44         LOLUNIT_ASSERT_EQUAL(b2, 10.0);
    45 
    46         double b3 = exp(real::R_LN2);
    47         LOLUNIT_ASSERT_EQUAL(b3, 2.0);
    48 
    49         double b4 = exp(real::R_LN10);
    50         LOLUNIT_ASSERT_EQUAL(b4, 10.0);
    51 
    52         double b5 = sin(real::R_PI);
    53         double b6 = cos(real::R_PI);
    54         LOLUNIT_ASSERT_DOUBLES_EQUAL(b5, 0.0, 1e-100);
    55         LOLUNIT_ASSERT_EQUAL(b6, -1.0);
    56 
    57         double b7 = sin(real::R_PI_2);
    58         double b8 = cos(real::R_PI_2);
    59         LOLUNIT_ASSERT_EQUAL(b7, 1.0);
    60         LOLUNIT_ASSERT_DOUBLES_EQUAL(b8, 0.0, 1e-100);
    61 
    62         double b9 = sin(real::R_PI_4) * sin(real::R_PI_4);
    63         double b10 = cos(real::R_PI_4) * cos(real::R_PI_4);
    64         LOLUNIT_ASSERT_EQUAL(b9, 0.5);
    65         LOLUNIT_ASSERT_EQUAL(b10, 0.5);
     37        double b1 = log(real::R_E);
     38        double b2 = log2(real::R_2);
     39        LOLUNIT_ASSERT_EQUAL(b1, 1.0);
     40        LOLUNIT_ASSERT_EQUAL(b2, 1.0);
     41
     42        double c1 = exp(re(real::R_LOG2E));
     43        double c2 = log(exp2(real::R_LOG2E));
     44        LOLUNIT_ASSERT_EQUAL(c1, 2.0);
     45        LOLUNIT_ASSERT_EQUAL(c2, 1.0);
     46
     47        double d1 = exp(re(real::R_LOG10E));
     48        LOLUNIT_ASSERT_EQUAL(d1, 10.0);
     49
     50        double e1 = exp(real::R_LN2);
     51        LOLUNIT_ASSERT_EQUAL(e1, 2.0);
     52
     53        double f1 = exp(real::R_LN10);
     54        LOLUNIT_ASSERT_EQUAL(f1, 10.0);
     55
     56        double g1 = sin(real::R_PI);
     57        double g2 = cos(real::R_PI);
     58        LOLUNIT_ASSERT_DOUBLES_EQUAL(g1, 0.0, 1e-100);
     59        LOLUNIT_ASSERT_EQUAL(g2, -1.0);
     60
     61        double h1 = sin(real::R_PI_2);
     62        double h2 = cos(real::R_PI_2);
     63        LOLUNIT_ASSERT_EQUAL(h1, 1.0);
     64        LOLUNIT_ASSERT_DOUBLES_EQUAL(h2, 0.0, 1e-100);
     65
     66        double i1 = sin(real::R_PI_4) * sin(real::R_PI_4);
     67        double i2 = cos(real::R_PI_4) * cos(real::R_PI_4);
     68        LOLUNIT_ASSERT_EQUAL(i1, 0.5);
     69        LOLUNIT_ASSERT_EQUAL(i2, 0.5);
    6670    }
    6771
Note: See TracChangeset for help on using the changeset viewer.