/* MIRANDA INTEGER PACKAGE */ /* package for unbounded precision integers */ /************************************************************************** * Copyright (C) Research Software Limited 1985-90. All rights reserved. * * The Miranda system is distributed as free software under the terms in * * the file "COPYING" which is included in the distribution. * * * * Revised to C11 standard and made 64bit compatible, January 2020 * *------------------------------------------------------------------------*/ #include "data.h" #include "lex.h" #include "big.h" #include /* 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 len(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(void) { logIBASE = log((double)IBASE); log10IBASE = log10((double)IBASE); big_one = make(INT, 1, 0); } int isnat(word x) { return (tag[x] == INT && poz(x)); } /* 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(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 */ } } /* 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); } /* 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); } } /* 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); } /* may assume y~=0 */ word bigdiv(word x, word 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; } /* 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); } /* 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; } /* 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); } /* NB - above have entier based handling of signed cases (as Miranda) in which remainder has sign of divisor. To get this:- if signs of divi(sor/dend) mixed negate quotient and if remainder non-zero take complement and add one to magnitude of quotient */ /* for alternative, truncate based handling of signed cases (usual in C):- magnitudes invariant under change of sign, remainder has sign of dividend, quotient negative if signs of divi(sor/dend) mixed */ /* 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 */ /* 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? */ } /* most significant 2 digits of big x (len>=2) */ word ms2d(word x) { word d = digit(x); x = rest(x); 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(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); } if (s) { return(-r); } return(r); } /* not compatible with std=c90, lib fns eg sqrtl broken */ /* entier */ word dbltobig(double x) { word s = (x < 0); word r = make(INT, 0, 0); word *p = &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 as factor so the low order bits are NOT JUNK -- 9.1.12 */ /* note on suppressed fix: choice of 1e9 arbitrary, chosen to prevent eg entier(100*0.29) = 28 but has undesirable effects, causing eg entier 1.9999999999 = 2 underlying problem is that computable floor on true Reals is _|_ at the exact integers. There are inherent deficiences in 64 bit fp, no point in trying to mask this */ /* 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 */ /* 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); } } return r; } /* 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); } return r; } word digitval(char c) { return isdigit(c) ? c - '0' : isupper(c) ? 10 + c - 'A' : 10 + c - 'a'; } /* 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; #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); } } } #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); } } } /* 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); } /* 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); } #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); } /* remove trailing zeros */ normalise(word x) { if (rest(x)) { rest(x) = norm1(rest(x)); } return (wff(x)); } norm1(word x) { if (rest(x)) { rest(x) = norm1(rest(x)); } return (!digit(x) && !rest(x) ? 0 : x); } #endif /* END OF MIRANDA INTEGER PACKAGE */