/*

 * Copyright (c) 2008-2024 Stefan Krah. All rights reserved.

 *

 * Redistribution and use in source and binary forms, with or without

 * modification, are permitted provided that the following conditions

 * are met:

 *

 * 1. Redistributions of source code must retain the above copyright

 *    notice, this list of conditions and the following disclaimer.

 * 2. Redistributions in binary form must reproduce the above copyright

 *    notice, this list of conditions and the following disclaimer in the

 *    documentation and/or other materials provided with the distribution.

 *

 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND

 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE

 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL

 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS

 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)

 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT

 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY

 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF

 * SUCH DAMAGE.

 */

#include <assert.h>

#include <stdio.h>

#include "basearith.h"

#include "constants.h"

#include "mpdecimal.h"

#include "typearith.h"

/*********************************************************************/

/*                   Calculations in base MPD_RADIX                  */

/*********************************************************************/

/*

 * Knuth, TAOCP, Volume 2, 4.3.1:

 *     w := sum of u (len m) and v (len n)

 *     n > 0 and m >= n

 * The calling function has to handle a possible final carry.

 */

mpd_uint_t

_mpd_baseadd(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,

             mpd_size_t m, mpd_size_t n)

{

    mpd_uint_t s;

    mpd_uint_t carry = 0;

    mpd_size_t i;

    assert(n > 0 && m >= n);

    /* add n members of u and v */

    for (i = 0; i < n; i++) {

        s = u[i] + (v[i] + carry);

        carry = (s < u[i]) | (s >= MPD_RADIX);

        w[i] = carry ? s-MPD_RADIX : s;

    }

    /* if there is a carry, propagate it */

    for (; carry && i < m; i++) {

        s = u[i] + carry;

        carry = (s == MPD_RADIX);

        w[i] = carry ? 0 : s;

    }

    /* copy the rest of u */

    for (; i < m; i++) {

        w[i] = u[i];

    }

    return carry;

}

/*

 * Add the contents of u to w. Carries are propagated further. The caller

 * has to make sure that w is big enough.

 */

void

_mpd_baseaddto(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n)

{

    mpd_uint_t s;

    mpd_uint_t carry = 0;

    mpd_size_t i;

    if (n == 0) return;

    /* add n members of u to w */

    for (i = 0; i < n; i++) {

        s = w[i] + (u[i] + carry);

        carry = (s < w[i]) | (s >= MPD_RADIX);

        w[i] = carry ? s-MPD_RADIX : s;

    }

    /* if there is a carry, propagate it */

    for (; carry; i++) {

        s = w[i] + carry;

        carry = (s == MPD_RADIX);

        w[i] = carry ? 0 : s;

    }

}

/*

 * Add v to w (len m). The calling function has to handle a possible

 * final carry. Assumption: m > 0.

 */

mpd_uint_t

_mpd_shortadd(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v)

{

    mpd_uint_t s;

    mpd_uint_t carry;

    mpd_size_t i;

    assert(m > 0);

    /* add v to w */

    s = w[0] + v;

    carry = (s < v) | (s >= MPD_RADIX);

    w[0] = carry ? s-MPD_RADIX : s;

    /* if there is a carry, propagate it */

    for (i = 1; carry && i < m; i++) {

        s = w[i] + carry;

        carry = (s == MPD_RADIX);

        w[i] = carry ? 0 : s;

    }

    return carry;

}

/* Increment u. The calling function has to handle a possible carry. */

mpd_uint_t

_mpd_baseincr(mpd_uint_t *u, mpd_size_t n)

{

    mpd_uint_t s;

    mpd_uint_t carry = 1;

    mpd_size_t i;

    assert(n > 0);

    /* if there is a carry, propagate it */

    for (i = 0; carry && i < n; i++) {

        s = u[i] + carry;

        carry = (s == MPD_RADIX);

        u[i] = carry ? 0 : s;

    }

    return carry;

}

/*

 * Knuth, TAOCP, Volume 2, 4.3.1:

 *     w := difference of u (len m) and v (len n).

 *     number in u >= number in v;

 */

void

_mpd_basesub(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,

             mpd_size_t m, mpd_size_t n)

{

    mpd_uint_t d;

    mpd_uint_t borrow = 0;

    mpd_size_t i;

    assert(m > 0 && n > 0);

    /* subtract n members of v from u */

    for (i = 0; i < n; i++) {

        d = u[i] - (v[i] + borrow);

        borrow = (u[i] < d);

        w[i] = borrow ? d + MPD_RADIX : d;

    }

    /* if there is a borrow, propagate it */

    for (; borrow && i < m; i++) {

        d = u[i] - borrow;

        borrow = (u[i] == 0);

        w[i] = borrow ? MPD_RADIX-1 : d;

    }

    /* copy the rest of u */

    for (; i < m; i++) {

        w[i] = u[i];

    }

}

/*

 * Subtract the contents of u from w. w is larger than u. Borrows are

 * propagated further, but eventually w can absorb the final borrow.

 */

void

_mpd_basesubfrom(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n)

{

    mpd_uint_t d;

    mpd_uint_t borrow = 0;

    mpd_size_t i;

    if (n == 0) return;

    /* subtract n members of u from w */

    for (i = 0; i < n; i++) {

        d = w[i] - (u[i] + borrow);

        borrow = (w[i] < d);

        w[i] = borrow ? d + MPD_RADIX : d;

    }

    /* if there is a borrow, propagate it */

    for (; borrow; i++) {

        d = w[i] - borrow;

        borrow = (w[i] == 0);

        w[i] = borrow ? MPD_RADIX-1 : d;

    }

}

/* w := product of u (len n) and v (single word) */

void

_mpd_shortmul(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)

{

    mpd_uint_t hi, lo;

    mpd_uint_t carry = 0;

    mpd_size_t i;

    assert(n > 0);

    for (i=0; i < n; i++) {

        _mpd_mul_words(&hi, &lo, u[i], v);

        lo = carry + lo;

        if (lo < carry) hi++;

        _mpd_div_words_r(&carry, &w[i], hi, lo);

    }

    w[i] = carry;

}

/*

 * Knuth, TAOCP, Volume 2, 4.3.1:

 *     w := product of u (len m) and v (len n)

 *     w must be initialized to zero

 */

void

_mpd_basemul(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,

             mpd_size_t m, mpd_size_t n)

{

    mpd_uint_t hi, lo;

    mpd_uint_t carry;

    mpd_size_t i, j;

    assert(m > 0 && n > 0);

    for (j=0; j < n; j++) {

        carry = 0;

        for (i=0; i < m; i++) {

            _mpd_mul_words(&hi, &lo, u[i], v[j]);

            lo = w[i+j] + lo;

            if (lo < w[i+j]) hi++;

            lo = carry + lo;

            if (lo < carry) hi++;

            _mpd_div_words_r(&carry, &w[i+j], hi, lo);

        }

        w[j+m] = carry;

    }

}

/*

 * Knuth, TAOCP Volume 2, 4.3.1, exercise 16:

 *     w := quotient of u (len n) divided by a single word v

 */

mpd_uint_t

_mpd_shortdiv(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)

{

    mpd_uint_t hi, lo;

    mpd_uint_t rem = 0;

    mpd_size_t i;

    assert(n > 0);

    for (i=n-1; i != MPD_SIZE_MAX; i--) {

        _mpd_mul_words(&hi, &lo, rem, MPD_RADIX);

        lo = u[i] + lo;

        if (lo < u[i]) hi++;

        _mpd_div_words(&w[i], &rem, hi, lo, v);

    }

    return rem;

}

/*

 * Knuth, TAOCP Volume 2, 4.3.1:

 *     q, r := quotient and remainder of uconst (len nplusm)

 *             divided by vconst (len n)

 *     nplusm >= n

 *

 * If r is not NULL, r will contain the remainder. If r is NULL, the

 * return value indicates if there is a remainder: 1 for true, 0 for

 * false.  A return value of -1 indicates an error.

 */

int

_mpd_basedivmod(mpd_uint_t *q, mpd_uint_t *r,

                const mpd_uint_t *uconst, const mpd_uint_t *vconst,

                mpd_size_t nplusm, mpd_size_t n)

{

    mpd_uint_t ustatic[MPD_MINALLOC_MAX];

    mpd_uint_t vstatic[MPD_MINALLOC_MAX];

    mpd_uint_t *u = ustatic;

    mpd_uint_t *v = vstatic;

    mpd_uint_t d, qhat, rhat, w2[2];

    mpd_uint_t hi, lo, x;

    mpd_uint_t carry;

    mpd_size_t i, j, m;

    int retval = 0;

    assert(n > 1 && nplusm >= n);

    m = sub_size_t(nplusm, n);

    /* D1: normalize */

    d = MPD_RADIX / (vconst[n-1] + 1);

    if (nplusm >= MPD_MINALLOC_MAX) {

        if ((u = mpd_alloc(nplusm+1, sizeof *u)) == NULL) {

            return -1;

        }

    }

    if (n >= MPD_MINALLOC_MAX) {

        if ((v = mpd_alloc(n+1, sizeof *v)) == NULL) {

            mpd_free(u);

            return -1;

        }

    }

    _mpd_shortmul(u, uconst, nplusm, d);

    _mpd_shortmul(v, vconst, n, d);

    /* D2: loop */

    for (j=m; j != MPD_SIZE_MAX; j--) {

        /* D3: calculate qhat and rhat */

        rhat = _mpd_shortdiv(w2, u+j+n-1, 2, v[n-1]);

        qhat = w2[1] * MPD_RADIX + w2[0];

        while (1) {

            if (qhat < MPD_RADIX) {

                _mpd_singlemul(w2, qhat, v[n-2]);

                if (w2[1] <= rhat) {

                    if (w2[1] != rhat || w2[0] <= u[j+n-2]) {

                        break;

                    }

                }

            }

            qhat -= 1;

            rhat += v[n-1];

            if (rhat < v[n-1] || rhat >= MPD_RADIX) {

                break;

            }

        }

        /* D4: multiply and subtract */

        carry = 0;

        for (i=0; i <= n; i++) {

            _mpd_mul_words(&hi, &lo, qhat, v[i]);

            lo = carry + lo;

            if (lo < carry) hi++;

            _mpd_div_words_r(&hi, &lo, hi, lo);

            x = u[i+j] - lo;

            carry = (u[i+j] < x);

            u[i+j] = carry ? x+MPD_RADIX : x;

            carry += hi;

        }

        q[j] = qhat;

        /* D5: test remainder */

        if (carry) {

            q[j] -= 1;

            /* D6: add back */

            (void)_mpd_baseadd(u+j, u+j, v, n+1, n);

        }

    }

    /* D8: unnormalize */

    if (r != NULL) {

        _mpd_shortdiv(r, u, n, d);

        /* we are not interested in the return value here */

        retval = 0;

    }

    else {

        retval = !_mpd_isallzero(u, n);

    }

if (u != ustatic) mpd_free(u);

if (v != vstatic) mpd_free(v);

return retval;

}

/*

 * Left shift of src by 'shift' digits; src may equal dest.

 *

 *  dest := area of n mpd_uint_t with space for srcdigits+shift digits.

 *  src  := coefficient with length m.

 *

 * The case splits in the function are non-obvious. The following

 * equations might help:

 *

 *  Let msdigits denote the number of digits in the most significant

 *  word of src. Then 1 <= msdigits <= rdigits.

 *

 *   1) shift = q * rdigits + r

 *   2) srcdigits = qsrc * rdigits + msdigits

 *   3) destdigits = shift + srcdigits

 *                 = q * rdigits + r + qsrc * rdigits + msdigits

 *                 = q * rdigits + (qsrc * rdigits + (r + msdigits))

 *

 *  The result has q zero words, followed by the coefficient that is left-

 *  shifted by r.  The case r == 0 is trivial. For r > 0, it is important

 *  to keep in mind that we always read m source words, but write m+1

 *  destination words if r + msdigits > rdigits, m words otherwise.

 */

void

_mpd_baseshiftl(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t n, mpd_size_t m,

                mpd_size_t shift)

{

#if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)

    /* spurious uninitialized warnings */

    mpd_uint_t l=l, lprev=lprev, h=h;

#else

    mpd_uint_t l, lprev, h;

#endif

    mpd_uint_t q, r;

    mpd_uint_t ph;

    assert(m > 0 && n >= m);

    _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);

    if (r != 0) {

        ph = mpd_pow10[r];

        --m; --n;

        _mpd_divmod_pow10(&h, &lprev, src[m--], MPD_RDIGITS-r);

        if (h != 0) { /* r + msdigits > rdigits <==> h != 0 */

            dest[n--] = h;

        }

        /* write m-1 shifted words */

        for (; m != MPD_SIZE_MAX; m--,n--) {

            _mpd_divmod_pow10(&h, &l, src[m], MPD_RDIGITS-r);

            dest[n] = ph * lprev + h;

            lprev = l;

        }

        /* write least significant word */

        dest[q] = ph * lprev;

    }

    else {

        while (--m != MPD_SIZE_MAX) {

            dest[m+q] = src[m];

        }

    }

    mpd_uint_zero(dest, q);

}

/*

 * Right shift of src by 'shift' digits; src may equal dest.

 * Assumption: srcdigits-shift > 0.

 *

 *  dest := area with space for srcdigits-shift digits.

 *  src  := coefficient with length 'slen'.

 *

 * The case splits in the function rely on the following equations:

 *

 *  Let msdigits denote the number of digits in the most significant

 *  word of src. Then 1 <= msdigits <= rdigits.

 *

 *  1) shift = q * rdigits + r

 *  2) srcdigits = qsrc * rdigits + msdigits

 *  3) destdigits = srcdigits - shift

 *                = qsrc * rdigits + msdigits - (q * rdigits + r)

 *                = (qsrc - q) * rdigits + msdigits - r

 *

 * Since destdigits > 0 and 1 <= msdigits <= rdigits:

 *

 *  4) qsrc >= q

 *  5) qsrc == q  ==>  msdigits > r

 *

 * The result has slen-q words if msdigits > r, slen-q-1 words otherwise.

 */

mpd_uint_t

_mpd_baseshiftr(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t slen,

                mpd_size_t shift)

{

#if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)

    /* spurious uninitialized warnings */

    mpd_uint_t l=l, h=h, hprev=hprev; /* low, high, previous high */

#else

    mpd_uint_t l, h, hprev; /* low, high, previous high */

#endif

    mpd_uint_t rnd, rest;   /* rounding digit, rest */

    mpd_uint_t q, r;

    mpd_size_t i, j;

    mpd_uint_t ph;

    assert(slen > 0);

    _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);

    rnd = rest = 0;

    if (r != 0) {

        ph = mpd_pow10[MPD_RDIGITS-r];

        _mpd_divmod_pow10(&hprev, &rest, src[q], r);

        _mpd_divmod_pow10(&rnd, &rest, rest, r-1);

        if (rest == 0 && q > 0) {

            rest = !_mpd_isallzero(src, q);

        }

        /* write slen-q-1 words */

        for (j=0,i=q+1; i<slen; i++,j++) {

            _mpd_divmod_pow10(&h, &l, src[i], r);

            dest[j] = ph * l + hprev;

            hprev = h;

        }

        /* write most significant word */

        if (hprev != 0) { /* always the case if slen==q-1 */

            dest[j] = hprev;

        }

    }

    else {

        if (q > 0) {

            _mpd_divmod_pow10(&rnd, &rest, src[q-1], MPD_RDIGITS-1);

            /* is there any non-zero digit below rnd? */

            if (rest == 0) rest = !_mpd_isallzero(src, q-1);

        }

        for (j = 0; j < slen-q; j++) {

            dest[j] = src[q+j];

        }

    }

    /* 0-4  ==> rnd+rest < 0.5   */

    /* 5    ==> rnd+rest == 0.5  */

    /* 6-9  ==> rnd+rest > 0.5   */

    return (rnd == 0 || rnd == 5) ? rnd + !!rest : rnd;

}

/*********************************************************************/

/*                      Calculations in base b                       */

/*********************************************************************/

/*

 * Add v to w (len m). The calling function has to handle a possible

 * final carry. Assumption: m > 0.

 */

mpd_uint_t

_mpd_shortadd_b(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v, mpd_uint_t b)

{

    mpd_uint_t s;

    mpd_uint_t carry;

    mpd_size_t i;

    assert(m > 0);

    /* add v to w */

    s = w[0] + v;

    carry = (s < v) | (s >= b);

    w[0] = carry ? s-b : s;

    /* if there is a carry, propagate it */

    for (i = 1; carry && i < m; i++) {

        s = w[i] + carry;

        carry = (s == b);

        w[i] = carry ? 0 : s;

    }

    return carry;

}

/* w := product of u (len n) and v (single word). Return carry. */

mpd_uint_t

_mpd_shortmul_c(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)

{

    mpd_uint_t hi, lo;

    mpd_uint_t carry = 0;

    mpd_size_t i;

    assert(n > 0);

    for (i=0; i < n; i++) {

        _mpd_mul_words(&hi, &lo, u[i], v);

        lo = carry + lo;

        if (lo < carry) hi++;

        _mpd_div_words_r(&carry, &w[i], hi, lo);

    }

    return carry;

}

/* w := product of u (len n) and v (single word) */

mpd_uint_t

_mpd_shortmul_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,

                mpd_uint_t v, mpd_uint_t b)

{

    mpd_uint_t hi, lo;

    mpd_uint_t carry = 0;

    mpd_size_t i;

    assert(n > 0);

    for (i=0; i < n; i++) {

        _mpd_mul_words(&hi, &lo, u[i], v);

        lo = carry + lo;

        if (lo < carry) hi++;

        _mpd_div_words(&carry, &w[i], hi, lo, b);

    }

    return carry;

}

/*

 * Knuth, TAOCP Volume 2, 4.3.1, exercise 16:

 *     w := quotient of u (len n) divided by a single word v

 */

mpd_uint_t

_mpd_shortdiv_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,

                mpd_uint_t v, mpd_uint_t b)

{

    mpd_uint_t hi, lo;

    mpd_uint_t rem = 0;

    mpd_size_t i;

    assert(n > 0);

    for (i=n-1; i != MPD_SIZE_MAX; i--) {

        _mpd_mul_words(&hi, &lo, rem, b);

        lo = u[i] + lo;

        if (lo < u[i]) hi++;

        _mpd_div_words(&w[i], &rem, hi, lo, v);

    }

    return rem;

}