squilu/mpdecimal/basearith.c

/*
 * Copyright (c) 2008-2010 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 "mpdecimal.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include "constants.h"
#include "memory.h"
#include "typearith.h"
#include "basearith.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.
 */
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 = 0;
	mpd_size_t i;

	/* add v to u */
	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_calloc(nplusm+1, sizeof *u)) == NULL) {
			return -1;
		}
	}
	if (n >= MPD_MINALLOC_MAX) {
		if ((v = mpd_calloc(n+1, sizeof *v)) == NULL) {
			mpd_free(u);
			return -1;
		}
	}

	_mpd_shortmul(u, uconst, nplusm, d);
	_mpd_shortmul(v, vconst, n, d);

	/* D2: loop */
	rhat = 0;
	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;
}

/* Leftshift of src by shift digits; src may equal dest. */
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) {
			dest[n--] = h;
		}
		for (; m != MPD_SIZE_MAX; m--,n--) {
			_mpd_divmod_pow10(&h, &l, src[m], MPD_RDIGITS-r);
			dest[n] = ph * lprev + h;
			lprev = l;
		}
		dest[q] = ph * lprev;
	}
	else {
		while (--m != MPD_SIZE_MAX) {
			dest[m+q] = src[m];
		}
	}

	mpd_uint_zero(dest, q);
}

/* Rightshift of src by shift digits; src may equal dest. */
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);
		}
		h = hprev;
		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;
		}
		if (hprev != 0) {
			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.
 */
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 = 0;
	mpd_size_t i;

	/* add v to u */
	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) */
void
_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);
	}
	w[i] = 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;
}