Changeset 1001


Ignore:
Timestamp:
Oct 3, 2011, 1:36:25 AM (10 years ago)
Author:
sam
Message:

core: implement asin() and acos() for real numbers and add unit tests for
these functions.

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r1000 r1001  
    9595
    9696    return u.d;
     97}
     98
     99real real::operator +() const
     100{
     101    return *this;
    97102}
    98103
     
    667672
    668673    return ret;
     674}
     675
     676static real asinacos(real const &x, bool is_asin, bool is_negative)
     677{
     678    /* Strategy for asin(): in [-0.5..0.5], use a Taylor series around
     679     * zero. In [0.5..1], use asin(x) = π/2 - 2*asin(sqrt((1-x)/2)), and
     680     * in [-1..-0.5] just revert the sign.
     681     * Strategy for acos(): use acos(x) = π/2 - asin(x) and try not to
     682     * lose the precision around x=1. */
     683    real absx = fabs(x);
     684    bool around_zero = (absx < (real::R_1 >> 1));
     685
     686    if (!around_zero)
     687        absx = sqrt((real::R_1 - absx) >> 1);
     688
     689    real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1;
     690    for (int i = 1; i < 280; i++)
     691    {
     692        xn *= x2;
     693        ret += (fact1 * xn / ((real)(2 * i + 1) * fact2)) >> (i * 2);
     694        fact1 *= (real)((2 * i + 1) * (2 * i + 2));
     695        fact2 *= (real)((i + 1) * (i + 1));
     696    }
     697
     698    if (is_negative)
     699        ret = -ret;
     700
     701    if (around_zero)
     702        ret = is_asin ? ret : real::R_PI_2 - ret;
     703    else
     704    {
     705        real adjust = is_negative ? real::R_PI : real::R_0;
     706        if (is_asin)
     707            ret = real::R_PI_2 - adjust - (ret << 1);
     708        else
     709            ret = adjust + (ret << 1);
     710    }
     711
     712    return ret;
     713}
     714
     715real asin(real const &x)
     716{
     717    return asinacos(x, true, x.m_signexp >> 31);
     718}
     719
     720real acos(real const &x)
     721{
     722    return asinacos(x, false, x.m_signexp >> 31);
    669723}
    670724
  • trunk/src/real.h

    r999 r1001  
    3737    operator unsigned int() const;
    3838
     39    real operator +() const;
    3940    real operator -() const;
    4041    real operator +(real const &x) const;
     
    7172    friend real sin(real const &x);
    7273    friend real cos(real const &x);
     74    friend real asin(real const &x);
     75    friend real acos(real const &x);
    7376    friend real atan(real const &x);
    7477
  • trunk/test/unit/real.cpp

    r998 r1001  
    237237        LOLUNIT_ASSERT(!!a);
    238238    }
     239
     240    LOLUNIT_TEST(AsinAcos)
     241    {
     242        double tests[14] =
     243        {
     244            -1024.0, -1023.0, -513.0, -512.0, -511.0, -1.0, -0.0,
     245            0.0, 1.0, 511.0, 512.0, 513.0, 1023.0, 1024.0
     246        };
     247
     248        for (int n = 0; n < 14; n++)
     249        {
     250            double a = tests[n] / 1024;
     251            double b = sin(asin((real)a));
     252            double c = cos(acos((real)a));
     253
     254            LOLUNIT_SET_CONTEXT(a);
     255            LOLUNIT_ASSERT_DOUBLES_EQUAL(b, a, 1e-100);
     256            LOLUNIT_ASSERT_DOUBLES_EQUAL(c, a, 1e-100);
     257        }
     258    }
    239259};
    240260
Note: See TracChangeset for help on using the changeset viewer.