apstruct.cpp

任意精度计算的实现

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


using namespace std;


// Round up to the nearest power of two or three times a power of two
size_t rnd23up (size_t x)
{
    size_t r = 1;

    if (!x) return 0;

    while (r < x)
    {
        if (r == 1)
            r = 2;
        else if (r == (r & -r))
            r = r / 2 * 3;
        else
            r = r / 3 * 4;
    }

    return r;
}

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);
}

// Simply move data from one mantissa to another
void moveblock (apstruct *d, apstruct *s, size_t dpos, size_t spos, size_t len)
{
    size_t l;
    modint *src, *dest;

    while (len)
    {
        l = min (len, Blocksize);
        len -= l;
        src = s->getdata (spos, l);
        dest = d->readydata (dpos, l);
        spos += l;
        dpos += l;

        moveraw (dest, src, l);

        d->putdata ();
        s->cleardata ();
    }
}

// Add data from one mantissa to another
// Return carry
rawtype addblock (apstruct *d, apstruct *s1, apstruct *s2, size_t dpos, size_t s1pos, size_t s2pos, size_t len, rawtype carry)
{
    size_t l;
    modint *src1, *src2, *dest;

    s1pos += len;
    s2pos += len;
    dpos += len;

    while (len)
    {
        l = min (len, Blocksize);
        len -= l;
        s1pos -= l;
        s2pos -= l;
        dpos -= l;
        src1 = s1->getdata (s1pos, l);
        src2 = s2->getdata (s2pos, l);
        dest = d->readydata (dpos, l);

        carry = baseadd ((rawtype *) dest, (rawtype *) src1, (rawtype *) src2, l, carry);

        d->putdata ();
        s2->cleardata ();
        s1->cleardata ();
    }

    return carry;
}

// Subtract data from one mantissa to another
// Return carry
rawtype subblock (apstruct *d, apstruct *s1, apstruct *s2, size_t dpos, size_t s1pos, size_t s2pos, size_t len, rawtype carry)
{
    size_t l;
    modint *src1, *src2, *dest;

    s1pos += len;
    s2pos += len;
    dpos += len;

    while (len)
    {
        l = min (len, Blocksize);
        len -= l;
        s1pos -= l;
        s2pos -= l;
        dpos -= l;
        src1 = s1->getdata (s1pos, l);
        src2 = s2->getdata (s2pos, l);
        dest = d->readydata (dpos, l);

        carry = basesub ((rawtype *) dest, (rawtype *) src1, (rawtype *) src2, l, carry);

        d->putdata ();
        s2->cleardata ();
        s1->cleardata ();
    }

    return carry;
}

// Negate data from one mantissa to another
// Return carry
rawtype negblock (apstruct *d, apstruct *s, size_t dpos, size_t spos, size_t len, rawtype carry)
{
    size_t l;
    modint *src, *dest;

    spos += len;
    dpos += len;

    while (len)
    {
        l = min (len, Blocksize);
        len -= l;
        spos -= l;
        dpos -= l;
        src = s->getdata (spos, l);
        dest = d->readydata (dpos, l);

        carry = basesub ((rawtype *) dest, 0, (rawtype *) src, l, carry);

        d->putdata ();
        s->cleardata ();
    }

    return carry;
}

// Shift data right one modint
void shrblock (apstruct *d, apstruct *s, size_t len, rawtype carry)
{
    modint *dest;

    dest = d->readydata (0, 1);
    dest[0] = carry;
    d->putdata ();

    moveblock (d, s, 1, 0, len);
}

// Shift data left one modint
void shlblock (apstruct *d, apstruct *s, size_t len, rawtype carry)
{
    modint *dest;

    moveblock (d, s, 0, 1, len);

    dest = d->readydata (len, 1);
    dest[0] = carry;
    d->putdata ();
}

// Propagate carries from subtract
// Return carry
rawtype subcarry (apstruct *d, apstruct *s, size_t dpos, size_t spos, size_t len, rawtype carry)
{
    size_t l;
    modint *src, *dest;

    spos += len;
    dpos += len;

    while (len && carry != 0)
    {
        l = min (len, Blocksize);
        len -= l;
        spos -= l;
        dpos -= l;
        src = s->getdata (spos, l);
        dest = d->readydata (dpos, l);

        carry = basesub ((rawtype *) dest, (rawtype *) src, 0, l, carry);

        d->putdata ();
        s->cleardata ();
    }

    if (len) moveblock (d, s, 0, 0, len);

    return carry;
}

// Propagate carries from negation
// Return carry
rawtype negcarry (apstruct *ds, size_t pos, size_t len, rawtype carry)
{
    size_t t, l;
    modint *src;

    pos += len;

    while (len && carry != 0)
    {
        l = min (len, Maxblocksize);
        len -= l;
        pos -= l;
        src = ds->getdata (pos, l);
        for (t = l; t--;)
        {
            if (bigsub ((rawtype *) src + t, &carry, 1) != 0)
            {
                bigadd ((rawtype *) src + t, &Base, 1);
                carry = 1;
            }
            else
            {
                carry = 0;
                break;
            }
        }
        ds->putdata ();
    }

    return carry;
}

// Propagate carries from addition
// Return carry
rawtype addcarry (apstruct *d, apstruct *s, size_t dpos, size_t spos, size_t len, rawtype carry)
{
    size_t l;
    modint *src, *dest;

    spos += len;
    dpos += len;

    while (len && carry != 0)
    {
        l = min (len, Blocksize);
        len -= l;
        spos -= l;
        dpos -= l;
        src = s->getdata (spos, l);
        dest = d->readydata (dpos, l);

        carry = baseadd ((rawtype *) dest, (rawtype *) src, 0, l, carry);

        d->putdata ();
        s->cleardata ();
    }

    if (len) moveblock (d, s, 0, 0, len);

    return carry;
}

// Propagate carries from adding zeros
// Return carry
rawtype poscarry (apstruct *ds, size_t pos, size_t len, rawtype carry)
{
    size_t t, l;
    modint *src;

    pos += len;

    while (len && carry != 0)
    {
        l = min (len, Maxblocksize);
        len -= l;
        pos -= l;
        src = ds->getdata (pos, l);
        for (t = l; t--;)
        {
            if (bigadd ((rawtype *) src + t, &carry, 1) != 0 ||
                bigcmp ((rawtype *) src + t, &Base, 1) >= 0)
            {
                bigsub ((rawtype *) src + t, &Base, 1);
                carry = 1;
            }
            else
            {
                carry = 0;
                break;
            }
        }
        ds->putdata ();
    }

    return carry;
}

// Returns number of trailing zeros
// len is optional number length
size_t lastzeros (apstruct *s, size_t len)
{
    size_t t, l, r = 0;
    modint *src;

    if (len == (size_t) DEFAULT) len = s->size;

    while (len)
    {
        l = min (len, Blocksize);
        len -= l;
        src = s->getdata (len, l);
        t = l;
        while (t-- && (rawtype) src[t] == 0) r++;
        s->cleardata ();
        if (t != (size_t) -1) break;
    }

    return r;
}

// The core of the apfloat addition and subtraction
// Return a new apstruct
apstruct *apaddsub (apstruct *a, apstruct *b, int sub)
{
    size_t t, p, l, r;
    apstruct *ap;
    modint *src, *dest;
    rawtype carry = 0;
    int truebsign, truefunc;

    assert (a);                 // Won't work on uninitialized apfloats
    assert (b);

    if (!b->sign)
    {
        a->nlinks++;
        return a;
    }

    if (!a->sign)
    {
        if (sub)
        {
            ap = new apstruct (*b);
            ap->sign = -ap->sign;
            return ap;
        }
        else
        {
            b->nlinks++;
            return b;
        }
    }

    // Now neither a or b is zero

    truebsign = (sub ? -b->sign : b->sign);
    truefunc = (a->sign * truebsign < 0 ? 1 : 0);

    if (a == b)
    {
        if (truefunc)
        {
            // Result is zero
            return new apstruct;
        }
        else
        {
            // Result is 2 * a
            size_t size = a->size;

            p = size;
            ap = new apstruct (a->sign, a->exp, a->prec, p, 0, (p > Maxblocksize ? DISK : MEMORY));

            while (p)
            {
                l = min (p, Blocksize);
                p -= l;
                src = a->getdata (p, l);
                dest = ap->readydata (p, l);

                carry = baseadd ((rawtype *) dest, (rawtype *) src, (rawtype *) src, l, carry);

                ap->putdata ();
                a->cleardata ();
            }

            size -= lastzeros (ap);

            if (carry != 0)
            {
                size++;                         // Precision increases always
                // Avoid unnecessary disk I/O
                apstruct *aptmp = new apstruct (ap->sign, ap->exp + 1, (ap->prec == INFINITE ? INFINITE : ap->prec + 1), size, 0, DEFAULT, 0);  // No fill

                // Shift mantissa
                shrblock (aptmp, ap, size - 1, carry);

                delete ap;
                ap = aptmp;
            }
            else ap->relocate (DEFAULT, size);

            return ap;
        }
    }
    else
    {
        int sign;
        size_t worklen, prec, done;
        long expdiff;
        apstruct *bigger, *smaller;

        if (truefunc)       // True subtract
        {
            int c;
            size_t shift;

            c = apcmp (a, b, 1);    // This *might* do a lot of unnecessary work

            switch (c)
            {