apstruct.cpp

任意精度计算的实现

                case 1: bigger = a; smaller = b; sign = a->sign; break;
                case -1: bigger = b; smaller = a; sign = truebsign; break;
                default: return new apstruct;           // Zero
            }

            expdiff = bigger->exp - smaller->exp;
            worklen = max (bigger->size, expdiff + smaller->size);
            prec = min (bigger->prec, (smaller->prec == INFINITE ? INFINITE : expdiff + smaller->prec));

            ap = new apstruct (sign, bigger->exp, prec, worklen, 0, (worklen > Maxblocksize ? DISK : MEMORY));

            done = worklen;

            if (done > (p = expdiff + smaller->size))       // Copy block
            {
                p = done - p;
                done -= p;
                moveblock (ap, bigger, done, done, p);
            }
            else if (done > (p = bigger->size))             // Copy negated block
            {
                p = done - p;
                if (p > done - expdiff)
                    p = done - expdiff;
                done -= p;
                carry = negblock (ap, smaller, done, done - expdiff, p, carry);
            }

            // Subtract
            if (done > (p = expdiff))
            {
                p = done - p;
                done -= p;
                carry = subblock (ap, bigger, smaller, done, done, done - expdiff, p, carry);
            }

            // Propagate last carries
            if (done > (p = bigger->size))
            {
                p = done - p;
                done -= p;
                carry = negcarry (ap, done, p, carry);
            }

            p = done;
            done -= p;
            subcarry (ap, bigger, done, done, p, carry);

            p = worklen;        // Normalize
            r = 0;
            shift = 0;

            while (p)           // Get denormalization
            {
                l = min (p, Blocksize);
                p -= l;
                src = ap->getdata (r, l);
                r += l;
                for (t = 0; t < l; t++, shift++)
                    if ((rawtype) src[t] != 0)
                    {
                        p = 0;
                        break;
                    }
                ap->cleardata ();
            }

            if (shift == worklen || shift >= prec)
            {
                delete ap;
                return new apstruct;            // Zero
            }

            worklen = min (worklen, prec);      // This is bad to do here

            worklen -= lastzeros (ap, worklen);

            if (shift)
            {
                // Avoid unnecessary disk I/O
                apstruct *aptmp = new apstruct (ap->sign, ap->exp - shift, (prec == INFINITE ? INFINITE : prec - shift), worklen - shift, 0, DEFAULT, 0);   // No fill

                // Shift mantissa
                moveblock (aptmp, ap, 0, shift, worklen - shift);

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

            return ap;
        }
        else        // Add the mantissas
        {
            if (a->exp >= b->exp)
            {
                bigger = a;
                smaller = b;
                sign = a->sign;
            }
            else
            {
                bigger = b;
                smaller = a;
                sign = truebsign;
            }

            expdiff = bigger->exp - smaller->exp;
            worklen = max (bigger->size, expdiff + smaller->size);
            prec = min (bigger->prec, (smaller->prec == INFINITE ? INFINITE : expdiff + smaller->prec));

            ap = new apstruct (sign, bigger->exp, prec, worklen, 0, (worklen > Maxblocksize ? DISK : MEMORY));

            done = worklen;

            if (done > (p = expdiff + smaller->size))       // Copy block
            {
                p = done - p;
                done -= p;
                moveblock (ap, bigger, done, done, p);
            }
            else if (done > (p = bigger->size))             // Copy block
            {
                p = done - p;
                if (p > done - expdiff)
                    p = done - expdiff;
                done -= p;
                moveblock (ap, smaller, done, done - expdiff, p);
            }

            // Add
            if (done > (p = expdiff))
            {
                p = done - p;
                done -= p;
                carry = addblock (ap, bigger, smaller, done, done, done - expdiff, p, carry);
            }

            // Propagate last carries
            if (done > (p = bigger->size))
            {
                p = done - p;
                done -= p;
                carry = poscarry (ap, done, p, carry);
            }

            p = done;
            done -= p;
            carry = addcarry (ap, bigger, done, done, p, carry);

            worklen = min (worklen, prec);      // This is bad to do here

            worklen -= lastzeros (ap, worklen);

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

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

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

            return ap;
        }
    }
}

// The core of the apfloat multiplication
// Return a new apstruct
apstruct *apmul (apstruct *a, apstruct *b)
{
    int i, sign;
    size_t n, prec, size, rsize;
    apstruct *ap;
    size_t trueasize, truebsize;

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

    sign = a->sign * b->sign;
    if (!sign) return new apstruct;             // Zero

    prec = min (a->prec, b->prec);
    size = min (prec, a->size + b->size);

    trueasize = min (a->size, prec);
    truebsize = min (b->size, prec);

    n = rnd23up (trueasize + truebsize);

    rsize = min (n, size + 2);

    if (a == b)
        ap = autoconvolution (a, rsize, trueasize, n, &i);
    else
        ap = convolution (a, b, rsize, trueasize, truebsize, n, &i);

    size -= lastzeros (ap, size);

    ap->relocate (DEFAULT, size);

    ap->sign = sign;
    ap->prec = prec;
    ap->exp = a->exp + b->exp - 1 + i;

    return ap;
}

// "Medium size" apfloat multiplication, size of b is <= USE_N2_MUL and <= size of a
// Uses the "schoolboy" multiplication scheme, which has n^2 complexity
// Return a new apstruct
apstruct *apmulmedium (apstruct *a, apstruct *b)
{
    int sign;
    size_t t, s, p, l, prec, size, ssize;
    apstruct *ap;
    modint *f, *src, *dest;
    rawtype carry, tmpdest[USE_N2_MUL];
    modint tmpf[USE_N2_MUL];

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

    sign = a->sign * b->sign;
    if (!sign) return new apstruct;             // Zero

    for (s = 0; s < b->size; s++)
        tmpdest[s] = 0;

    if (b->location == MEMORY)
    {
        f = b->getdata (0, b->size);            // Pointer to a persistent memory address
        b->cleardata ();                        // Now we know that f is not invalidated
    }
    else
    {
        src = b->getdata (0, b->size);

        for (s = 0; s < b->size; s++)           // Create a copy of b in case a and b are the same
            tmpf[s] = src[s];

        f = tmpf;

        b->cleardata ();
    }

    prec = min (a->prec, b->prec);
    size = min (prec, a->size + b->size - 1);

    ssize = min (prec + b->size - 1, a->size + b->size - 1);

    ap = new apstruct (sign, a->exp + b->exp - 1, prec, ssize, 0, (ssize > Maxblocksize ? DISK : MEMORY));

    p = ssize - b->size + 1;

    while (p)
    {
        l = min (p, Blocksize);
        p -= l;
        src = a->getdata (p, l);
        dest = ap->readydata (p + b->size - 1, l);
        for (t = l; t--;)
        {
            carry = basemuladd (tmpdest, (rawtype *) f, tmpdest, (rawtype) src[t], b->size, 0);

            dest[t] = tmpdest[b->size - 1];
            for (s = b->size - 1; s--;)
                tmpdest[s + 1] = tmpdest[s];
            tmpdest[0] = carry;
        }
        ap->putdata ();
        a->cleardata ();
    }

    dest = ap->readydata (0, b->size - 1);

    for (s = b->size - 1; s--;)
        dest[s] = tmpdest[s + 1];

    ap->putdata ();

    size -= lastzeros (ap, size);

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

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

        delete ap;
        ap = aptmp;

        size -= lastzeros (ap);
    }

    ap->relocate (DEFAULT, size);

    return ap;
}

// "Short" apfloat multiplication, size of b is 1
// Return a new apstruct
apstruct *apmulshort (apstruct *a, apstruct *b)
{
    int sign;
    size_t p, l, prec, size;
    rawtype f;
    apstruct *ap;
    modint *src, *dest;
    rawtype carry = 0;

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

    sign = a->sign * b->sign;
    if (!sign) return new apstruct;             // Zero

    assert (b->size == 1);                      // For trapping bugs

    src = b->getdata (0, 1);
    f = src[0];
    b->cleardata ();

    prec = min (a->prec, b->prec);
    size = min (prec, a->size);

    p = size;

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

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

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

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

        size -= lastzeros (ap);

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

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

            delete ap;
            ap = aptmp;

            size -= lastzeros (ap);
        }
    }
    else
    {
        moveblock (ap, a, 0, 0, p);
    }

    ap->relocate (DEFAULT, size);

    return ap;
}

// "Short" apfloat division, size of b is 1
// Return a new apstruct
apstruct *apdivshort (apstruct *a, apstruct *b)
{
    int sign;
    size_t t, p, l, r, prec, getsize, ressize;
    rawtype f;
    apstruct *ap;
    modint *src, *dest;
    rawtype tmp, tmp1[2] = {0}, tmp2[2] = {0}, carry = 0;

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

    assert (b->sign);                           // Infinity

    sign = a->sign * b->sign;
    if (!sign) return new apstruct;             // Zero

    assert (b->size == 1);

    src = b->getdata (0, 1);
    f = src[0];
    b->cleardata ();

    prec = min (a->prec, b->prec);
    getsize = min (prec, a->size);

    if (f != 1)
    {
        // Check for finite or infinite result sequence
        tmp1[0] = f;
        for (t = 0; t < NBasefactors; t++)
            while (bigdiv (tmp2, tmp1, Basefactors[t], 1) == 0) tmp1[0] = tmp2[0];
        if (tmp1[0] != 1)
        {
            // Infinite nonzero sequence
            ressize = prec;
        }
        else
        {
            // Calculate maximum sequence length
            carry = 1;
            for (t = 0; carry != 0; t++)
            {
                tmp1[1] = bigmul (tmp1, &carry, Base, 1);
                carry = bigdiv (tmp1, tmp1, f, 2);
            }
            ressize = min (prec, getsize + t);