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

static double logIBASE;
static double 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 && 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 = big_is_negative(y)) {
    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 */

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 big_is_positive */
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);
}

double bigtodbl(x)
word x;
{
  word s = big_is_negative(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=big_is_negative(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);
  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);
  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);
}

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

double biglog(x)		/* logarithm of big x */
word x;
{
  word n = 0;
  double r = digit(x);
  if (big_is_negative(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 (big_is_negative(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);
}

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

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 */
  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(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);
}

extern char *dicp;

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