zmath.c

Calc Software Package for Number Calc

/*
 * zmath - extended precision integral arithmetic primitives
 *
 * Copyright (C) 1999  David I. Bell and Ernest Bowen
 *
 * Primary author:  David I. Bell
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.3 $
 * @(#) $Id: zmath.c,v 29.3 2006/12/15 16:20:04 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/RCS/zmath.c,v $
 *
 * Under source code control:	1990/02/15 01:48:28
 * File existed as early as:	before 1990
 *
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */


#include <stdio.h>
#include "zmath.h"

HALF _zeroval_[] = { 0 };
HALF _oneval_[] = { 1 };
HALF _twoval_[] = { 2 };
HALF _threeval_[] = { 3 };
HALF _fourval_[] = { 4 };
HALF _fiveval_[] = { 5 };
HALF _sixval_[] = { 6 };
HALF _sevenval_[] = { 7 };
HALF _eightval_[] = { 8 };
HALF _nineval_[] = { 9 };
HALF _tenval_[] = { 10 };
HALF _elevenval_[] = { 11 };
HALF _twelveval_[] = { 12 };
HALF _thirteenval_[] = { 13 };
HALF _fourteenval_[] = { 14 };
HALF _fifteenval_[] = { 15 };
HALF _sixteenval_[] = { 16 };
HALF _seventeenval_[] = { 17 };
HALF _eightteenval_[] = { 18 };
HALF _nineteenval_[] = { 19 };
HALF _twentyval_[] = { 20 };
HALF _sqbaseval_[] = { 0, 1 };
HALF _pow4baseval_[] = { 0, 0, 1 };
HALF _pow8baseval_[] = { 0, 0, 0, 0, 1 };

ZVALUE zconst[] = {
    { _zeroval_, 1, 0}, { _oneval_, 1, 0}, { _twoval_, 1, 0},
    { _threeval_, 1, 0}, { _fourval_, 1, 0}, { _fiveval_, 1, 0},
    { _sixval_, 1, 0}, { _sevenval_, 1, 0}, { _eightval_, 1, 0},
    { _nineval_, 1, 0}, { _tenval_, 1, 0}, { _elevenval_, 1, 0},
    { _twelveval_, 1, 0}, { _thirteenval_, 1, 0}, { _fourteenval_, 1, 0},
    { _fifteenval_, 1, 0}, { _sixteenval_, 1, 0}, { _seventeenval_, 1, 0},
    { _eightteenval_, 1, 0}, { _nineteenval_, 1, 0}, { _twentyval_, 1, 0}
};

ZVALUE _zero_ = { _zeroval_, 1, 0};
ZVALUE _one_ = { _oneval_, 1, 0 };
ZVALUE _two_ = { _twoval_, 1, 0 };
ZVALUE _ten_ = { _tenval_, 1, 0 };
ZVALUE _sqbase_ = { _sqbaseval_, 2, 0 };
ZVALUE _pow4base_ = { _pow4baseval_, 4, 0 };
ZVALUE _pow8base_ = { _pow8baseval_, 4, 0 };
ZVALUE _neg_one_ = { _oneval_, 1, 1 };

/*
 * 2^64 as a ZVALUE
 */
#if BASEB == 32
ZVALUE _b32_ = { _sqbaseval_, 2, 0 };
ZVALUE _b64_ = { _pow4baseval_, 3, 0 };
#elif BASEB == 16
ZVALUE _b32_ = { _pow4baseval_, 3, 0 };
ZVALUE _b64_ = { _pow8baseval_, 5, 0 };
#else
	-=@=- BASEB not 16 or 32 -=@=-
#endif


/*
 * highhalf[i] - masks off the upper i bits of a HALF
 * rhighhalf[i] - masks off the upper BASEB-i bits of a HALF
 * lowhalf[i] - masks off the upper i bits of a HALF
 * rlowhalf[i] - masks off the upper BASEB-i bits of a HALF
 * bitmask[i] - (1 << i) for  0 <= i <= BASEB*2
 */
HALF highhalf[BASEB+1] = {
#if BASEB == 32
	0x00000000,
	0x80000000, 0xC0000000, 0xE0000000, 0xF0000000,
	0xF8000000, 0xFC000000, 0xFE000000, 0xFF000000,
	0xFF800000, 0xFFC00000, 0xFFE00000, 0xFFF00000,
	0xFFF80000, 0xFFFC0000, 0xFFFE0000, 0xFFFF0000,
	0xFFFF8000, 0xFFFFC000, 0xFFFFE000, 0xFFFFF000,
	0xFFFFF800, 0xFFFFFC00, 0xFFFFFE00, 0xFFFFFF00,
	0xFFFFFF80, 0xFFFFFFC0, 0xFFFFFFE0, 0xFFFFFFF0,
	0xFFFFFFF8, 0xFFFFFFFC, 0xFFFFFFFE, 0xFFFFFFFF
#elif BASEB == 16
	0x0000,
	0x8000, 0xC000, 0xE000, 0xF000,
	0xF800, 0xFC00, 0xFE00, 0xFF00,
	0xFF80, 0xFFC0, 0xFFE0, 0xFFF0,
	0xFFF8, 0xFFFC, 0xFFFE, 0xFFFF
#else
	-=@=- BASEB not 16 or 32 -=@=-
#endif
};
HALF rhighhalf[BASEB+1] = {
#if BASEB == 32
	0xFFFFFFFF, 0xFFFFFFFE, 0xFFFFFFFC, 0xFFFFFFF8,
	0xFFFFFFF0, 0xFFFFFFE0, 0xFFFFFFC0, 0xFFFFFF80,
	0xFFFFFF00, 0xFFFFFE00, 0xFFFFFC00, 0xFFFFF800,
	0xFFFFF000, 0xFFFFE000, 0xFFFFC000, 0xFFFF8000,
	0xFFFF0000, 0xFFFE0000, 0xFFFC0000, 0xFFF80000,
	0xFFF00000, 0xFFE00000, 0xFFC00000, 0xFF800000,
	0xFF000000, 0xFE000000, 0xFC000000, 0xF8000000,
	0xF0000000, 0xE0000000, 0xC0000000, 0x80000000,
	0x00000000
#elif BASEB == 16
	0xFFFF, 0xFFFE, 0xFFFC, 0xFFF8,
	0xFFF0, 0xFFE0, 0xFFC0, 0xFF80,
	0xFF00, 0xFE00, 0xFC00, 0xF800,
	0xF000, 0xE000, 0xC000, 0x8000,
	0x0000
#else
	-=@=- BASEB not 16 or 32 -=@=-
#endif
};
HALF lowhalf[BASEB+1] = {
	0x0,
	0x1,  0x3, 0x7, 0xF,
	0x1F,  0x3F, 0x7F, 0xFF,
	0x1FF,	0x3FF, 0x7FF, 0xFFF,
	0x1FFF,	 0x3FFF, 0x7FFF, 0xFFFF
#if BASEB == 32
       ,0x1FFFF,  0x3FFFF, 0x7FFFF, 0xFFFFF,
	0x1FFFFF,  0x3FFFFF, 0x7FFFFF, 0xFFFFFF,
	0x1FFFFFF,  0x3FFFFFF, 0x7FFFFFF, 0xFFFFFFF,
	0x1FFFFFFF,  0x3FFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF
#endif
};
HALF rlowhalf[BASEB+1] = {
#if BASEB == 32
	0xFFFFFFFF, 0x7FFFFFFF, 0x3FFFFFFF, 0x1FFFFFFF,
	0xFFFFFFF, 0x7FFFFFF, 0x3FFFFFF, 0x1FFFFFF,
	0xFFFFFF, 0x7FFFFF, 0x3FFFFF, 0x1FFFFF,
	0xFFFFF, 0x7FFFF, 0x3FFFF, 0x1FFFF,
#endif
	0xFFFF, 0x7FFF, 0x3FFF, 0x1FFF,
	0xFFF, 0x7FF, 0x3FF, 0x1FF,
	0xFF, 0x7F, 0x3F, 0x1F,
	0xF, 0x7, 0x3, 0x1,
	0x0
};
HALF bitmask[(2*BASEB)+1] = {
#if BASEB == 32
	0x00000001, 0x00000002, 0x00000004, 0x00000008,
	0x00000010, 0x00000020, 0x00000040, 0x00000080,
	0x00000100, 0x00000200, 0x00000400, 0x00000800,
	0x00001000, 0x00002000, 0x00004000, 0x00008000,
	0x00010000, 0x00020000, 0x00040000, 0x00080000,
	0x00100000, 0x00200000, 0x00400000, 0x00800000,
	0x01000000, 0x02000000, 0x04000000, 0x08000000,
	0x10000000, 0x20000000, 0x40000000, 0x80000000,
	0x00000001, 0x00000002, 0x00000004, 0x00000008,
	0x00000010, 0x00000020, 0x00000040, 0x00000080,
	0x00000100, 0x00000200, 0x00000400, 0x00000800,
	0x00001000, 0x00002000, 0x00004000, 0x00008000,
	0x00010000, 0x00020000, 0x00040000, 0x00080000,
	0x00100000, 0x00200000, 0x00400000, 0x00800000,
	0x01000000, 0x02000000, 0x04000000, 0x08000000,
	0x10000000, 0x20000000, 0x40000000, 0x80000000,
	0x00000001
#elif BASEB == 16
	0x0001, 0x0002, 0x0004, 0x0008,
	0x0010, 0x0020, 0x0040, 0x0080,
	0x0100, 0x0200, 0x0400, 0x0800,
	0x1000, 0x2000, 0x4000, 0x8000,
	0x0001, 0x0002, 0x0004, 0x0008,
	0x0010, 0x0020, 0x0040, 0x0080,
	0x0100, 0x0200, 0x0400, 0x0800,
	0x1000, 0x2000, 0x4000, 0x8000,
	0x0001
#else
	-=@=- BASEB not 16 or 32 -=@=-
#endif
};		/* actual rotation thru 8 cycles */

BOOL _math_abort_;		/* nonzero to abort calculations */


/*
 * popcnt - popcnt[x] number of 1 bits in 0 <= x < 256
 */
char popcnt[256] = {
    0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
};



#ifdef ALLOCTEST
static long nalloc = 0;
static long nfree = 0;
#endif


HALF *
alloc(LEN len)
{
	HALF *hp;

	if (_math_abort_) {
		math_error("Calculation aborted");
		/*NOTREACHED*/
	}
	hp = (HALF *) malloc((len+1) * sizeof(HALF));
	if (hp == 0) {
		math_error("Not enough memory");
		/*NOTREACHED*/
	}
#ifdef ALLOCTEST
	++nalloc;
#endif
	return hp;
}


#ifdef ALLOCTEST
void
freeh(HALF *h)
{
	if ((h != _zeroval_) && (h != _oneval_)) {
		free(h);
		++nfree;
	}
}


void
allocStat(void)
{
	fprintf(stderr, "nalloc: %ld nfree: %ld kept: %ld\n",
		nalloc, nfree, nalloc - nfree);
}
#endif


/*
 * Convert a normal integer to a number.
 */
void
itoz(long i, ZVALUE *res)
{
	long diddle, len;

	res->len = 1;
	res->sign = 0;
	diddle = 0;
	if (i == 0) {
		res->v = _zeroval_;
		return;
	}
	if (i < 0) {
		res->sign = 1;
		i = -i;
		if (i < 0) {	/* fix most negative number */
			diddle = 1;
			i--;
		}
	}
	if (i == 1) {
		res->v = _oneval_;
		return;
	}
	len = 1 + (((FULL) i) >= BASE);
	res->len = (LEN)len;
	res->v = alloc((LEN)len);
	res->v[0] = (HALF) (i + diddle);
	if (len == 2)
		res->v[1] = (HALF) (i / BASE);
}


/*
 * Convert a number to a normal integer, as far as possible.
 * If the number is out of range, the largest number is returned.
 */
long
ztoi(ZVALUE z)
{
	long i;

	if (zgtmaxlong(z)) {
		i = MAXLONG;
		return (z.sign ? -i : i);
	}
	i = ztolong(z);
	return (z.sign ? -i : i);
}


/*
 * Convert a normal unsigned integer to a number.
 */
void
utoz(FULL i, ZVALUE *res)
{
	long len;

	res->len = 1;
	res->sign = 0;
	if (i == 0) {
		res->v = _zeroval_;
		return;
	}
	if (i == 1) {
		res->v = _oneval_;
		return;
	}
	len = 1 + (((FULL) i) >= BASE);
	res->len = (LEN)len;
	res->v = alloc((LEN)len);
	res->v[0] = (HALF)i;
	if (len == 2)
		res->v[1] = (HALF) (i / BASE);
}


/*
 * Convert a normal signed integer to a number.
 */
void
stoz(SFULL i, ZVALUE *res)
{
	long diddle, len;

	res->len = 1;
	res->sign = 0;
	diddle = 0;
	if (i == 0) {
		res->v = _zeroval_;
		return;
	}
	if (i < 0) {
		res->sign = 1;
		i = -i;
		if (i < 0) {	/* fix most negative number */
			diddle = 1;
			i--;
		}
	}
	if (i == 1) {
		res->v = _oneval_;
		return;
	}
	len = 1 + (((FULL) i) >= BASE);
	res->len = (LEN)len;
	res->v = alloc((LEN)len);
	res->v[0] = (HALF) (i + diddle);
	if (len == 2)
		res->v[1] = (HALF) (i / BASE);
}


/*
 * Convert a number to a unsigned normal integer, as far as possible.
 * If the number is out of range, the largest number is returned.
 * The absolute value of z is converted.
 */
FULL
ztou(ZVALUE z)
{
	if (z.len > 2) {
		return MAXUFULL;
	}
	return ztofull(z);
}


/*
 * Convert a number to a signed normal integer, as far as possible.
 *
 * If the number is too large to fit into an integer, than the largest
 * positive or largest negative integer is returned depending on the sign.
 */
SFULL
ztos(ZVALUE z)
{
	FULL val;	/* absolute value of the return value */

	/* negative value processing */
	if (z.sign) {
	    if (z.len > 2) {
		return MINSFULL;
	    }
	    val = ztofull(z);
	    if (val > TOPFULL) {
		return MINSFULL;
	    }
	    return (SFULL)(-val);
	}

	/* positive value processing */
	if (z.len > 2) {
	    return (SFULL)MAXFULL;
	}
	val = ztofull(z);
	if (val > MAXFULL) {
	    return (SFULL)MAXFULL;
	}
	return (SFULL)val;
}


/*
 * Make a copy of an integer value
 */
void
zcopy(ZVALUE z, ZVALUE *res)
{
	res->sign = z.sign;
	res->len = z.len;
	if (zisabsleone(z)) {	/* zero or plus or minus one are easy */
		res->v = (z.v[0] ? _oneval_ : _zeroval_);
		return;
	}
	res->v = alloc(z.len);
	zcopyval(z, *res);
}


/*
 * Add together two integers.
 */
void
zadd(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
	ZVALUE dest;
	HALF *p1, *p2, *pd;
	long len;
	FULL carry;
	SIUNION sival;

	if (z1.sign && !z2.sign) {
		z1.sign = 0;
		zsub(z2, z1, res);
		return;
	}
	if (z2.sign && !z1.sign) {
		z2.sign = 0;
		zsub(z1, z2, res);
		return;
	}
	if (z2.len > z1.len) {
		pd = z1.v; z1.v = z2.v; z2.v = pd;
		len = z1.len; z1.len = z2.len; z2.len = (LEN)len;
	}
	dest.len = z1.len + 1;
	dest.v = alloc(dest.len);
	dest.sign = z1.sign;
	carry = 0;
	pd = dest.v;
	p1 = z1.v;
	p2 = z2.v;
	len = z2.len;
	while (len--) {
		sival.ivalue = ((FULL) *p1++) + ((FULL) *p2++) + carry;
		/* ignore Saber-C warning #112 - get ushort from uint */
		/*	  ok to ignore on name zadd`sival */
		*pd++ = sival.silow;
		carry = sival.sihigh;
	}
	len = z1.len - z2.len;
	while (len--) {
		sival.ivalue = ((FULL) *p1++) + carry;
		*pd++ = sival.silow;
		carry = sival.sihigh;
	}
	*pd = (HALF)carry;
	zquicktrim(dest);
	*res = dest;
}


/*
 * Subtract two integers.
 */
void
zsub(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
	register HALF *h1, *h2, *hd;
	long len1, len2;
	FULL carry;
	SIUNION sival;
	ZVALUE dest;

	if (z1.sign != z2.sign) {
		z2.sign = z1.sign;
		zadd(z1, z2, res);
		return;
	}
	len1 = z1.len;
	len2 = z2.len;
	if (len1 == len2) {
		h1 = z1.v + len1;
		h2 = z2.v + len2;
		while ((len1 > 0) && ((FULL)*--h1 == (FULL)*--h2)) {
			len1--;
		}
		if (len1 == 0) {
			*res = _zero_;
			return;
		}
		len2 = len1;
		carry = ((FULL)*h1 < (FULL)*h2);
	} else {
		carry = (len1 < len2);
	}
	dest.sign = z1.sign;
	h1 = z1.v;
	h2 = z2.v;
	if (carry) {
		carry = len1;
		len1 = len2;
		len2 = (long)carry;
		h1 = z2.v;
		h2 = z1.v;
		dest.sign = !dest.sign;
	}
	hd = alloc((LEN)len1);
	dest.v = hd;
	dest.len = (LEN)len1;
	len1 -= len2;
	carry = 0;
	while (--len2 >= 0) {
		/* ignore Saber-C warning #112 - get ushort from uint */
		/*	  ok to ignore on name zsub`sival */
		sival.ivalue = (BASE1 - ((FULL) *h1++)) + *h2++ + carry;
		*hd++ = (HALF)(BASE1 - sival.silow);
		carry = sival.sihigh;
	}
	while (--len1 >= 0) {
		sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
		*hd++ = (HALF)(BASE1 - sival.silow);
		carry = sival.sihigh;
	}
	if (hd[-1] == 0)
		ztrim(&dest);
	*res = dest;
}


/*
 * Multiply an integer by a small number.
 */
void
zmuli(ZVALUE z, long n, ZVALUE *res)
{
	register HALF *h1, *sd;
	FULL low;
	FULL high;
	FULL carry;
	long len;
	SIUNION sival;
	ZVALUE dest;

	if ((n == 0) || ziszero(z)) {
		*res = _zero_;
		return;
	}
	if (n < 0) {
		n = -n;
		z.sign = !z.sign;
	}
	if (n == 1) {
		zcopy(z, res);
		return;
	}
#if LONG_BITS > BASEB
	low = ((FULL) n) & BASE1;
	high = ((FULL) n) >> BASEB;
#else
	low = (FULL)n;
	high = 0;
#endif
	dest.len = z.len + 2;
	dest.v = alloc(dest.len);
	dest.sign = z.sign;
	/*
	 * Multiply by the low digit.
	 */
	h1 = z.v;
	sd = dest.v;
	len = z.len;
	carry = 0;
	while (len--) {
		/* ignore Saber-C warning #112 - get ushort from uint */
		/*	  ok to ignore on name zmuli`sival */
		sival.ivalue = ((FULL) *h1++) * low + carry;
		*sd++ = sival.silow;
		carry = sival.sihigh;
	}
	*sd = (HALF)carry;
	/*
	 * If there was only one digit, then we are all done except
	 * for trimming the number if there was no last carry.
	 */
	if (high == 0) {
		dest.len--;
		if (carry == 0)
			dest.len--;
		*res = dest;
		return;
	}
	/*
	 * Need to multiply by the high digit and add it into the
	 * previous value.  Clear the final word of rubbish first.
	 */
	*(++sd) = 0;
	h1 = z.v;
	sd = dest.v + 1;
	len = z.len;
	carry = 0;
	while (len--) {
		sival.ivalue = ((FULL) *h1++) * high + ((FULL) *sd) + carry;
		*sd++ = sival.silow;
		carry = sival.sihigh;
	}
	*sd = (HALF)carry;
	zquicktrim(dest);
	*res = dest;
}


/*
 * Divide two numbers by their greatest common divisor.
 * This is useful for reducing the numerator and denominator of
 * a fraction to its lowest terms.