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

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

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.