summaryrefslogtreecommitdiff
path: root/big.c
diff options
context:
space:
mode:
Diffstat (limited to 'big.c')
-rw-r--r--big.c1377
1 files changed, 829 insertions, 548 deletions
diff --git a/big.c b/big.c
index 05f346b..4b20d04 100644
--- a/big.c
+++ b/big.c
@@ -14,245 +14,371 @@
#include "big.h"
#include <errno.h>
-static double logIBASE,log10IBASE;
+/* contains remainder from last call to longdiv or shortdiv */
+word b_rem;
+
+extern char *dicp;
+
+#define maxval (1ll<<60)
+
+static double logIBASE, log10IBASE;
word big_one;
-static word big_plus(word,word,int);
-static word big_sub(word,word);
+static word big_plus(word, word, int);
+static word big_sub(word, word);
static word len(word);
-static word longdiv(word,word);
+static word longdiv(word, word);
static word ms2d(word);
static word msd(word);
-static word shift(word,word);
-static word shortdiv(word,word);
-static word stimes(word,word);
-
-void bigsetup()
-{ logIBASE=log((double)IBASE);
- log10IBASE=log10((double)IBASE);
- big_one=make(INT,1,0);
-}
-
-int isnat(x)
-word x;
-{ return(tag[x]==INT&&poz(x));
-}
-
-word sto_int(i) /* store C long long as mira bigint */
-long long i;
-{ word s,x;
- if(i<0)s=SIGNBIT,i= -i; else s=0;
- x=make(INT,s|i&MAXDIGIT,0);
- if(i>>=DIGITWIDTH)
- { word *p = &rest(x);
- *p=make(INT,i&MAXDIGIT,0),p= &rest(*p);
- while(i>>=DIGITWIDTH)
- *p=make(INT,i&MAXDIGIT,0),p= &rest(*p); }
- return(x);
-} /* change to long long, DT Oct 2019 */
+static word shift(word, word);
+static word shortdiv(word, word);
+static word stimes(word, word);
+
+void bigsetup(void)
+{
+ logIBASE = log((double)IBASE);
+ log10IBASE = log10((double)IBASE);
+ big_one = make(INT, 1, 0);
+}
-#define maxval (1ll<<60)
+int isnat(word x)
+{
+ return (tag[x] == INT && poz(x));
+}
-long long get_int(x) /* mira bigint to C long long */
-word x;
-{ long long n=digit0(x);
- word sign=neg(x);
- if(!(x=rest(x)))return(sign?-n:n);
-{ word w=DIGITWIDTH;
- while(x&&w<60)n+=(long long)digit(x)<<w,w+=DIGITWIDTH,x=rest(x);
- if(x)n=maxval; /* overflow, return large value */
- return(sign?-n:n);
-}} /* change to long long, DT Oct 2019 */
-
-word bignegate(x)
-word x;
-{ if(bigzero(x))return(x);
- return(make(INT,hd[x]&SIGNBIT?hd[x]&MAXDIGIT:SIGNBIT|hd[x],tl[x]));
-}
-
-word bigplus(x,y)
-word x,y;
-{ if(poz(x))
- if(poz(y))return(big_plus(x,y,0));
- else return(big_sub(x,y));
- else
- if(poz(y))return(big_sub(y,x));
- else return(big_plus(x,y,SIGNBIT)); /* both negative */
-}
-
-word big_plus(x,y,signbit) /* ignore input signs, treat x,y as positive */
-word x,y; int signbit;
-{ word d=digit0(x)+digit0(y);
- word carry = ((d&IBASE)!=0);
- word r = make(INT,signbit|d&MAXDIGIT,0); /* result */
- word *z = &rest(r); /* pointer to rest of result */
- x = rest(x); y = rest(y);
- while(x&&y) /* this loop has been unwrapped once, see above */
- { d = carry+digit(x)+digit(y);
- carry = ((d&IBASE)!=0);
- *z = make(INT,d&MAXDIGIT,0);
- x = rest(x); y = rest(y); z = &rest(*z); }
- if(y)x=y; /* by convention x is the longer one */
- while(x)
- { d = carry+digit(x);
- carry = ((d&IBASE)!=0);
- *z = make(INT,d&MAXDIGIT,0);
- x = rest(x); z = &rest(*z); }
- if(carry)*z=make(INT,1,0);
- return(r);
+/* store C long long as mira bigint */
+word sto_int(long long i)
+{
+ word s, x;
+ if (i < 0) {
+ s = SIGNBIT, i = -i;
+ } else {
+ s = 0;
+ }
+ x = make(INT, (s | i) & MAXDIGIT, 0);
+ if (i >>= DIGITWIDTH) {
+ word *p = &rest(x);
+ *p = make(INT, i & MAXDIGIT, 0), p = &rest(*p);
+ while (i >>= DIGITWIDTH) {
+ *p = make(INT, i & MAXDIGIT, 0), p = &rest(*p);
+ }
+ }
+ return (x);
+}
+
+/* mira bigint to C long long */
+long long get_int(word x)
+{
+ long long n = digit0(x);
+ word sign = neg(x);
+ if (!(x = rest(x))) {
+ return (sign ? -n : n);
+ }
+ word w = DIGITWIDTH;
+ while (x && w < 60) {
+ n += (long long)digit(x) << w, w += DIGITWIDTH, x = rest(x);
+ }
+ if (x) {
+ n = maxval; /* overflow, return large value */
+ }
+ return (sign ? -n : n);
+}
+
+word bignegate(word x)
+{
+ if (bigzero(x)) {
+ return (x);
+ }
+ return (make
+ (INT, hd[x] & SIGNBIT ? hd[x] & MAXDIGIT : SIGNBIT | hd[x], tl[x]));
+}
+
+word bigplus(word x, word y)
+{
+ if (poz(x)) {
+ if (poz(y)) {
+ return (big_plus(x, y, 0));
+ } else {
+ return (big_sub(x, y));
+ }
+ } else if (poz(y)) {
+ return (big_sub(y, x));
+ } else {
+ return (big_plus(x, y, SIGNBIT)); /* both negative */
+ }
+}
+
+/* ignore input signs, treat x,y as positive */
+word big_plus(word x, word y, int signbit)
+{
+ word d = digit0(x) + digit0(y);
+ word carry = ((d & IBASE) != 0);
+ word r = make(INT, (signbit | d) & MAXDIGIT, 0); /* result */
+ word *z = &rest(r); /* pointer to rest of result */
+ x = rest(x);
+ y = rest(y);
+ while (x && y) { /* this loop has been unwrapped once, see above */
+ d = carry + digit(x) + digit(y);
+ carry = ((d & IBASE) != 0);
+ *z = make(INT, d & MAXDIGIT, 0);
+ x = rest(x);
+ y = rest(y);
+ z = &rest(*z);
+ }
+ if (y) {
+ x = y; /* by convention x is the longer one */
+ }
+ while (x) {
+ d = carry + digit(x);
+ carry = ((d & IBASE) != 0);
+ *z = make(INT, d & MAXDIGIT, 0);
+ x = rest(x);
+ z = &rest(*z);
+ }
+ if (carry) {
+ *z = make(INT, 1, 0);
+ }
+ return (r);
}
-word bigsub(x,y)
-word x,y;
-{ if(poz(x))
- if(poz(y))return(big_sub(x,y));
- else return(big_plus(x,y,0)); /* poz x, negative y */
- else
- if(poz(y))return(big_plus(x,y,SIGNBIT)); /* negative x, poz y */
- else return(big_sub(y,x)); /* both negative */
+word bigsub(word x, word y)
+{
+ if (poz(x)) {
+ if (poz(y)) {
+ return (big_sub(x, y));
+ } else {
+ return (big_plus(x, y, 0)); /* poz x, negative y */
+ }
+ } else if (poz(y)) {
+ return (big_plus(x, y, SIGNBIT)); /* negative x, poz y */
+ } else {
+ return (big_sub(y, x)); /* both negative */
+ }
}
-word big_sub(x,y) /* ignore input signs, treat x,y as positive */
-word x,y;
-{ word d = digit0(x)-digit0(y);
- word borrow = (d&IBASE)!=0;
- word r=make(INT,d&MAXDIGIT,0); /* result */
+/* ignore input signs, treat x,y as positive */
+word big_sub(word x, word y)
+{
+ word d = digit0(x) - digit0(y);
+ word borrow = (d & IBASE) != 0;
+ word r = make(INT, d & MAXDIGIT, 0); /* result */
word *z = &rest(r);
- word *p=NULL; /* pointer to trailing zeros, if any */
- x = rest(x); y = rest(y);
- while(x&&y) /* this loop has been unwrapped once, see above */
- { d = digit(x)-digit(y)-borrow;
- borrow = (d&IBASE)!=0;
- d = d&MAXDIGIT;
- *z = make(INT,d,0);
- if(d)p=NULL; else if(!p)p=z;
- x = rest(x); y = rest(y); z = &rest(*z); }
- while(y) /* at most one of these two loops will be invoked */
- { d = -digit(y)-borrow;
- borrow = ((d&IBASE)!=0);
- d = d&MAXDIGIT;
- *z = make(INT,d,0);
- if(d)p=NULL; else if(!p)p=z;
- y = rest(y); z = &rest(*z); }
- while(x) /* alternative loop */
- { d = digit(x)-borrow;
- borrow = ((d&IBASE)!=0);
- d = d&MAXDIGIT;
- *z = make(INT,d,0);
- if(d)p=NULL; else if(!p)p=z;
- x = rest(x); z = &rest(*z); }
- if(borrow) /* result is negative - take complement and add 1 */
- { p=NULL;
- d = (digit(r)^MAXDIGIT) + 1;
- borrow = ((d&IBASE)!=0); /* borrow now means `carry' (sorry) */
- digit(r) = SIGNBIT|d; /* set sign bit of result */
- z = &rest(r);
- while(*z)
- { d = (digit(*z)^MAXDIGIT)+borrow;
- borrow = ((d&IBASE)!=0);
- digit(*z) = d = d&MAXDIGIT;
- if(d)p=NULL; else if(!p)p=z;
- z = &rest(*z); }
- }
- if(p)*p=0; /* remove redundant (ie trailing) zeros */
- return(r);
+ word *p = NULL; /* pointer to trailing zeros, if any */
+ x = rest(x);
+ y = rest(y);
+ while (x && y) { /* this loop has been unwrapped once, see above */
+ d = digit(x) - digit(y) - borrow;
+ borrow = (d & IBASE) != 0;
+ d = d & MAXDIGIT;
+ *z = make(INT, d, 0);
+ if (d)
+ p = NULL;
+ else if (!p)
+ p = z;
+ x = rest(x);
+ y = rest(y);
+ z = &rest(*z);
+ }
+ while (y) { /* at most one of these two loops will be invoked */
+ d = -digit(y) - borrow;
+ borrow = ((d & IBASE) != 0);
+ d = d & MAXDIGIT;
+ *z = make(INT, d, 0);
+ if (d)
+ p = NULL;
+ else if (!p)
+ p = z;
+ y = rest(y);
+ z = &rest(*z);
+ }
+ while (x) { /* alternative loop */
+ d = digit(x) - borrow;
+ borrow = ((d & IBASE) != 0);
+ d = d & MAXDIGIT;
+ *z = make(INT, d, 0);
+ if (d)
+ p = NULL;
+ else if (!p)
+ p = z;
+ x = rest(x);
+ z = &rest(*z);
+ }
+ if (borrow) { /* result is negative - take complement and add 1 */
+ p = NULL;
+ d = (digit(r) ^ MAXDIGIT) + 1;
+ borrow = ((d & IBASE) != 0); /* borrow now means `carry' (sorry) */
+ digit(r) = SIGNBIT | d; /* set sign bit of result */
+ z = &rest(r);
+ while (*z) {
+ d = (digit(*z) ^ MAXDIGIT) + borrow;
+ borrow = ((d & IBASE) != 0);
+ digit(*z) = d = d & MAXDIGIT;
+ if (d)
+ p = NULL;
+ else if (!p)
+ p = z;
+ z = &rest(*z);
+ }
+ }
+ if (p)
+ *p = 0; /* remove redundant (ie trailing) zeros */
+ return (r);
}
-int bigcmp(x,y) /* returns +ve,0,-ve as x greater than, equal, less than y */
-word x,y;
-{ word d,r,s=neg(x);
- if(neg(y)!=s)return(s?-1:1);
- r=digit0(x)-digit0(y);
- for(;;)
- { x=rest(x); y=rest(y);
- if(!x)if(y)return(s?1:-1);
- else return(s?-r:r);
- if(!y)return(s?-1:1);
- d=digit(x)-digit(y);
- if(d)r=d; }
-}
-
-word bigtimes(x,y) /* naive multiply - quadratic */
-word x,y;
-{ if(len(x)<len(y))
- { word hold=x; x=y; y=hold; } /* important optimisation */
- word r=make(INT,0,0);
+/* returns +ve,0,-ve as x greater than, equal, less than y */
+int bigcmp(word x, word y)
+{
+ word d, r, s = neg(x);
+ if (neg(y) != s) {
+ return (s ? -1 : 1);
+ }
+ r = digit0(x) - digit0(y);
+ for (;;) {
+ x = rest(x);
+ y = rest(y);
+ if (!x) {
+ if (y) {
+ return (s ? 1 : -1);
+ } else {
+ return (s ? -r : r);
+ }
+ }
+ if (!y) {
+ return (s ? -1 : 1);
+ }
+ d = digit(x) - digit(y);
+ if (d) {
+ r = d;
+ }
+ }
+}
+
+/* naive multiply - quadratic */
+word bigtimes(word x, word y)
+{
+ if (len(x) < len(y)) {
+ word hold = x;
+ x = y;
+ y = hold;
+ } /* important optimisation */
+ word r = make(INT, 0, 0);
word d = digit0(y);
- word s=neg(y);
- word n=0;
- if(bigzero(x))return(r); /* short cut */
- for(;;)
- { if(d)r = bigplus(r,shift(n,stimes(x,d)));
- n++;
- y = rest(y);
- if(!y)
- return(s!=neg(x)?bignegate(r):r);
- d=digit(y); }
-}
-
-
-word shift(n,x) /* multiply big x by n'th power of IBASE */
-word n,x;
-{ while(n--)x=make(INT,0,x);
- return(x);
-} /* NB - we assume x non-zero, else unnormalised result */
-
-word stimes(x,n) /* multiply big x (>=0) by digit n (>0) */
-word x,n;
-{ unsigned d= n*digit0(x); /* ignore sign of x */
- word carry=d>>DIGITWIDTH;
- word r = make(INT,d&MAXDIGIT,0);
+ word s = neg(y);
+ word n = 0;
+ if (bigzero(x)) {
+ return (r); /* short cut */
+ }
+ for (;;) {
+ if (d) {
+ r = bigplus(r, shift(n, stimes(x, d)));
+ }
+ n++;
+ y = rest(y);
+ if (!y) {
+ return (s != neg(x) ? bignegate(r) : r);
+ }
+ d = digit(y);
+ }
+}
+
+/* multiply big x by n'th power of IBASE */
+/* NB - we assume x non-zero, else unnormalised result */
+word shift(word n, word x)
+{
+ while (n--) {
+ x = make(INT, 0, x);
+ }
+ return (x);
+}
+
+/* multiply big x (>=0) by digit n (>0) */
+word stimes(word x, word n)
+{
+ unsigned d = n * digit0(x); /* ignore sign of x */
+ word carry = d >> DIGITWIDTH;
+ word r = make(INT, d & MAXDIGIT, 0);
word *y = &rest(r);
- while(x=rest(x))
- d=n*digit(x)+carry,
- *y=make(INT,d&MAXDIGIT,0),
- y = &rest(*y),
- carry=d>>DIGITWIDTH;
- if(carry)*y=make(INT,carry,0);
- return(r);
+ while ((x = rest(x))) {
+ d = n * digit(x) + carry,
+ *y = make(INT, d & MAXDIGIT, 0), y = &rest(*y), carry = d >> DIGITWIDTH;
+ }
+ if (carry) {
+ *y = make(INT, carry, 0);
+ }
+ return (r);
}
-word b_rem; /* contains remainder from last call to longdiv or shortdiv */
+/* may assume y~=0 */
+word bigdiv(word x, word y)
+{
+ word s1, s2, q;
-word bigdiv(x,y) /* may assume y~=0 */
-word x,y;
-{ word s1,s2,q;
/* make x,y positive and remember signs */
- if(s1=neg(y))y=make(INT,digit0(y),rest(y));
- if(neg(x))
- x=make(INT,digit0(x),rest(x)),s2=!s1;
- else s2=s1;
+ if ((s1 = neg(y))) {
+ y = make(INT, digit0(y), rest(y));
+ }
+ if (neg(x)) {
+ x = make(INT, digit0(x), rest(x)), s2 = !s1;
+ } else {
+ s2 = s1;
+ }
+
/* effect: s1 set iff y negative, s2 set iff signs mixed */
- if(rest(y))q=longdiv(x,y);
- else q=shortdiv(x,digit(y));
- if(s2){ if(!bigzero(b_rem))
- { x=q;
- while((digit(x)+=1)==IBASE) /* add 1 to q in situ */
- { digit(x)=0;
- if(!rest(x)){ rest(x)=make(INT,1,0); break; }
- else x=rest(x);
- }
- }
- if(!bigzero(q))digit(q)=SIGNBIT|digit(q);
- }
- return(q);
-}
-
-word bigmod(x,y) /* may assume y~=0 */
-word x,y;
-{ word s1,s2;
+ if (rest(y)) {
+ q = longdiv(x, y);
+ } else {
+ q = shortdiv(x, digit(y));
+ }
+
+ if (s2) {
+ if (!bigzero(b_rem)) {
+ x = q;
+ while ((digit(x) += 1) == IBASE) { /* add 1 to q in situ */
+ digit(x) = 0;
+ if (!rest(x)) {
+ rest(x) = make(INT, 1, 0);
+ break;
+ } else {
+ x = rest(x);
+ }
+ }
+ }
+ if (!bigzero(q)) {
+ digit(q) = SIGNBIT | digit(q);
+ }
+ }
+ return (q);
+}
+
+/* may assume y~=0 */
+word bigmod(word x, word y)
+{
+ word s1, s2;
+
/* make x,y positive and remember signs */
- if(s1=neg(y))y=make(INT,digit0(y),rest(y));
- if(neg(x))
- x=make(INT,digit0(x),rest(x)),s2=!s1;
- else s2=s1;
+ if ((s1 = neg(y))) {
+ y = make(INT, digit0(y), rest(y));
+ }
+
+ if (neg(x)) {
+ x = make(INT, digit0(x), rest(x)), s2 = !s1;
+ } else {
+ s2 = s1;
+ }
+
/* effect: s1 set iff y negative, s2 set iff signs mixed */
- if(rest(y))longdiv(x,y);
- else shortdiv(x,digit(y));
- if(s2){ if(!bigzero(b_rem))
- b_rem = bigsub(y,b_rem);
- }
- return(s1?bignegate(b_rem):b_rem);
+ if (rest(y)) {
+ longdiv(x, y);
+ } else {
+ shortdiv(x, digit(y));
+ }
+
+ if (s2) {
+ if (!bigzero(b_rem)) {
+ b_rem = bigsub(y, b_rem);
+ }
+ }
+ return (s1 ? bignegate(b_rem) : b_rem);
}
/* NB - above have entier based handling of signed cases (as Miranda) in
@@ -264,141 +390,213 @@ word x,y;
magnitudes invariant under change of sign, remainder has sign of
dividend, quotient negative if signs of divi(sor/dend) mixed */
-word shortdiv(x,n) /* divide big x by single digit n returning big quotient
- and setting external b_rem as side effect */
- /* may assume - x>=0,n>0 */
-word x,n;
-{ word d=digit(x),s_rem,q=0;
- while(x=rest(x)) /* reverse rest(x) into q */
- q=make(INT,d,q),d=digit(x); /* leaving most sig. digit in d */
- { word tmp;
- x=q; s_rem=d%n; d=d/n;
- if(d||!q)q=make(INT,d,0); /* put back first digit (if not leading 0) */
- else q=0;
- while(x) /* in situ division of q by n AND destructive reversal */
- d=s_rem*IBASE+digit(x),digit(x)=d/n,s_rem=d%n,
- tmp=x,x=rest(x),rest(tmp)=q,q=tmp;
- }
- b_rem=make(INT,s_rem,0);
- return(q);
-}
-
-word longdiv(x,y) /* divide big x by big y returning quotient, leaving
- remainder in extern variable b_rem */
- /* may assume - x>=0,y>0 */
-word x,y;
-{ word n,q,ly,y1,scale;
- if(bigcmp(x,y)<0){ b_rem=x; return(make(INT,0,0)); }
- y1=msd(y);
- if((scale=IBASE/(y1+1))>1) /* rescale if necessary */
- x=stimes(x,scale),y=stimes(y,scale),y1=msd(y);
- n=q=0;ly=len(y);
- while(bigcmp(x,y=make(INT,0,y))>=0)n++;
- y=rest(y); /* want largest y not exceeding x */
+/* divide big x by single digit n returning big quotient
+ and setting external b_rem as side effect */
+
+/* may assume - x>=0,n>0 */
+word shortdiv(word x, word n)
+{
+ word d = digit(x), s_rem, q = 0;
+ while ((x = rest(x))) { /* reverse rest(x) into q */
+ q = make(INT, d, q), d = digit(x); /* leaving most sig. digit in d */
+ }
+
+ word tmp;
+ x = q;
+ s_rem = d % n;
+ d = d / n;
+
+ if (d || !q) {
+ q = make(INT, d, 0); /* put back first digit (if not leading 0) */
+ } else {
+ q = 0;
+ }
+
+ while (x) { /* in situ division of q by n AND destructive reversal */
+ d = s_rem * IBASE + digit(x), digit(x) = d / n, s_rem = d % n,
+ tmp = x, x = rest(x), rest(tmp) = q, q = tmp;
+ }
+
+ b_rem = make(INT, s_rem, 0);
+
+ return (q);
+}
+
+/* divide big x by big y returning quotient, leaving
+ remainder in extern variable b_rem */
+/* may assume - x>=0,y>0 */
+word longdiv(word x, word y)
+{
+ word n, q, ly, y1, scale;
+ if (bigcmp(x, y) < 0) {
+ b_rem = x;
+ return (make(INT, 0, 0));
+ }
+ y1 = msd(y);
+ if ((scale = IBASE / (y1 + 1)) > 1) { /* rescale if necessary */
+ x = stimes(x, scale), y = stimes(y, scale), y1 = msd(y);
+ }
+
+ n = q = 0;
+ ly = len(y);
+ while (bigcmp(x, y = make(INT, 0, y)) >= 0) {
+ n++;
+ }
+
+ y = rest(y); /* want largest y not exceeding x */
ly += n;
- for(;;)
- { word d,lx=len(x);
- if(lx<ly)d=0; else
- if(lx==ly)
- if(bigcmp(x,y)>=0)x=bigsub(x,y),d=1;
- else d=0;
- else{ d=ms2d(x)/y1;
- if(d>MAXDIGIT)d=MAXDIGIT;
- if((d -= 2)>0)x=bigsub(x,stimes(y,d));
- else d=0;
- if(bigcmp(x,y)>=0)
- { x=bigsub(x,y),d++;
- if(bigcmp(x,y)>=0)
- x=bigsub(x,y),d++; }
- }
- q = make(INT,d,q);
- if(n-- ==0)
- { b_rem = scale==1?x:shortdiv(x,scale); return(q); }
- ly-- ; y = rest(y); }
-} /* see Bird & Wadler p82 for explanation */
-
-word len(x) /* no of digits in big x */
-word x;
-{ word n=1;
- while(x=rest(x))n++;
- return(n);
-}
-
-word msd(x) /* most significant digit of big x */
-word x;
-{ while(rest(x))x=rest(x);
- return(digit(x)); /* sign? */
-}
-
-word ms2d(x) /* most significant 2 digits of big x (len>=2) */
-word x;
-{ word d=digit(x);
- x=rest(x);
- while(rest(x))d=digit(x),x=rest(x);
- return(digit(x)*IBASE+d);
-}
-
-word bigpow(x,y) /* assumes y poz */
-word x,y;
-{ word d,r=make(INT,1,0);
- while(rest(y)) /* this loop has been unwrapped once, see below */
- { word i=DIGITWIDTH;
- d=digit(y);
- while(i--)
- { if(d&1)r=bigtimes(r,x);
- x = bigtimes(x,x);
- d >>= 1; }
- y=rest(y);
- }
- d=digit(y);
- if(d&1)r=bigtimes(r,x);
- while(d>>=1)
- { x = bigtimes(x,x);
- if(d&1)r=bigtimes(r,x); }
- return(r);
+ for (;;) {
+ word d, lx = len(x);
+ if (lx < ly) {
+ d = 0;
+ } else if (lx == ly) {
+ if (bigcmp(x, y) >= 0) {
+ x = bigsub(x, y), d = 1;
+ } else {
+ d = 0;
+ }
+ } else {
+ d = ms2d(x) / y1;
+ if (d > MAXDIGIT) {
+ d = MAXDIGIT;
+ }
+ if ((d -= 2) > 0) {
+ x = bigsub(x, stimes(y, d));
+ } else {
+ d = 0;
+ }
+ if (bigcmp(x, y) >= 0) {
+ x = bigsub(x, y), d++;
+ if (bigcmp(x, y) >= 0) {
+ x = bigsub(x, y), d++;
+ }
+ }
+ }
+ q = make(INT, d, q);
+ if (n-- == 0) {
+ b_rem = scale == 1 ? x : shortdiv(x, scale);
+ return (q);
+ }
+ ly--;
+ y = rest(y);
+ }
+} /* see Bird & Wadler p82 for explanation */
+
+/* no of digits in big x */
+word len(word x)
+{
+ word n = 1;
+ while ((x = rest(x))) {
+ n++;
+ }
+ return (n);
+}
+
+/* most significant digit of big x */
+word msd(word x)
+{
+ while (rest(x)) {
+ x = rest(x);
+ }
+ return (digit(x)); /* sign? */
}
-double bigtodbl(x)
-word x;
-{ word s=neg(x);
- double b=1.0, r=(double)digit0(x);
+/* most significant 2 digits of big x (len>=2) */
+word ms2d(word x)
+{
+ word d = digit(x);
x = rest(x);
- while(x)b=b*IBASE,r=r+b*digit(x),x=rest(x);
- if(s)return(-r);
- return(r);
-} /* small end first */
+ while (rest(x)) {
+ d = digit(x), x = rest(x);
+ }
+ return (digit(x) * IBASE + d);
+}
+
+/* assumes y poz */
+word bigpow(word x, word y)
+{
+ word d, r = make(INT, 1, 0);
+ while (rest(y)) { /* this loop has been unwrapped once, see below */
+ word i = DIGITWIDTH;
+ d = digit(y);
+ while (i--) {
+ if (d & 1) {
+ r = bigtimes(r, x);
+ }
+ x = bigtimes(x, x);
+ d >>= 1;
+ }
+ y = rest(y);
+ }
+ d = digit(y);
+ if (d & 1) {
+ r = bigtimes(r, x);
+ }
+ while (d >>= 1) {
+ x = bigtimes(x, x);
+ if (d & 1) {
+ r = bigtimes(r, x);
+ }
+ }
+ return (r);
+}
+
+double bigtodbl(word x)
+{
+ word s = neg(x);
+ double b = 1.0, r = (double)digit0(x);
+ x = rest(x);
+ while (x) {
+ b = b * IBASE, r = r + b * digit(x), x = rest(x);
+ }
+ if (s) {
+ return (-r);
+ }
+ return (r);
+} /* small end first */
+
/* note: can return oo, -oo
but is used without surrounding sto_/set)dbl() only in compare() */
-/* not currently used
-long double bigtoldbl(x)
-word x;
-{ int s=neg(x);
+/* not currently used */
+long double bigtoldbl(word x)
+{
+ int s=neg(x);
long double b=1.0L, r=digit0(x);
x = rest(x);
- while(x)b=b*IBASE,r=r+b*digit(x),x=rest(x);
-/*printf("bigtoldbl returns %Le\n",s?-r:r); /* DEBUG
- if(s)return(-r);
+ while (x) {
+ b=b*IBASE,r=r+b*digit(x),x=rest(x);
+ }
+ if (s) {
+ return(-r);
+ }
return(r);
} /* not compatible with std=c90, lib fns eg sqrtl broken */
-word dbltobig(x) /* entier */
-double x;
-{ word s= (x<0);
- word r=make(INT,0,0);
+/* entier */
+word dbltobig(double x)
+{
+ word s = (x < 0);
+ word r = make(INT, 0, 0);
word *p = &r;
- double y= floor(x);
-/*if(fabs(y-x+1.0)<1e-9)y += 1.0; /* trick due to Peter Bartke, see note */
- for(y=fabs(y);;)
- { double n = fmod(y,(double)IBASE);
- digit(*p) = (word)n;
- y = (y-n)/(double)IBASE;
- if(y>0.0)rest(*p)=make(INT,0,0),p=&rest(*p);
- else break;
- }
- if(s)digit(r)=SIGNBIT|digit(r);
- return(r);
+ double y = floor(x);
+
+ for (y = fabs(y);;) {
+ double n = fmod(y, (double)IBASE);
+ digit(*p) = (word) n;
+ y = (y - n) / (double)IBASE;
+ if (y > 0.0) {
+ rest(*p) = make(INT, 0, 0), p = &rest(*p);
+ } else {
+ break;
+ }
+ }
+ if (s) {
+ digit(r) = SIGNBIT | digit(r);
+ }
+ return (r);
}
+
/* produces junk in low order digits if x exceeds range in which integer
can be held without error as a double -- NO, see next comment */
/* hugs, ghci, mira produce same integer for floor/entier hugenum, has 2^971
@@ -411,246 +609,329 @@ double x;
the exact integers. There are inherent deficiences in 64 bit fp,
no point in trying to mask this */
-double biglog(x) /* logarithm of big x */
-word x;
-{ word n=0;
- double r=digit(x);
- if(neg(x)||bigzero(x))errno=EDOM,math_error("log");
- while(x=rest(x))n++,r=digit(x)+r/IBASE;
- return(log(r)+n*logIBASE);
-}
-
-double biglog10(x) /* logarithm of big x */
-word x;
-{ word n=0;
- double r=digit(x);
- if(neg(x)||bigzero(x))errno=EDOM,math_error("log10");
- while(x=rest(x))n++,r=digit(x)+r/IBASE;
- return(log10(r)+n*log10IBASE);
-}
-
-word bigscan(p) /* read a big number (in decimal) */
- /* NB does NOT check for malformed number, assumes already done */
-char *p; /* p is a pointer to a null terminated string of digits */
-{ word s=0,r=make(INT,0,0);
- if(*p=='-')s=1,p++; /* optional leading `-' (for NUMVAL) */
- while(*p)
- { word d= *p-'0',f=10;
- p++;
- while(*p&&f<PTEN)d=10*d+*p-'0',f=10*f,p++;
- /* rest of loop does r=f*r+d; (in situ) */
- d= f*digit(r)+d;
- { word carry=d>>DIGITWIDTH;
- word *x = &rest(r);
- digit(r)=d&MAXDIGIT;
- while(*x)
- d=f*digit(*x)+carry,
- digit(*x)=d&MAXDIGIT,
- carry=d>>DIGITWIDTH,
- x = &rest(*x);
- if(carry)*x=make(INT,carry,0);
- }}
-/*if(*p=='e')
- { int s=bigscan(p+1);
- r = bigtimes(r,bigpow(make(INT,10,0),s); } */
- if(s&&!bigzero(r))digit(r)=digit(r)|SIGNBIT;
- return(r);
+/* logarithm of big x */
+double biglog(word x)
+{
+ word n = 0;
+ double r = digit(x);
+ if (neg(x) || bigzero(x)) {
+ errno = EDOM, math_error("log");
+ }
+ while ((x = rest(x))) {
+ n++, r = digit(x) + r / IBASE;
+ }
+ return (log(r) + n * logIBASE);
}
+
+/* logarithm of big x */
+double biglog10(word x)
+{
+ word n = 0;
+ double r = digit(x);
+ if (neg(x) || bigzero(x)) {
+ errno = EDOM, math_error("log10");
+ }
+ while ((x = rest(x))) {
+ n++, r = digit(x) + r / IBASE;
+ }
+ return (log10(r) + n * log10IBASE);
+}
+
+/* NB does NOT check for malformed number, assumes already done */
+/* p is a pointer to a null terminated string of digits */
+/* read a big number (in decimal) */
+word bigscan(char *p)
+{
+ word s = 0, r = make(INT, 0, 0);
+ if (*p == '-') {
+ s = 1, p++; /* optional leading `-' (for NUMVAL) */
+ }
+
+ while (*p) {
+ word d = *p - '0', f = 10;
+ p++;
+ while (*p && f < PTEN) {
+ d = 10 * d + *p - '0', f = 10 * f, p++;
+ }
+
+ /* rest of loop does r=f*r+d; (in situ) */
+ d = f * digit(r) + d;
+ word carry = d >> DIGITWIDTH;
+ word *x = &rest(r);
+ digit(r) = d & MAXDIGIT;
+ while (*x) {
+ d = f * digit(*x) + carry,
+ digit(*x) = d & MAXDIGIT, carry = d >> DIGITWIDTH, x = &rest(*x);
+ }
+ if (carry) {
+ *x = make(INT, carry, 0);
+ }
+ }
+ if (s && !bigzero(r)) {
+ digit(r) = digit(r) | SIGNBIT;
+ }
+ return (r);
+}
+
/* code to handle (unsigned) exponent commented out */
-word bigxscan(p,q) /* read unsigned hex number in '\0'-terminated string p to q */
- /* assumes redundant leading zeros removed */
-char *p, *q;
-{ word r; /* will hold result */
+/* assumes redundant leading zeros removed */
+/* read unsigned hex number in '\0'-terminated string p to q */
+word bigxscan(char *p, char *q)
+{
+ word r; /* will hold result */
word *x = &r;
- if(*p=='0'&&!p[1])return make(INT,0,0);
- while(q>p)
- { unsigned long long hold;
- q = q-p<15 ? p : q-15; /* read upto 15 hex digits from small end */
- sscanf(q,"%llx",&hold);
- *q = '\0';
- word count=4; /* 15 hex digits => 4 bignum digits */
- while(count-- && !(hold==0 && q==p))
- *x = make(INT,hold&MAXDIGIT,0),
- hold >>= DIGITWIDTH,
- x = &rest(*x);
- }
+ if (*p == '0' && !p[1]) {
+ return make(INT, 0, 0);
+ }
+ while (q > p) {
+ unsigned long long hold;
+ q = q - p < 15 ? p : q - 15; /* read upto 15 hex digits from small end */
+ sscanf(q, "%llx", &hold);
+ *q = '\0';
+ word count = 4; /* 15 hex digits => 4 bignum digits */
+ while (count-- && !(hold == 0 && q == p)) {
+ *x = make(INT, hold & MAXDIGIT, 0), hold >>= DIGITWIDTH, x = &rest(*x);
+ }
+ }
return r;
}
-word bigoscan(p,q) /* read unsigned octal number in '\0'-terminated string p to q */
- /* assumes redundant leading zeros removed */
-char *p, *q;
-{ word r; /* will hold result */
+/* read unsigned octal number in '\0'-terminated string p to q */
+word bigoscan(char *p, char *q)
+ /* assumes redundant leading zeros removed */
+{
+ word r; /* will hold result */
word *x = &r;
- while(q>p)
- { unsigned hold;
- q = q-p<5 ? p : q-5; /* read (upto) 5 octal digits from small end */
- sscanf(q,"%o",&hold);
- *q = '\0';
- *x = make(INT,hold,0),
- x = &rest(*x);
- }
+ while (q > p) {
+ unsigned hold;
+ q = q - p < 5 ? p : q - 5; /* read (upto) 5 octal digits from small end */
+ sscanf(q, "%o", &hold);
+ *q = '\0';
+ *x = make(INT, hold, 0), x = &rest(*x);
+ }
return r;
}
-word digitval(c)
-char c;
-{ return isdigit(c)?c-'0':
- isupper(c)?10+c-'A':
- 10+c-'a'; }
-
-word strtobig(z,base) /* numeral (as Miranda string) to big number */
- /* does NOT check for malformed numeral, assumes
- done and that z fully evaluated */
-word z; int base;
-{ word s=0,r=make(INT,0,0),PBASE=PTEN;
- if(base==16)PBASE=PSIXTEEN; else
- if(base==8)PBASE=PEIGHT;
- if(z!=NIL&&hd[z]=='-')s=1,z=tl[z]; /* optional leading `-' (for NUMVAL) */
- if(base!=10)z=tl[tl[z]]; /* remove "0x" or "0o" */
- while(z!=NIL)
- { word d=digitval(hd[z]),f=base;
- z=tl[z];
- while(z!=NIL&&f<PBASE)d=base*d+digitval(hd[z]),f=base*f,z=tl[z];
- /* rest of loop does r=f*r+d; (in situ) */
- d= f*digit(r)+d;
- { word carry=d>>DIGITWIDTH;
- word *x = &rest(r);
- digit(r)=d&MAXDIGIT;
- while(*x)
- d=f*digit(*x)+carry,
- digit(*x)=d&MAXDIGIT,
- carry=d>>DIGITWIDTH,
- x = &rest(*x);
- if(carry)*x=make(INT,carry,0);
- }}
- if(s&&!bigzero(r))digit(r)=digit(r)|SIGNBIT;
- return(r);
+word digitval(char c)
+{
+ return isdigit(c) ? c - '0' : isupper(c) ? 10 + c - 'A' : 10 + c - 'a';
}
-extern char *dicp;
+/* numeral (as Miranda string) to big number */
+/* does NOT check for malformed numeral, assumes
+ done and that z fully evaluated */
+word strtobig(word z, int base)
+{
+ word s = 0, r = make(INT, 0, 0), PBASE = PTEN;
+
+ if (base == 16) {
+ PBASE = PSIXTEEN;
+ } else if (base == 8) {
+ PBASE = PEIGHT;
+ }
+
+ if (z != NIL && hd[z] == '-') {
+ s = 1, z = tl[z]; /* optional leading `-' (for NUMVAL) */
+ }
+
+ if (base != 10) {
+ z = tl[tl[z]]; /* remove "0x" or "0o" */
+ }
+
+ while (z != NIL) {
+ word d = digitval(hd[z]), f = base;
+ z = tl[z];
+ while (z != NIL && f < PBASE) {
+ d = base * d + digitval(hd[z]), f = base * f, z = tl[z];
+ }
+
+ /* rest of loop does r=f*r+d; (in situ) */
+ d = f * digit(r) + d;
+ word carry = d >> DIGITWIDTH;
+ word *x = &rest(r);
+ digit(r) = d & MAXDIGIT;
+ while (*x) {
+ d = f * digit(*x) + carry,
+ digit(*x) = d & MAXDIGIT, carry = d >> DIGITWIDTH, x = &rest(*x);
+ }
+ if (carry) {
+ *x = make(INT, carry, 0);
+ }
+ }
+
+ if (s && !bigzero(r)) {
+ digit(r) = digit(r) | SIGNBIT;
+ }
+
+ return (r);
+}
+
+
+word bigtostr(word x) /* number to decimal string (as Miranda list) */
+{
+ word x1, sign, s = NIL;
-word bigtostr(x) /* number to decimal string (as Miranda list) */
-word x;
-{ word x1,sign,s=NIL;
#ifdef DEBUG
extern int debug;
- if(debug&04) /* print octally */
- { word d=digit0(x);
- sign=neg(x);
- for(;;)
- { word i=OCTW;
- while(i--||d)s=cons('0'+(d&07),s),d >>= 3;
- x=rest(x);
- if(x)s=cons(' ',s),d=digit(x);
- else return(sign?cons('-',s):s); }
+ if (debug & 04) { /* print octally */
+ word d = digit0(x);
+ sign = neg(x);
+ for (;;) {
+ word i = OCTW;
+ while (i-- || d) {
+ s = cons('0' + (d & 07), s), d >>= 3;
+ }
+ x = rest(x);
+ if (x) {
+ s = cons(' ', s), d = digit(x);
+ } else {
+ return (sign ? cons('-', s) : s);
+ }
}
+ }
#endif
- if(rest(x)==0)
- { extern char *dicp;
- sprintf(dicp,"%ld",getsmallint(x));
- return(str_conv(dicp)); }
- sign=neg(x);
- x1=make(INT,digit0(x),0); /* reverse x into x1 */
- while(x=rest(x))x1=make(INT,digit(x),x1);
- x=x1;
- for(;;)
- { /* in situ division of (reversed order) x by PTEN */
- word d=digit(x),rem=d%PTEN;
- d=d/PTEN; x1=rest(x);
- if(d)digit(x)=d;
- else x=x1; /* remove leading zero from result */
- while(x1)
- d=rem*IBASE+digit(x1),
- digit(x1)=d/PTEN,
- rem=d%PTEN,
- x1=rest(x1);
- /* end of in situ division (also uses x1 as temporary) */
- if(x)
- { word i=TENW;
- while(i--)s=cons('0'+rem%10,s),rem=rem/10; }
- else
- { while(rem)s=cons('0'+rem%10,s),rem=rem/10;
- return(sign?cons('-',s):s); }
- }
-}
-
-word bigtostrx(x) /* integer to hexadecimal string (as Miranda list) */
-word x;
-{ word r=NIL, s=neg(x);
- while(x)
- { word count=4; /* 60 bits => 20 octal digits => 4 bignum digits */
- unsigned long long factor=1;
- unsigned long long hold=0;
- while(count-- && x) /* calculate value of (upto) 4 bignum digits */
- hold=hold+factor*digit0(x),
- /* printf("::%llx\n",hold), /* DEBUG */
- factor<<=15,
- x=rest(x);
- sprintf(dicp,"%.15llx",hold); /* 15 hex digits = 60 bits */
- /* printf(":::%s\n",dicp); /* DEBUG */
- char *q=dicp+15;
- while(--q>=dicp)r = cons(*q,r);
- }
- while(digit(r)=='0'&&rest(r)!=NIL)r=rest(r); /* remove redundant leading 0's */
- r = cons('0',cons('x',r));
- if(s)r = cons('-',r);
- return(r);
-}
-word bigtostr8(x) /* integer to octal string (as Miranda list) */
-word x;
-{ word r=NIL, s=neg(x);
- while(x)
- { char *q = dicp+5;
- sprintf(dicp,"%.5lo",digit0(x));
- while(--q>=dicp)r = cons(*q,r);
- x = rest(x); }
- while(digit(r)=='0'&&rest(r)!=NIL)r=rest(r); /* remove redundant leading 0's */
- r = cons('0',cons('o',r));
- if(s)r = cons('-',r);
- return(r);
+ if (rest(x) == 0) {
+ extern char *dicp;
+ sprintf(dicp, "%ld", getsmallint(x));
+ return (str_conv(dicp));
+ }
+
+ sign = neg(x);
+ x1 = make(INT, digit0(x), 0); /* reverse x into x1 */
+ while ((x = rest(x))) {
+ x1 = make(INT, digit(x), x1);
+ }
+
+ x = x1;
+ for (;;) { /* in situ division of (reversed order) x by PTEN */
+ word d = digit(x), rem = d % PTEN;
+ d = d / PTEN;
+ x1 = rest(x);
+
+ if (d) {
+ digit(x) = d;
+ } else {
+ x = x1; /* remove leading zero from result */
+ }
+
+ while (x1) {
+ d = rem * IBASE + digit(x1),
+ digit(x1) = d / PTEN, rem = d % PTEN, x1 = rest(x1);
+ }
+
+ /* end of in situ division (also uses x1 as temporary) */
+ if (x) {
+ word i = TENW;
+ while (i--) {
+ s = cons('0' + rem % 10, s), rem = rem / 10;
+ }
+ } else {
+ while (rem) {
+ s = cons('0' + rem % 10, s), rem = rem / 10;
+ }
+ return (sign ? cons('-', s) : s);
+ }
+ }
}
-#ifdef DEBUG
-wff(x) /* check for well-formation of integer */
-word x;
-{ word y=x;
- if(tag[x]!=INT)printf("BAD TAG %d\n",tag[x]);
- if(neg(x)&&!digit0(x)&&!rest(x))printf("NEGATIVE ZERO!\n");
- if(digit0(x)&(~MAXDIGIT))printf("OVERSIZED DIGIT!\n");
- while(x=rest(x))
- if(tag[x]!=INT)printf("BAD INTERNAL TAG %d\n",tag[x]); else
- if(digit(x)&(~MAXDIGIT))
- printf("OVERSIZED DIGIT!\n"); else
- if(!digit(x)&&!rest(x))
- printf("TRAILING ZERO!\n");
- return(y);
+/* integer to hexadecimal string (as Miranda list) */
+word bigtostrx(word x)
+{
+ word r = NIL, s = neg(x);
+ while (x) {
+ word count = 4; /* 60 bits => 20 octal digits => 4 bignum digits */
+ unsigned long long factor = 1;
+ unsigned long long hold = 0;
+ while (count-- && x) { /* calculate value of (upto) 4 bignum digits */
+ hold = hold + factor * digit0(x),
+ factor <<= 15, x = rest(x);
+ }
+ sprintf(dicp, "%.15llx", hold); /* 15 hex digits = 60 bits */
+ char *q = dicp + 15;
+ while (--q >= dicp) {
+ r = cons(*q, r);
+ }
+ }
+ while (digit(r) == '0' && rest(r) != NIL) {
+ r = rest(r); /* remove redundant leading 0's */
+ }
+
+ r = cons('0', cons('x', r));
+ if (s) {
+ r = cons('-', r);
+ }
+
+ return (r);
}
-normalise(x) /* remove trailing zeros */
-word x;
-{ if(rest(x))rest(x)=norm1(rest(x));
- return(wff(x));
+/* integer to octal string (as Miranda list) */
+word bigtostr8(word x)
+{
+ word r = NIL, s = neg(x);
+ while (x) {
+ char *q = dicp + 5;
+ sprintf(dicp, "%.5lo", digit0(x));
+ while (--q >= dicp) {
+ r = cons(*q, r);
+ }
+ x = rest(x);
+ }
+ while (digit(r) == '0' && rest(r) != NIL) {
+ r = rest(r); /* remove redundant leading 0's */
+ }
+ r = cons('0', cons('o', r));
+ if (s) {
+ r = cons('-', r);
+ }
+ return (r);
}
-norm1(x)
-word x;
-{ if(rest(x))rest(x)=norm1(rest(x));
- return(!digit(x)&&!rest(x)?0:x);
+#ifdef DEBUG
+/* check for well-formation of integer */
+wff(word x)
+{
+ word y = x;
+ if (tag[x] != INT) {
+ printf("BAD TAG %d\n", tag[x]);
+ }
+ if (neg(x) && !digit0(x) && !rest(x)) {
+ printf("NEGATIVE ZERO!\n");
+ }
+ if (digit0(x) & (~MAXDIGIT)) {
+ printf("OVERSIZED DIGIT!\n");
+ }
+
+ while (x = rest(x)) {
+ if (tag[x] != INT) {
+ printf("BAD INTERNAL TAG %d\n", tag[x]);
+ } else if (digit(x) & (~MAXDIGIT)) {
+ printf("OVERSIZED DIGIT!\n");
+ } else if (!digit(x) && !rest(x)) {
+ printf("TRAILING ZERO!\n");
+ }
+ }
+ return (y);
}
-#endif
+/* remove trailing zeros */
+normalise(word x)
+{
+ if (rest(x)) {
+ rest(x) = norm1(rest(x));
+ }
+ return (wff(x));
+}
-/* stall(s)
-char *s;
-{ fprintf(stderr,"big integer %s not yet implemented\n",s);
- exit(0);
+norm1(word x)
+{
+ if (rest(x)) {
+ rest(x) = norm1(rest(x));
+ }
+ return (!digit(x) && !rest(x) ? 0 : x);
}
-#define destrev(x,y,z) while(x)z=x,x=rest(x),rest(z)=y,y=z;
-/* destructively reverse x into y using z as temp */
+#endif
/* END OF MIRANDA INTEGER PACKAGE */
-
generated by cgit v1.2.1 (git 2.18.0) at 2024-12-23 00:48:40 +0000