apstruct.cpp

任意精度计算的实现

        }

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

        p = getsize;
        r = 0;

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

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

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

        p = ressize - getsize;

        src = a->getdata (0, 1);
        tmp = src[0];
        a->cleardata ();

        if (tmp < f)
        {
            // Avoid unnecessary disk I/O
            apstruct *aptmp = new apstruct (sign, ap->exp - 1, prec, ressize, 0, DEFAULT, 0);   // No fill

            tmp1[1] = bigmul (tmp1, &carry, Base, 1);
            carry = bigdiv (tmp1, tmp1, f, 2);

            // Shift mantissa
            shlblock (aptmp, ap, getsize - 1, tmp1[0]);

            delete ap;
            ap = aptmp;
        }

        while (p)
        {
            l = min (p, Maxblocksize);
            p -= l;
            dest = ap->readydata (r, l);
            r += l;

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

            ap->putdata ();
        }
    }
    else
    {
        ap = new apstruct (sign, a->exp - b->exp + 1, prec, getsize, 0, (getsize > Maxblocksize ? DISK : MEMORY), 0);   // No fill

        moveblock (ap, a, 0, 0, getsize);
    }

    ap->relocate (DEFAULT, ap->size - lastzeros (ap));

    return ap;
}

// Apfloat comparison
// Return 1 if a > b
//       -1 if a < b
//        0 if a == b
int apcmp (apstruct *a, apstruct *b, int absolute)
{
    size_t m, t, p, l, r;
    modint *da, *db;
    int trueasign, truebsign;

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

    if (a == b) return 0;       // Same number

    if (absolute)
    {
        trueasign = abs (a->sign);
        truebsign = abs (b->sign);
    }
    else
    {
        trueasign = a->sign;
        truebsign = b->sign;
    }

    if (trueasign < truebsign) return -1;
    else if (trueasign > truebsign) return 1;

    if (trueasign < 0)
    {
        if (a->exp < b->exp) return 1;
        else if (a->exp > b->exp) return -1;
    }
    else if (trueasign > 0)
    {
        if (a->exp < b->exp) return -1;
        else if (a->exp > b->exp) return 1;
    }

    p = m = min (a->size, b->size);
    r = 0;

    while (p)
    {
        l = min (p, Blocksize);
        p -= l;
        da = a->getdata (r, l);
        db = b->getdata (r, l);
        r += l;
        for (t = 0; t < l; t++)
        {
            if ((rawtype) da[t] < (rawtype) db[t])
            {
                a->cleardata ();
                b->cleardata ();
                return -trueasign;
            }
            else if ((rawtype) da[t] > (rawtype) db[t])
            {
                a->cleardata ();
                b->cleardata ();
                return trueasign;
            }
        }
        a->cleardata ();
        b->cleardata ();
    }

    if (a->size > m)
        return trueasign;
    else if (b->size > m)
        return -trueasign;

    return 0;
}

// Convert apfloat to double
double ap2double (apstruct *a)
{
    size_t doubledigits = (size_t) (log (MAX_PRECISE_DOUBLE) / log ((double) Basedigit)),     // Approximate accuracy of a double
           size = (doubledigits + 2 * Basedigits - 2) / Basedigits;
    modint *data;
    double i = 0.0;

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

    if (!a->sign) return 0.0;

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

    data = a->getdata (0, size);

    while (size--)
    {
        i += (double) (rawtype) data[size];
        i /= Base;
    }

    a->cleardata ();

    return a->sign * i * pow ((double) Base, (double) a->exp);
}

// Convert apfloat to rawtype
// The apfloat must be integer with size one
rawtype ap2rawtype (apstruct *a)
{
    rawtype r;

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

    if (a->sign)
    {
        assert (a->size == 1);
        assert (a->exp == 1);

        modint *data = a->getdata (0, 1);
        r = (rawtype) data[0];
        a->cleardata ();
    }
    else
    {
        r = 0;
    }

    return r;
}

// Rounds down in _absolute value_ an apfloat to the nearest integer
// Returns a new apstruct
apstruct *apabsfloor (apstruct *a)
{
    size_t size;
    apstruct *ap;

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

    if (!a->sign || a->exp <= 0)
        return new apstruct;                            // Zero

    size = min (a->exp, a->size);                       // Size of integer part

    size -= lastzeros (a, size);

    ap = new apstruct (a->sign, a->exp, INFINITE, size, 0, DEFAULT, 0);         // No fill

    moveblock (ap, a, 0, 0, size);

    return ap;
}

// Rounds up in _absolute value_ an apfloat to the nearest integer
// Assumes that the input is normalized and has nonzero last digit (=base unit)
// Returns a new apstruct
apstruct *apabsceil (apstruct *a)
{
    size_t size;
    rawtype carry;
    apstruct *ap;

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

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

    if (a->exp <= 0)                                    // Result is 1 (or -1)
    {
        modint *data = new modint[1];

        data[0] = (rawtype) 1;

        return new apstruct (a->sign, 1, INFINITE, 1, data, DEFAULT);
    }

    // Check if the fractional part is nonzero
    if (a->size > a->exp)
        carry = 1;
    else
        carry = 0;

    size = min (a->exp, a->size);                       // Size of integer part

    if (carry == 0)                                     // Simply truncate
    {
        size -= lastzeros (a, size);

        ap = new apstruct (a->sign, a->exp, INFINITE, size, 0, DEFAULT, 0);     // No fill

        moveblock (ap, a, 0, 0, size);

        return ap;
    }

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

    // Add 1 to the integer part
    carry = addcarry (ap, a, 0, 0, size, carry);

    size -= lastzeros (ap);

    if (carry != 0)
    {
        size++;

        // Avoid unnecessary disk I/O
        apstruct *aptmp = new apstruct (ap->sign, ap->exp + 1, ap->prec, 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;
}

// Apfloat comparison (for AGM)
// Return the number of matching _base units_
size_t apeq (apstruct *a, apstruct *b)
{
    size_t t, p, l, r;
    modint *da, *db;

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

    if (a->sign != b->sign) return 0;           // No match

    if (!a->sign) return INFINITE;              // Both zero

    if (a->exp != b->exp) return 0;             // No match

    if (a == b) return a->prec;                 // Same number

    p = min (a->size, b->size);
    r = 0;

    while (p)
    {
        l = min (p, Blocksize);
        p -= l;
        da = a->getdata (r, l);
        db = b->getdata (r, l);

        for (t = 0; t < l; t++, r++)
            if ((rawtype) da[t] != (rawtype) db[t])
            {
                a->cleardata ();
                b->cleardata ();
                return r;
            }

        a->cleardata ();
        b->cleardata ();
    }

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

    if (a->size > r)
    {
        p = min (a->size, p) - r;

        while (p)
        {
            l = min (p, Blocksize);
            p -= l;
            da = a->getdata (r, l);

            for (t = 0; t < l; t++, r++)
                if ((rawtype) da[t] != 0)
                {
                    a->cleardata ();
                    return r;
                }

            a->cleardata ();
        }
    }
    else if (b->size > r)
    {
        p = min (b->size, p) - r;

        while (p)
        {
            l = min (p, Blocksize);
            p -= l;
            db = b->getdata (r, l);

            for (t = 0; t < l; t++, r++)
                if ((rawtype) db[t] != 0)
                {
                    b->cleardata ();
                    return r;
                }

            b->cleardata ();
        }
    }

    return min (a->prec, b->prec);
}

typedef struct
{
    size_t size;
    size_t prec;
    long exp;
    int sign;
} apdiskheader;

// Read apstruct from file
// Returns a new apstruct
apstruct *apswapfrom (const char *filename)
{
    apstruct *ap = new apstruct ();
    apdiskheader h;

    ap->capture (filename);

    fstream &fs = ap->openstream ();
    fs.seekg (- (int) sizeof (h), ios::end);
    fs.read ((char *) &h, sizeof (h));
    assert (fs.good ());
    ap->closestream ();

    ap->prec = h.prec;
    ap->exp = h.exp;
    ap->sign = h.sign;
    ap->relocate (DEFAULT, h.size);

    return ap;
}

// Writes apstruct to file
// Deletes a
void apswapto (apstruct *a, const char *filename)
{
    apdiskheader h;

    assert (a);

    h.prec = a->prec;
    h.exp = a->exp;
    h.sign = a->sign;
    h.size = a->size;

    a->relocate (DISK);

    fstream &fs = a->openstream ();
    fs.seekp (0, ios::end);
    fs.write ((char *) &h, sizeof (h));
    assert (fs.good ());
    a->closestream ();

    a->release (filename);

    delete a;
}