Changeset 999


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

core: implement atan() for real numbers.

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r998 r999  
    668668}
    669669
     670real atan(real const &x)
     671{
     672    /* Computing atan(x): we choose a different Taylor series depending on
     673     * the value of x to help with convergence.
     674     *
     675     * If |x| < 0.5 we evaluate atan(y) near 0:
     676     *  atan(y) = y - y^3/3 + y^5/5 - y^7/7 + y^9/9 ...
     677     *
     678     * If 0.5 <= |x| < 1.5 we evaluate atan(1+y) near 0:
     679     *  atan(1+y) = π/4 + y/(1*2^1) - y^2/(2*2^1) + y^3/(3*2^2)
     680     *                  - y^5/(5*2^3) + y^6/(6*2^3) - y^7/(7*2^4)
     681     *                  + y^9/(9*2^5) - y^10/(10*2^5) + y^11/(11*2^6) ...
     682     *
     683     * If 1.5 <= |x| < 2 we evaluate atan(sqrt(3)+2y) near 0:
     684     *  atan(sqrt(3)+2y) = π/3 + 1/2 y - sqrt(3)/2 y^2/2
     685     *                         + y^3/3 - sqrt(3)/2 y^4/4 + 1/2 y^5/5
     686     *                         - 1/2 y^7/7 + sqrt(3)/2 y^8/8
     687     *                         - y^9/9 + sqrt(3)/2 y^10/10 - 1/2 y^11/11
     688     *                         + 1/2 y^13/13 - sqrt(3)/2 y^14/14
     689     *                         + y^15/15 - sqrt(3)/2 y^16/16 + 1/2 y^17/17 ...
     690     *
     691     * If |x| >= 2 we evaluate atan(y) near +∞:
     692     *  atan(y) = π/2 - y^-1 + y^-3/3 - y^-5/5 + y^-7/7 - y^-9/9 ...
     693     */
     694    real absx = fabs(x);
     695
     696    if (absx < (real::R_1 >> 1))
     697    {
     698        real ret = x, xn = x, mx2 = -x * x;
     699        for (int i = 3; i < 100; i += 2)
     700        {
     701            xn *= mx2;
     702            ret += xn / (real)i;
     703        }
     704        return ret;
     705    }
     706
     707    real ret = 0;
     708
     709    if (absx < (real::R_3 >> 1))
     710    {
     711        real y = real::R_1 - absx;
     712        real yn = y, my2 = -y * y;
     713        for (int i = 0; i < 200; i += 2)
     714        {
     715            ret += (yn / (real)(2 * i + 1)) >> (i + 1);
     716            yn *= y;
     717            ret += (yn / (real)(2 * i + 2)) >> (i + 1);
     718            yn *= y;
     719            ret += (yn / (real)(2 * i + 3)) >> (i + 2);
     720            yn *= my2;
     721        }
     722        ret = real::R_PI_4 - ret;
     723    }
     724    else if (absx < real::R_2)
     725    {
     726        real y = (absx - real::R_SQRT3) >> 1;
     727        real yn = y, my2 = -y * y;
     728        for (int i = 1; i < 200; i += 6)
     729        {
     730            ret += (yn / (real)i) >> 1;
     731            yn *= y;
     732            ret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 1);
     733            yn *= y;
     734            ret += yn / (real)(i + 2);
     735            yn *= y;
     736            ret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 3);
     737            yn *= y;
     738            ret += (yn / (real)(i + 4)) >> 1;
     739            yn *= my2;
     740        }
     741        ret = real::R_PI_3 + ret;
     742    }
     743    else
     744    {
     745        real y = re(absx);
     746        real yn = y, my2 = -y * y;
     747        ret = y;
     748        for (int i = 3; i < 120; i += 2)
     749        {
     750            yn *= my2;
     751            ret += yn / (real)i;
     752        }
     753        ret = real::R_PI_2 - ret;
     754    }
     755
     756    /* Propagate sign */
     757    ret.m_signexp |= (x.m_signexp & 0x80000000u);
     758    return ret;
     759}
     760
    670761void real::print(int ndigits) const
    671762{
     
    745836real const real::R_1        = (real)1.0;
    746837real const real::R_2        = (real)2.0;
     838real const real::R_3        = (real)3.0;
    747839real const real::R_10       = (real)10.0;
    748840
     
    754846real const real::R_PI       = fast_pi();
    755847real const real::R_PI_2     = R_PI >> 1;
     848real const real::R_PI_3     = R_PI / R_3;
    756849real const real::R_PI_4     = R_PI >> 2;
    757850real const real::R_1_PI     = re(R_PI);
     
    759852real const real::R_2_SQRTPI = re(sqrt(R_PI)) << 1;
    760853real const real::R_SQRT2    = sqrt(R_2);
     854real const real::R_SQRT3    = sqrt(R_3);
    761855real const real::R_SQRT1_2  = R_SQRT2 >> 1;
    762856
  • trunk/src/real.h

    r998 r999  
    7171    friend real sin(real const &x);
    7272    friend real cos(real const &x);
     73    friend real atan(real const &x);
    7374
    7475    void print(int ndigits = 150) const;
     
    7778    static real const R_1;
    7879    static real const R_2;
     80    static real const R_3;
    7981    static real const R_10;
    8082
     
    8688    static real const R_PI;
    8789    static real const R_PI_2;
     90    static real const R_PI_3;
    8891    static real const R_PI_4;
    8992    static real const R_1_PI;
     
    9194    static real const R_2_SQRTPI;
    9295    static real const R_SQRT2;
     96    static real const R_SQRT3;
    9397    static real const R_SQRT1_2;
    9498
Note: See TracChangeset for help on using the changeset viewer.