path: root/big.c

/* 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 <errno.h>

#define MAX_INT_BITS		60
#define LONG_LONG_OVERFLOW	(1ll<<MAX_INT_BITS)

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)
{
  big_one = make(INT, 1, 0);
}

int isnat(word x)
{
  return (tag[x] == INT && big_is_positive(x));
}

/* store C long long as mira bigint */
word sto_int(long long i)
{
  word sign = 0;

  if (i < 0) {
    sign = SIGNBIT;
    i = -i;
  }

  word x = make(INT, (sign | i) & MAXDIGIT, 0);
  if (i >>= DIGITWIDTH) {
    word *p = &rest(x);
    do {
      *p = make(INT, i & MAXDIGIT, 0);
      p = &rest(*p);
    } while (i >>= DIGITWIDTH);
  }

  return x;
}

/* mira bigint to C long long */
long long get_int(word x)
{
  long long n = digit0(x);
  word sign = big_is_negative(x);
  if (!(x = rest(x))) {
    return (sign ? -n : n);
  }

  word w = DIGITWIDTH;
  while (x && w < MAX_INT_BITS) {
    n += (long long)digit(x) << w;
    w += DIGITWIDTH;
    x = rest(x);
  }

  if (x) {
    n = LONG_LONG_OVERFLOW;
  }

  return (sign ? -n : n);
}

word bignegate(word x)
{
  if (bigzero(x)) {
    return x;
  }

  if (hd[x] & SIGNBIT) {
    return make(INT, hd[x] & MAXDIGIT, tl[x]);
  }
  return make(INT, SIGNBIT | hd[x], tl[x]);
}

word bigplus(word x, word y)
{
  if (big_is_positive(x)) {
    if (big_is_positive(y)) {
      return (big_plus(x, y, 0));
    }
    return (big_sub(x, y));
  }

  if (big_is_positive(y)) {
    return (big_sub(y, x));
  }

  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 r = make(INT, signbit | d & MAXDIGIT, 0);	/* result */
  word *z = &rest(r);		/* pointer to rest of result */

  word carry = ((d & IBASE) != 0);
  x = rest(x);
  y = rest(y);

  while (x && y) {
    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 (big_is_positive(x)) {
    if (big_is_positive(y)) {
      return (big_sub(x, y));		/* both positive */
    }
    return (big_plus(x, y, 0));		/* positive x, negative y */
  }

  if (big_is_positive(y)) {
    return (big_plus(x, y, SIGNBIT));	/* negative x, positive y */
  }

  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);

  /* this loop has been unwrapped once, see above */
  while (x && y) {
    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);
  }

  /* at most one of these two loops will be invoked */
  while (y) {
    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);
  }

  /* alternative loop */
  while (x) {
    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);
  }

  /* result is negative  - take complement and add 1 */
  if (borrow) {
    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);
    }
  }

  /* remove redundant (ie trailing) zeros */
  if (p) {
    *p = 0;
  }

  return r;
}

/* returns +ve,0,-ve as x greater than, equal, less than y */
int bigcmp(word x, word y)
{
  word s = big_is_negative(x);
  if (big_is_negative(y) != s) {
    return (s ? -1 : 1);
  }

  word r = digit0(x) - digit0(y);
  for (;;) {
    x = rest(x);
    y = rest(y);

    if (!x) {
      if (y) {
	return (s ? 1 : -1);
      }
      return (s ? -r : r);
    }

    if (!y) {
      return (s ? -1 : 1);
    }

    word d = digit(x) - digit(y);
    if (d) {
      r = d;
    }
  }
}

/* naive multiply - quadratic */
word bigtimes(word x, word y)
{
  /* important optimisation */
  if (len(x) < len(y)) {
    word hold = x;
    x = y;
    y = hold;
  }

  word r = make(INT, 0, 0);

  /* short cut */
  if (bigzero(x)) {
    return (r);
  }

  word s = big_is_negative(y);
  for (word n = 0, d = digit0(y); y; n++, d = digit(y)) {
    if (d) {
      r = bigplus(r, shift(n, stimes(x, d)));
    }

    y = rest(y);
  }
  return (s != big_is_negative(x) ? bignegate(r) : r);
}

/* multiply big x by n'th power of IBASE */
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;
}

/* TODO: make bigdiv()/bigmod() behave like div()/ldiv()/lldiv() */
word b_rem;			/* contains remainder from last call to longdiv or shortdiv */

/* may assume y~=0 */
word bigdiv(word x, word y)
{
  word s1 = big_is_negative(y);
  word s2 = s1;

  /* make x,y positive and remember signs */
  if (s1) {
    y = make(INT, digit0(y), rest(y));
  }

  /* effect: s1 set iff y negative, s2 set iff signs mixed */
  if (big_is_negative(x)) {
    x = make(INT, digit0(x), rest(x));
    s2 = !s1;
  }

  word q = rest(y) ? longdiv(x, y) : 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;
	}
        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 = big_is_negative(y);
  word s2 = s1;

  /* make x,y positive and remember signs */
  if (s1) {
    y = make(INT, digit0(y), rest(y));
  }

  if (big_is_negative(x)) {
    x = make(INT, digit0(x), rest(x));
    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);
  word q = 0;

  /* reverse rest(x) into q */
  while (x = rest(x)) {
    q = make(INT, d, q);
    d = digit(x);	/* leaving most sig. digit in d */
  }

  x = q;
  word s_rem = d % n;
  d = d / n;

  /* put back first digit (if not leading 0) */
  if (d || !q) {
    q = make(INT, d, 0);
  } else {
    q = 0;
  }

  /* in situ division of q by n AND destructive reversal */
  while (x) {
    d = s_rem * IBASE + digit(x);
    digit(x) = d / n;
    s_rem = d % n;
    word 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)
{
  if (bigcmp(x, y) < 0) {
    b_rem = x;
    return (make(INT, 0, 0));
  }

  word y1 = msd(y);

  /* rescale if necessary */
  word scale = IBASE / (y1 + 1);
  if (scale > 1) {
    x = stimes(x, scale);
    y = stimes(y, scale);
    y1 = msd(y);
  }

  word n = 0;
  word q = 0;
  word 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 big_is_positive */
word bigpow(word x, word y)
{
  word d;
  word r = make(INT, 1, 0);

  /* this loop has been unwrapped once, see below */
  while (rest(y)) {	
    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 = big_is_negative(x);
  double b = 1.0;
  double r = (double)digit0(x);
  x = rest(x);
  while (x) {
    b = b * IBASE;
    r = r + b * digit(x);
    x = rest(x);
  }

  return s ? -r : r;
}				/* small end first */

/* note: can return oo, -oo
   but is used without surrounding sto_/set)dbl() only in compare() */

/* 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 */

static double big__log(word x, double logbase, char *name, double (*fn)(double))
{
  word n = 0;
  double r = digit(x);

  if (big_is_negative(x) || bigzero(x)) {
    errno = EDOM;
    math_error(name);
  }

  while ((x = rest(x)) != 0) {
    n++;
    r = digit(x) + r / IBASE;
  }

  return (fn(r) + n * logbase);
}

/* logarithm of big x */
double biglog(word x)
{
  static double logIBASE = 0.0;
  if (logIBASE == 0.0) {
    logIBASE = log((double)IBASE);
  }

  return big__log(x, logIBASE, "log", log);
}

/* logarithm of big x */
double biglog10(word x)
{
  static double log10IBASE = 0.0;
  if (log10IBASE == 0.0) {
    log10IBASE = log10((double)IBASE);
  }

  return big__log(x, log10IBASE, "log10", log10);
}

/* read a big number (in decimal) */
/* NB does NOT check for malformed number, assumes already done */
/* p is a pointer to a null terminated string of digits */
word bigscan(char *p)
{
  word s = 0;
  word r = make(INT, 0, 0);

  /* optional leading `-' (for NUMVAL) */
  if (*p == '-') {
    s = 1;
    p++;
  }

  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 */

/* read unsigned hex number in '\0'-terminated string p to q */
/* assumes redundant leading zeros removed */
word bigxscan(char *p, char *q)
{
  if (*p == '0' && !p[1]) {
    return make(INT, 0, 0);
  }

  word r;			/* will hold result */
  word *x = &r;

  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 */
/* assumes redundant leading zeros removed */
word bigoscan(char *p, char *q)
{
  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;
}

/* TODO: this assumes ASCII/UTF-8 */
word digitval(char c)
{
  return isdigit(c) ? c - '0' : 10 + tolower(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;
  word r = make(INT, 0, 0);
 
  word PBASE = PTEN;
  if (base == 16) {
    PBASE = PSIXTEEN;
  } else if (base == 8) {
    PBASE = PEIGHT;
  }

  /* optional leading `-' (for NUMVAL) */
  if (z != NIL && hd[z] == '-') {
    s = 1;
    z = tl[z];
  }

  /* remove "0x" or "0o" */
  if (base != 10) {
    z = tl[tl[z]];
  }

  while (z != NIL) {
    word d = digitval(hd[z]);
    word 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;
}


/* number to decimal string (as Miranda list) */
word bigtostr(word x)
{
  word x1, sign, s = NIL;

#ifdef DEBUG
  extern int debug;
  if (debug & 04) {		/* print octally */
    word d = digit0(x);
    sign = big_is_negative(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 = big_is_negative(x);

  x1 = make(INT, digit0(x), 0);	/* reverse x into x1 */
  while (x = rest(x)) {
    x1 = make(INT, digit(x), x1);
  }
  x = x1;

  /* in situ division of (reversed order) x by PTEN */
  for (;;) {
    word d = digit(x);
    word 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)
{
  extern char *dicp;
  word r = NIL;
  word s = big_is_negative(x);

  while (x) {
    word count = 4;		/* 60 bits => 20 octal digits => 4 bignum digits */
    unsigned long long factor = 1;
    unsigned long long hold = 0;

    /* calculate value of (upto) 4 bignum digits */
    while (count-- && x) {
      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));

  return s ? cons('-', r) : r;
}

word bigtostr8(x)		/* integer to octal string (as Miranda list) */
word x;
{
  extern char *dicp;
  word r = NIL, s = big_is_negative(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);
}