apcfunc.cpp

任意精度计算的实现

    if (M_PI - fabs (dx) <= 1.0 / sqrt (MAX_PRECISE_DOUBLE) &&
        M_PI - fabs (dy) <= 1.0 / sqrt (MAX_PRECISE_DOUBLE) &&
        dx * dy < 0.0)
    {
        // Both imaginary parts are near +-pi and of opposite sign
        if (y->im.sign () < 0)
            y->im += twopi;
        else
            y->im -= twopi;
    }
}

// Calculate needed precision for Readypi
size_t piprec (apcomplex *z)
{
    size_t prec;

    prec = z->im.prec ();

    if (prec != INFINITE)
    {
        prec -= z->im.exp () - Basedigits;
        if ((long) prec < Basedigits)
            prec = Basedigits;
    }

    if (prec == INFINITE)
        prec = z->prec ();

    return prec;
}

// Exponent function, calculated using Newton's iteration for the inverse of log
apcomplex exp (apcomplex u)
{
    size_t prec, minprec, maxprec, destprec = u.prec (), pprec = piprec (&u), doubledigits, fprec;
    int k, f;
    double d;
    apfloat x, y, twopi;
    apcomplex z, t;

    if (!u.re.sign () && !u.im.sign ()) return apcomplex (1);

    checkpi (pprec);

    twopi = Readypi;
    twopi.prec (pprec);
    twopi += twopi;

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

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

    // First handle the real part
    if (!u.re.sign ())
    {
        x = 1;
    }
    else if (u.re.exp () < -Basedigits)
    {
        // Taylor series: exp(x) = 1 + x + x^2/2 + ...

        x = u.re;
        x.prec (doubledigits);
        x += 1;
    }
    else
    {
        // Approximate starting value for iteration
        double i, f;

        // If u.re is too big, an overflow will occur (somewhere)
        d = ap2double (u.re.ap);
        i = floor (d);
        f = d - i;

        d = i / log ((double) Base);

        x = apfloat (exp (f) * pow ((double) Base, d - floor (d)), fprec);
        x.exp (x.exp () + Basedigits * (long) floor (d));
    }

    // Then handle the imaginary part
    u.im = fmod (u.im, twopi);
    if (u.im > Readypi)
        u.im -= twopi;
    else if (u.im <= -Readypi)
        u.im += twopi;

    // Re-calculate destination precision as u.im may have changed considerably
    destprec = min (destprec, piprec (&u));
    destprec = min (destprec, pprec);
    checklogconst (destprec);

    d = ap2double (u.im.ap);

    if (fabs (d) <= 1.0 / sqrt (MAX_PRECISE_DOUBLE))
    {
        // Taylor series: exp(z) = 1 + z + z^2/2 + ...

        y = u.im;
        y.prec (doubledigits);
        z = apcomplex (1, y);
    }
    else if (M_PI - fabs (d) <= 1.0 / sqrt (MAX_PRECISE_DOUBLE))
    {
        // exp(z + i*pi) = exp(z)*exp(i*pi) = -exp(z)
        // exp(z - i*pi) = exp(z)/exp(i*pi) = -exp(z)
        // Taylor series: exp(z) = 1 + z + z^2/2 + ...

        y = u.im;

        if (y.sign () < 0)
            y = y + Readypi;
        else
            y = y - Readypi;

        y.sign (-y.sign ());
        y.prec (doubledigits);

        z = apcomplex (-1, y);
    }
    else
    {
        // Approximate starting value for iteration

        z = apcomplex (apfloat (cos (d), fprec), apfloat (sin (d), fprec));
    }

    z *= x;

    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)
    {
        z.re.prec (minprec + 3 * Basedigits);
        z.im.prec (minprec + 3 * Basedigits);
        while (prec < minprec)
        {
            t = log (z);
            checkimsign (u, &t, twopi);

            z += z * (u - t);
            prec *= 2;
        }
        prec = minprec;
    }

    // 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 >= 3 * Basedigits && (minprec - 3 * Basedigits) << f >= destprec)
            break;

    // Newton's iteration
    while (k--)
    {
        prec *= 2;
        // Complex log needs a bit extra precision for convergence
        z.re.prec (max (4 * Basedigits, min (prec, destprec)));
        z.im.prec (max (4 * Basedigits, min (prec, destprec)));

        t = log (z);
        checkimsign (u, &t, twopi);

        t = u - t;
        if (k < f)
        {
            t.re.prec (prec / 2);
            t.im.prec (prec / 2);
        }

        z += z * t;

        // Precising iteration
        if (k == f)
        {
            t = log (z);
            checkimsign (u, &t, twopi);

            z += z * (u - t);
        }
    }

    z.re.prec (destprec);
    z.im.prec (destprec);

    return z;
}

// Arbitrary power, calculated using log and exp
apcomplex pow (apcomplex z, apcomplex w)
{
    size_t destprec = min (z.prec (), w.prec ());

    checklogconst (destprec);

    z.re.prec (destprec);
    z.im.prec (destprec);

    return exp (w * log (z));
}

apcomplex pow (apcomplex z, apfloat y)
{
    size_t destprec = min (z.prec (), y.prec ());

    checklogconst (destprec);

    z.re.prec (destprec);
    z.im.prec (destprec);

    return exp (y * log (z));
}

apcomplex pow (apfloat x, apcomplex w)
{
    size_t destprec = min (x.prec (), w.prec ());

    checklogconst (destprec);

    x.prec (destprec);

    return exp (w * log (x));
}


// Trigonometric and hyperbolic functions and their inverses

apcomplex acos (apcomplex z)
{
    apcomplex i = apcomplex (0, 1), w;

    if (z.re.sign () >= 0)
        w = i * log (z + sqrt (z * z - 1));
    else
        w = -i * log (z - sqrt (z * z - 1));

    if (z.re.sign () * z.im.sign () >= 0)
        return -w;
    else
        return w;
}

apcomplex acosh (apcomplex z)
{
    apcomplex w;

    if (z.re.sign () >= 0)
        return log (z + sqrt (z * z - 1));
    else
        return log (z - sqrt (z * z - 1));
}

apcomplex asin (apcomplex z)
{
    apcomplex i = apcomplex (0, 1);

    if (z.im.sign () >= 0)
        return i * log (sqrt (1 - z * z) - i * z);
    else
        return -i * log (i * z + sqrt (1 - z * z));
}

apcomplex asinh (apcomplex z)
{
    if (z.re.sign () >= 0)
        return log (sqrt (z * z + 1) + z);
    else
        return -log (sqrt (z * z + 1) - z);
}

apcomplex atan (apcomplex z)
{
    apcomplex i = apcomplex (0, 1);

    return log ((i + z) / (i - z)) * i / 2;
}

apcomplex atanh (apcomplex z)
{
    return log ((1 + z) / (1 - z)) / 2;
}

apcomplex cos (apcomplex z)
{
    apcomplex i = apcomplex (0, 1);
    apcomplex w = exp (i * z);

    return (w + 1 / w) / 2;
}

apcomplex cosh (apcomplex z)
{
    apcomplex w = exp (z);

    return (w + 1 / w) / 2;
}

apcomplex sin (apcomplex z)
{
    apcomplex i = apcomplex (0, 1);
    apcomplex w = exp (i * z);

    return (1 / w - w) * i / 2;
}

apcomplex sinh (apcomplex z)
{
    apcomplex w = exp (z);

    return (w - 1 / w) / 2;
}

apcomplex tan (apcomplex z)
{
    apcomplex i = apcomplex (0, 1);
    apcomplex w = exp (2 * i * z);

    return i * (1 - w) / (1 + w);
}

apcomplex tanh (apcomplex z)
{
    apcomplex w = exp (2 * z);

    return (w - 1) / (w + 1);
}


// Real trigonometric and hyperbolic functions and their inverses
// use complex functions

apfloat acos (apfloat x)
{
    apcomplex i = apcomplex (0, 1);

    return imag (log (x + i * sqrt (1 - x * x)));
}

apfloat asin (apfloat x)
{
    apcomplex i = apcomplex (0, 1);

    return -imag (log (sqrt (1 - x * x) - i * x));
}

apfloat atan (apfloat x)
{
    apcomplex i = apcomplex (0, 1);

    return imag (log ((i - x) / (i + x))) / 2;
}

apfloat atan2 (apfloat x, apfloat y)
{
    long tmpexp;
    apfloat t;

    if (!x.sign ())
    {
        assert (y.sign ());

        checkpi (y.prec ());

        t = Readypi;
        t.prec (y.prec ());

        return y.sign () * t / 2;
    }
    else if (!y.sign ())
    {
        if (x.sign () > 0) return 0;

        checkpi (x.prec ());

        t = Readypi;
        t.prec (x.prec ());

        return t;
    }
    else
    {
        tmpexp = x.exp ();
        if (y.exp () > tmpexp) tmpexp = y.exp ();
    }

    x.exp (x.exp () - tmpexp);
    y.exp (y.exp () - tmpexp);

    return imag (rawlog (apcomplex (x, y)));
}

apfloat cos (apfloat x)
{
    return real (exp (apcomplex (0, x)));
}

apfloat sin (apfloat x)
{
    return imag (exp (apcomplex (0, x)));
}

apfloat tan (apfloat x)
{
    apcomplex w = exp (apcomplex (0, x));

    return imag (w) / real (w);
}