apfloat.cpp

任意精度计算的实现

#include <iomanip>
#include <sstream>
#include <cstring>
#include <cmath>
#include "ap.h"


using namespace std;


const int Maxint = (unsigned) -1 >> 1;
const unsigned Maxunsigned = (unsigned) -1;
const long Maxlong = (unsigned long) -1 >> 1;
const unsigned long Maxunsignedlong = (unsigned long) -1;

rawtype INVALID = 256;

// Numeric value to character, up to base 36
char chartable[] = {'0', '1', '2', '3', '4', '5', '6', '7', '8',
                    '9', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H',
                    'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q',
                    'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z'};

// Character to numeric value, up to base 36
rawtype valuetable[] = {INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID,
                        0, 1, 2, 3, 4, 5, 6, 7, 8, 9, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
                        25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
                        25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID,
                        INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID, INVALID};

// Characters used in exponent notation (do not mix with e.g. hex digits)
const char Expchar[2] = {'e', 'E'};

// Construct an apfloat from an integer
apfloat::apfloat (int value, size_t prec, int location)
{
    int sign;
    size_t t, size;
    long exp;
    modint tmpdata[3], *data;

    if (value > 0)
    {
        sign = 1;
    }
    else if (value < 0)
    {
        sign = -1;
        value = -value;
    }
    else
    {
        ap = new apstruct;
        return;
    }

    if (Base <= Maxint)
    {
        for (size = 0; size < 3 && value > 0; size++)
        {
            tmpdata[2 - size] = (rawtype) (value % (int) Base);
            value /= (int) Base;
        }
    }
    else
    {
        size = 1;                               // Nonzero
        tmpdata[2] = (rawtype) value;
    }

    exp = size;

    if (prec != INFINITE)
    {
        prec = (prec + Basedigits - 1) / Basedigits;
        if (prec < size)
            size = prec;
    }

    while ((rawtype) tmpdata[3 - exp + size - 1] == 0) size--;

    data = new modint[size];

    for (t = 0; t < size; t++)
        data[t] = tmpdata[t + 3 - exp];

    ap = new apstruct (sign, exp, prec, size, data, location);
}

// Construct an apfloat from an unsigned integer
apfloat::apfloat (unsigned value, size_t prec, int location)
{
    int sign;
    size_t t, size;
    long exp;
    modint tmpdata[3], *data;

    if (value > 0)
    {
        sign = 1;
    }
    else
    {
        ap = new apstruct;
        return;
    }

    if (Base <= Maxunsigned)
    {
        for (size = 0; size < 3 && value > 0; size++)
        {
            tmpdata[2 - size] = (rawtype) (value % (unsigned) Base);
            value /= (unsigned) Base;
        }
    }
    else
    {
        size = 1;                               // Nonzero
        tmpdata[2] = (rawtype) value;
    }

    exp = size;

    if (prec != INFINITE)
    {
        prec = (prec + Basedigits - 1) / Basedigits;
        if (prec < size)
            size = prec;
    }

    while ((rawtype) tmpdata[3 - exp + size - 1] == 0) size--;

    data = new modint[size];

    for (t = 0; t < size; t++)
        data[t] = tmpdata[t + 3 - exp];

    ap = new apstruct (sign, exp, prec, size, data, location);
}

// Construct an apfloat from a long integer
apfloat::apfloat (long value, size_t prec, int location)
{
    int sign;
    size_t t, size;
    long exp;
    modint tmpdata[3], *data;

    if (value > 0)
    {
        sign = 1;
    }
    else if (value < 0)
    {
        sign = -1;
        value = -value;
    }
    else
    {
        ap = new apstruct;
        return;
    }

    if (Base <= Maxlong)
    {
        for (size = 0; size < 3 && value > 0; size++)
        {
            tmpdata[2 - size] = (rawtype) (value % (long) Base);
            value /= (long) Base;
        }
    }
    else
    {
        size = 1;                               // Nonzero
        tmpdata[2] = (rawtype) value;
    }

    exp = size;

    if (prec != INFINITE)
    {
        prec = (prec + Basedigits - 1) / Basedigits;
        if (prec < size)
            size = prec;
    }

    while ((rawtype) tmpdata[3 - exp + size - 1] == 0) size--;

    data = new modint[size];

    for (t = 0; t < size; t++)
        data[t] = tmpdata[t + 3 - exp];

    ap = new apstruct (sign, exp, prec, size, data, location);
}

// Construct an apfloat from an unsigned long integer
apfloat::apfloat (unsigned long value, size_t prec, int location)
{
    int sign;
    size_t t, size;
    long exp;
    modint tmpdata[3], *data;

    if (value > 0)
    {
        sign = 1;
    }
    else
    {
        ap = new apstruct;
        return;
    }

    if (Base <= Maxunsignedlong)
    {
        for (size = 0; size < 3 && value > 0; size++)
        {
            tmpdata[2 - size] = (rawtype) (value % (unsigned long) Base);
            value /= (unsigned long) Base;
        }
    }
    else
    {
        size = 1;                               // Nonzero
        tmpdata[2] = (rawtype) value;
    }

    exp = size;

    if (prec != INFINITE)
    {
        prec = (prec + Basedigits - 1) / Basedigits;
        if (prec < size)
            size = prec;
    }

    while ((rawtype) tmpdata[3 - exp + size - 1] == 0) size--;

    data = new modint[size];

    for (t = 0; t < size; t++)
        data[t] = tmpdata[t + 3 - exp];

    ap = new apstruct (sign, exp, prec, size, data, location);
}

// Construct an apfloat from a double
apfloat::apfloat (double value, size_t prec, int location)
{
    int sign;
    size_t t, size;
    long exp;
    double tmp;
    modint tmpdata[4], *data;

    if (value > 0.0)
    {
        sign = 1;
    }
    else if (value < 0.0)
    {
        sign = -1;
        value = -value;
    }
    else
    {
        ap = new apstruct;
        return;
    }

    exp = (long) floor (log (value) / log ((double) Base));

    // Avoid overflow in intermediate value
    if (exp > 0)
    {
        value *= pow ((double) Base, (double) -exp);
    }
    else if (exp < 0)
    {
        value *= pow ((double) Base, (double) (-exp - 4));
        value *= pow ((double) Base, (double) 4);
    }

    exp++;

    for (size = 0; size < 4 && value > 0.0; size++)
    {
        tmp = floor (value);
        tmpdata[size] = (rawtype) tmp;
        value -= tmp;
        value *= Base;
    }

    if (prec == (size_t) DEFAULT)
    {
        // Approximate accuracy of a double
        size_t doubledigits = (size_t) (log (MAX_PRECISE_DOUBLE) / log ((double) Basedigit));

        prec = (doubledigits - (size_t) (log ((double) (rawtype) tmpdata[0]) / log ((double) Basedigit)) + 2 * Basedigits - 1) / Basedigits;
    }
    else
    {
        if (prec != INFINITE)
            prec = (prec + Basedigits - 1) / Basedigits;
    }

    if (prec < size)
        size = prec;

    while ((rawtype) tmpdata[size - 1] == 0) size--;

    data = new modint[size];

    for (t = 0; t < size; t++)
        data[t] = tmpdata[t];

    ap = new apstruct (sign, exp, prec, size, data, location);
}

// Construct an apfloat from a character string
apfloat::apfloat (char *valuestring, size_t prec, int location)
{
    size_t t, r, l, e, d;
    int sign = 1, dot = 0;
    bool nonzero = false;
    size_t intsize = 0, fracsize = 0, leadzeros = 0, size = 0, lastzeros = 0;
    long exp = 0;
    unsigned long tmpexp;
    rawtype val, digit;
    modint *data;
    char c;

    // First scan the size of the number, integer and fractional parts
    // and exponent

    // Sign & whitespace

    for (t = 0; (c = valuestring[t]) != '\0'; t++)
        if (c == '-')
            sign = -sign;
        else if (c != '+' && c != ' ')
            break;

    valuestring += t;

    // Integer part

    for (intsize = 0; (c = valuestring[intsize]) != '\0'; intsize++)
        if (c == '.' || c == Expchar[0] || c == Expchar[1])
            break;
        else if (c == chartable[0])
        {
            if (!nonzero) leadzeros++;
            lastzeros++;
        }
        else
        {
            lastzeros = 0;
            nonzero = true;
        }

    // Fractional part

    if (c == '.')
    {
        dot = 1;
        for (fracsize = 1 + intsize; (c = valuestring[fracsize]) != '\0'; fracsize++)
            if (c == Expchar[0] || c == Expchar[1])
                break;
            else if (c == chartable[0])
            {
                if (!nonzero) leadzeros++;
                lastzeros++;
            }
            else
            {
                lastzeros = 0;
                nonzero = true;
            }
        fracsize -= 1 + intsize;
    }

    // Exponent, always in base 10

    if (c == Expchar[0] || c == Expchar[1])
    {
        istringstream s (valuestring + intsize + dot + fracsize + 1);
        s >> dec >> exp;
    }

    if (!nonzero)
    {
        ap = new apstruct;              // It's zero
        return;
    }

    exp += intsize - leadzeros;
    tmpexp = exp;
    tmpexp += (Maxlong / Basedigits) * Basedigits;

    // Align exponent

    e = tmpexp % Basedigits;
    t = (e ? Basedigits - e : e);

    // Fracional part only or both integer and fractional part

    if (leadzeros >= intsize)
    {
        valuestring += leadzeros + 1;
        l = intsize + fracsize - leadzeros - lastzeros;
        size = (l + t + Basedigits - 1) / Basedigits;
    }
    else
    {
        valuestring += leadzeros;
        l = intsize + fracsize - leadzeros - lastzeros + dot;
        size = (l - dot + t + Basedigits - 1) / Basedigits;
    }

    if (prec == (size_t) DEFAULT)
        prec = size;
    else
    {
        if (prec != INFINITE)
            prec = (prec + Basedigits - 1) / Basedigits;
        if (prec < size)
            size = prec;
    }

    data = new modint[size];

    if (!e) e = Basedigits;

    // Get the actual data

    val = 0;
    for (t = 0, r = 0, d = 0; t < l && d < size; t++)
    {
        if ((c = valuestring[t]) == '.') continue;
        digit = valuetable[(unsigned char) c];
        if (digit == INVALID) digit = 0;        // Ignore invalid characters
        val = val * Basedigit + digit;
        if (++r == e)
        {
            r = 0;
            data[d++] = val;
            val = 0;
            e = Basedigits;
        }
    }

    // Last base unit

    if (r && d < size)
    {
        for (t = r; t < e; t++) val *= Basedigit;
        data[d++] = val;
    }

    exp = (exp >= 0 ? (exp + Basedigits - 1) / Basedigits : exp / Basedigits);

    ap = new apstruct (sign, exp, prec, size, data, location);
}

apfloat::apfloat (const apfloat &d)
{
    ap = d.ap;

    if (!ap) return;

    ap->nlinks++;
}

apfloat::~apfloat ()
{
    if (!ap) return;

    if (!(--(ap->nlinks)))
        delete ap;
}

int apfloat::sign (void) const
{
    if (ap) return ap->sign;
    else return 0;
}

void apfloat::sign (int newsign)
{
    if (!ap) return;

    if (ap->sign == newsign) return;

    unique ();

    ap->sign = newsign;
}

long apfloat::exp (void) const
{
    if (ap) return ap->exp * Basedigits;
    else return 0;
}

void apfloat::exp (long newexp)
{
    if (!ap) return;
    if (!ap->sign) return;              // Zero is zero

    newexp /= Basedigits;

    if (ap->exp == newexp) return;

    unique ();

    ap->exp = newexp;
}

size_t apfloat::prec (void) const
{
    if (ap)
    {
        if (ap->prec == INFINITE)
            return INFINITE;
        else
            return ap->prec * Basedigits;
    }
    else return 0;
}

void apfloat::prec (size_t newprec)
{
    if (!ap) return;
    if (!ap->sign) return;              // Zero is zero

    if (newprec != INFINITE)