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.
1504 lines
36 KiB
C
1504 lines
36 KiB
C
/*
|
|
* qmath - extended precision rational arithmetic primitive routines
|
|
*
|
|
* Copyright (C) 1999-2007,2014,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:21
|
|
* File existed as early as: before 1990
|
|
*
|
|
* Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
|
|
*/
|
|
|
|
|
|
#include <stdio.h>
|
|
#include "qmath.h"
|
|
#include "config.h"
|
|
|
|
|
|
#include "errtbl.h"
|
|
#include "banned.h" /* include after system header <> includes */
|
|
|
|
|
|
NUMBER _qzero_ = { { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
|
|
NUMBER _qone_ = { { _oneval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
|
|
NUMBER _qtwo_ = { { _twoval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
|
|
NUMBER _qten_ = { { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
|
|
NUMBER _qnegone_ = { { _oneval_, 1, 1 }, { _oneval_, 1, 0 }, 1, NULL };
|
|
NUMBER _qonehalf_ = { { _oneval_, 1, 0 }, { _twoval_, 1, 0 }, 1, NULL };
|
|
NUMBER _qneghalf_ = { { _oneval_, 1, 1 }, { _twoval_, 1, 0 }, 1, NULL };
|
|
NUMBER _qonesqbase_ = { { _oneval_, 1, 0 }, { _sqbaseval_, 2, 0 }, 1, NULL };
|
|
|
|
NUMBER * initnumbs[] = {&_qzero_, &_qone_, &_qtwo_,
|
|
&_qten_, &_qnegone_, &_qonehalf_, &_qneghalf_,
|
|
&_qonesqbase_,
|
|
NULL /* must be last */
|
|
};
|
|
|
|
|
|
/*
|
|
* Create another copy of a number.
|
|
* q2 = qcopy(q1);
|
|
*/
|
|
NUMBER *
|
|
qcopy(NUMBER *q)
|
|
{
|
|
register NUMBER *r;
|
|
|
|
r = qalloc();
|
|
r->num.sign = q->num.sign;
|
|
if (!zisunit(q->num)) {
|
|
r->num.len = q->num.len;
|
|
r->num.v = alloc(r->num.len);
|
|
zcopyval(q->num, r->num);
|
|
}
|
|
if (!zisunit(q->den)) {
|
|
r->den.len = q->den.len;
|
|
r->den.v = alloc(r->den.len);
|
|
zcopyval(q->den, r->den);
|
|
}
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Convert a number to a normal integer.
|
|
* i = qtoi(q);
|
|
*/
|
|
long
|
|
qtoi(NUMBER *q)
|
|
{
|
|
long i;
|
|
ZVALUE res;
|
|
|
|
if (qisint(q))
|
|
return ztoi(q->num);
|
|
zquo(q->num, q->den, &res, 0);
|
|
i = ztoi(res);
|
|
zfree(res);
|
|
return i;
|
|
}
|
|
|
|
|
|
/*
|
|
* Convert a normal integer into a number.
|
|
* q = itoq(i);
|
|
*/
|
|
NUMBER *
|
|
itoq(long i)
|
|
{
|
|
register NUMBER *q;
|
|
|
|
if ((i >= -1) && (i <= 10)) {
|
|
switch ((int) i) {
|
|
case 0: q = &_qzero_; break;
|
|
case 1: q = &_qone_; break;
|
|
case 2: q = &_qtwo_; break;
|
|
case 10: q = &_qten_; break;
|
|
case -1: q = &_qnegone_; break;
|
|
default: q = NULL;
|
|
}
|
|
if (q)
|
|
return qlink(q);
|
|
}
|
|
q = qalloc();
|
|
itoz(i, &q->num);
|
|
return q;
|
|
}
|
|
|
|
|
|
/*
|
|
* Convert a number to a normal unsigned integer.
|
|
* u = qtou(q);
|
|
*/
|
|
FULL
|
|
qtou(NUMBER *q)
|
|
{
|
|
FULL i;
|
|
ZVALUE res;
|
|
|
|
if (qisint(q))
|
|
return ztou(q->num);
|
|
zquo(q->num, q->den, &res, 0);
|
|
i = ztou(res);
|
|
zfree(res);
|
|
return i;
|
|
}
|
|
|
|
|
|
/*
|
|
* Convert a number to a normal signed integer.
|
|
* s = qtos(q);
|
|
*/
|
|
SFULL
|
|
qtos(NUMBER *q)
|
|
{
|
|
SFULL i;
|
|
ZVALUE res;
|
|
|
|
if (qisint(q))
|
|
return ztos(q->num);
|
|
zquo(q->num, q->den, &res, 0);
|
|
i = ztos(res);
|
|
zfree(res);
|
|
return i;
|
|
}
|
|
|
|
|
|
/*
|
|
* Convert a normal unsigned integer into a number.
|
|
* q = utoq(i);
|
|
*/
|
|
NUMBER *
|
|
utoq(FULL i)
|
|
{
|
|
register NUMBER *q;
|
|
|
|
if (i <= 10) {
|
|
switch ((int) i) {
|
|
case 0: q = &_qzero_; break;
|
|
case 1: q = &_qone_; break;
|
|
case 2: q = &_qtwo_; break;
|
|
case 10: q = &_qten_; break;
|
|
default: q = NULL;
|
|
}
|
|
if (q)
|
|
return qlink(q);
|
|
}
|
|
q = qalloc();
|
|
utoz(i, &q->num);
|
|
return q;
|
|
}
|
|
|
|
|
|
/*
|
|
* Convert a normal signed integer into a number.
|
|
* q = stoq(s);
|
|
*/
|
|
NUMBER *
|
|
stoq(SFULL i)
|
|
{
|
|
register NUMBER *q;
|
|
|
|
if (i <= 10) {
|
|
switch ((int) i) {
|
|
case 0: q = &_qzero_; break;
|
|
case 1: q = &_qone_; break;
|
|
case 2: q = &_qtwo_; break;
|
|
case 10: q = &_qten_; break;
|
|
default: q = NULL;
|
|
}
|
|
if (q)
|
|
return qlink(q);
|
|
}
|
|
q = qalloc();
|
|
stoz(i, &q->num);
|
|
return q;
|
|
}
|
|
|
|
|
|
/*
|
|
* Create a number from the given FULL numerator and denominator.
|
|
* q = uutoq(inum, iden);
|
|
*/
|
|
NUMBER *
|
|
uutoq(FULL inum, FULL iden)
|
|
{
|
|
register NUMBER *q;
|
|
FULL d;
|
|
bool sign;
|
|
|
|
if (iden == 0) {
|
|
math_error("Division by zero");
|
|
not_reached();
|
|
}
|
|
if (inum == 0)
|
|
return qlink(&_qzero_);
|
|
sign = 0;
|
|
d = uugcd(inum, iden);
|
|
inum /= d;
|
|
iden /= d;
|
|
if (iden == 1)
|
|
return utoq(inum);
|
|
q = qalloc();
|
|
if (inum != 1)
|
|
utoz(inum, &q->num);
|
|
utoz(iden, &q->den);
|
|
q->num.sign = sign;
|
|
return q;
|
|
}
|
|
|
|
|
|
/*
|
|
* Create a number from the given integral numerator and denominator.
|
|
* q = iitoq(inum, iden);
|
|
*/
|
|
NUMBER *
|
|
iitoq(long inum, long iden)
|
|
{
|
|
register NUMBER *q;
|
|
long d;
|
|
bool sign;
|
|
|
|
if (iden == 0) {
|
|
math_error("Division by zero");
|
|
not_reached();
|
|
}
|
|
if (inum == 0)
|
|
return qlink(&_qzero_);
|
|
sign = 0;
|
|
if (inum < 0) {
|
|
sign = 1;
|
|
inum = -inum;
|
|
}
|
|
if (iden < 0) {
|
|
sign = 1 - sign;
|
|
iden = -iden;
|
|
}
|
|
d = iigcd(inum, iden);
|
|
inum /= d;
|
|
iden /= d;
|
|
if (iden == 1)
|
|
return itoq(sign ? -inum : inum);
|
|
q = qalloc();
|
|
if (inum != 1)
|
|
itoz(inum, &q->num);
|
|
itoz(iden, &q->den);
|
|
q->num.sign = sign;
|
|
return q;
|
|
}
|
|
|
|
|
|
/*
|
|
* Add two numbers to each other.
|
|
* q3 = qqadd(q1, q2);
|
|
*/
|
|
NUMBER *
|
|
qqadd(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
NUMBER *r;
|
|
ZVALUE t1, t2, temp, d1, d2, vpd1, upd1;
|
|
|
|
if (qiszero(q1))
|
|
return qlink(q2);
|
|
if (qiszero(q2))
|
|
return qlink(q1);
|
|
r = qalloc();
|
|
/*
|
|
* If either number is an integer, then the result is easy.
|
|
*/
|
|
if (qisint(q1) && qisint(q2)) {
|
|
zadd(q1->num, q2->num, &r->num);
|
|
return r;
|
|
}
|
|
if (qisint(q2)) {
|
|
zmul(q1->den, q2->num, &temp);
|
|
zadd(q1->num, temp, &r->num);
|
|
zfree(temp);
|
|
zcopy(q1->den, &r->den);
|
|
return r;
|
|
}
|
|
if (qisint(q1)) {
|
|
zmul(q2->den, q1->num, &temp);
|
|
zadd(q2->num, temp, &r->num);
|
|
zfree(temp);
|
|
zcopy(q2->den, &r->den);
|
|
return r;
|
|
}
|
|
/*
|
|
* Both arguments are true fractions, so we need more work.
|
|
* If the denominators are relatively prime, then the answer is the
|
|
* straightforward cross product result with no need for reduction.
|
|
*/
|
|
zgcd(q1->den, q2->den, &d1);
|
|
if (zisunit(d1)) {
|
|
zfree(d1);
|
|
zmul(q1->num, q2->den, &t1);
|
|
zmul(q1->den, q2->num, &t2);
|
|
zadd(t1, t2, &r->num);
|
|
zfree(t1);
|
|
zfree(t2);
|
|
zmul(q1->den, q2->den, &r->den);
|
|
return r;
|
|
}
|
|
/*
|
|
* The calculation is now more complicated.
|
|
* See Knuth Vol 2 for details.
|
|
*/
|
|
zquo(q2->den, d1, &vpd1, 0);
|
|
zquo(q1->den, d1, &upd1, 0);
|
|
zmul(q1->num, vpd1, &t1);
|
|
zmul(q2->num, upd1, &t2);
|
|
zadd(t1, t2, &temp);
|
|
zfree(t1);
|
|
zfree(t2);
|
|
zfree(vpd1);
|
|
zgcd(temp, d1, &d2);
|
|
zfree(d1);
|
|
if (zisunit(d2)) {
|
|
zfree(d2);
|
|
r->num = temp;
|
|
zmul(upd1, q2->den, &r->den);
|
|
zfree(upd1);
|
|
return r;
|
|
}
|
|
zquo(temp, d2, &r->num, 0);
|
|
zfree(temp);
|
|
zquo(q2->den, d2, &temp, 0);
|
|
zfree(d2);
|
|
zmul(temp, upd1, &r->den);
|
|
zfree(temp);
|
|
zfree(upd1);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Subtract one number from another.
|
|
* q3 = qsub(q1, q2);
|
|
*/
|
|
NUMBER *
|
|
qsub(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
NUMBER *r;
|
|
|
|
if (q1 == q2)
|
|
return qlink(&_qzero_);
|
|
if (qiszero(q2))
|
|
return qlink(q1);
|
|
if (qisint(q1) && qisint(q2)) {
|
|
r = qalloc();
|
|
zsub(q1->num, q2->num, &r->num);
|
|
return r;
|
|
}
|
|
q2 = qneg(q2);
|
|
if (qiszero(q1))
|
|
return q2;
|
|
r = qqadd(q1, q2);
|
|
qfree(q2);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Increment a number by one.
|
|
*/
|
|
NUMBER *
|
|
qinc(NUMBER *q)
|
|
{
|
|
NUMBER *r;
|
|
|
|
r = qalloc();
|
|
if (qisint(q)) {
|
|
zadd(q->num, _one_, &r->num);
|
|
return r;
|
|
}
|
|
zadd(q->num, q->den, &r->num);
|
|
zcopy(q->den, &r->den);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Decrement a number by one.
|
|
*/
|
|
NUMBER *
|
|
qdec(NUMBER *q)
|
|
{
|
|
NUMBER *r;
|
|
|
|
r = qalloc();
|
|
if (qisint(q)) {
|
|
zsub(q->num, _one_, &r->num);
|
|
return r;
|
|
}
|
|
zsub(q->num, q->den, &r->num);
|
|
zcopy(q->den, &r->den);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Add a normal small integer value to an arbitrary number.
|
|
*/
|
|
NUMBER *
|
|
qaddi(NUMBER *q1, long n)
|
|
{
|
|
NUMBER addnum; /* temporary number */
|
|
HALF addval[2]; /* value of small number */
|
|
bool neg; /* true if number is neg */
|
|
#if LONG_BITS > BASEB
|
|
FULL nf;
|
|
#endif
|
|
|
|
if (n == 0)
|
|
return qlink(q1);
|
|
if (n == 1)
|
|
return qinc(q1);
|
|
if (n == -1)
|
|
return qdec(q1);
|
|
if (qiszero(q1))
|
|
return itoq(n);
|
|
addnum.num.sign = 0;
|
|
addnum.num.v = addval;
|
|
addnum.den = _one_;
|
|
neg = (n < 0);
|
|
if (neg)
|
|
n = -n;
|
|
addval[0] = (HALF) n;
|
|
#if LONG_BITS > BASEB
|
|
nf = (((FULL) n) >> BASEB);
|
|
if (nf) {
|
|
addval[1] = (HALF) nf;
|
|
addnum.num.len = 2;
|
|
}
|
|
#else
|
|
addnum.num.len = 1;
|
|
#endif
|
|
if (neg)
|
|
return qsub(q1, &addnum);
|
|
else
|
|
return qqadd(q1, &addnum);
|
|
}
|
|
|
|
|
|
/*
|
|
* Multiply two numbers.
|
|
* q3 = qmul(q1, q2);
|
|
*/
|
|
NUMBER *
|
|
qmul(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
NUMBER *r; /* returned value */
|
|
ZVALUE n1, n2, d1, d2; /* numerators and denominators */
|
|
ZVALUE tmp;
|
|
|
|
if (qiszero(q1) || qiszero(q2))
|
|
return qlink(&_qzero_);
|
|
if (qisone(q1))
|
|
return qlink(q2);
|
|
if (qisone(q2))
|
|
return qlink(q1);
|
|
if (qisint(q1) && qisint(q2)) { /* easy results if integers */
|
|
r = qalloc();
|
|
zmul(q1->num, q2->num, &r->num);
|
|
return r;
|
|
}
|
|
n1 = q1->num;
|
|
n2 = q2->num;
|
|
d1 = q1->den;
|
|
d2 = q2->den;
|
|
if (ziszero(d1) || ziszero(d2)) {
|
|
math_error("Division by zero");
|
|
not_reached();
|
|
}
|
|
if (ziszero(n1) || ziszero(n2))
|
|
return qlink(&_qzero_);
|
|
if (!zisunit(n1) && !zisunit(d2)) { /* possibly reduce */
|
|
zgcd(n1, d2, &tmp);
|
|
if (!zisunit(tmp)) {
|
|
zequo(q1->num, tmp, &n1);
|
|
zequo(q2->den, tmp, &d2);
|
|
}
|
|
zfree(tmp);
|
|
}
|
|
if (!zisunit(n2) && !zisunit(d1)) { /* again possibly reduce */
|
|
zgcd(n2, d1, &tmp);
|
|
if (!zisunit(tmp)) {
|
|
zequo(q2->num, tmp, &n2);
|
|
zequo(q1->den, tmp, &d1);
|
|
}
|
|
zfree(tmp);
|
|
}
|
|
r = qalloc();
|
|
zmul(n1, n2, &r->num);
|
|
zmul(d1, d2, &r->den);
|
|
if (q1->num.v != n1.v)
|
|
zfree(n1);
|
|
if (q1->den.v != d1.v)
|
|
zfree(d1);
|
|
if (q2->num.v != n2.v)
|
|
zfree(n2);
|
|
if (q2->den.v != d2.v)
|
|
zfree(d2);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Multiply a number by a small integer.
|
|
* q2 = qmuli(q1, n);
|
|
*/
|
|
NUMBER *
|
|
qmuli(NUMBER *q, long n)
|
|
{
|
|
NUMBER *r;
|
|
long d; /* gcd of multiplier and denominator */
|
|
int sign;
|
|
|
|
if ((n == 0) || qiszero(q))
|
|
return qlink(&_qzero_);
|
|
if (n == 1)
|
|
return qlink(q);
|
|
r = qalloc();
|
|
if (qisint(q)) {
|
|
zmuli(q->num, n, &r->num);
|
|
return r;
|
|
}
|
|
sign = 1;
|
|
if (n < 0) {
|
|
n = -n;
|
|
sign = -1;
|
|
}
|
|
d = zmodi(q->den, n);
|
|
d = iigcd(d, n);
|
|
zmuli(q->num, (n * sign) / d, &r->num);
|
|
(void) zdivi(q->den, d, &r->den);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Divide two numbers (as fractions).
|
|
* q3 = qqdiv(q1, q2);
|
|
*/
|
|
NUMBER *
|
|
qqdiv(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
NUMBER temp;
|
|
|
|
if (qiszero(q2)) {
|
|
math_error("Division by zero");
|
|
not_reached();
|
|
}
|
|
if ((q1 == q2) || !qcmp(q1, q2))
|
|
return qlink(&_qone_);
|
|
if (qisone(q1))
|
|
return qinv(q2);
|
|
temp.num = q2->den;
|
|
temp.den = q2->num;
|
|
temp.num.sign = temp.den.sign;
|
|
temp.den.sign = 0;
|
|
temp.links = 1;
|
|
return qmul(q1, &temp);
|
|
}
|
|
|
|
|
|
/*
|
|
* Divide a number by a small integer.
|
|
* q2 = qdivi(q1, n);
|
|
*/
|
|
NUMBER *
|
|
qdivi(NUMBER *q, long n)
|
|
{
|
|
NUMBER *r;
|
|
long d; /* gcd of divisor and numerator */
|
|
int sign;
|
|
|
|
if (n == 0) {
|
|
math_error("Division by zero");
|
|
not_reached();
|
|
}
|
|
if ((n == 1) || qiszero(q))
|
|
return qlink(q);
|
|
sign = 1;
|
|
if (n < 0) {
|
|
n = -n;
|
|
sign = -1;
|
|
}
|
|
r = qalloc();
|
|
d = zmodi(q->num, n);
|
|
d = iigcd(d, n);
|
|
(void) zdivi(q->num, d * sign, &r->num);
|
|
zmuli(q->den, n / d, &r->den);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the integer quotient of a pair of numbers
|
|
* If q1/q2 is an integer qquo(q1, q2) returns this integer
|
|
* If q2 is zero, zero is returned
|
|
* In other cases whether rounding is down, up, towards zero, etc.
|
|
* is determined by rnd.
|
|
*/
|
|
NUMBER *
|
|
qquo(NUMBER *q1, NUMBER *q2, long rnd)
|
|
{
|
|
ZVALUE tmp, tmp1, tmp2;
|
|
NUMBER *q;
|
|
|
|
if (qiszero(q1) || qiszero(q2))
|
|
return qlink(&_qzero_);
|
|
if (qisint(q1) && qisint(q2)) {
|
|
zquo(q1->num, q2->num, &tmp, rnd);
|
|
} else {
|
|
zmul(q1->num, q2->den, &tmp1);
|
|
zmul(q2->num, q1->den, &tmp2);
|
|
zquo(tmp1, tmp2, &tmp, rnd);
|
|
zfree(tmp1);
|
|
zfree(tmp2);
|
|
}
|
|
if (ziszero(tmp)) {
|
|
zfree(tmp);
|
|
return qlink(&_qzero_);
|
|
}
|
|
q = qalloc();
|
|
q->num = tmp;
|
|
return q;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the absolute value of a number.
|
|
* q2 = qqabs(q1);
|
|
*/
|
|
NUMBER *
|
|
qqabs(NUMBER *q)
|
|
{
|
|
register NUMBER *r;
|
|
|
|
if (q->num.sign == 0)
|
|
return qlink(q);
|
|
r = qalloc();
|
|
if (!zisunit(q->num))
|
|
zcopy(q->num, &r->num);
|
|
if (!zisunit(q->den))
|
|
zcopy(q->den, &r->den);
|
|
r->num.sign = 0;
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Negate a number.
|
|
* q2 = qneg(q1);
|
|
*/
|
|
NUMBER *
|
|
qneg(NUMBER *q)
|
|
{
|
|
register NUMBER *r;
|
|
|
|
if (qiszero(q))
|
|
return qlink(&_qzero_);
|
|
r = qalloc();
|
|
if (!zisunit(q->num))
|
|
zcopy(q->num, &r->num);
|
|
if (!zisunit(q->den))
|
|
zcopy(q->den, &r->den);
|
|
r->num.sign = !q->num.sign;
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the sign of a number (-1, 0 or 1)
|
|
*/
|
|
NUMBER *
|
|
qsign(NUMBER *q)
|
|
{
|
|
if (qiszero(q))
|
|
return qlink(&_qzero_);
|
|
if (qisneg(q))
|
|
return qlink(&_qnegone_);
|
|
return qlink(&_qone_);
|
|
}
|
|
|
|
|
|
/*
|
|
* Invert a number.
|
|
* q2 = qinv(q1);
|
|
*/
|
|
NUMBER *
|
|
qinv(NUMBER *q)
|
|
{
|
|
register NUMBER *r;
|
|
|
|
if (qisunit(q)) {
|
|
r = (qisneg(q) ? &_qnegone_ : &_qone_);
|
|
return qlink(r);
|
|
}
|
|
if (qiszero(q)) {
|
|
math_error("Division by zero");
|
|
not_reached();
|
|
}
|
|
r = qalloc();
|
|
if (!zisunit(q->num))
|
|
zcopy(q->num, &r->den);
|
|
if (!zisunit(q->den))
|
|
zcopy(q->den, &r->num);
|
|
r->num.sign = q->num.sign;
|
|
r->den.sign = 0;
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return just the numerator of a number.
|
|
* q2 = qnum(q1);
|
|
*/
|
|
NUMBER *
|
|
qnum(NUMBER *q)
|
|
{
|
|
register NUMBER *r;
|
|
|
|
if (qisint(q))
|
|
return qlink(q);
|
|
if (zisunit(q->num)) {
|
|
r = (qisneg(q) ? &_qnegone_ : &_qone_);
|
|
return qlink(r);
|
|
}
|
|
r = qalloc();
|
|
zcopy(q->num, &r->num);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return just the denominator of a number.
|
|
* q2 = qden(q1);
|
|
*/
|
|
NUMBER *
|
|
qden(NUMBER *q)
|
|
{
|
|
register NUMBER *r;
|
|
|
|
if (qisint(q))
|
|
return qlink(&_qone_);
|
|
r = qalloc();
|
|
zcopy(q->den, &r->num);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the fractional part of a number.
|
|
* q2 = qfrac(q1);
|
|
*/
|
|
NUMBER *
|
|
qfrac(NUMBER *q)
|
|
{
|
|
register NUMBER *r;
|
|
|
|
if (qisint(q))
|
|
return qlink(&_qzero_);
|
|
if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
|
|
(q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
|
|
return qlink(q);
|
|
r = qalloc();
|
|
zmod(q->num, q->den, &r->num, 2);
|
|
zcopy(q->den, &r->den);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the integral part of a number.
|
|
* q2 = qint(q1);
|
|
*/
|
|
NUMBER *
|
|
qint(NUMBER *q)
|
|
{
|
|
register NUMBER *r;
|
|
|
|
if (qisint(q))
|
|
return qlink(q);
|
|
if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
|
|
(q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
|
|
return qlink(&_qzero_);
|
|
r = qalloc();
|
|
zquo(q->num, q->den, &r->num, 2);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compute the square of a number.
|
|
*/
|
|
NUMBER *
|
|
qsquare(NUMBER *q)
|
|
{
|
|
ZVALUE num, zden;
|
|
|
|
if (qiszero(q))
|
|
return qlink(&_qzero_);
|
|
if (qisunit(q))
|
|
return qlink(&_qone_);
|
|
num = q->num;
|
|
zden = q->den;
|
|
q = qalloc();
|
|
if (!zisunit(num))
|
|
zsquare(num, &q->num);
|
|
if (!zisunit(zden))
|
|
zsquare(zden, &q->den);
|
|
return q;
|
|
}
|
|
|
|
|
|
/*
|
|
* Shift an integer by a given number of bits. This multiplies the number
|
|
* by the appropriate power of two. Positive numbers shift left, negative
|
|
* ones shift right. Low bits are truncated when shifting right.
|
|
*/
|
|
NUMBER *
|
|
qshift(NUMBER *q, long n)
|
|
{
|
|
register NUMBER *r;
|
|
|
|
if (qisfrac(q)) {
|
|
math_error("Shift of non-integer");
|
|
not_reached();
|
|
}
|
|
if (qiszero(q) || (n == 0))
|
|
return qlink(q);
|
|
if (n <= -(q->num.len * BASEB))
|
|
return qlink(&_qzero_);
|
|
r = qalloc();
|
|
zshift(q->num, n, &r->num);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Scale a number by a power of two, as in:
|
|
* ans = q * 2^n.
|
|
* This is similar to shifting, except that fractions work.
|
|
*/
|
|
NUMBER *
|
|
qscale(NUMBER *q, long pow)
|
|
{
|
|
long numshift, denshift, tmp;
|
|
NUMBER *r;
|
|
|
|
if (qiszero(q) || (pow == 0))
|
|
return qlink(q);
|
|
numshift = zisodd(q->num) ? 0 : zlowbit(q->num);
|
|
denshift = zisodd(q->den) ? 0 : zlowbit(q->den);
|
|
if (pow > 0) {
|
|
tmp = pow;
|
|
if (tmp > denshift)
|
|
tmp = denshift;
|
|
denshift = -tmp;
|
|
numshift = (pow - tmp);
|
|
} else {
|
|
pow = -pow;
|
|
tmp = pow;
|
|
if (tmp > numshift)
|
|
tmp = numshift;
|
|
numshift = -tmp;
|
|
denshift = (pow - tmp);
|
|
}
|
|
r = qalloc();
|
|
if (numshift)
|
|
zshift(q->num, numshift, &r->num);
|
|
else
|
|
zcopy(q->num, &r->num);
|
|
if (denshift)
|
|
zshift(q->den, denshift, &r->den);
|
|
else
|
|
zcopy(q->den, &r->den);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the minimum of two numbers.
|
|
*/
|
|
NUMBER *
|
|
qmin(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
if (q1 == q2)
|
|
return qlink(q1);
|
|
if (qrel(q1, q2) > 0)
|
|
q1 = q2;
|
|
return qlink(q1);
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the maximum of two numbers.
|
|
*/
|
|
NUMBER *
|
|
qmax(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
if (q1 == q2)
|
|
return qlink(q1);
|
|
if (qrel(q1, q2) < 0)
|
|
q1 = q2;
|
|
return qlink(q1);
|
|
}
|
|
|
|
|
|
/*
|
|
* Perform the bitwise OR of two integers.
|
|
*/
|
|
NUMBER *
|
|
qor(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
register NUMBER *r;
|
|
NUMBER *q1tmp, *q2tmp, *q;
|
|
|
|
if (qisfrac(q1) || qisfrac(q2)) {
|
|
math_error("Non-integers for bitwise or");
|
|
not_reached();
|
|
}
|
|
if (qcmp(q1,q2) == 0 || qiszero(q2))
|
|
return qlink(q1);
|
|
if (qiszero(q1))
|
|
return qlink(q2);
|
|
if (qisneg(q1)) {
|
|
q1tmp = qcomp(q1);
|
|
if (qisneg(q2)) {
|
|
q2tmp = qcomp(q2);
|
|
q = qand(q1tmp,q2tmp);
|
|
r = qcomp(q);
|
|
qfree(q1tmp);
|
|
qfree(q2tmp);
|
|
qfree(q);
|
|
return r;
|
|
}
|
|
q = qandnot(q1tmp, q2);
|
|
qfree(q1tmp);
|
|
r = qcomp(q);
|
|
qfree(q);
|
|
return r;
|
|
}
|
|
if (qisneg(q2)) {
|
|
q2tmp = qcomp(q2);
|
|
q = qandnot(q2tmp, q1);
|
|
qfree(q2tmp);
|
|
r = qcomp(q);
|
|
qfree(q);
|
|
return r;
|
|
}
|
|
r = qalloc();
|
|
zor(q1->num, q2->num, &r->num);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Perform the bitwise AND of two integers.
|
|
*/
|
|
NUMBER *
|
|
qand(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
register NUMBER *r;
|
|
NUMBER *q1tmp, *q2tmp, *q;
|
|
ZVALUE res;
|
|
|
|
if (qisfrac(q1) || qisfrac(q2)) {
|
|
math_error("Non-integers for bitwise and");
|
|
not_reached();
|
|
}
|
|
if (qcmp(q1, q2) == 0)
|
|
return qlink(q1);
|
|
if (qiszero(q1) || qiszero(q2))
|
|
return qlink(&_qzero_);
|
|
if (qisneg(q1)) {
|
|
q1tmp = qcomp(q1);
|
|
if (qisneg(q2)) {
|
|
q2tmp = qcomp(q2);
|
|
q = qor(q1tmp, q2tmp);
|
|
qfree(q1tmp);
|
|
qfree(q2tmp);
|
|
r = qcomp(q);
|
|
qfree(q);
|
|
return r;
|
|
}
|
|
r = qandnot(q2, q1tmp);
|
|
qfree(q1tmp);
|
|
return r;
|
|
}
|
|
if (qisneg(q2)) {
|
|
q2tmp = qcomp(q2);
|
|
r = qandnot(q1, q2tmp);
|
|
qfree(q2tmp);
|
|
return r;
|
|
}
|
|
zand(q1->num, q2->num, &res);
|
|
if (ziszero(res)) {
|
|
zfree(res);
|
|
return qlink(&_qzero_);
|
|
}
|
|
r = qalloc();
|
|
r->num = res;
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Perform the bitwise XOR of two integers.
|
|
*/
|
|
NUMBER *
|
|
qxor(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
register NUMBER *r;
|
|
NUMBER *q1tmp, *q2tmp, *q;
|
|
ZVALUE res;
|
|
|
|
if (qisfrac(q1) || qisfrac(q2)) {
|
|
math_error("Non-integers for bitwise xor");
|
|
not_reached();
|
|
}
|
|
if (qcmp(q1,q2) == 0)
|
|
return qlink(&_qzero_);
|
|
if (qiszero(q1))
|
|
return qlink(q2);
|
|
if (qiszero(q2))
|
|
return qlink(q1);
|
|
if (qisneg(q1)) {
|
|
q1tmp = qcomp(q1);
|
|
if (qisneg(q2)) {
|
|
q2tmp = qcomp(q2);
|
|
r = qxor(q1tmp, q2tmp);
|
|
qfree(q1tmp);
|
|
qfree(q2tmp);
|
|
return r;
|
|
}
|
|
q = qxor(q1tmp, q2);
|
|
qfree(q1tmp);
|
|
r = qcomp(q);
|
|
qfree(q);
|
|
return r;
|
|
}
|
|
if (qisneg(q2)) {
|
|
q2tmp = qcomp(q2);
|
|
q = qxor(q1, q2tmp);
|
|
qfree(q2tmp);
|
|
r = qcomp(q);
|
|
qfree(q);
|
|
return r;
|
|
}
|
|
zxor(q1->num, q2->num, &res);
|
|
if (ziszero(res)) {
|
|
zfree(res);
|
|
return qlink(&_qzero_);
|
|
}
|
|
r = qalloc();
|
|
r->num = res;
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Perform the bitwise AND-NOT of two integers.
|
|
*/
|
|
NUMBER *
|
|
qandnot(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
register NUMBER *r;
|
|
NUMBER *q1tmp, *q2tmp, *q;
|
|
|
|
if (qisfrac(q1) || qisfrac(q2)) {
|
|
math_error("Non-integers for bitwise xor");
|
|
not_reached();
|
|
}
|
|
if (qcmp(q1,q2) == 0 || qiszero(q1))
|
|
return qlink(&_qzero_);
|
|
if (qiszero(q2))
|
|
return qlink(q1);
|
|
if (qisneg(q1)) {
|
|
q1tmp = qcomp(q1);
|
|
if (qisneg(q2)) {
|
|
q2tmp = qcomp(q2);
|
|
r = qandnot(q2tmp, q1tmp);
|
|
qfree(q1tmp);
|
|
qfree(q2tmp);
|
|
return r;
|
|
}
|
|
q = qor(q1tmp, q2);
|
|
qfree(q1tmp);
|
|
r = qcomp(q);
|
|
qfree(q);
|
|
return r;
|
|
}
|
|
if (qisneg(q2)) {
|
|
q2tmp = qcomp(q2);
|
|
r = qand(q1, q2tmp);
|
|
qfree(q2tmp);
|
|
return r;
|
|
}
|
|
r = qalloc();
|
|
zandnot(q1->num, q2->num, &r->num);
|
|
return r;
|
|
}
|
|
|
|
/*
|
|
* Return the bitwise "complement" of a number. This is - q -1 if q is an
|
|
* integer, - q otherwise.
|
|
*/
|
|
NUMBER *
|
|
qcomp(NUMBER *q)
|
|
{
|
|
NUMBER *qtmp;
|
|
NUMBER *res;
|
|
|
|
if (qiszero(q))
|
|
return qlink(&_qnegone_);
|
|
if (qisnegone(q))
|
|
return qlink(&_qzero_);
|
|
qtmp = qneg(q);
|
|
if (qisfrac(q))
|
|
return qtmp;
|
|
res = qdec(qtmp);
|
|
qfree(qtmp);
|
|
return res;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the number whose binary representation only has the specified
|
|
* bit set (counting from zero). This thus produces a given power of two.
|
|
*/
|
|
NUMBER *
|
|
qbitvalue(long n)
|
|
{
|
|
register NUMBER *r;
|
|
|
|
if (n == 0)
|
|
return qlink(&_qone_);
|
|
r = qalloc();
|
|
if (n > 0)
|
|
zbitvalue(n, &r->num);
|
|
else
|
|
zbitvalue(-n, &r->den);
|
|
return r;
|
|
}
|
|
|
|
/*
|
|
* Return 10^n
|
|
*/
|
|
NUMBER *
|
|
qtenpow(long n)
|
|
{
|
|
register NUMBER *r;
|
|
|
|
if (n == 0)
|
|
return qlink(&_qone_);
|
|
r = qalloc();
|
|
if (n > 0)
|
|
ztenpow(n, &r->num);
|
|
else
|
|
ztenpow(-n, &r->den);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* Return the precision of a number (usually for examining an epsilon value).
|
|
* The precision of a number e less than 1 is the positive
|
|
* integer p for which e = 2^-p * f, where 1 <= f < 2.
|
|
* Numbers greater than or equal to one have a precision of zero.
|
|
* For example, the precision of e is 6 if 1/64 <= e < 1/32.
|
|
*/
|
|
long
|
|
qprecision(NUMBER *q)
|
|
{
|
|
long r;
|
|
|
|
if (qiszero(q) || qisneg(q)) {
|
|
math_error("Non-positive number for precision");
|
|
not_reached();
|
|
}
|
|
r = - qilog2(q);
|
|
return (r < 0 ? 0 : r);
|
|
}
|
|
|
|
|
|
/*
|
|
* Determine whether or not one number exactly divides another one.
|
|
* Returns true if the first number is an integer multiple of the second one.
|
|
*/
|
|
bool
|
|
qdivides(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
if (qiszero(q1))
|
|
return true;
|
|
if (qisint(q1) && qisint(q2)) {
|
|
if (qisunit(q2))
|
|
return true;
|
|
return zdivides(q1->num, q2->num);
|
|
}
|
|
return zdivides(q1->num, q2->num) && zdivides(q2->den, q1->den);
|
|
}
|
|
|
|
|
|
/*
|
|
* Compare two numbers and return an integer indicating their relative size.
|
|
* i = qrel(q1, q2);
|
|
*/
|
|
FLAG
|
|
qrel(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
ZVALUE z1, z2;
|
|
long wc1, wc2;
|
|
int sign;
|
|
int z1f = 0, z2f = 0;
|
|
|
|
if (q1 == q2)
|
|
return 0;
|
|
sign = q2->num.sign - q1->num.sign;
|
|
if (sign)
|
|
return sign;
|
|
if (qiszero(q2))
|
|
return !qiszero(q1);
|
|
if (qiszero(q1))
|
|
return -1;
|
|
/*
|
|
* Make a quick comparison by calculating the number of words
|
|
* resulting as if we multiplied through by the denominators,
|
|
* and then comparing the word counts.
|
|
*/
|
|
sign = 1;
|
|
if (qisneg(q1))
|
|
sign = -1;
|
|
wc1 = q1->num.len + q2->den.len;
|
|
wc2 = q2->num.len + q1->den.len;
|
|
if (wc1 < wc2 - 1)
|
|
return -sign;
|
|
if (wc2 < wc1 - 1)
|
|
return sign;
|
|
/*
|
|
* Quick check failed, must actually do the full comparison.
|
|
*/
|
|
if (zisunit(q2->den)) {
|
|
z1 = q1->num;
|
|
} else if (zisone(q1->num)) {
|
|
z1 = q2->den;
|
|
} else {
|
|
z1f = 1;
|
|
zmul(q1->num, q2->den, &z1);
|
|
}
|
|
if (zisunit(q1->den)) {
|
|
z2 = q2->num;
|
|
} else if (zisone(q2->num)) {
|
|
z2 = q1->den;
|
|
} else {
|
|
z2f = 1;
|
|
zmul(q2->num, q1->den, &z2);
|
|
}
|
|
sign = zrel(z1, z2);
|
|
if (z1f)
|
|
zfree(z1);
|
|
if (z2f)
|
|
zfree(z2);
|
|
return sign;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compare two numbers to see if they are equal.
|
|
* This differs from qrel in that the numbers are not ordered.
|
|
* Returns true if they differ.
|
|
*/
|
|
bool
|
|
qcmp(NUMBER *q1, NUMBER *q2)
|
|
{
|
|
if (q1 == q2)
|
|
return false;
|
|
if ((q1->num.sign != q2->num.sign) || (q1->num.len != q2->num.len) ||
|
|
(q1->den.len != q2->den.len) || (*q1->num.v != *q2->num.v) ||
|
|
(*q1->den.v != *q2->den.v))
|
|
return true;
|
|
if (zcmp(q1->num, q2->num))
|
|
return true;
|
|
if (qisint(q1))
|
|
return false;
|
|
return zcmp(q1->den, q2->den);
|
|
}
|
|
|
|
|
|
/*
|
|
* Compare a number against a normal small integer.
|
|
* Returns 1, 0, or -1, according to whether the first number is greater,
|
|
* equal, or less than the second number.
|
|
* res = qreli(q, n);
|
|
*/
|
|
FLAG
|
|
qreli(NUMBER *q, long n)
|
|
{
|
|
ZVALUE z1, z2;
|
|
FLAG res;
|
|
|
|
if (qiszero(q))
|
|
return ((n > 0) ? -1 : (n < 0));
|
|
|
|
if (n == 0)
|
|
return (q->num.sign ? -1 : 0);
|
|
|
|
if (q->num.sign != (n < 0))
|
|
return ((n < 0) ? 1 : -1);
|
|
|
|
itoz(n, &z1);
|
|
|
|
if (qisfrac(q)) {
|
|
zmul(q->den, z1, &z2);
|
|
zfree(z1);
|
|
z1 = z2;
|
|
}
|
|
|
|
res = zrel(q->num, z1);
|
|
zfree(z1);
|
|
return res;
|
|
}
|
|
|
|
|
|
/*
|
|
* Compare a number against a small integer to see if they are equal.
|
|
* Returns true if they differ.
|
|
*/
|
|
bool
|
|
qcmpi(NUMBER *q, long n)
|
|
{
|
|
long len;
|
|
#if LONG_BITS > BASEB
|
|
FULL nf;
|
|
#endif
|
|
|
|
len = q->num.len;
|
|
if (qisfrac(q) || (q->num.sign != (n < 0)))
|
|
return true;
|
|
if (n < 0)
|
|
n = -n;
|
|
if (((HALF)(n)) != q->num.v[0])
|
|
return true;
|
|
#if LONG_BITS > BASEB
|
|
nf = ((FULL) n) >> BASEB;
|
|
return ((nf != 0 || len > 1) && (len != 2 || nf != q->num.v[1]));
|
|
#else
|
|
return (len > 1);
|
|
#endif
|
|
}
|
|
|
|
|
|
/*
|
|
* Number node allocation routines
|
|
*/
|
|
|
|
#define NNALLOC 1000
|
|
|
|
|
|
STATIC NUMBER *freeNum = NULL;
|
|
STATIC NUMBER **firstNums = NULL;
|
|
STATIC long blockcount = 0;
|
|
|
|
|
|
NUMBER *
|
|
qalloc(void)
|
|
{
|
|
NUMBER *temp;
|
|
NUMBER ** newfn;
|
|
|
|
if (freeNum == NULL) {
|
|
freeNum = (NUMBER *) malloc(sizeof (NUMBER) * NNALLOC);
|
|
if (freeNum == NULL) {
|
|
math_error("Not enough memory");
|
|
not_reached();
|
|
}
|
|
freeNum[NNALLOC - 1].next = NULL;
|
|
freeNum[NNALLOC - 1].links = 0;
|
|
|
|
/*
|
|
* We prevent the temp pointer from walking behind freeNum
|
|
* by stopping one short of the end and running the loop one
|
|
* more time.
|
|
*
|
|
* We would stop the loop with just temp >= freeNum, but
|
|
* doing this helps make code checkers such as insure happy.
|
|
*/
|
|
for (temp = freeNum + NNALLOC - 2; temp > freeNum; --temp) {
|
|
temp->next = temp + 1;
|
|
temp->links = 0;
|
|
}
|
|
/* run the loop manually one last time */
|
|
temp->next = temp + 1;
|
|
temp->links = 0;
|
|
|
|
blockcount++;
|
|
if (firstNums == NULL) {
|
|
newfn = (NUMBER **) malloc(blockcount * sizeof(NUMBER *));
|
|
} else {
|
|
newfn = (NUMBER **)
|
|
realloc(firstNums, blockcount * sizeof(NUMBER *));
|
|
}
|
|
if (newfn == NULL) {
|
|
math_error("Cannot allocate new number block");
|
|
not_reached();
|
|
}
|
|
firstNums = newfn;
|
|
firstNums[blockcount - 1] = freeNum;
|
|
}
|
|
temp = freeNum;
|
|
freeNum = temp->next;
|
|
temp->links = 1;
|
|
temp->num = _one_;
|
|
temp->den = _one_;
|
|
return temp;
|
|
}
|
|
|
|
|
|
void
|
|
qfreenum(NUMBER *q)
|
|
{
|
|
if (q == NULL) {
|
|
math_error("Calling qfreenum with null argument!!!");
|
|
not_reached();
|
|
}
|
|
if (q->links != 0) {
|
|
math_error("Calling qfreenum with non-zero links!!!");
|
|
not_reached();
|
|
}
|
|
zfree(q->num);
|
|
zfree(q->den);
|
|
q->next = freeNum;
|
|
freeNum = q;
|
|
}
|
|
|
|
void
|
|
shownumbers(void)
|
|
{
|
|
NUMBER *vp;
|
|
long i, j, k;
|
|
long count = 0;
|
|
|
|
printf("Index Links Digits Value\n");
|
|
printf("----- ----- ------ -----\n");
|
|
|
|
for (i = 0, k = 0; initnumbs[i] != NULL; i++) {
|
|
count++;
|
|
vp = initnumbs[i];
|
|
printf("%6ld %4ld ", k++, vp->links);
|
|
fitprint(vp, 40);
|
|
printf("\n");
|
|
}
|
|
|
|
for (i = 0; i < blockcount; i++) {
|
|
vp = firstNums[i];
|
|
for (j = 0; j < NNALLOC; j++, k++, vp++) {
|
|
if (vp->links > 0) {
|
|
count++;
|
|
printf("%6ld %4ld ", k, vp->links);
|
|
fitprint(vp, 40);
|
|
printf("\n");
|
|
}
|
|
}
|
|
}
|
|
printf("\nNumber: %ld\n", count);
|
|
}
|