home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
ftp.cs.arizona.edu
/
ftp.cs.arizona.edu.tar
/
ftp.cs.arizona.edu
/
icon
/
historic
/
v941a.tgz
/
v941a.tar
/
icon.v941src
/
src
/
runtime
/
rlrgint.r
< prev
next >
Wrap
Text File
|
2002-01-18
|
53KB
|
2,318 lines
/*
* File: rlrgint.r
* Large integer arithmetic
*/
#ifdef LargeInts
extern int over_flow;
/*
* Conventions:
*
* Lrgints entering this module and leaving it are too large to
* be represented with T_Integer. So, externally, a given value
* is always T_Integer or always T_Lrgint.
*
* Routines outside this module operate on bignums by calling
* a routine like
*
* bigadd(da, db, dx)
*
* where da, db, and dx are pointers to tended descriptors.
* For the common case where one argument is a T_Integer, these
* call routines like
*
* bigaddi(da, IntVal(*db), dx).
*
* The bigxxxi routines can convert an integer to bignum form;
* they use itobig.
*
* The routines that actually do the work take (length, address)
* pairs specifying unsigned base-B digit strings. The sign handling
* is done in the bigxxx routines.
*/
/*
* Type for doing arithmetic on (2 * NB)-bit nonnegative numbers.
* Normally unsigned but may be signed (with NB reduced appropriately)
* if unsigned arithmetic is slow.
*/
/* The bignum radix, B */
#define B ((word)1 << NB)
/* Lrgint digits in a word */
#define WORDLEN (WordBits / NB + (WordBits % NB != 0))
/* size of a bignum block that will hold an integer */
#define INTBIGBLK sizeof(struct b_bignum) + sizeof(DIGIT) * WORDLEN
/* lo(uword d) : the low digit of a uword
hi(uword d) : the rest, d is unsigned
signed_hi(uword d) : the rest, d is signed
dbl(DIGIT a, DIGIT b) : the two-digit uword [a,b] */
#define lo(d) ((d) & (B - 1))
#define hi(d) ((uword)(d) >> NB)
#define dbl(a,b) (((uword)(a) << NB) + (b))
#if ((-1) >> 1) < 0
#define signed_hi(d) ((word)(d) >> NB)
#else
#define signbit ((uword)1 << (WordBits - NB - 1))
#define signed_hi(d) ((word)((((uword)(d) >> NB) ^ signbit) - signbit))
#endif
/* LrgInt(dptr dp) : the struct b_bignum pointed to by dp */
#define LrgInt(dp) ((struct b_bignum *)&BlkLoc(*dp)->bignumblk)
/* LEN(struct b_bignum *b) : number of significant digits */
#define LEN(b) ((b)->lsd - (b)->msd + 1)
/* DIG(struct b_bignum *b, word i): pointer to ith most significant digit */
/* (NOTE: This macro expansion often results in a very long string,
* that has been known to cause problems with the C compiler under VMS.
* So when DIG is used, keep it to one use per line.)
*/
#define DIG(b,i) (&(b)->digits[(b)->msd+(i)])
/* ceil, ln: ceil may be 1 too high in case ln is inaccurate */
#undef ceil
#define ceil(x) ((word)((x) + 1.01))
#define ln(n) (log((double)n))
/* determine the number of words needed for a bignum block with n digits */
#define LrgNeed(n) ( ((sizeof(struct b_bignum) + ((n) - 1) * sizeof(DIGIT)) \
+ WordSize - 1) & -WordSize )
/* copied from rconv.c */
#define tonum(c) (isdigit(c) ? (c)-'0' : 10+(((c)|(040))-'a'))
/* copied from oref.c */
#define RandVal (RanScale*(k_random=(RandA*(long)k_random+RandC)&0x7fffffffL))
/*
* Prototypes.
*/
static int mkdesc (struct b_bignum *x, dptr dx);
static void itobig (word i, struct b_bignum *x, dptr dx);
static void decout (FILE *f, DIGIT *n, word l);
static int bigaddi (dptr da, word i, dptr dx);
static int bigsubi (dptr da, word i, dptr dx);
static int bigmuli (dptr da, word i, dptr dx);
static int bigdivi (dptr da, word i, dptr dx);
static int bigmodi (dptr da, word i, dptr dx);
static int bigpowi (dptr da, word i, dptr dx);
static int bigpowii (word a, word i, dptr dx);
static word bigcmpi (dptr da, word i);
static DIGIT add1 (DIGIT *u, DIGIT *v, DIGIT *w, word n);
static word sub1 (DIGIT *u, DIGIT *v, DIGIT *w, word n);
static void mul1 (DIGIT *u, DIGIT *v, DIGIT *w, word n, word m);
static int div1
(DIGIT *a, DIGIT *b, DIGIT *q, DIGIT *r, word m, word n, struct b_bignum *b1, struct b_bignum *b2);
static void compl1 (DIGIT *u, DIGIT *w, word n);
static word cmp1 (DIGIT *u, DIGIT *v, word n);
static DIGIT addi1 (DIGIT *u, word k, DIGIT *w, word n);
static void subi1 (DIGIT *u, word k, DIGIT *w, word n);
static DIGIT muli1 (DIGIT *u, word k, int c, DIGIT *w, word n);
static DIGIT divi1 (DIGIT *u, word k, DIGIT *w, word n);
static DIGIT shifti1 (DIGIT *u, word k, DIGIT c, DIGIT *w, word n);
static word cmpi1 (DIGIT *u, word k, word n);
#define bdzero(dest,l) memset(dest, '\0', (l) * sizeof(DIGIT))
#define bdcopy(src, dest, l) memcpy(dest, src, (l) * sizeof(DIGIT))
/*
* mkdesc -- put value into a descriptor
*/
static int mkdesc(x, dx)
struct b_bignum *x;
dptr dx;
{
word xlen, cmp;
static DIGIT maxword[WORDLEN] = { 1 << ((WordBits - 1) % NB) };
/* suppress leading zero digits */
while (x->msd != x->lsd &&
*DIG(x,0) == 0)
x->msd++;
/* put it into a word if it fits, otherwise return the bignum */
xlen = LEN(x);
if (xlen < WORDLEN ||
(xlen == WORDLEN &&
((cmp = cmp1(DIG(x,0), maxword, (word)WORDLEN)) < 0 ||
(cmp == (word)0 && x->sign)))) {
word val = -(word)*DIG(x,0);
word i;
for (i = x->msd; ++i <= x->lsd; )
val = (val << NB) - x->digits[i];
if (!x->sign)
val = -val;
dx->dword = D_Integer;
IntVal(*dx) = val;
}
else {
dx->dword = D_Lrgint;
BlkLoc(*dx) = (union block *)x;
}
return Succeeded;
}
/*
* i -> big
*/
static void itobig(i, x, dx)
word i;
struct b_bignum *x;
dptr dx;
{
x->lsd = WORDLEN - 1;
x->msd = WORDLEN;
x->sign = 0;
if (i == 0) {
x->msd--;
*DIG(x,0) = 0;
}
else if (i < 0) {
word d = lo(i);
if (d != 0) {
d = B - d;
i += B;
}
i = - signed_hi(i);
x->msd--;
*DIG(x,0) = d;
x->sign = 1;
}
while (i != 0) {
x->msd--;
*DIG(x,0) = lo(i);
i = hi(i);
}
dx->dword = D_Lrgint;
BlkLoc(*dx) = (union block *)x;
}
/*
* string -> bignum
*/
word bigradix(sign, r, s, end_s, result)
int sign; /* '-' or not */
int r; /* radix 2 .. 36 */
char *s, *end_s; /* input string */
union numeric *result; /* output T_Integer or T_Lrgint */
{
struct b_bignum *b;
DIGIT *bd;
word len;
int c;
if (r == 0)
return CvtFail;
len = ceil((end_s - s) * ln(r) / ln(B));
Protect(b = alcbignum(len), return Error);
bd = DIG(b,0);
bdzero(bd, len);
if (r < 2 || r > 36)
return CvtFail;
for (c = ((s < end_s) ? *s++ : ' '); isalnum(c);
c = ((s < end_s) ? *s++ : ' ')) {
c = tonum(c);
if (c >= r)
return CvtFail;
muli1(bd, (word)r, c, bd, len);
}
/*
* Skip trailing white space and make sure there is nothing else left
* in the string. Note, if we have already reached end-of-string,
* c has been set to a space.
*/
while (isspace(c) && s < end_s)
c = *s++;
if (!isspace(c))
return CvtFail;
if (sign == '-')
b->sign = 1;
/* put value into dx and return the type */
{ struct descrip dx;
(void)mkdesc(b, &dx);
if (Type(dx) == T_Lrgint)
result->big = (struct b_bignum *)BlkLoc(dx);
else
result->integer = IntVal(dx);
return Type(dx);
}
}
/*
* bignum -> real
*/
double bigtoreal(da)
dptr da;
{
word i;
double r = 0;
struct b_bignum *b = &BlkLoc(*da)->bignumblk;
for (i = b->msd; i <= b->lsd; i++)
r = r * B + b->digits[i];
return (b->sign ? -r : r);
}
/*
* real -> bignum
*/
int realtobig(da, dx)
dptr da, dx;
{
#ifdef Double
double x;
#else /* Double */
double x = BlkLoc(*da)->realblk.realval;
#endif /* Double */
struct b_bignum *b;
word i, blen;
word d;
int sgn;
#ifdef Double
{
int *rp, *rq;
rp = (int *) &(BlkLoc(*da)->realblk.realval);
rq = (int *) &x;
*rq++ = *rp++;
*rq = *rp;
}
#endif /* Double */
if (x > 0.9999 * MinLong && x < 0.9999 * MaxLong) {
MakeInt((word)x, dx);
return Succeeded; /* got lucky; a simple integer suffices */
}
if (sgn = x < 0)
x = -x;
blen = ln(x) / ln(B) + 0.99;
for (i = 0; i < blen; i++)
x /= B;
if (x >= 1.0) {
x /= B;
blen += 1;
}
Protect(b = alcbignum(blen), return Error);
for (i = 0; i < blen; i++) {
d = (x *= B);
*DIG(b,i) = d;
x -= d;
}
b->sign = sgn;
return mkdesc(b, dx);
}
/*
* bignum -> string
*/
int bigtos(da, dx)
dptr da, dx;
{
tended struct b_bignum *a, *temp;
word alen = LEN(LrgInt(da));
word slen = ceil(alen * ln(B) / ln(10));
char *p, *q;
a = LrgInt(da);
Protect(temp = alcbignum(alen), fatalerr(0,NULL));
if (a->sign)
slen++;
Protect(q = alcstr(NULL,slen), fatalerr(0,NULL));
bdcopy(DIG(a,0),
DIG(temp,0),
alen);
p = q += slen;
#ifdef VMS
{
DIGIT *tempdg;
for (;;) {
tempdg = DIG(temp,0);
if (!cmpi1(tempdg, (word)0, alen))
break;
*--p = '0' + divi1(tempdg, (word)10, tempdg, alen);
}
}
#else
while (cmpi1(DIG(temp,0),
(word)0, alen))
*--p = '0' + divi1(DIG(temp,0),
(word)10,
DIG(temp,0),
alen);
#endif /* VMS */
if (a->sign)
*--p = '-';
StrLen(*dx) = q - p;
StrLoc(*dx) = p;
return NoCvt; /* The mnemonic is wrong, but the signal means */
/* that the string is allocated and not null- */
/* terminated. */
}
/*
* bignum -> file
*/
void bigprint(f, da)
FILE *f;
dptr da;
{
struct b_bignum *a, *temp;
word alen = LEN(LrgInt(da));
word slen, dlen;
struct b_bignum *blk = &BlkLoc(*da)->bignumblk;
slen = blk->lsd - blk->msd;
dlen = slen * NB * 0.3010299956639812 /* 1 / log2(10) */
+ log((double)blk->digits[blk->msd]) * 0.4342944819032518 + 0.5;
/* 1 / ln(10) */
if (dlen >= MaxDigits) {
fprintf(f, "integer(~10^%ld)",(long)dlen);
return;
}
/* not worth passing this one back */
Protect(temp = alcbignum(alen), fatalerr(0, NULL));
a = LrgInt(da);
bdcopy(DIG(a,0),
DIG(temp,0),
alen);
if (a->sign)
putc('-', f);
decout(f,
DIG(temp,0),
alen);
}
/*
* decout - given a base B digit string, print the number in base 10.
*/
static void decout(f, n, l)
FILE *f;
DIGIT *n;
word l;
{
DIGIT i = divi1(n, (word)10, n, l);
if (cmpi1(n, (word)0, l))
decout(f, n, l);
putc('0' + i, f);
}
/*
* da -> dx
*/
int cpbignum(da, dx)
dptr da, dx;
{
struct b_bignum *a, *x;
word alen = LEN(LrgInt(da));
Protect(x = alcbignum(alen), return Error);
a = LrgInt(da);
bdcopy(DIG(a,0),
DIG(x,0),
alen);
x->sign = a->sign;
return mkdesc(x, dx);
}
/*
* da + db -> dx
*/
int bigadd(da, db, dx)
dptr da, db;
dptr dx;
{
tended struct b_bignum *a, *b;
struct b_bignum *x;
word alen, blen;
word c;
if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
alen = LEN(LrgInt(da));
blen = LEN(LrgInt(db));
a = LrgInt(da);
b = LrgInt(db);
if (a->sign == b->sign) {
if (alen > blen) {
Protect(x = alcbignum(alen + 1), return Error);
c = add1(DIG(a,alen-blen),
DIG(b,0),
DIG(x,alen-blen+1),
blen);
*DIG(x,0) =
addi1(DIG(a,0),
c,
DIG(x,1),
alen-blen);
}
else if (alen == blen) {
Protect(x = alcbignum(alen + 1), return Error);
*DIG(x,0) =
add1(DIG(a,0),
DIG(b,0),
DIG(x,1),
alen);
}
else {
Protect(x = alcbignum(blen + 1), return Error);
c = add1(DIG(b,blen-alen),
DIG(a,0),
DIG(x,blen-alen+1),
alen);
*DIG(x,0) =
addi1(DIG(b,0),
c,
DIG(x,1),
blen-alen);
}
x->sign = a->sign;
}
else {
if (alen > blen) {
Protect(x = alcbignum(alen), return Error);
c = sub1(DIG(a,alen-blen),
DIG(b,0),
DIG(x,alen-blen),
blen);
subi1(DIG(a,0),
-c,
DIG(x,0),
alen-blen);
x->sign = a->sign;
}
else if (alen == blen) {
Protect(x = alcbignum(alen), return Error);
if (cmp1(DIG(a,0),
DIG(b,0),
alen) > 0) {
(void)sub1(DIG(a,0),
DIG(b,0),
DIG(x,0),
alen);
x->sign = a->sign;
}
else {
(void)sub1(DIG(b,0),
DIG(a,0),
DIG(x,0),
alen);
x->sign = b->sign;
}
}
else {
Protect(x = alcbignum(blen), return Error);
c = sub1(DIG(b,blen-alen),
DIG(a,0),
DIG(x,blen-alen),
alen);
subi1(DIG(b,0),
-c,
DIG(x,0),
blen-alen);
x->sign = b->sign;
}
}
return mkdesc(x, dx);
}
else if (Type(*da) == T_Lrgint) /* bignum + integer */
return bigaddi(da, IntVal(*db), dx);
else if (Type(*db) == T_Lrgint) /* integer + bignum */
return bigaddi(db, IntVal(*da), dx);
else { /* integer + integer */
struct descrip td;
char tdigits[INTBIGBLK];
itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
return bigaddi(&td, IntVal(*db), dx);
}
}
/*
* da - db -> dx
*/
int bigsub(da, db, dx)
dptr da, db, dx;
{
struct descrip td;
char tdigits[INTBIGBLK];
tended struct b_bignum *a, *b;
struct b_bignum *x;
word alen, blen;
word c;
if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
alen = LEN(LrgInt(da));
blen = LEN(LrgInt(db));
a = LrgInt(da);
b = LrgInt(db);
if (a->sign != b->sign) {
if (alen > blen) {
Protect(x = alcbignum(alen + 1), return Error);
c = add1(DIG(a,alen-blen),
DIG(b,0),
DIG(x,alen-blen+1),
blen);
*DIG(x,0) =
addi1(DIG(a,0),
c,
DIG(x,1),
alen-blen);
}
else if (alen == blen) {
Protect(x = alcbignum(alen + 1), return Error);
*DIG(x,0) =
add1(DIG(a,0),
DIG(b,0),
DIG(x,1),
alen);
}
else {
Protect(x = alcbignum(blen + 1), return Error);
c = add1(DIG(b,blen-alen),
DIG(a,0),
DIG(x,blen-alen+1),
alen);
*DIG(x,0) =
addi1(DIG(b,0),
c,
DIG(x,1),
blen-alen);
}
x->sign = a->sign;
}
else {
if (alen > blen) {
Protect(x = alcbignum(alen), return Error);
c = sub1(DIG(a,alen-blen),
DIG(b,0),
DIG(x,alen-blen),
blen);
subi1(DIG(a,0),
-c,
DIG(x,0),
alen-blen);
x->sign = a->sign;
}
else if (alen == blen) {
Protect(x = alcbignum(alen), return Error);
if (cmp1(DIG(a,0),
DIG(b,0),
alen) > 0) {
(void)sub1(DIG(a,0),
DIG(b,0),
DIG(x,0),
alen);
x->sign = a->sign;
}
else {
(void)sub1(DIG(b,0),
DIG(a,0),
DIG(x,0),
alen);
x->sign = 1 ^ b->sign;
}
}
else {
Protect(x = alcbignum(blen), return Error);
c = sub1(DIG(b,blen-alen),
DIG(a,0),
DIG(x,blen-alen),
alen);
subi1(DIG(b,0),
-c,
DIG(x,0),
blen-alen);
x->sign = 1 ^ b->sign;
}
}
return mkdesc(x, dx);
}
else if (Type(*da) == T_Lrgint) /* bignum - integer */
return bigsubi(da, IntVal(*db), dx);
else if (Type(*db) == T_Lrgint) { /* integer - bignum */
itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
alen = LEN(LrgInt(&td));
blen = LEN(LrgInt(db));
a = LrgInt(&td);
b = LrgInt(db);
if (a->sign != b->sign) {
if (alen == blen) {
Protect(x = alcbignum(alen + 1), return Error);
*DIG(x,0) =
add1(DIG(a,0),
DIG(b,0),
DIG(x,1),
alen);
}
else {
Protect(x = alcbignum(blen + 1), return Error);
c = add1(DIG(b,blen-alen),
DIG(a,0),
DIG(x,blen-alen+1),
alen);
*DIG(x,0) =
addi1(DIG(b,0),
c,
DIG(x,1),
blen-alen);
}
x->sign = a->sign;
}
else {
if (alen == blen) {
Protect(x = alcbignum(alen), return Error);
if (cmp1(DIG(a,0),
DIG(b,0),
alen) > 0) {
(void)sub1(DIG(a,0),
DIG(b,0),
DIG(x,0),
alen);
x->sign = a->sign;
}
else {
(void)sub1(DIG(b,0),
DIG(a,0),
DIG(x,0),
alen);
x->sign = 1 ^ b->sign;
}
}
else {
Protect(x = alcbignum(blen), return Error);
c = sub1(DIG(b,blen-alen),
DIG(a,0),
DIG(x,blen-alen),
alen);
subi1(DIG(b,0),
-c,
DIG(x,0),
blen-alen);
x->sign = 1 ^ b->sign;
}
}
return mkdesc(x, dx);
}
else { /* integer - integer */
itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
return bigsubi(&td, IntVal(*db), dx);
}
}
/*
* da * db -> dx
*/
int bigmul(da, db, dx)
dptr da, db, dx;
{
tended struct b_bignum *a, *b;
struct b_bignum *x;
word alen, blen;
if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
alen = LEN(LrgInt(da));
blen = LEN(LrgInt(db));
a = LrgInt(da);
b = LrgInt(db);
Protect(x = alcbignum(alen + blen), return Error);
mul1(DIG(a,0),
DIG(b,0),
DIG(x,0),
alen, blen);
x->sign = a->sign ^ b->sign;
return mkdesc(x, dx);
}
else if (Type(*da) == T_Lrgint) /* bignum * integer */
return bigmuli(da, IntVal(*db), dx);
else if (Type(*db) == T_Lrgint) /* integer * bignum */
return bigmuli(db, IntVal(*da), dx);
else { /* integer * integer */
struct descrip td;
char tdigits[INTBIGBLK];
itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
return bigmuli(&td, IntVal(*db), dx);
}
}
/*
* da / db -> dx
*/
int bigdiv(da, db, dx)
dptr da, db, dx;
{
tended struct b_bignum *a, *b, *x, *tu, *tv;
word alen, blen;
struct descrip td;
char tdigits[INTBIGBLK];
/* Put *da into large integer format. */
if (Type(*da) != T_Lrgint) {
itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
da = &td;
}
if (Type(*db) == T_Lrgint) { /* bignum / bignum */
alen = LEN(LrgInt(da));
blen = LEN(LrgInt(db));
if (alen < blen) {
MakeInt(0, dx);
return Succeeded;
}
a = LrgInt(da);
b = LrgInt(db);
Protect(x = alcbignum(alen - blen + 1), return Error);
if (blen == 1)
divi1(DIG(a,0),
(word)*DIG(b,0),
DIG(x,0),
alen);
else {
Protect(tu = alcbignum(alen + 1), return Error);
Protect(tv = alcbignum(blen), return Error);
if (div1(DIG(a,0),
DIG(b,0),
DIG(x,0),
NULL, alen-blen, blen, tu, tv) == Error)
return Error;
}
x->sign = a->sign ^ b->sign;
return mkdesc(x, dx);
}
else /* bignum / integer */
return bigdivi(da, IntVal(*db), dx);
}
/*
* da % db -> dx
*/
int bigmod(da, db, dx)
dptr da, db, dx;
{
tended struct b_bignum *a, *b, *x, *temp, *tu, *tv;
word alen, blen;
struct descrip td;
char tdigits[INTBIGBLK];
/* Put *da into large integer format. */
if (Type(*da) != T_Lrgint) {
itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
da = &td;
}
if (Type(*db) == T_Lrgint) { /* bignum % bignum */
alen = LEN(LrgInt(da));
blen = LEN(LrgInt(db));
if (alen < blen) {
cpbignum(da, dx);
return Succeeded;
}
a = LrgInt(da);
b = LrgInt(db);
Protect(x = alcbignum(blen), return Error);
if (blen == 1) {
Protect(temp = alcbignum(alen), return Error);
*DIG(x,0) =
divi1(DIG(a,0),
(word)*DIG(b,0),
DIG(temp,0),
alen);
}
else {
Protect(tu = alcbignum(alen + 1), return Error);
Protect(tv = alcbignum(blen), return Error);
if (div1(DIG(a,0),
DIG(b,0),
NULL,
DIG(x,0),
alen-blen, blen, tu, tv) == Error)
return Error;
}
x->sign = a->sign;
return mkdesc(x, dx);
}
else /* bignum % integer */
return bigmodi(da, IntVal(*db), dx);
}
/*
* -i -> dx
*/
int bigneg(da, dx)
dptr da, dx;
{
struct descrip td;
char tdigits[INTBIGBLK];
int cpstat;
/* Put *da into large integer format. */
if (Type(*da) != T_Lrgint) {
itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
da = &td;
}
LrgInt(da)->sign ^= 1; /* Temporarily change the sign */
cpstat = cpbignum(da, dx);
LrgInt(da)->sign ^= 1; /* Change it back */
return cpstat;
}
/*
* da ^ db -> dx
*/
int bigpow(da, db, dx)
dptr da, db, dx;
{
if (Type(*db) == T_Lrgint) {
struct b_bignum *b;
b = LrgInt ( db );
if (Type(*da) == T_Lrgint) {
if ( b->sign ) {
/* bignum ^ -bignum = 0 */
MakeInt ( 0, dx );
return Succeeded;
}
else
/* bignum ^ +bignum = guaranteed overflow */
ReturnErrNum(307, Error);
}
else if ( b->sign )
/* integer ^ -bignum */
switch ( IntVal ( *da ) ) {
case 1:
MakeInt ( 1, dx );
return Succeeded;
case -1:
/* Result is +1 / -1, depending on whether *b is even or odd. */
if ( ( b->digits[ b->lsd ] ) & 01 )
MakeInt ( -1, dx );
else
MakeInt ( 1, dx );
return Succeeded;
case 0:
ReturnErrNum(204,Error);
default:
/* da ^ (negative int) = 0 for all non-special cases */
MakeInt(0, dx);
return Succeeded;
}
else {
/* integer ^ +bignum */
word n, blen;
register DIGIT nth_dig, mask;
b = LrgInt ( db );
blen = LEN ( b );
/* We scan the bits of b from the most to least significant.
* The bit position in b is represented by the pair ( n, mask )
* where n is the DIGIT number (0 = most sig.) and mask is the
* the bit mask for the current bit.
*
* For each bit (most sig to least) in b,
* for each zero, square the partial result;
* for each one, square it and multiply it by a */
MakeInt ( 1, dx );
for ( n = 0; n < blen; ++n ) {
nth_dig = *DIG ( b, n );
for ( mask = 1 << ( NB - 1 ); mask; mask >>= 1 ) {
if ( bigmul ( dx, dx, dx ) == Error )
return Error;
if ( nth_dig & mask )
if ( bigmul ( dx, da, dx ) == Error )
return Error;
}
}
}
return Succeeded;
}
else if (Type(*da) == T_Lrgint) /* bignum ^ integer */
return bigpowi(da, IntVal(*db), dx);
else /* integer ^ integer */
return bigpowii(IntVal(*da), IntVal(*db), dx);
}
int bigpowri( a, db, drslt )
double a;
dptr db, drslt;
{
register double retval;
register word n;
register DIGIT nth_dig, mask;
struct b_bignum *b;
word blen;
b = LrgInt ( db );
blen = LEN ( b );
if ( b->sign ) {
if ( a == 0.0 )
ReturnErrNum(204, Error);
else
a = 1.0 / a;
}
/* We scan the bits of b from the most to least significant.
* The bit position in b is represented by the pair ( n, mask )
* where n is the DIGIT number (0 = most sig.) and mask is the
* the bit mask for the current bit.
*
* For each bit (most sig to least) in b,
* for each zero, square the partial result;
* for each one, square it and multiply it by a */
retval = 1.0;
for ( n = 0; n < blen; ++n ) {
nth_dig = *DIG ( b, n );
for ( mask = 1 << ( NB - 1 ); mask; mask >>= 1 ) {
retval *= retval;
if ( nth_dig & mask )
retval *= a;
}
}
Protect(BlkLoc(*drslt) = (union block *)alcreal(retval), return Error);
drslt->dword = D_Real;
return Succeeded;
}
/*
* iand(da, db) -> dx
*/
int bigand(da, db, dx)
dptr da, db, dx;
{
tended struct b_bignum *a, *b, *x, *tad, *tbd;
word alen, blen, xlen;
word i;
DIGIT *ad, *bd;
struct descrip td;
char tdigits[INTBIGBLK];
if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
alen = LEN(LrgInt(da));
blen = LEN(LrgInt(db));
xlen = alen > blen ? alen : blen;
a = LrgInt(da);
b = LrgInt(db);
Protect(x = alcbignum(xlen), return Error);
if (alen == xlen && !a->sign)
ad = DIG(a,0);
else {
Protect(tad = alcbignum(xlen), return Error);
ad = DIG(tad,0);
bdzero(ad, xlen - alen);
bdcopy(DIG(a,0),
&ad[xlen-alen], alen);
if (a->sign)
compl1(ad, ad, xlen);
}
if (blen == xlen && !b->sign)
bd = DIG(b,0);
else {
Protect(tbd = alcbignum(xlen), return Error);
bd = DIG(tbd,0);
bdzero(bd, xlen - blen);
bdcopy(DIG(b,0),
&bd[xlen-blen], blen);
if (b->sign)
compl1(bd, bd, xlen);
}
for (i = 0; i < xlen; i++)
*DIG(x,i) =
ad[i] & bd[i];
if (a->sign & b->sign) {
x->sign = 1;
compl1(DIG(x,0),
DIG(x,0),
xlen);
}
}
else if (Type(*da) == T_Lrgint) { /* iand(bignum,integer) */
itobig(IntVal(*db), (struct b_bignum *)tdigits, &td);
alen = LEN(LrgInt(da));
blen = LEN(LrgInt(&td));
xlen = alen > blen ? alen : blen;
a = LrgInt(da);
b = LrgInt(&td);
Protect(x = alcbignum(alen), return Error);
if (alen == xlen && !a->sign)
ad = DIG(a,0);
else {
Protect(tad = alcbignum(xlen), return Error);
ad = DIG(tad,0);
bdzero(ad, xlen - alen);
bdcopy(DIG(a,0),
&ad[xlen-alen], alen);
if (a->sign)
compl1(ad, ad, xlen);
}
if (blen == xlen && !b->sign)
bd = DIG(b,0);
else {
Protect(tbd = alcbignum(xlen), return Error);
bd = DIG(tbd,0);
bdzero(bd, xlen - blen);
bdcopy(DIG(b,0),
&bd[xlen-blen], blen);
if (b->sign)
compl1(bd, bd, xlen);
}
for (i = 0; i < xlen; i++)
*DIG(x,i) =
ad[i] & bd[i];
if (a->sign & b->sign) {
x->sign = 1;
compl1(DIG(x,0),
DIG(x,0),
xlen);
}
}
else if (Type(*db) == T_Lrgint) { /* iand(integer,bignum) */
itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
alen = LEN(LrgInt(&td));
blen = LEN(LrgInt(db));
xlen = alen > blen ? alen : blen;
a = LrgInt(&td);
b = LrgInt(db);
Protect(x = alcbignum(blen), return Error);
if (alen == xlen && !a->sign)
ad = DIG(a,0);
else {
Protect(tad = alcbignum(xlen), return Error);
ad = DIG(tad,0);
bdzero(ad, xlen - alen);
bdcopy(DIG(a,0),
&ad[xlen-alen], alen);
if (a->sign)
compl1(ad, ad, xlen);
}
if (blen == xlen && !b->sign)
bd = DIG(b,0);
else {
Protect(tbd = alcbignum(xlen), return Error);
bd = DIG(tbd,0);
bdzero(bd, xlen - blen);
bdcopy(DIG(b,0),
&bd[xlen-blen], blen);
if (b->sign)
compl1(bd, bd, xlen);
}
for (i = 0; i < xlen; i++)
*DIG(x,i) =
ad[i] & bd[i];
if (a->sign & b->sign) {
x->sign = 1;
compl1(DIG(x,0),
DIG(x,0),
xlen);
}
}
/* not called for iand(integer,integer) */
return mkdesc(x, dx);
}
/*
* ior(da, db) -> dx
*/
int bigor(da, db, dx)
dptr da, db, dx;
{
tended struct b_bignum *a, *b, *x, *tad, *tbd;
word alen, blen, xlen;
word i;
DIGIT *ad, *bd;
struct descrip td;
char tdigits[INTBIGBLK];
if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
alen = LEN(LrgInt(da));
blen = LEN(LrgInt(db));
xlen = alen > blen ? alen : blen;
a = LrgInt(da);
b = LrgInt(db);
Protect(x = alcbignum(xlen), return Error);
if (alen == xlen && !a->sign)
ad = DIG(a,0);
else {
Protect(tad = alcbignum(xlen), return Error);
ad = DIG(tad,0);
bdzero(ad, xlen - alen);
bdcopy(DIG(a,0),
&ad[xlen-alen], alen);
if (a->sign)
compl1(ad, ad, xlen);
}
if (blen == xlen && !b->sign)
bd = DIG(b,0);
else {
Protect(tbd = alcbignum(xlen), return Error);
bd = DIG(tbd,0);
bdzero(bd, xlen - blen);
bdcopy(DIG(b,0),
&bd[xlen-blen], blen);
if (b->sign)
compl1(bd, bd, xlen);
}
for (i = 0; i < xlen; i++)
*DIG(x,i) =
ad[i] | bd[i];
if (a->sign | b->sign) {
x->sign = 1;
compl1(DIG(x,0),
DIG(x,0),
xlen);
}
}
else if (Type(*da) == T_Lrgint) { /* ior(bignum,integer) */
itobig(IntVal(*db), (struct b_bignum *)tdigits, &td);
alen = LEN(LrgInt(da));
blen = LEN(LrgInt(&td));
xlen = alen > blen ? alen : blen;
a = LrgInt(da);
b = LrgInt(&td);
Protect(x = alcbignum(alen), return Error);
if (alen == xlen && !a->sign)
ad = DIG(a,0);
else {
Protect(tad = alcbignum(xlen), return Error);
ad = DIG(tad,0);
bdzero(ad, xlen - alen);
bdcopy(DIG(a,0),
&ad[xlen-alen], alen);
if (a->sign)
compl1(ad, ad, xlen);
}
if (blen == xlen && !b->sign)
bd = DIG(b,0);
else {
Protect(tbd = alcbignum(xlen), return Error);
bd = DIG(tbd,0);
bdzero(bd, xlen - blen);
bdcopy(DIG(b,0),
&bd[xlen-blen], blen);
if (b->sign)
compl1(bd, bd, xlen);
}
for (i = 0; i < xlen; i++)
*DIG(x,i) =
ad[i] | bd[i];
if (a->sign | b->sign) {
x->sign = 1;
compl1(DIG(x,0),
DIG(x,0),
xlen);
}
}
else if (Type(*db) == T_Lrgint) { /* ior(integer,bignym) */
itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
alen = LEN(LrgInt(&td));
blen = LEN(LrgInt(db));
xlen = alen > blen ? alen : blen;
a = LrgInt(&td);
b = LrgInt(db);
Protect(x = alcbignum(blen), return Error);
if (alen == xlen && !a->sign)
ad = DIG(a,0);
else {
Protect(tad = alcbignum(xlen), return Error);
ad = DIG(tad,0);
bdzero(ad, xlen - alen);
bdcopy(DIG(a,0),
&ad[xlen-alen], alen);
if (a->sign)
compl1(ad, ad, xlen);
}
if (blen == xlen && !b->sign)
bd = DIG(b,0);
else {
Protect(tbd = alcbignum(xlen), return Error);
bd = DIG(tbd,0);
bdzero(bd, xlen - blen);
bdcopy(DIG(b,0),
&bd[xlen-blen], blen);
if (b->sign)
compl1(bd, bd, xlen);
}
for (i = 0; i < xlen; i++)
*DIG(x,i) =
ad[i] | bd[i];
if (a->sign | b->sign) {
x->sign = 1;
compl1(DIG(x,0),
DIG(x,0),
xlen);
}
}
/* not called for ior(integer,integer) */
return mkdesc(x, dx);
}
/*
* xor(da, db) -> dx
*/
int bigxor(da, db, dx)
dptr da, db, dx;
{
tended struct b_bignum *a, *b, *x, *tad, *tbd;
word alen, blen, xlen;
word i;
DIGIT *ad, *bd;
struct descrip td;
char tdigits[INTBIGBLK];
if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
alen = LEN(LrgInt(da));
blen = LEN(LrgInt(db));
xlen = alen > blen ? alen : blen;
a = LrgInt(da);
b = LrgInt(db);
Protect(x = alcbignum(xlen), return Error);
if (alen == xlen && !a->sign)
ad = DIG(a,0);
else {
Protect(tad = alcbignum(xlen), return Error);
ad = DIG(tad,0);
bdzero(ad, xlen - alen);
bdcopy(DIG(a,0),
&ad[xlen-alen], alen);
if (a->sign)
compl1(ad, ad, xlen);
}
if (blen == xlen && !b->sign)
bd = DIG(b,0);
else {
Protect(tbd = alcbignum(xlen), return Error);
bd = DIG(tbd,0);
bdzero(bd, xlen - blen);
bdcopy(DIG(b,0),
&bd[xlen-blen], blen);
if (b->sign)
compl1(bd, bd, xlen);
}
for (i = 0; i < xlen; i++)
*DIG(x,i) =
ad[i] ^ bd[i];
if (a->sign ^ b->sign) {
x->sign = 1;
compl1(DIG(x,0),
DIG(x,0),
xlen);
}
}
else if (Type(*da) == T_Lrgint) { /* ixor(bignum,integer) */
itobig(IntVal(*db), (struct b_bignum *)tdigits, &td);
alen = LEN(LrgInt(da));
blen = LEN(LrgInt(&td));
xlen = alen > blen ? alen : blen;
a = LrgInt(da);
b = LrgInt(&td);
Protect(x = alcbignum(alen), return Error);
if (alen == xlen && !a->sign)
ad = DIG(a,0);
else {
Protect(tad = alcbignum(xlen), return Error);
ad = DIG(tad,0);
bdzero(ad, xlen - alen);
bdcopy(DIG(a,0),
&ad[xlen-alen], alen);
if (a->sign)
compl1(ad, ad, xlen);
}
if (blen == xlen && !b->sign)
bd = DIG(b,0);
else {
Protect(tbd = alcbignum(xlen), return Error);
bd = DIG(tbd,0);
bdzero(bd, xlen - blen);
bdcopy(DIG(b,0),
&bd[xlen-blen], blen);
if (b->sign)
compl1(bd, bd, xlen);
}
for (i = 0; i < xlen; i++)
*DIG(x,i) =
ad[i] ^ bd[i];
if (a->sign ^ b->sign) {
x->sign = 1;
compl1(DIG(x,0),
DIG(x,0),
xlen);
}
}
else if (Type(*db) == T_Lrgint) { /* ixor(integer,bignum) */
itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
alen = LEN(LrgInt(&td));
blen = LEN(LrgInt(db));
xlen = alen > blen ? alen : blen;
a = LrgInt(&td);
b = LrgInt(db);
Protect(x = alcbignum(blen), return Error);
if (alen == xlen && !a->sign)
ad = DIG(a,0);
else {
Protect(tad = alcbignum(xlen), return Error);
ad = DIG(tad,0);
bdzero(ad, xlen - alen);
bdcopy(DIG(a,0),
&ad[xlen-alen], alen);
if (a->sign)
compl1(ad, ad, xlen);
}
if (blen == xlen && !b->sign)
bd = DIG(b,0);
else {
Protect(tbd = alcbignum(xlen), return Error);
bd = DIG(tbd,0);
bdzero(bd, xlen - blen);
bdcopy(DIG(b,0),
&bd[xlen-blen], blen);
if (b->sign)
compl1(bd, bd, xlen);
}
for (i = 0; i < xlen; i++)
*DIG(x,i) =
ad[i] ^ bd[i];
if (a->sign ^ b->sign) {
x->sign = 1;
compl1(DIG(x,0),
DIG(x,0),
xlen);
}
}
/* not called for ixor(integer,integer) */
return mkdesc(x, dx);
}
/*
* bigshift(da, db) -> dx
*/
int bigshift(da, db, dx)
dptr da, db, dx;
{
tended struct b_bignum *a, *x, *tad;
word alen;
word r = IntVal(*db) % NB;
word q = (r >= 0 ? IntVal(*db) : (IntVal(*db) - (r += NB))) / NB;
word xlen;
DIGIT *ad;
struct descrip td;
char tdigits[INTBIGBLK];
if (Type(*da) == T_Integer) {
itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
da = &td;
}
alen = LEN(LrgInt(da));
xlen = alen + q + 1;
if (xlen <= 0) {
MakeInt(-LrgInt(da)->sign, dx);
return Succeeded;
}
else {
a = LrgInt(da);
Protect(x = alcbignum(xlen), return Error);
if (a->sign) {
Protect(tad = alcbignum(alen), return Error);
ad = DIG(tad,0);
bdcopy(DIG(a,0),
ad, alen);
compl1(ad, ad, alen);
}
else
ad = DIG(a,0);
if (q >= 0) {
*DIG(x,0) =
shifti1(ad, r, (DIGIT)0,
DIG(x,1),
alen);
bdzero(DIG(x,alen+1),
q);
}
else
*DIG(x,0) =
shifti1(ad, r, ad[alen+q] >> (NB-r),
DIG(x,1), alen+q);
if (a->sign) {
x->sign = 1;
*DIG(x,0) |=
B - (1 << r);
compl1(DIG(x,0),
DIG(x,0),
xlen);
}
return mkdesc(x, dx);
}
}
/*
* negative if da < db
* zero if da == db
* positive if da > db
*/
word bigcmp(da, db)
dptr da, db;
{
struct b_bignum *a = LrgInt(da);
struct b_bignum *b = LrgInt(db);
word alen, blen;
if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
if (a->sign != b->sign)
return (b->sign - a->sign);
alen = LEN(a);
blen = LEN(b);
if (alen != blen)
return (a->sign ? blen - alen : alen - blen);
if (a->sign)
return cmp1(DIG(b,0),
DIG(a,0),
alen);
else
return cmp1(DIG(a,0),
DIG(b,0),
alen);
}
else if (Type(*da) == T_Lrgint) /* cmp(bignum, integer) */
return bigcmpi(da, IntVal(*db));
else /* cmp(integer, bignum) */
return -bigcmpi(db, IntVal(*da));
}
/*
* ?da -> dx
*/
int bigrand(da, dx)
dptr da, dx;
{
tended struct b_bignum *x, *a, *td, *tu, *tv;
word alen = LEN(LrgInt(da));
DIGIT *d;
word i;
double rval;
Protect(x = alcbignum(alen), return Error);
Protect(td = alcbignum(alen + 1), return Error);
d = DIG(td,0);
a = LrgInt(da);
for (i = alen; i >= 0; i--) {
rval = RandVal;
d[i] = rval * B;
}
Protect(tu = alcbignum(alen + 2), return Error);
Protect(tv = alcbignum(alen), return Error);
if (div1(d, DIG(a,0),
NULL,
DIG(x,0),
(word)1, alen, tu, tv) == Error)
return Error;
addi1(DIG(x,0),
(word)1,
DIG(x,0),
alen);
return mkdesc(x, dx);
}
/*
* da + i -> dx
*/
static int bigaddi(da, i, dx)
dptr da, dx;
word i;
{
tended struct b_bignum *a;
struct b_bignum *x;
word alen;
if (i < 0 && i > MinLong)
return bigsubi(da, -i, dx);
else if (i < 0 || i >= B ) {
struct descrip td;
char tdigits[INTBIGBLK];
itobig(i, (struct b_bignum *)tdigits, &td);
return bigadd(da, &td, dx);
}
else {
alen = LEN(LrgInt(da));
a = LrgInt(da);
if (a->sign) {
Protect(x = alcbignum(alen), return Error);
subi1(DIG(a,0),
i,
DIG(x,0),
alen);
}
else {
Protect(x = alcbignum(alen + 1), return Error);
*DIG(x,0) =
addi1(DIG(a,0),
i,
DIG(x,1),
alen);
}
x->sign = a->sign;
return mkdesc(x, dx);
}
}
/*
* da - i -> dx
*/
static int bigsubi(da, i, dx)
dptr da, dx;
word i;
{
tended struct b_bignum *a;
struct b_bignum *x;
word alen;
if (i < 0 && i > MinLong)
return bigaddi(da, -i, dx);
else if (i < 0 || i >= B) {
struct descrip td;
char tdigits[INTBIGBLK];
itobig(i, (struct b_bignum *)tdigits, &td);
return bigsub(da, &td, dx);
}
else {
alen = LEN(LrgInt(da));
a = LrgInt(da);
if (a->sign) {
Protect(x = alcbignum(alen + 1), return Error);
*DIG(x,0) =
addi1(DIG(a,0),
i,
DIG(x,1),
alen);
}
else {
Protect(x = alcbignum(alen), return Error);
subi1(DIG(a,0),
i,
DIG(x,0),
alen);
}
x->sign = a->sign;
return mkdesc(x, dx);
}
}
/*
* da * i -> dx
*/
static int bigmuli(da, i, dx)
dptr da, dx;
word i;
{
tended struct b_bignum *a;
struct b_bignum *x;
word alen;
if (i <= -B || i >= B) {
struct descrip td;
char tdigits[INTBIGBLK];
itobig(i, (struct b_bignum *)tdigits, &td);
return bigmul(da, &td, dx);
}
else {
alen = LEN(LrgInt(da));
a = LrgInt(da);
Protect(x = alcbignum(alen + 1), return Error);
if (i >= 0)
x->sign = a->sign;
else {
x->sign = 1 ^ a->sign;
i = -i;
}
*DIG(x,0) =
muli1(DIG(a,0),
i, 0,
DIG(x,1),
alen);
return mkdesc(x, dx);
}
}
/*
* da / i -> dx
*/
static int bigdivi(da, i, dx)
dptr da, dx;
word i;
{
tended struct b_bignum *a;
struct b_bignum *x;
word alen;
if (i <= -B || i >= B) {
struct descrip td;
char tdigits[INTBIGBLK];
itobig(i, (struct b_bignum *)tdigits, &td);
return bigdiv(da, &td, dx);
}
else {
alen = LEN(LrgInt(da));
a = LrgInt(da);
Protect(x = alcbignum(alen), return Error);
if (i >= 0)
x->sign = a->sign;
else {
x->sign = 1 ^ a->sign;
i = -i;
}
divi1(DIG(a,0),
i,
DIG(x,0),
alen);
return mkdesc(x, dx);
}
}
/*
* da % i -> dx
*/
static int bigmodi(da, i, dx)
dptr da, dx;
word i;
{
tended struct b_bignum *a, *temp;
word alen;
word x;
if (i <= -B || i >= B) {
struct descrip td;
char tdigits[INTBIGBLK];
itobig(i, (struct b_bignum *)tdigits, &td);
return bigmod(da, &td, dx);
}
else {
alen = LEN(LrgInt(da));
a = LrgInt(da);
temp = a; /* avoid trash pointer */
Protect(temp = alcbignum(alen), return Error);
x = divi1(DIG(a,0),
Abs(i),
DIG(temp,0),
alen);
if (a->sign)
x = -x;
MakeInt(x, dx);
return Succeeded;
}
}
/*
* da ^ i -> dx
*/
static int bigpowi(da, i, dx)
dptr da, dx;
word i;
{
int n = WordBits;
if (i > 0) {
/* scan bits left to right. skip leading 1. */
while (--n >= 0)
if (i & ((word)1 << n))
break;
/* then, for each zero, square the partial result;
for each one, square it and multiply it by a */
*dx = *da;
while (--n >= 0) {
if (bigmul(dx, dx, dx) == Error)
return Error;
if (i & ((word)1 << n))
if (bigmul(dx, da, dx) == Error)
return Error;
}
}
else if (i == 0) {
MakeInt(1, dx);
}
else {
MakeInt(0, dx);
}
return Succeeded;
}
/*
* a ^ i -> dx
*/
static int bigpowii(a, i, dx)
word a, i;
dptr dx;
{
word x, y;
int n = WordBits;
int isbig = 0;
if (a == 0 || i <= 0) { /* special cases */
if (a == 0 && i <= 0) /* 0 ^ negative -> error */
ReturnErrNum(204,Error);
if (i == 0) {
MakeInt(1, dx);
return Succeeded;
}
if (a == -1) { /* -1 ^ [odd,even] -> [-1,+1] */
if (!(i & 1))
a = 1;
}
else if (a != 1) { /* 1 ^ any -> 1 */
a = 0;
} /* others ^ negative -> 0 */
MakeInt(a, dx);
}
else {
struct descrip td;
char tdigits[INTBIGBLK];
/* scan bits left to right. skip leading 1. */
while (--n >= 0)
if (i & ((word)1 << n))
break;
/* then, for each zero, square the partial result;
for each one, square it and multiply it by a */
x = a;
while (--n >= 0) {
if (isbig) {
if (bigmul(dx, dx, dx) == Error)
return Error;
}
else {
y = mul(x, x);
if (!over_flow)
x = y;
else {
itobig(x, (struct b_bignum *)tdigits, &td);
if (bigmul(&td, &td, dx) == Error)
return Error;
isbig = (Type(*dx) == T_Lrgint);
}
}
if (i & ((word)1 << n)) {
if (isbig) {
if (bigmuli(dx, a, dx) == Error)
return Error;
}
else {
y = mul(x, a);
if (!over_flow)
x = y;
else {
itobig(x, (struct b_bignum *)tdigits, &td);
if (bigmuli(&td, a, dx) == Error)
return Error;
isbig = (Type(*dx) == T_Lrgint);
}
}
}
}
if (!isbig) {
MakeInt(x, dx);
}
}
return Succeeded;
}
/*
* negative if da < i
* zero if da == i
* positive if da > i
*/
static word bigcmpi(da, i)
dptr da;
word i;
{
struct b_bignum *a = LrgInt(da);
word alen = LEN(a);
if (i > -B && i < B) {
if (i >= 0)
if (a->sign)
return -1;
else
return cmpi1(DIG(a,0),
i, alen);
else
if (a->sign)
return -cmpi1(DIG(a,0),
-i, alen);
else
return 1;
}
else {
struct descrip td;
char tdigits[INTBIGBLK];
itobig(i, (struct b_bignum *)tdigits, &td);
return bigcmp(da, &td);
}
}
/* These are all straight out of Knuth vol. 2, Sec. 4.3.1. */
/*
* (u,n) + (v,n) -> (w,n)
*
* returns carry, 0 or 1
*/
static DIGIT add1(u, v, w, n)
DIGIT *u, *v, *w;
word n;
{
uword dig, carry;
word i;
carry = 0;
for (i = n; --i >= 0; ) {
dig = (uword)u[i] + v[i] + carry;
w[i] = lo(dig);
carry = hi(dig);
}
return carry;
}
/*
* (u,n) - (v,n) -> (w,n)
*
* returns carry, 0 or -1
*/
static word sub1(u, v, w, n)
DIGIT *u, *v, *w;
word n;
{
uword dig, carry;
word i;
carry = 0;
for (i = n; --i >= 0; ) {
dig = (uword)u[i] - v[i] + carry;
w[i] = lo(dig);
carry = signed_hi(dig);
}
return carry;
}
/*
* (u,n) * (v,m) -> (w,m+n)
*/
static void mul1(u, v, w, n, m)
DIGIT *u, *v, *w;
word n, m;
{
word i, j;
uword dig, carry;
bdzero(&w[m], n);
for (j = m; --j >= 0; ) {
carry = 0;
for (i = n; --i >= 0; ) {
dig = (uword)u[i] * v[j] + w[i+j+1] + carry;
w[i+j+1] = lo(dig);
carry = hi(dig);
}
w[j] = carry;
}
}
/*
* (a,m+n) / (b,n) -> (q,m+1) (r,n)
*
* if q or r is NULL, the quotient or remainder is discarded
*/
static int div1(a, b, q, r, m, n, tu, tv)
DIGIT *a, *b, *q, *r;
word m, n;
struct b_bignum *tu, *tv;
{
uword qhat, rhat;
uword dig, carry;
DIGIT *u, *v;
word d;
word i, j;
u = DIG(tu,0);
v = DIG(tv,0);
/* D1 */
for (d = 0; d < NB; d++)
if (b[0] & (1 << (NB - 1 - d)))
break;
u[0] = shifti1(a, d, (DIGIT)0, &u[1], m+n);
shifti1(b, d, (DIGIT)0, v, n);
/* D2, D7 */
for (j = 0; j <= m; j++) {
/* D3 */
if (u[j] == v[0]) {
qhat = B - 1;
rhat = (uword)v[0] + u[j+1];
}
else {
uword numerator = dbl(u[j], u[j+1]);
qhat = numerator / (uword)v[0];
rhat = numerator % (uword)v[0];
}
while (rhat < (uword)B && qhat * (uword)v[1] > (uword)dbl(rhat, u[j+2])) {
qhat -= 1;
rhat += v[0];
}
/* D4 */
carry = 0;
for (i = n; i > 0; i--) {
dig = u[i+j] - v[i-1] * qhat + carry; /* -BSQ+B .. B-1 */
u[i+j] = lo(dig);
if ((uword)dig < (uword)B)
carry = hi(dig);
else carry = hi(dig) | -B;
}
carry = (word)(carry + u[j]) < 0;
/* D5 */
if (q)
q[j] = qhat;
/* D6 */
if (carry) {
if (q)
q[j] -= 1;
carry = 0;
for (i = n; i > 0; i--) {
dig = (uword)u[i+j] + v[i-1] + carry;
u[i+j] = lo(dig);
carry = hi(dig);
}
}
}
if (r) {
if (d == 0)
shifti1(&u[m+1], (word)d, (DIGIT)0, r, n);
else
r[0] = shifti1(&u[m+1], (word)(NB - d), u[m+n]>>d, &r[1], n - 1);
}
return Succeeded;
}
/*
* - (u,n) -> (w,n)
*
*/
static void compl1(u, w, n)
DIGIT *u, *w;
word n;
{
uword dig, carry = 0;
word i;
for (i = n; --i >= 0; ) {
dig = carry - u[i];
w[i] = lo(dig);
carry = signed_hi(dig);
}
}
/*
* (u,n) : (v,n)
*/
static word cmp1(u, v, n)
DIGIT *u, *v;
word n;
{
word i;
for (i = 0; i < n; i++)
if (u[i] != v[i])
return u[i] > v[i] ? 1 : -1;
return 0;
}
/*
* (u,n) + k -> (w,n)
*
* k in 0 .. B-1
* returns carry, 0 or 1
*/
static DIGIT addi1(u, k, w, n)
DIGIT *u, *w;
word k;
word n;
{
uword dig, carry;
word i;
carry = k;
for (i = n; --i >= 0; ) {
dig = (uword)u[i] + carry;
w[i] = lo(dig);
carry = hi(dig);
}
return carry;
}
/*
* (u,n) - k -> (w,n)
*
* k in 0 .. B-1
* u must be greater than k
*/
static void subi1(u, k, w, n)
DIGIT *u, *w;
word k;
word n;
{
uword dig, carry;
word i;
carry = -k;
for (i = n; --i >= 0; ) {
dig = (uword)u[i] + carry;
w[i] = lo(dig);
carry = signed_hi(dig);
}
}
/*
* (u,n) * k + c -> (w,n)
*
* k in 0 .. B-1
* returns carry, 0 .. B-1
*/
static DIGIT muli1(u, k, c, w, n)
DIGIT *u, *w;
word k;
int c;
word n;
{
uword dig, carry;
word i;
carry = c;
for (i = n; --i >= 0; ) {
dig = (uword)k * u[i] + carry;
w[i] = lo(dig);
carry = hi(dig);
}
return carry;
}
/*
* (u,n) / k -> (w,n)
*
* k in 0 .. B-1
* returns remainder, 0 .. B-1
*/
static DIGIT divi1(u, k, w, n)
DIGIT *u, *w;
word k;
word n;
{
uword dig, remain;
word i;
remain = 0;
for (i = 0; i < n; i++) {
dig = dbl(remain, u[i]);
w[i] = dig / k;
remain = dig % k;
}
return remain;
}
/*
* ((u,n) << k) + c -> (w,n)
*
* k in 0 .. NB-1
* c in 0 .. B-1
* returns carry, 0 .. B-1
*/
static DIGIT shifti1(u, k, c, w, n)
DIGIT *u, c, *w;
word k;
word n;
{
uword dig;
word i;
if (k == 0) {
bdcopy(u, w, n);
return 0;
}
for (i = n; --i >= 0; ) {
dig = ((uword)u[i] << k) + c;
w[i] = lo(dig);
c = hi(dig);
}
return c;
}
/*
* (u,n) : k
*
* k in 0 .. B-1
*/
static word cmpi1(u, k, n)
DIGIT *u;
word k;
word n;
{
word i;
for (i = 0; i < n-1; i++)
if (u[i])
return 1;
if (u[n - 1] == (DIGIT)k)
return 0;
return u[n - 1] > (DIGIT)k ? 1 : -1;
}
#else /* LargeInts */
static char x; /* prevent empty module */
#endif /* LargeInts */