mirror of
https://github.com/lcn2/calc.git
synced 2025-08-16 01:03:29 +03:00
db77e29a23
Converted all ASCII tabs to ASCII spaces using a 8 character tab stop, for all files, except for all Makefiles (plus rpm.mk). The `git diff -w` reports no changes.
2166 lines
56 KiB
C
2166 lines
56 KiB
C
/*
|
|
* zmath - extended precision integral arithmetic primitives
|
|
*
|
|
* Copyright (C) 1999-2007,2021-2023 David I. Bell, Landon Curt Noll 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.
|
|
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
*
|
|
* 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 "int.h"
|
|
#include "alloc.h"
|
|
#include "zmath.h"
|
|
|
|
|
|
#include "errtbl.h"
|
|
#include "banned.h" /* include after system header <> includes */
|
|
|
|
|
|
HALF _zeroval_[] = { 0 };
|
|
ZVALUE _zero_ = { _zeroval_, 1, 0};
|
|
|
|
HALF _oneval_[] = { 1 };
|
|
ZVALUE _one_ = { _oneval_, 1, 0 };
|
|
ZVALUE _neg_one_ = { _oneval_, 1, 1 };
|
|
|
|
HALF _twoval_[] = { 2 };
|
|
ZVALUE _two_ = { _twoval_, 1, 0 };
|
|
|
|
HALF _tenval_[] = { 10 };
|
|
ZVALUE _ten_ = { _tenval_, 1, 0 };
|
|
|
|
HALF _sqbaseval_[] = { 0, 1 };
|
|
ZVALUE _sqbase_ = { _sqbaseval_, 2, 0 };
|
|
|
|
HALF _pow4baseval_[] = { 0, 0, 1 };
|
|
ZVALUE _pow4base_ = { _pow4baseval_, 4, 0 };
|
|
|
|
HALF _pow8baseval_[] = { 0, 0, 0, 0, 1 };
|
|
ZVALUE _pow8base_ = { _pow8baseval_, 4, 0 };
|
|
|
|
/*
|
|
* 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
|
|
|
|
/*
|
|
* ZVALUE - values that should not be freed
|
|
*/
|
|
HALF *half_tbl[] = {
|
|
_zeroval_,
|
|
_oneval_,
|
|
_twoval_,
|
|
_tenval_,
|
|
_sqbaseval_,
|
|
_pow4baseval_,
|
|
_pow8baseval_,
|
|
NULL /* must be last */
|
|
};
|
|
|
|
|
|
/*
|
|
* 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 lower i bits of a HALF
|
|
* rlowhalf[i] - masks off the lower BASEB-i bits of a HALF
|
|
* bitmask[i] - (1 << i) for 0 <= i <= BASEB*2
|
|
*
|
|
* NOTE: In all cases 0 <= i <= BASEB
|
|
*/
|
|
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");
|
|
not_reached();
|
|
}
|
|
hp = (HALF *) malloc((len+1) * sizeof(HALF));
|
|
if (hp == 0) {
|
|
math_error("Not enough memory");
|
|
not_reached();
|
|
}
|
|
#ifdef ALLOCTEST
|
|
++nalloc;
|
|
#endif
|
|
return hp;
|
|
}
|
|
|
|
|
|
/*
|
|
* is_const - determine if a HALF array is an pre-allocated array
|
|
*
|
|
* given:
|
|
* h pointer to the beginning of the HALF array
|
|
*
|
|
* returns:
|
|
* true - h is found in the half_tbl array
|
|
* false - is is not found in the half_tbl array
|
|
*/
|
|
int
|
|
is_const(HALF *h)
|
|
{
|
|
HALF **h_p; /* half_tbl array pointer */
|
|
|
|
/* firewall */
|
|
if (h == NULL) {
|
|
math_error("%s: h NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
/* search the half_tbl for h */
|
|
for (h_p = &half_tbl[0]; *h_p != NULL; ++h_p) {
|
|
if (h == *h_p) {
|
|
return true; /* found in the half_tbl array */
|
|
}
|
|
}
|
|
|
|
/* not found in the half_tbl array */
|
|
return false;
|
|
}
|
|
|
|
|
|
#ifdef ALLOCTEST
|
|
void
|
|
freeh(HALF *h)
|
|
{
|
|
/* firewall */
|
|
if (h == NULL) {
|
|
math_error("%s: h NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
/* free h is not a constant */
|
|
if (!is_const(h)) {
|
|
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;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
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;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
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;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
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)
|
|
{
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
/* make copy */
|
|
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;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
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;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
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;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
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.
|
|
*/
|
|
void
|
|
zreduce(ZVALUE z1, ZVALUE z2, ZVALUE *z1res, ZVALUE *z2res)
|
|
{
|
|
ZVALUE tmp;
|
|
|
|
/* firewall */
|
|
if (z1res == NULL) {
|
|
math_error("%s: z1res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
if (z2res == NULL) {
|
|
math_error("%s: z2res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (zisabsleone(z1) || zisabsleone(z2))
|
|
tmp = _one_;
|
|
else
|
|
zgcd(z1, z2, &tmp);
|
|
if (zisunit(tmp)) {
|
|
zcopy(z1, z1res);
|
|
zcopy(z2, z2res);
|
|
} else {
|
|
zequo(z1, tmp, z1res);
|
|
zequo(z2, tmp, z2res);
|
|
}
|
|
zfree(tmp);
|
|
}
|
|
|
|
|
|
|
|
/*
|
|
* Compute the quotient and remainder for division of an integer by an
|
|
* integer, rounding criteria determined by rnd. Returns the sign of
|
|
* the remainder.
|
|
*/
|
|
long
|
|
zdiv(ZVALUE z1, ZVALUE z2, ZVALUE *quo, ZVALUE *rem, long rnd)
|
|
{
|
|
register HALF *a, *b;
|
|
HALF s, u;
|
|
HALF *A, *B, *a1, *b0;
|
|
FULL f, g, h, x;
|
|
bool adjust, onebit;
|
|
LEN m, n, len, i, p, j1, j2, k;
|
|
long t, val;
|
|
|
|
/* firewall */
|
|
if (quo == NULL) {
|
|
math_error("%s: quo NULL", __func__);
|
|
not_reached();
|
|
}
|
|
if (rem == NULL) {
|
|
math_error("%s: rem NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (ziszero(z2)) {
|
|
math_error("Division by zero in zdiv");
|
|
not_reached();
|
|
}
|
|
m = z1.len;
|
|
n = z2.len;
|
|
B = z2.v;
|
|
s = 0;
|
|
if (m < n) {
|
|
A = alloc(n + 1);
|
|
memcpy(A, z1.v, m * sizeof(HALF));
|
|
for (i = m; i <= n; i++)
|
|
A[i] = 0;
|
|
a1 = A + n;
|
|
len = 1;
|
|
goto done;
|
|
}
|
|
A = alloc(m + 2);
|
|
memcpy(A, z1.v, m * sizeof(HALF));
|
|
A[m] = 0;
|
|
A[m + 1] = 0;
|
|
len = m - n + 1; /* quotient length will be len or len +/- 1 */
|
|
a1 = A + n; /* start of digits for quotient */
|
|
b0 = B;
|
|
p = n;
|
|
while (!*b0) { /* b0: working start for divisor */
|
|
++b0;
|
|
--p;
|
|
}
|
|
if (p == 1) {
|
|
u = *b0;
|
|
if (u == 1) {
|
|
for (; m >= n; m--)
|
|
A[m] = A[m - 1];
|
|
A[m] = 0;
|
|
goto done;
|
|
}
|
|
f = 0;
|
|
a = A + m;
|
|
i = len;
|
|
while (i--) {
|
|
f = f << BASEB | *--a;
|
|
a[1] = (HALF)(f / u);
|
|
f = f % u;
|
|
}
|
|
*a = (HALF)f;
|
|
m = n;
|
|
goto done;
|
|
}
|
|
f = B[n - 1];
|
|
k = 1;
|
|
while (f >>= 1)
|
|
k++; /* k: number of bits in top divisor digit */
|
|
j1 = BASEB - k;
|
|
j2 = BASEB + j1;
|
|
h = j1 ? ((FULL) B[n - 1] << j1 | B[n - 2] >> k) : B[n-1];
|
|
onebit = (bool)((B[n - 2] >> (k - 1)) & 1);
|
|
m++;
|
|
while (m > n) {
|
|
m--;
|
|
f = (FULL) A[m] << j2 | (FULL) A[m - 1] << j1;
|
|
if (j1) f |= A[m - 2] >> k;
|
|
if (s) f = ~f;
|
|
x = f / h;
|
|
if (x) {
|
|
if (onebit && x > TOPHALF + f % h)
|
|
x--;
|
|
a = A + m - p;
|
|
b = b0;
|
|
u = 0;
|
|
i = p;
|
|
if (s) {
|
|
while (i--) {
|
|
f = (FULL) *a + u + x * *b++;
|
|
*a++ = (HALF) f;
|
|
u = (HALF) (f >> BASEB);
|
|
}
|
|
s = *a + u;
|
|
A[m] = (HALF) (~x + !s);
|
|
} else {
|
|
while (i--) {
|
|
f = (FULL) *a - u - x * *b++;
|
|
*a++ = (HALF) f;
|
|
u = -(HALF)(f >> BASEB);
|
|
}
|
|
s = *a - u;
|
|
A[m] = (HALF)(x + s);
|
|
}
|
|
}
|
|
else
|
|
A[m] = s;
|
|
}
|
|
done: while (m > 0 && A[m - 1] == 0)
|
|
m--;
|
|
if (m == 0 && s == 0) {
|
|
*rem = _zero_;
|
|
val = 0;
|
|
if (a1[len - 1] == 0)
|
|
len--;
|
|
if (len == 0) {
|
|
*quo = _zero_;
|
|
} else {
|
|
quo->len = len;
|
|
quo->v = alloc(len);
|
|
memcpy(quo->v, a1, len * sizeof(HALF));
|
|
quo->sign = z1.sign ^ z2.sign;
|
|
}
|
|
freeh(A);
|
|
return val;
|
|
}
|
|
if (rnd & 8)
|
|
adjust = (((*a1 ^ rnd) & 1) ? true : false);
|
|
else
|
|
adjust = (((rnd & 1) ^ z1.sign ^ z2.sign) ? true : false);
|
|
if (rnd & 2)
|
|
adjust ^= z1.sign ^ z2.sign;
|
|
if (rnd & 4)
|
|
adjust ^= z2.sign;
|
|
if (rnd & 16) { /* round-off case */
|
|
a = A + n;
|
|
b = B + n;
|
|
i = n + 1;
|
|
f = g = 0;
|
|
t = -1;
|
|
if (s) {
|
|
while (--i > 0) {
|
|
g = (FULL) *--a + (*--b >> 1 | f);
|
|
f = *b & 1 ? TOPHALF : 0;
|
|
if (g != BASE1)
|
|
break;
|
|
}
|
|
if (g == BASE && f == 0) {
|
|
while ((--i > 0) && ((*--a | *--b) == 0));
|
|
t = (i > 0);
|
|
} else if (g >= BASE) {
|
|
t = 1;
|
|
}
|
|
} else {
|
|
while (--i > 0) {
|
|
g = (FULL) *--a - (*--b >> 1 | f);
|
|
f = *b & 1 ? TOPHALF : 0;
|
|
if (g != 0)
|
|
break;
|
|
}
|
|
if (g > 0 && g < BASE)
|
|
t = 1;
|
|
else if (g == 0 && f == 0)
|
|
t = 0;
|
|
}
|
|
if (t)
|
|
adjust = (t > 0);
|
|
}
|
|
if (adjust) {
|
|
i = len;
|
|
a = a1;
|
|
while (i > 0 && *a == BASE1) {
|
|
i--;
|
|
*a++ = 0;
|
|
}
|
|
(*a)++;
|
|
if (i == 0)
|
|
len++;
|
|
}
|
|
if (s && adjust) {
|
|
i = 0;
|
|
while (A[i] == 0)
|
|
i++;
|
|
A[i] = -A[i];
|
|
while (++i < n)
|
|
A[i] = ~A[i];
|
|
m = n;
|
|
while (A[m - 1] == 0)
|
|
m--;
|
|
}
|
|
if (!s && adjust) {
|
|
a = A;
|
|
b = B;
|
|
i = n;
|
|
u = 0;
|
|
while (i--) {
|
|
f = (FULL) *b++ - *a - (FULL) u;
|
|
*a++ = (HALF) f;
|
|
u = -(HALF)(f >> BASEB);
|
|
}
|
|
m = n;
|
|
while (A[m - 1] == 0)
|
|
m--;
|
|
}
|
|
if (s && !adjust) {
|
|
a = A;
|
|
b = B;
|
|
i = n;
|
|
f = 0;
|
|
while (i--) {
|
|
f = (FULL) *b++ + *a + (f >> BASEB);
|
|
*a++ = (HALF) f;
|
|
}
|
|
m = n;
|
|
while (A[m-1] == 0)
|
|
m--;
|
|
}
|
|
rem->len = m;
|
|
rem->v = alloc(m);
|
|
memcpy(rem->v, A, m * sizeof(HALF));
|
|
rem->sign = z1.sign ^ adjust;
|
|
val = rem->sign ? -1 : 1;
|
|
if (a1[len - 1] == 0)
|
|
len--;
|
|
if (len == 0) {
|
|
*quo = _zero_;
|
|
} else {
|
|
quo->len = len;
|
|
quo->v = alloc(len);
|
|
memcpy(quo->v, a1, len * sizeof(HALF));
|
|
quo->sign = z1.sign ^ z2.sign;
|
|
}
|
|
freeh(A);
|
|
return val;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute and store at a specified location the integer quotient
|
|
* of two integers, the type of rounding being determined by rnd.
|
|
* Returns the sign of z1/z2 - calculated quotient.
|
|
*/
|
|
long
|
|
zquo(ZVALUE z1, ZVALUE z2, ZVALUE *res, long rnd)
|
|
{
|
|
ZVALUE tmp;
|
|
long val;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
val = zdiv(z1, z2, res, &tmp, rnd);
|
|
if (z2.sign)
|
|
val = -val;
|
|
zfree(tmp);
|
|
return val;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute and store at a specified location the remainder for
|
|
* division of an integer by an integer, the type of rounding
|
|
* used being determined by rnd. Returns the sign of the remainder.
|
|
*/
|
|
long
|
|
zmod(ZVALUE z1, ZVALUE z2, ZVALUE *res, long rnd)
|
|
{
|
|
ZVALUE tmp;
|
|
long val;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
val = zdiv(z1, z2, &tmp, res, rnd);
|
|
zfree(tmp);
|
|
return val;
|
|
}
|
|
|
|
|
|
/*
|
|
* Computes the quotient z1/z2 on the assumption that this is exact.
|
|
* There is no thorough check on the exactness of the division
|
|
* so this should not be called unless z1/z2 is an integer.
|
|
*/
|
|
void
|
|
zequo(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
LEN i, j, k, len, m, n, o, p;
|
|
HALF *a, *a0, *A, *b, *B, u, v, w, x;
|
|
FULL f, g;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (ziszero(z1)) {
|
|
*res = _zero_;
|
|
return;
|
|
}
|
|
if (ziszero(z2)) {
|
|
math_error("Division by zero");
|
|
not_reached();
|
|
}
|
|
if (zisone(z2)) {
|
|
zcopy(z1, res);
|
|
return;
|
|
}
|
|
if (zhighbit(z1) < zhighbit(z2)) {
|
|
math_error("Bad call to zequo");
|
|
not_reached();
|
|
}
|
|
B = z2.v;
|
|
o = 0;
|
|
while (!*B) {
|
|
++B;
|
|
++o;
|
|
}
|
|
m = z1.len - o;
|
|
n = z2.len - o;
|
|
len = m - n + 1; /* Maximum length of quotient */
|
|
v = *B;
|
|
A = alloc(len+1);
|
|
memcpy(A, z1.v + o, len * sizeof(HALF));
|
|
A[len] = 0;
|
|
if (n == 1) {
|
|
if (v > 1) {
|
|
a = A + len;
|
|
i = len;
|
|
f = 0;
|
|
while (i--) {
|
|
f = f << BASEB | *--a;
|
|
*a = (HALF)(f / v);
|
|
f %= v;
|
|
}
|
|
}
|
|
} else {
|
|
k = 0;
|
|
while (!(v & 1)) {
|
|
k++;
|
|
v >>= 1;
|
|
}
|
|
j = BASEB - k;
|
|
if (n > 1 && k > 0)
|
|
v |= B[1] << j;
|
|
u = v - 1;
|
|
w = x = 1;
|
|
while (u) { /* To find w = inverse of v modulo BASE */
|
|
do {
|
|
v <<= 1;
|
|
x <<= 1;
|
|
}
|
|
while (!(u & x));
|
|
u += v;
|
|
w |= x;
|
|
}
|
|
a0 = A;
|
|
p = len;
|
|
while (p > 1) {
|
|
if (!*a0) {
|
|
while (!*++a0 && p > 1)
|
|
p--;
|
|
--a0;
|
|
}
|
|
if (p == 1)
|
|
break;
|
|
x = k ? w * (*a0 >> k | a0[1] << j) : w * *a0;
|
|
g = x;
|
|
if (x) {
|
|
a = a0;
|
|
b = B;
|
|
u = 0;
|
|
i = n;
|
|
if (i > p)
|
|
i = p;
|
|
while (i--) {
|
|
f = (FULL) *a - g * *b++ - (FULL) u;
|
|
*a++ = (HALF)f;
|
|
u = -(HALF)(f >> BASEB);
|
|
}
|
|
if (u && p > n) {
|
|
i = p - n;
|
|
while (u && i--) {
|
|
f = (FULL) *a - u;
|
|
*a++ = (HALF) f;
|
|
u = -(HALF)(f >> BASEB);
|
|
}
|
|
}
|
|
}
|
|
*a0++ = x;
|
|
p--;
|
|
}
|
|
if (k == 0) {
|
|
*a0 = w * *a0;
|
|
} else {
|
|
u = (HALF)(w * *a0) >> k;
|
|
x = (HALF)(((FULL) z1.v[z1.len - 1] << BASEB
|
|
| z1.v[z1.len - 2])
|
|
/((FULL) B[n-1] << BASEB | B[n-2]));
|
|
if ((x ^ u) & 1) x--;
|
|
*a0 = x;
|
|
}
|
|
}
|
|
if (!A[len - 1]) len--;
|
|
res->v = A;
|
|
res->len = len;
|
|
res->sign = z1.sign != z2.sign;
|
|
}
|
|
|
|
|
|
|
|
/*
|
|
* Return the quotient and remainder of an integer divided by a small
|
|
* number. A nonzero remainder is only meaningful when both numbers
|
|
* are positive.
|
|
*/
|
|
long
|
|
zdivi(ZVALUE z, long n, ZVALUE *res)
|
|
{
|
|
HALF *h1, *sd;
|
|
FULL val;
|
|
HALF divval[2];
|
|
ZVALUE div;
|
|
ZVALUE dest;
|
|
LEN len;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (n == 0) {
|
|
math_error("Division by zero");
|
|
not_reached();
|
|
}
|
|
if (ziszero(z)) {
|
|
*res = _zero_;
|
|
return 0;
|
|
}
|
|
if (n < 0) {
|
|
n = -n;
|
|
z.sign = !z.sign;
|
|
}
|
|
if (n == 1) {
|
|
zcopy(z, res);
|
|
return 0;
|
|
}
|
|
/*
|
|
* If the division is by a large number, then call the normal
|
|
* divide routine.
|
|
*/
|
|
if (n & ~BASE1) {
|
|
div.sign = 0;
|
|
div.v = divval;
|
|
divval[0] = (HALF) n;
|
|
#if LONG_BITS > BASEB
|
|
divval[1] = (HALF)(((FULL) n) >> BASEB);
|
|
div.len = 2;
|
|
#else
|
|
div.len = 1;
|
|
#endif
|
|
zdiv(z, div, res, &dest, 0);
|
|
n = ztolong(dest);
|
|
zfree(dest);
|
|
return n;
|
|
}
|
|
/*
|
|
* Division is by a small number, so we can be quick about it.
|
|
*/
|
|
len = z.len;
|
|
dest.sign = z.sign;
|
|
dest.len = len;
|
|
dest.v = alloc(len);
|
|
h1 = z.v + len;
|
|
sd = dest.v + len;
|
|
val = 0;
|
|
while (len--) {
|
|
val = ((val << BASEB) + ((FULL) *--h1));
|
|
*--sd = (HALF)(val / n);
|
|
val %= n;
|
|
}
|
|
zquicktrim(dest);
|
|
*res = dest;
|
|
return (long)val;
|
|
}
|
|
|
|
|
|
|
|
/*
|
|
* Calculate the mod of an integer by a small number.
|
|
* This is only defined for positive moduli.
|
|
*/
|
|
long
|
|
zmodi(ZVALUE z, long n)
|
|
{
|
|
register HALF *h1;
|
|
FULL val;
|
|
HALF divval[2];
|
|
ZVALUE div;
|
|
ZVALUE temp;
|
|
long len;
|
|
|
|
if (n == 0) {
|
|
math_error("Division by zero");
|
|
not_reached();
|
|
}
|
|
if (n < 0) {
|
|
math_error("Non-positive modulus");
|
|
not_reached();
|
|
}
|
|
if (ziszero(z) || (n == 1))
|
|
return 0;
|
|
if (zisone(z))
|
|
return 1;
|
|
/*
|
|
* If the modulus is by a large number, then call the normal
|
|
* modulo routine.
|
|
*/
|
|
if (n & ~BASE1) {
|
|
div.sign = 0;
|
|
div.v = divval;
|
|
divval[0] = (HALF) n;
|
|
#if LONG_BITS > BASEB
|
|
divval[1] = (HALF)(((FULL) n) >> BASEB);
|
|
div.len = 2;
|
|
#else
|
|
div.len = 1;
|
|
#endif
|
|
zmod(z, div, &temp, 0);
|
|
n = ztolong(temp);
|
|
zfree(temp);
|
|
return n;
|
|
}
|
|
/*
|
|
* The modulus is by a small number, so we can do this quickly.
|
|
*/
|
|
len = z.len;
|
|
h1 = z.v + len;
|
|
val = 0;
|
|
while (len-- > 0)
|
|
val = ((val << BASEB) + ((FULL) *--h1)) % n;
|
|
if (val && z.sign)
|
|
val = n - val;
|
|
return (long)val;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return whether or not one number exactly divides another one.
|
|
* Returns true if division occurs with no remainder.
|
|
* z1 is the number to be divided by z2.
|
|
*/
|
|
|
|
bool
|
|
zdivides(ZVALUE z1, ZVALUE z2)
|
|
{
|
|
LEN i, j, k, m, n;
|
|
HALF u, v, w, x;
|
|
HALF *a, *a0, *A, *b, *B, *c, *d;
|
|
FULL f;
|
|
bool ans;
|
|
|
|
if (zisunit(z2) || ziszero(z1)) return true;
|
|
if (ziszero(z2)) return false;
|
|
|
|
m = z1.len;
|
|
n = z2.len;
|
|
if (m < n) return false;
|
|
|
|
c = z1.v;
|
|
d = z2.v;
|
|
|
|
while (!*d) {
|
|
if (*c++) return false;
|
|
d++;
|
|
m--;
|
|
n--;
|
|
}
|
|
|
|
j = 0;
|
|
u = *c;
|
|
v = *d;
|
|
while (!(v & 1)) { /* Counting and checking zero bits */
|
|
if (u & 1) return false;
|
|
u >>= 1;
|
|
v >>= 1;
|
|
j++;
|
|
}
|
|
|
|
if (n == 1 && v == 1) return true;
|
|
if (!*c) { /* Removing any further zeros */
|
|
while(!*++c)
|
|
m--;
|
|
c--;
|
|
}
|
|
|
|
if (m < n) return false;
|
|
|
|
if (j) {
|
|
B = alloc((LEN)n); /* Array for shifted z2 */
|
|
d += n;
|
|
b = B + n;
|
|
i = n;
|
|
f = 0;
|
|
while(i--) {
|
|
f = f << BASEB | *--d;
|
|
*--b = (HALF)(f >> j);
|
|
}
|
|
if (!B[n - 1]) n--;
|
|
}
|
|
else B = d;
|
|
u = *B;
|
|
v = x = 1;
|
|
w = 0;
|
|
while (x) { /* Finding minv(*B, BASE) */
|
|
if (v & x) {
|
|
v -= x * u;
|
|
w |= x;
|
|
}
|
|
x <<= 1;
|
|
}
|
|
|
|
A = alloc((LEN)(m + 2)); /* Main working array */
|
|
memcpy(A, c, m * sizeof(HALF));
|
|
A[m + 1] = 1;
|
|
A[m] = 0;
|
|
|
|
k = m - n + 1; /* Length of presumed quotient */
|
|
|
|
a0 = A;
|
|
|
|
while (k--) {
|
|
if (*a0) {
|
|
x = w * *a0;
|
|
a = a0;
|
|
b = B;
|
|
i = n;
|
|
u = 0;
|
|
while (i--) {
|
|
f = (FULL) *a - (FULL) x * *b++ - u;
|
|
*a++ = (HALF)f;
|
|
u = -(HALF)(f >> BASEB);
|
|
}
|
|
f = (FULL) *a - u;
|
|
*a++ = (HALF)f;
|
|
u = -(HALF)(f >> BASEB);
|
|
if (u) {
|
|
while (*a == 0) *a++ = BASE1;
|
|
(*a)--;
|
|
}
|
|
}
|
|
a0++;
|
|
}
|
|
ans = false;
|
|
if (A[m + 1]) {
|
|
a = A + m;
|
|
while (--n && !*--a);
|
|
if (!n) ans = true;
|
|
}
|
|
freeh(A);
|
|
if (j) freeh(B);
|
|
return ans;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute the bitwise OR of two integers
|
|
*/
|
|
void
|
|
zor(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
register HALF *sp, *dp;
|
|
long len;
|
|
ZVALUE bz, lz, dest;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (z1.len >= z2.len) {
|
|
bz = z1;
|
|
lz = z2;
|
|
} else {
|
|
bz = z2;
|
|
lz = z1;
|
|
}
|
|
dest.len = bz.len;
|
|
dest.v = alloc(dest.len);
|
|
dest.sign = 0;
|
|
zcopyval(bz, dest);
|
|
len = lz.len;
|
|
sp = lz.v;
|
|
dp = dest.v;
|
|
while (len--)
|
|
*dp++ |= *sp++;
|
|
*res = dest;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute the bitwise AND of two integers
|
|
*/
|
|
void
|
|
zand(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
HALF *h1, *h2, *hd;
|
|
LEN len;
|
|
ZVALUE dest;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
len = ((z1.len <= z2.len) ? z1.len : z2.len);
|
|
h1 = &z1.v[len-1];
|
|
h2 = &z2.v[len-1];
|
|
while ((len > 1) && ((*h1 & *h2) == 0)) {
|
|
h1--;
|
|
h2--;
|
|
len--;
|
|
}
|
|
dest.len = len;
|
|
dest.v = alloc(len);
|
|
dest.sign = 0;
|
|
h1 = z1.v;
|
|
h2 = z2.v;
|
|
hd = dest.v;
|
|
while (len--)
|
|
*hd++ = (*h1++ & *h2++);
|
|
*res = dest;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute the bitwise XOR of two integers.
|
|
*/
|
|
void
|
|
zxor(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
HALF *dp, *h1, *h2;
|
|
LEN len, j, k;
|
|
ZVALUE dest;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
h1 = z1.v;
|
|
h2 = z2.v;
|
|
len = z1.len;
|
|
j = z2.len;
|
|
if (z1.len < z2.len) {
|
|
len = z2.len;
|
|
j = z1.len;
|
|
h1 = z2.v;
|
|
h2 = z1.v;
|
|
} else if (z1.len == z2.len) {
|
|
while (len > 1 && z1.v[len-1] == z2.v[len-1])
|
|
len--;
|
|
j = len;
|
|
}
|
|
k = len - j;
|
|
dest.len = len;
|
|
dest.v = alloc(len);
|
|
dest.sign = 0;
|
|
dp = dest.v;
|
|
while (j-- > 0)
|
|
*dp++ = *h1++ ^ *h2++;
|
|
while (k-- > 0)
|
|
*dp++ = *h1++;
|
|
*res = dest;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute the bitwise AND NOT of two integers.
|
|
*/
|
|
void
|
|
zandnot(ZVALUE z1, ZVALUE z2, ZVALUE *res)
|
|
{
|
|
HALF *dp, *h1, *h2;
|
|
LEN len, j, k;
|
|
ZVALUE dest;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
len = z1.len;
|
|
if (z2.len >= len) {
|
|
while (len > 1 && (z1.v[len-1] & ~z2.v[len-1]) == 0)
|
|
len--;
|
|
j = len;
|
|
k = 0;
|
|
} else {
|
|
j = z2.len;
|
|
k = len - z2.len;
|
|
}
|
|
dest.len = len;
|
|
dest.v = alloc(len);
|
|
dest.sign = 0;
|
|
dp = dest.v;
|
|
h1 = z1.v;
|
|
h2 = z2.v;
|
|
while (j-- > 0)
|
|
*dp++ = *h1++ & ~*h2++;
|
|
while (k-- > 0)
|
|
*dp++ = *h1++;
|
|
*res = dest;
|
|
}
|
|
|
|
|
|
/*
|
|
* Shift a number left (or right) by the specified number of bits.
|
|
* Positive shift means to the left. When shifting right, rightmost
|
|
* bits are lost. The sign of the number is preserved.
|
|
*/
|
|
void
|
|
zshift(ZVALUE z, long n, ZVALUE *res)
|
|
{
|
|
ZVALUE ans;
|
|
LEN hc; /* number of halfwords shift is by */
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (ziszero(z)) {
|
|
*res = _zero_;
|
|
return;
|
|
}
|
|
if (n == 0) {
|
|
zcopy(z, res);
|
|
return;
|
|
}
|
|
/*
|
|
* If shift value is negative, then shift right.
|
|
* Check for large shifts, and handle word-sized shifts quickly.
|
|
*/
|
|
if (n < 0) {
|
|
n = -n;
|
|
if ((n < 0) || (n >= (z.len * BASEB))) {
|
|
*res = _zero_;
|
|
return;
|
|
}
|
|
hc = (LEN)(n / BASEB);
|
|
n %= BASEB;
|
|
z.v += hc;
|
|
z.len -= hc;
|
|
ans.len = z.len;
|
|
ans.v = alloc(ans.len);
|
|
ans.sign = z.sign;
|
|
zcopyval(z, ans);
|
|
if (n > 0) {
|
|
zshiftr(ans, n);
|
|
ztrim(&ans);
|
|
}
|
|
if (ziszero(ans)) {
|
|
zfree(ans);
|
|
ans = _zero_;
|
|
}
|
|
*res = ans;
|
|
return;
|
|
}
|
|
/*
|
|
* Shift value is positive, so shift leftwards.
|
|
* Check specially for a shift of the value 1, since this is common.
|
|
* Also handle word-sized shifts quickly.
|
|
*/
|
|
if (zisunit(z)) {
|
|
zbitvalue(n, res);
|
|
res->sign = z.sign;
|
|
return;
|
|
}
|
|
hc = (LEN)(n / BASEB);
|
|
n %= BASEB;
|
|
ans.len = z.len + hc + 1;
|
|
ans.v = alloc(ans.len);
|
|
ans.sign = z.sign;
|
|
if (hc > 0)
|
|
memset((char *) ans.v, 0, hc * sizeof(HALF));
|
|
memcpy((char *) (ans.v + hc),
|
|
(char *) z.v, z.len * sizeof(HALF));
|
|
ans.v[ans.len - 1] = 0;
|
|
if (n > 0) {
|
|
ans.v += hc;
|
|
ans.len -= hc;
|
|
zshiftl(ans, n);
|
|
ans.v -= hc;
|
|
ans.len += hc;
|
|
}
|
|
ztrim(&ans);
|
|
*res = ans;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the position of the lowest bit which is set in the binary
|
|
* representation of a number (counting from zero). This is the highest
|
|
* power of two which evenly divides the number.
|
|
*/
|
|
long
|
|
zlowbit(ZVALUE z)
|
|
{
|
|
register HALF *zp;
|
|
long n;
|
|
HALF dataval;
|
|
HALF *bitval;
|
|
|
|
n = 0;
|
|
for (zp = z.v; *zp == 0; zp++)
|
|
if (++n >= z.len)
|
|
return 0;
|
|
dataval = *zp;
|
|
bitval = bitmask;
|
|
/* ignore Saber-C warning #530 about empty while statement */
|
|
/* OK to ignore in proc zlowbit */
|
|
while ((*(bitval++) & dataval) == 0) {
|
|
}
|
|
return (n*BASEB)+(bitval-bitmask-1);
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the position of the highest bit which is set in the binary
|
|
* representation of a number (counting from zero). This is the highest power
|
|
* of two which is less than or equal to the number (which is assumed nonzero).
|
|
*/
|
|
LEN
|
|
zhighbit(ZVALUE z)
|
|
{
|
|
HALF dataval;
|
|
HALF *bitval;
|
|
|
|
dataval = z.v[z.len-1];
|
|
if (dataval == 0) {
|
|
return 0;
|
|
}
|
|
bitval = bitmask+BASEB;
|
|
if (dataval) {
|
|
/* ignore Saber-C warning #530 about empty while statement */
|
|
/* OK to ignore in proc zhighbit */
|
|
while ((*(--bitval) & dataval) == 0) {
|
|
}
|
|
}
|
|
return (z.len*BASEB)+(LEN)(bitval-bitmask-BASEB);
|
|
}
|
|
|
|
|
|
/*
|
|
* Return whether or not the specified bit number is set in a number.
|
|
* Rightmost bit of a number is bit 0.
|
|
*/
|
|
bool
|
|
zisset(ZVALUE z, long n)
|
|
{
|
|
if ((n < 0) || ((n / BASEB) >= z.len))
|
|
return false;
|
|
return ((z.v[n / BASEB] & (((HALF) 1) << (n % BASEB))) != 0);
|
|
}
|
|
|
|
|
|
/*
|
|
* Check whether or not a number has exactly one bit set, and
|
|
* thus is an exact power of two. Returns true if so.
|
|
*/
|
|
bool
|
|
zisonebit(ZVALUE z)
|
|
{
|
|
register HALF *hp;
|
|
register LEN len;
|
|
|
|
if (ziszero(z) || zisneg(z))
|
|
return false;
|
|
hp = z.v;
|
|
len = z.len;
|
|
while (len > 4) {
|
|
len -= 4;
|
|
if (*hp++ || *hp++ || *hp++ || *hp++)
|
|
return false;
|
|
}
|
|
while (--len > 0) {
|
|
if (*hp++)
|
|
return false;
|
|
}
|
|
return ((*hp & -*hp) == *hp); /* NEEDS 2'S COMPLEMENT */
|
|
}
|
|
|
|
|
|
/*
|
|
* Check whether or not a number has all of its bits set below some
|
|
* bit position, and thus is one less than an exact power of two.
|
|
* Returns true if so.
|
|
*/
|
|
bool
|
|
zisallbits(ZVALUE z)
|
|
{
|
|
register HALF *hp;
|
|
register LEN len;
|
|
HALF digit;
|
|
|
|
if (ziszero(z) || zisneg(z))
|
|
return false;
|
|
hp = z.v;
|
|
len = z.len;
|
|
while (len > 4) {
|
|
len -= 4;
|
|
if ((*hp++ != BASE1) || (*hp++ != BASE1) ||
|
|
(*hp++ != BASE1) || (*hp++ != BASE1))
|
|
return false;
|
|
}
|
|
while (--len > 0) {
|
|
if (*hp++ != BASE1)
|
|
return false;
|
|
}
|
|
digit = (HALF)(*hp + 1);
|
|
return ((digit & -digit) == digit); /* NEEDS 2'S COMPLEMENT */
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the number whose binary representation contains only one bit which
|
|
* is in the specified position (counting from zero). This is equivalent
|
|
* to raising two to the given power.
|
|
*/
|
|
void
|
|
zbitvalue(long n, ZVALUE *res)
|
|
{
|
|
ZVALUE z;
|
|
|
|
/* firewall */
|
|
if (res == NULL) {
|
|
math_error("%s: res NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
if (n < 0) n = 0;
|
|
z.sign = 0;
|
|
z.len = (LEN)((n / BASEB) + 1);
|
|
z.v = alloc(z.len);
|
|
zclearval(z);
|
|
z.v[z.len-1] = (((HALF) 1) << (n % BASEB));
|
|
*res = z;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compare a number against zero.
|
|
* Returns the sgn function of the number (-1, 0, or 1).
|
|
*/
|
|
FLAG
|
|
ztest(ZVALUE z)
|
|
{
|
|
register int sign;
|
|
register HALF *h;
|
|
register long len;
|
|
|
|
sign = 1;
|
|
if (z.sign)
|
|
sign = -sign;
|
|
h = z.v;
|
|
len = z.len;
|
|
while (len--)
|
|
if (*h++)
|
|
return sign;
|
|
return 0;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the sign of z1 - z2, i.e. 1 if the first integer is greater,
|
|
* 0 if they are equal, -1 otherwise.
|
|
*/
|
|
FLAG
|
|
zrel(ZVALUE z1, ZVALUE z2)
|
|
{
|
|
HALF *h1, *h2;
|
|
LEN len;
|
|
int sign;
|
|
|
|
sign = 1;
|
|
if (z1.sign < z2.sign)
|
|
return 1;
|
|
if (z2.sign < z1.sign)
|
|
return -1;
|
|
if (z2.sign)
|
|
sign = -1;
|
|
if (z1.len != z2.len)
|
|
return (z1.len > z2.len) ? sign : -sign;
|
|
len = z1.len;
|
|
h1 = z1.v + len;
|
|
h2 = z2.v + len;
|
|
|
|
while (len > 0) {
|
|
if (*--h1 != *--h2)
|
|
break;
|
|
len--;
|
|
}
|
|
if (len > 0)
|
|
return (*h1 > *h2) ? sign : -sign;
|
|
return 0;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the sign of abs(z1) - abs(z2), i.e. 1 if the first integer
|
|
* has greater absolute value, 0 is they have equal absolute value,
|
|
* -1 otherwise.
|
|
*/
|
|
FLAG
|
|
zabsrel(ZVALUE z1, ZVALUE z2)
|
|
{
|
|
HALF *h1, *h2;
|
|
LEN len;
|
|
|
|
if (z1.len != z2.len)
|
|
return (z1.len > z2.len) ? 1 : -1;
|
|
|
|
len = z1.len;
|
|
h1 = z1.v + len;
|
|
h2 = z2.v + len;
|
|
while (len > 0) {
|
|
if (*--h1 != *--h2)
|
|
break;
|
|
len--;
|
|
}
|
|
if (len > 0)
|
|
return (*h1 > *h2) ? 1 : -1;
|
|
return 0;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compare two numbers to see if they are equal or not.
|
|
* Returns true if they differ.
|
|
*/
|
|
bool
|
|
zcmp(ZVALUE z1, ZVALUE z2)
|
|
{
|
|
register HALF *h1, *h2;
|
|
register long len;
|
|
|
|
if ((z1.sign != z2.sign) || (z1.len != z2.len) || (*z1.v != *z2.v))
|
|
return true;
|
|
len = z1.len;
|
|
h1 = z1.v;
|
|
h2 = z2.v;
|
|
while (--len > 0) {
|
|
if (*++h1 != *++h2)
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
|
|
/*
|
|
* Utility to calculate the gcd of two FULL integers.
|
|
*/
|
|
FULL
|
|
uugcd(FULL f1, FULL f2)
|
|
{
|
|
FULL temp;
|
|
|
|
while (f1) {
|
|
temp = f2 % f1;
|
|
f2 = f1;
|
|
f1 = temp;
|
|
}
|
|
return (FULL) f2;
|
|
}
|
|
|
|
|
|
/*
|
|
* Utility to calculate the gcd of two small integers.
|
|
*/
|
|
long
|
|
iigcd(long i1, long i2)
|
|
{
|
|
FULL f1, f2, temp;
|
|
|
|
f1 = (FULL) ((i1 >= 0) ? i1 : -i1);
|
|
f2 = (FULL) ((i2 >= 0) ? i2 : -i2);
|
|
while (f1) {
|
|
temp = f2 % f1;
|
|
f2 = f1;
|
|
f1 = temp;
|
|
}
|
|
return (long) f2;
|
|
}
|
|
|
|
|
|
void
|
|
ztrim(ZVALUE *z)
|
|
{
|
|
HALF *h;
|
|
LEN len;
|
|
|
|
/* firewall */
|
|
if (z == NULL) {
|
|
math_error("%s: z NULL", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
h = z->v + z->len - 1;
|
|
len = z->len;
|
|
while (*h == 0 && len > 1) {
|
|
--h;
|
|
--len;
|
|
}
|
|
z->len = len;
|
|
}
|
|
|
|
|
|
/*
|
|
* Utility routine to shift right.
|
|
*
|
|
* NOTE: The ZVALUE length is not adjusted instead, the value is
|
|
* zero padded from the left. One may need to call ztrim()
|
|
* or use zshift() instead.
|
|
*/
|
|
void
|
|
zshiftr(ZVALUE z, long n)
|
|
{
|
|
register HALF *h, *lim;
|
|
FULL mask, maskt;
|
|
long len;
|
|
|
|
if (n >= BASEB) {
|
|
len = n / BASEB;
|
|
h = z.v;
|
|
lim = z.v + z.len - len;
|
|
while (h < lim) {
|
|
h[0] = h[len];
|
|
++h;
|
|
}
|
|
n -= BASEB * len;
|
|
lim = z.v + z.len;
|
|
while (h < lim)
|
|
*h++ = 0;
|
|
}
|
|
if (n) {
|
|
len = z.len;
|
|
h = z.v + len;
|
|
mask = 0;
|
|
while (len--) {
|
|
maskt = (((FULL) *--h) << (BASEB - n)) & BASE1;
|
|
*h = ((*h >> n) | (HALF)mask);
|
|
mask = maskt;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
/*
|
|
* Utility routine to shift left.
|
|
*
|
|
* NOTE: The ZVALUE length is not adjusted. The bits in the upper
|
|
* HALF are simply tossed. You may want to use zshift() instead.
|
|
*/
|
|
void
|
|
zshiftl(ZVALUE z, long n)
|
|
{
|
|
register HALF *h;
|
|
FULL mask, i;
|
|
long len;
|
|
|
|
if (n >= BASEB) {
|
|
len = n / BASEB;
|
|
h = z.v + z.len - 1;
|
|
while (!*h)
|
|
--h;
|
|
while (h >= z.v) {
|
|
h[len] = h[0];
|
|
--h;
|
|
}
|
|
n -= BASEB * len;
|
|
while (len)
|
|
h[len--] = 0;
|
|
}
|
|
if (n > 0) {
|
|
len = z.len;
|
|
h = z.v;
|
|
mask = 0;
|
|
while (len--) {
|
|
i = (((FULL) *h) << n) | mask;
|
|
if (i > BASE1) {
|
|
mask = i >> BASEB;
|
|
i &= BASE1;
|
|
} else {
|
|
mask = 0;
|
|
}
|
|
*h = (HALF) i;
|
|
++h;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
/*
|
|
* popcnt - count the number of 0 or 1 bits in an integer
|
|
*
|
|
* We ignore all 0 bits above the highest bit.
|
|
*/
|
|
long
|
|
zpopcnt(ZVALUE z, int bitval)
|
|
{
|
|
long cnt = 0; /* number of times found */
|
|
HALF h; /* HALF to count */
|
|
int i;
|
|
|
|
/*
|
|
* count 1's
|
|
*/
|
|
if (bitval) {
|
|
|
|
/*
|
|
* count each HALF
|
|
*/
|
|
for (i=0; i < z.len; ++i) {
|
|
/* count each octet */
|
|
for (h = z.v[i]; h; h >>= 8) {
|
|
cnt += (long)popcnt[h & 0xff];
|
|
}
|
|
}
|
|
|
|
/*
|
|
* count 0's
|
|
*/
|
|
} else {
|
|
|
|
/*
|
|
* count each HALF up until the last
|
|
*/
|
|
for (i=0; i < z.len-1; ++i) {
|
|
|
|
/* count each octet */
|
|
cnt += BASEB;
|
|
for (h = z.v[i]; h; h >>= 8) {
|
|
cnt -= (long)popcnt[h & 0xff];
|
|
}
|
|
}
|
|
|
|
/*
|
|
* count the last octet up until the highest 1 bit
|
|
*/
|
|
for (h = z.v[z.len-1]; h; h>>=1) {
|
|
/* count each 0 bit */
|
|
if ((h & 0x1) == 0) {
|
|
++cnt;
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* return count
|
|
*/
|
|
return cnt;
|
|
}
|