apfloat.cpp

任意精度计算的实现

        newprec = (newprec + Basedigits - 1) / Basedigits;

    if (ap->prec == newprec) return;

    if (ap->nlinks > 1)
    {
        ap->nlinks--;

        if (newprec < ap->size)
            ap = new apstruct (*ap, newprec - lastzeros (ap, newprec));
        else
            ap = new apstruct (*ap);
    }
    else
    {
        if (newprec < ap->size)
            ap->resize (newprec - lastzeros (ap, newprec));
    }

    ap->prec = newprec;
}

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

void apfloat::location (int newlocation)
{
    if (!ap) return;

    if (ap->location == newlocation) return;

    unique ();

    ap->relocate (newlocation);
}

void apfloat::unique (void)
{
    if (!ap) return;

    if (ap->nlinks > 1)
    {
        ap->nlinks--;
        ap = new apstruct (*ap);
    }
}

apfloat apfloat::operator+ () const
{
    return *this;
}

apfloat apfloat::operator- () const
{
    apstruct *m = new apstruct (*ap);

    m->sign = -m->sign;

    return apfloat (m);
}


// Definitions of compound-assignment operator member functions

apfloat &apfloat::operator++ ()
{
    *this = *this + apfloat (1);

    return *this;
}

apfloat &apfloat::operator-- ()
{
    *this = *this - apfloat (1);

    return *this;
}

// Post-increment operators by Adam Pawlowski

apfloat apfloat::operator++ (int d)
{
    apfloat tmp = *this;
    ++(*this);

    return tmp;
}

apfloat apfloat::operator-- (int d)
{
    apfloat tmp = *this;
    --(*this);

    return tmp;
}

apfloat &apfloat::operator= (const apfloat &d)
{
    if (ap != d.ap)
    {
        if (ap)
            if (!(--(ap->nlinks)))
                delete ap;

        ap = d.ap;
        if (ap)
            ap->nlinks++;
    }
    return *this;
}


// Definitions of non-member binary operator functions

apfloat operator+ (const apfloat &d1, const apfloat &d2)
{
    return apfloat (apaddsub (d1.ap, d2.ap, 0));
}

apfloat operator- (const apfloat &d1, const apfloat &d2)
{
    return apfloat (apaddsub (d1.ap, d2.ap, 1));
}

apfloat operator* (const apfloat &d1, const apfloat &d2)
{
    assert (d1.ap);
    assert (d2.ap);

    if (d1.ap->size <= 1) return apfloat (apmulshort (d2.ap, d1.ap));
    else if (d2.ap->size <= 1) return apfloat (apmulshort (d1.ap, d2.ap));
    else if (d1.ap->size <= USE_N2_MUL && d1.ap->size <= d2.ap->size) return apfloat (apmulmedium (d2.ap, d1.ap));
    else if (d2.ap->size <= USE_N2_MUL) return apfloat (apmulmedium (d1.ap, d2.ap));
    else return apfloat (apmul (d1.ap, d2.ap));
}

apfloat operator/ (const apfloat &d1, const apfloat &d2)
{
    assert (d1.ap);
    assert (d2.ap);

    if (d2.ap->size <= 1) return apfloat (apdivshort (d1.ap, d2.ap));
    else return d1 * invroot (d2, 1, (d1.prec () < d2.prec () ? d1.prec () : d2.prec ()));
}

bool operator== (const apfloat &d1, const apfloat &d2)
{
    return apcmp (d1.ap, d2.ap) == 0;
}

bool operator!= (const apfloat &d1, const apfloat &d2)
{
    return apcmp (d1.ap, d2.ap) != 0;
}

bool operator>= (const apfloat &d1, const apfloat &d2)
{
    return apcmp (d1.ap, d2.ap) >= 0;
}

bool operator<= (const apfloat &d1, const apfloat &d2)
{
    return apcmp (d1.ap, d2.ap) <= 0;
}

bool operator> (const apfloat &d1, const apfloat &d2)
{
    return apcmp (d1.ap, d2.ap) > 0;
}

bool operator< (const apfloat &d1, const apfloat &d2)
{
    return apcmp (d1.ap, d2.ap) < 0;
}

ostream &pretty (ostream &str)
{
    apfloat::prettyprint = true;

    return str;
}

// Output a base unit in arbitrary base
void baseout (ostream &str, rawtype x, int zerofill = 0)
{
    size_t t;
    char tmp[65];

    tmp[Basedigits] = '\0';

    for (t = Basedigits; t && x > 0;)
        tmp[--t] = chartable[(int) bigdiv (&x, &x, Basedigit, 1)];

    if (zerofill)
    {
        while (t--)
            tmp[t] = chartable[0];
        str << tmp;
    }
    else
    {
        str << (tmp + t);
    }
}

// Output an apfloat to a stream
ostream &operator<< (ostream &str, const apfloat &d)
{
    int f;
    size_t t, p, l, r;
    long e = d.ap->exp;
    modint *data;

    assert (d.ap);

    f = d.ap->sign;
    if (f != 1)
    {
        if (f == -1)
            str << "-";
        else
        {
            str << chartable[0];
            apfloat::prettyprint = false;
            return str;
        }
    }
    f = 0;

    if (apfloat::prettyprint)
    {
        if (e <= 0)
        {
            str << chartable[0] << '.';
            for (; e; e++)
                for (t = Basedigits; t--;)
                    str << chartable[0];
            e--;
            f = 1;
        }
    }
    else
    {
        str << chartable[0] << '.';
        e = -1;
        f = 1;
    }

    p = (d.ap->prec < d.ap->size ? d.ap->prec : d.ap->size);
    r = 0;

    while (p)
    {
        l = (p < Maxblocksize ? p : Maxblocksize);
        p -= l;
        data = d.ap->getdata (r, l);
        r += l;
        for (t = 0; t < l; t++)
        {
            if (!e) str << '.';
            e--;
            baseout (str, (rawtype) data[t], f);
            f = 1;
        }
        d.ap->cleardata ();
    }

    for (; e > 0; e--)
        for (t = Basedigits; t--;)
            str << chartable[0];

    if (!apfloat::prettyprint)
        str << Expchar[0] << d.ap->exp * Basedigits;

    apfloat::prettyprint = false;

    return str;
}

// Flush output buffer if full
void buffercheck (size_t *t, apstruct *ap, modint **data)
{
    if (*t == Memorytreshold)
    {
        if (ap->location == MEMORY)
        {
            // Memory buffer full, relocate to disk since we don't know input size
            ap->relocate (DISK);                    // This does delete[] data
            *data = new modint[Memorytreshold];     // New disk write buffer
        }
        else
        {
            // Flush buffer to disk
            fstream &fs = ap->openstream ();
            fs.seekp (0, ios::end);
            fs.write ((char *) *data, sizeof (modint) * *t);
            assert (fs.good ());
            ap->closestream ();
        }

        *t = 0;
    }
}

// Input an apfloat from a stream
istream &operator>> (istream &str, apfloat &d)
{
    size_t t = 0, r = 0, l, e;
    int sign = 1;
    bool nonzero = false, gotsign = false, gotdot = false, gotexp = false;
    size_t count = 0, intcount = 0, leadzeros = 0, size = 0, zerosize = 0;
    long exp = 0;
    unsigned long tmpexp;
    rawtype val = 0, digit;
    modint *data;
    apstruct *ap;
    char c;

    data = new modint[Memorytreshold];      // Allocate max memory since we don't know input size
    ap = new apstruct (sign, exp, INFINITE, Memorytreshold, data, MEMORY);

    // Whitespace

    str >> ws;

    while ((c = str.get ()) != EOF)
    {
        // Sign

        if (c == '-' && !gotsign)
        {
            sign = -sign;
            gotsign = true;

            continue;
        }
        else if (c == '+' && !gotsign)
        {
            gotsign = true;

            continue;
        }

        // Dot

        if (c == '.' && !gotdot)
        {
            gotdot = true;
            intcount = count;

            if (!size && nonzero)
            {
                // Optimized one base unit integer part
                data[t++] = val;
                size++;
                r = 0;
                val = 0;
            }

            continue;
        }

        // Exponent, always in base 10

        if ((c == Expchar[0] || c == Expchar[1]) && !gotexp)
        {
            gotexp = true;
            if (!(str >> dec >> exp))
                size = 0;               // Error inputting exponent

            break;
        }

        // Valid digit?

        digit = valuetable[(unsigned char) c];
        if (digit >= Basedigit || digit == INVALID)
        {
            str.putback (c);

            break;            // Stop at first invalid character
        }

        // Integer or fractional part

        if (c == chartable[0])
        {
            if (!nonzero)
            {
                leadzeros++;

                if (!gotdot)
                {
                    count++;
                    continue;           // Skip leading zeros before decimal point
                }
            }
        }
        else
        {
            nonzero = true;
        }

        // Get the actual data

        val = val * Basedigit + valuetable[(unsigned char) c];
        if (++r == Basedigits)
        {
            if (nonzero)
            {
                buffercheck (&t, ap, &data);

                data[t++] = val;
                size++;

                if (val == 0)
                    zerosize++;
                else
                    zerosize = 0;
            }

            r = 0;
            val = 0;
        }

        count++;
    }

    if (!count)
    {
        delete ap;                              // No valid input
        str.clear (ios::failbit);
        return str;
    }

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

    if (!gotdot)
        intcount = count;

    // Last base unit

    if (r && val != 0)
    {
        buffercheck (&t, ap, &data);

        if (size || intcount <= leadzeros)      // Not first integer part base unit
            for (; r < Basedigits; r++)
                val *= Basedigit;

        data[t++] = val;
        size++;
        zerosize = 0;
    }

    if (ap->location != MEMORY)
    {
        // Flush buffer to disk
        fstream &fs = ap->openstream ();
        fs.seekp (0, ios::end);
        fs.write ((char *) data, sizeof (modint) * t);
        assert (fs.good ());
        ap->closestream ();

        // Deallocate disk write buffer
        delete[] data;

        // Size is now > Memorytreshold
        // Set the size to the actual size including trailing zeros
        ap->size = size;
    }

    // Align exponent

    t = Maxlong / Basedigits;
    l = t * Basedigits;

    if (intcount <= leadzeros)
    {
        // Fractional part only
        tmpexp = (leadzeros - intcount) / Basedigits;
        e = (exp + l) % Basedigits;
        exp = (exp + l) / Basedigits - t - tmpexp;
    }
    else if (intcount < leadzeros + Basedigits)
    {
        // One integer base unit
        tmpexp = 1;
        e = (exp + l) % Basedigits;
        exp = (exp + l) / Basedigits - t + tmpexp;
    }
    else
    {
        // More than 1 integer base unit
        tmpexp = intcount - leadzeros;
        e = tmpexp % Basedigits;
        tmpexp = tmpexp / Basedigits;
        e += (exp + l) % Basedigits;
        exp = (exp + l) / Basedigits - t + tmpexp;
        if (e >= Basedigits)
        {
            e -= Basedigits;
            exp++;
        }
    }

    ap->sign = sign;
    ap->exp = exp;

    // Calculate actual size (ignoring leading and trailing zeros)

    size -= zerosize;
    ap->resize (size);                  // Set actual size, removing trailing zeros

    if (e)
        d = apfloat (ap) * pow (apfloat ((int) Basedigit), e);
    else
        d = apfloat (ap);

    d.prec (d.ap->size * Basedigits);   // Set precision to the same as size

    return str;
}

// Retrieve an apfloat from disk, essentially deletes the file
void apfloat::swapfrom (const char *filename)
{
    *this = apfloat (apswapfrom (filename));
}

// Store an apfloat to disk, deleting the object from memory
void apfloat::swapto (const char *newfilename)
{
    unique ();

    apswapto (ap, newfilename);

    ap = 0;
}