apfunc.cpp

任意精度计算的实现

#include <cmath>
#include "ap.h"


using namespace std;


#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif


// Overloaded apfloat functions


// Constants needed for log & exp
apfloat Readypi;
apfloat PiA, PiB, PiJ, Pione, Pitwo, Pifive, Pisix;
apfloat ReadyT;
apfloat ReadyQ;
apfloat ReadyP;
apfloat Readyinvroot;
size_t  Readypiterms;
rawtype ReadypiBase;
apfloat Logbase;
rawtype LogbaseBase;

inline size_t min (size_t a, size_t b)
{
    return (a < b ? a : b);
}

inline size_t max (size_t a, size_t b)
{
    return (a > b ? a : b);
}


// Positive integer power
// Algorithm improvements by Bernd Kellner
apfloat pow (apfloat base, unsigned long exp)
{
    apfloat r;
    int b2pow = 0;

    if (!exp) return 1;

    while (!(exp & 1))
    {
        b2pow++;
        exp >>= 1;
    }
    r = base;

    while (exp >>= 1)
    {
        base *= base;
        if (exp & 1) r *= base;
    }

    while (b2pow--)
        r *= r;

    return r;
}

apfloat pow (apfloat base, unsigned exp)
{
    return pow (base, (unsigned long) exp);
}

// Integer power
apfloat pow (apfloat base, long exp)
{
    if (exp < 0)
        return invroot (pow (base, (unsigned long) -exp), 1);
    else
        return pow (base, (unsigned long) exp);
}

apfloat pow (apfloat base, int exp)
{
    if (exp < 0)
        return invroot (pow (base, (unsigned long) -exp), 1);
    else
        return pow (base, (unsigned long) exp);
}

// Square root
apfloat sqrt (apfloat x)
{
    if (!x.sign ())
        return apfloat (new apstruct);          // Avoid division by zero
    else
        return x * invroot (x, 2);
}

// Cube root
apfloat cbrt (apfloat x)
{
    if (!x.sign ())
        return apfloat (new apstruct);          // Avoid division by zero
    else
        return x * invroot (x * x, 3);
}

// Positive integer root
apfloat root (apfloat x, unsigned n)
{
    if (!x.sign ())
        return apfloat (new apstruct);          // Avoid division by zero
    else
        return 1 / invroot (x, n);
}

// Inverse positive integer root
apfloat invroot (apfloat u, unsigned n, size_t destprec, apfloat initguess, size_t initprec)
{
    size_t prec, minprec, maxprec, doubledigits, fprec;
    int k, f;
    long p, r;
    apfloat one = 1, d = n, x;

    if (!n) return 1;

    assert (u.sign ());                         // Infinity

    if (!(n & 1))
        assert (u.sign () > 0);                 // Complex result

    if (destprec == (size_t) DEFAULT) destprec = u.prec ();

    // Approximate accuracy of a double
    doubledigits = (size_t) (log (MAX_PRECISE_DOUBLE) / log ((double) Basedigit));

    // Initial guess accuracy
    fprec = max (doubledigits, 2 * Basedigits);

    if (!initguess.ap)
    {
        p = u.ap->exp / (long) n;
        r = u.ap->exp - p * (long) n;

        x = u;
        x.prec (fprec);
        x.exp (Basedigits * r);         // Allow exponents in exess of doubles'

        // No initial guess, calculate it from u
        x = apfloat (pow (ap2double (x.ap), -1.0 / (double) n), fprec);
        x.exp (x.exp () - p * Basedigits);

        prec = min (doubledigits, Basedigits);

        // Check if a factor of 3 should be used in length
        maxprec = rnd23up (destprec / Basedigits);
        if (maxprec != (maxprec & -maxprec))
            minprec = 3 * Basedigits;
        else
            minprec = Basedigits;

        // Highly ineffective unless precision is 2^n * Basedigits (or 3*2^n)
        if (prec < minprec)
        {
            x.prec (minprec + 3 * Basedigits);
            while (prec < minprec)
            {
                x += x * (one - u * pow (x, n)) / d;
                prec *= 2;
            }
            prec = minprec;
        }
    }
    else
    {
        // Take initial guess
        x = initguess;

        if (initprec == (size_t) DEFAULT)
            prec = initguess.prec ();
        else
            prec = initprec;
    }

    // Check where the precising iteration should be done
    for (k = 0, maxprec = prec; maxprec < destprec; k++, maxprec <<= 1);
    for (f = k, minprec = prec; f; f--, minprec <<= 1)
        if (minprec >= 2 * Basedigits && (minprec - 2 * Basedigits) << f >= destprec)
            break;

    // Newton's iteration
    while (k--)
    {
        apfloat t;

        prec *= 2;
        x.prec (min (prec, destprec));

        t = one - u * pow (x, n);
        if (k < f) t.prec (prec / 2);

        if (n > 1)
            x += x * t / d;
        else
            x += x * t;

        // Precising iteration
        if (k == f)
        {
            if (n > 1)
                x += x * (one - u * pow (x, n)) / d;
            else
                x += x * (one - u * pow (x, n));
        }
    }

    x.prec (destprec);

    return x;
}

// Floor function
apfloat floor (apfloat x)
{
    if (x.sign () >= 0)
        return apfloat (apabsfloor (x.ap));
    else
        return apfloat (apabsceil (x.ap));
}

// Ceiling function
apfloat ceil (apfloat x)
{
    if (x.sign () >= 0)
        return apfloat (apabsceil (x.ap));
    else
        return apfloat (apabsfloor (x.ap));
}

// Splits x to integer and fractional parts
apfloat modf (apfloat x, apfloat *ipart)
{
    *ipart = apfloat (apabsfloor (x.ap));

    return x - *ipart;
}

// Returns x modulo y, where 0 <= result < y
apfloat fmod (apfloat x, apfloat y)
{
    size_t s, destprec;
    apfloat t, a, b, tx, ty;

    if (!y.sign ())
        return 0;                   // By definition

    a = abs (x);
    b = abs (y);

    if (a < b)
        return x;                   // abs (x) < abs (y)
    else if (x.prec () <= x.exp () - y.exp ())
        return 0;                   // Degenerate case; not enough precision to make any sense
    else
        s = x.exp () - y.exp () + 3 * Basedigits;   // Some extra precision

    tx = x;
    ty = y;

    tx.prec (s);
    ty.prec (s);

    t = tx / ty;                    // Approximate division
    t = apfloat (apabsfloor (t.ap));

    tx = x;
    ty = y;

    destprec = min (y.prec () == INFINITE ? INFINITE : y.prec () + x.exp () - y.exp (), x.prec ());

    tx.prec (destprec);
    ty.prec (destprec);

    t = tx - t * ty;

    a = abs (t);

    if (a >= b)                     // Fix division round-off error
        t = x.sign () * (a - b);

    return t;
}

// Absolute value
apfloat abs (apfloat x)
{
    if (x.sign () >= 0)
        return x;
    else
        return -x;
}

// Calculates pi to destprec digits using the Chudnovskys' binary splitting algorithm
// Uses previously calculated terms (if such exist) to improve the precision of the calculation
apfloat pi_a (size_t n)
{
    apfloat v = PiA + PiB * n;

    v.sign (1 - 2 * (n & 1));

    return v;
}

apfloat pi_p (size_t n)
{
    apfloat f = n, sixf = Pisix * f;

    if (!n)
        return Pione;
    else
        return apfloat (sixf - Pione) * (Pitwo * f - Pione) * (sixf - Pifive);
}

apfloat pi_q (size_t n)
{
    apfloat f = n;

    if (!n)
        return Pione;
    else
        return PiJ * f * f * f;
}

void pi_r (size_t N1, size_t N2, apfloat *T, apfloat *Q, apfloat *P)
{
    switch (N2 - N1)
    {
        case 0:
        {
            assert (N1 != N2);

            break;
        }
        case 1:
        {
            apfloat p0 = pi_p (N1);

            *T = pi_a (N1) * p0;
            *Q = pi_q (N1);
            *P = p0;

            break;
        }
        case 2:
        {
            apfloat p0 = pi_p (N1), p01 = p0 * pi_p (N1 + 1),
                    q1 = pi_q (N1 + 1);

            *T = q1 * pi_a (N1) * p0 +
                 pi_a (N1 + 1) * p01;
            *Q = pi_q (N1) * q1;
            *P = p01;

            break;
        }
        case 3:
        {
            apfloat p0 = pi_p (N1), p01 = p0 * pi_p (N1 + 1), p012 = p01 * pi_p (N1 + 2),
                    q2 = pi_q (N1 + 2), q12 = pi_q (N1 + 1) * q2;

            *T = q12 * pi_a (N1) * p0 +
                 q2 * pi_a (N1 + 1) * p01 +
                 pi_a (N1 + 2) * p012;
            *Q = pi_q (N1) * q12;
            *P = p012;