basex.cpp

// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.

//-----------------------------------------------------------------------------
//  Package Title  ratpak
//  File           basex.c
//  Copyright      (C) 1995-97 Microsoft
//  Date           03-14-97
//
//
//  Description
//
//     Contains number routines for internal base computations, these assume
//  internal base is a power of 2.
//
//-----------------------------------------------------------------------------
#include "ratpak.h"
#include <cstring> // for memmove

void _mulnumx(PNUMBER* pa, PNUMBER b);

//----------------------------------------------------------------------------
//
//    FUNCTION: mulnumx
//
//    ARGUMENTS: pointer to a number and a second number, the
//               base is always BASEX.
//
//    RETURN: None, changes first pointer.
//
//    DESCRIPTION: Does the number equivalent of *pa *= b.
//    This is a stub which prevents multiplication by 1, this is a big speed
//    improvement.
//
//----------------------------------------------------------------------------

void mulnumx(_Inout_ PNUMBER* pa, _In_ PNUMBER b)

{
    if (b->cdigit > 1 || b->mant[0] != 1 || b->exp != 0)
    {
        // If b is not one we multiply
        if ((*pa)->cdigit > 1 || (*pa)->mant[0] != 1 || (*pa)->exp != 0)
        {
            // pa and b are both non-one.
            _mulnumx(pa, b);
        }
        else
        {
            // if pa is one and b isn't just copy b. and adjust the sign.
            int32_t sign = (*pa)->sign;
            DUPNUM(*pa, b);
            (*pa)->sign *= sign;
        }
    }
    else
    {
        // B is +/- 1, But we do have to set the sign.
        (*pa)->sign *= b->sign;
    }
}

//----------------------------------------------------------------------------
//
//    FUNCTION: _mulnumx
//
//    ARGUMENTS: pointer to a number and a second number, the
//               base is always BASEX.
//
//    RETURN: None, changes first pointer.
//
//    DESCRIPTION: Does the number equivalent of *pa *= b.
//    Assumes the base is BASEX of both numbers.  This algorithm is the
//    same one you learned in grade school, except the base isn't 10 it's
//    BASEX.
//
//----------------------------------------------------------------------------

void _mulnumx(PNUMBER* pa, PNUMBER b)

{
    PNUMBER c = nullptr;  // c will contain the result.
    PNUMBER a = nullptr;  // a is the dereferenced number pointer from *pa
    MANTTYPE* ptra;       // ptra is a pointer to the mantissa of a.
    MANTTYPE* ptrb;       // ptrb is a pointer to the mantissa of b.
    MANTTYPE* ptrc;       // ptrc is a pointer to the mantissa of c.
    MANTTYPE* ptrcoffset; // ptrcoffset, is the anchor location of the next
                          // single digit multiply partial result.
    int32_t iadigit = 0;  // Index of digit being used in the first number.
    int32_t ibdigit = 0;  // Index of digit being used in the second number.
    MANTTYPE da = 0;      // da is the digit from the fist number.
    TWO_MANTTYPE cy = 0;  // cy is the carry resulting from the addition of
                          // a multiplied row into the result.
    TWO_MANTTYPE mcy = 0; // mcy is the resultant from a single
                          // multiply, AND the carry of that multiply.
    int32_t icdigit = 0;  // Index of digit being calculated in final result.

    a = *pa;

    ibdigit = a->cdigit + b->cdigit - 1;
    createnum(c, ibdigit + 1);
    c->cdigit = ibdigit;
    c->sign = a->sign * b->sign;

    c->exp = a->exp + b->exp;
    ptra = a->mant;
    ptrcoffset = c->mant;

    for (iadigit = a->cdigit; iadigit > 0; iadigit--)
    {
        da = *ptra++;
        ptrb = b->mant;

        // Shift ptrc, and ptrcoffset, one for each digit
        ptrc = ptrcoffset++;

        for (ibdigit = b->cdigit; ibdigit > 0; ibdigit--)
        {
            cy = 0;
            mcy = (uint64_t)da * (*ptrb);
            if (mcy)
            {
                icdigit = 0;
                if (ibdigit == 1 && iadigit == 1)
                {
                    c->cdigit++;
                }
            }

            // If result is nonzero, or while result of carry is nonzero...
            while (mcy || cy)
            {
                // update carry from addition(s) and multiply.
                cy += (TWO_MANTTYPE)ptrc[icdigit] + ((uint32_t)mcy & ((uint32_t)~BASEX));

                // update result digit from
                ptrc[icdigit++] = (MANTTYPE)((uint32_t)cy & ((uint32_t)~BASEX));

                // update carries from
                mcy >>= BASEXPWR;
                cy >>= BASEXPWR;
            }

            ptrb++;
            ptrc++;
        }
    }

    // prevent different kinds of zeros, by stripping leading duplicate zeros.
    // digits are in order of increasing significance.
    while (c->cdigit > 1 && c->mant[c->cdigit - 1] == 0)
    {
        c->cdigit--;
    }

    destroynum(*pa);
    *pa = c;
}
//-----------------------------------------------------------------------------
//
//    FUNCTION: numpowi32x
//
//    ARGUMENTS: root as number power as int32_t
//               number.
//
//    RETURN: None root is changed.
//
//    DESCRIPTION: changes numeric representation of root to
//    root ** power. Assumes base BASEX
//    decomposes the exponent into it's sums of powers of 2, so on average
//    it will take n+n/2 multiplies where n is the highest on bit.
//
//-----------------------------------------------------------------------------

void numpowi32x(_Inout_ PNUMBER* proot, int32_t power)

{
    PNUMBER lret = i32tonum(1, BASEX);

    // Once the power remaining is zero we are done.
    while (power > 0)
    {
        // If this bit in the power decomposition is on, multiply the result
        // by the root number.
        if (power & 1)
        {
            mulnumx(&lret, *proot);
        }

        // multiply the root number by itself to scale for the next bit (i.e.
        // square it.
        mulnumx(proot, *proot);

        // move the next bit of the power into place.
        power >>= 1;
    }
    destroynum(*proot);
    *proot = lret;
}

void _divnumx(PNUMBER* pa, PNUMBER b, int32_t precision);

//----------------------------------------------------------------------------
//
//    FUNCTION: divnumx
//
//    ARGUMENTS: pointer to a number, a second number and precision.
//
//    RETURN: None, changes first pointer.
//
//    DESCRIPTION: Does the number equivalent of *pa /= b.
//    Assumes radix is the internal radix representation.
//    This is a stub which prevents division by 1, this is a big speed
//    improvement.
//
//----------------------------------------------------------------------------

void divnumx(_Inout_ PNUMBER* pa, _In_ PNUMBER b, int32_t precision)

{
    if (b->cdigit > 1 || b->mant[0] != 1 || b->exp != 0)
    {
        // b is not one.
        if ((*pa)->cdigit > 1 || (*pa)->mant[0] != 1 || (*pa)->exp != 0)
        {
            // pa and b are both not one.
            _divnumx(pa, b, precision);
        }
        else
        {
            // if pa is one and b is not one, just copy b, and adjust the sign.
            int32_t sign = (*pa)->sign;
            DUPNUM(*pa, b);
            (*pa)->sign *= sign;
        }
    }
    else
    {
        // b is one so don't divide, but set the sign.
        (*pa)->sign *= b->sign;
    }
}

//----------------------------------------------------------------------------
//
//    FUNCTION: _divnumx
//
//    ARGUMENTS: pointer to a number, a second number and precision.
//
//    RETURN: None, changes first pointer.
//
//    DESCRIPTION: Does the number equivalent of *pa /= b.
//    Assumes radix is the internal radix representation.
//
//----------------------------------------------------------------------------

void _divnumx(PNUMBER* pa, PNUMBER b, int32_t precision)

{
    PNUMBER a = nullptr;       // a is the dereferenced number pointer from *pa
    PNUMBER c = nullptr;       // c will contain the result.
    PNUMBER lasttmp = nullptr; // lasttmp allows a backup when the algorithm
                               // guesses one bit too far.
    PNUMBER tmp = nullptr;     // current guess being worked on for divide.
    PNUMBER rem = nullptr;     // remainder after applying guess.
    int32_t cdigits;           // count of digits for answer.
    MANTTYPE* ptrc;            // ptrc is a pointer to the mantissa of c.

    int32_t thismax = precision + g_ratio; // set a maximum number of internal digits
                                           // to shoot for in the divide.

    a = *pa;
    if (thismax < a->cdigit)
    {
        // a has more digits than precision specified, bump up digits to shoot
        // for.
        thismax = a->cdigit;
    }

    if (thismax < b->cdigit)
    {
        // b has more digits than precision specified, bump up digits to shoot
        // for.
        thismax = b->cdigit;
    }

    // Create c (the divide answer) and set up exponent and sign.
    createnum(c, thismax + 1);
    c->exp = (a->cdigit + a->exp) - (b->cdigit + b->exp) + 1;
    c->sign = a->sign * b->sign;

    ptrc = c->mant + thismax;
    cdigits = 0;

    DUPNUM(rem, a);
    rem->sign = b->sign;
    rem->exp = b->cdigit + b->exp - rem->cdigit;

    while (cdigits++ < thismax && !zernum(rem))
    {
        int32_t digit = 0;
        *ptrc = 0;
        while (!lessnum(rem, b))
        {
            digit = 1;
            DUPNUM(tmp, b);
            destroynum(lasttmp);
            lasttmp = i32tonum(0, BASEX);
            while (lessnum(tmp, rem))
            {
                destroynum(lasttmp);
                DUPNUM(lasttmp, tmp);
                addnum(&tmp, tmp, BASEX);
                digit *= 2;
            }
            if (lessnum(rem, tmp))
            {
                // too far, back up...
                destroynum(tmp);
                digit /= 2;
                tmp = lasttmp;
                lasttmp = nullptr;
            }

            tmp->sign *= -1;
            addnum(&rem, tmp, BASEX);
            destroynum(tmp);
            destroynum(lasttmp);
            *ptrc |= digit;
        }
        rem->exp++;
        ptrc--;
    }
    cdigits--;
    if (c->mant != ++ptrc)
    {
        memmove(c->mant, ptrc, (int)(cdigits * sizeof(MANTTYPE)));
    }

    if (!cdigits)
    {
        // A zero, make sure no weird exponents creep in
        c->exp = 0;
        c->cdigit = 1;
    }
    else
    {
        c->cdigit = cdigits;
        c->exp -= cdigits;
        // prevent different kinds of zeros, by stripping leading duplicate
        // zeros. digits are in order of increasing significance.
        while (c->cdigit > 1 && c->mant[c->cdigit - 1] == 0)
        {
            c->cdigit--;
        }
    }

    destroynum(rem);

    destroynum(*pa);
    *pa = c;
}