/* Bignum arithmetic
 */

#include <math.h>
#include <ctype.h>

#include "scheme.h"

Object Make_Uninitialized_Bignum (size) {
    register char *p;
    Object big;
    
    p = Get_Bytes ((sizeof (struct S_Bignum) - sizeof (gran_t)) +
		   (size * sizeof (gran_t)));
    SET(big, T_Bignum, (struct S_Bignum *)p);
    BIGNUM(big)->minusp = False;
    BIGNUM(big)->size = size;
    BIGNUM(big)->usize = 0;
    return big;
}

Object Copy_Bignum (x) Object x; {
    Object big;
    register size;
    GC_Node;

    GC_Link (x);
    big = Make_Uninitialized_Bignum (size = BIGNUM(x)->usize);
    BIGNUM(big)->minusp = BIGNUM(x)->minusp;
    BIGNUM(big)->usize = size;
    bcopy ((char *)BIGNUM(x)->data, (char *)BIGNUM(big)->data,
	  size * sizeof (gran_t));
    GC_Unlink;
    return big;
}

Object Copy_S_Bignum (s) struct S_Bignum *s; {
    Object big;
    register size;

    big = Make_Uninitialized_Bignum (size = s->usize);
    BIGNUM(big)->minusp = s->minusp;
    BIGNUM(big)->usize = size;
    bcopy ((char *)s->data, (char *)BIGNUM(big)->data,
	  size * sizeof (gran_t));
    return big;
}

Object Make_Bignum (buf, neg, base) char *buf; {
    Object big;
    register char *p;
    register c;
    register size = (strlen (buf) + 4) / 4;

    big = Make_Uninitialized_Bignum (size);
    BIGNUM(big)->minusp = neg ? True : False;
    p = buf;
    while (c = *p++) {
	Bignum_Mult_In_Place (BIGNUM(big), base);
	if (base == 16) {
	    if (isupper (c))
		c = tolower (c);
	    if (c >= 'a')
		c = '9' + c - 'a' + 1;
	}
	Bignum_Add_In_Place (BIGNUM(big), c - '0');
    }
    Bignum_Normalize_In_Place (BIGNUM(big)); /* to avoid -0 */
    return big;
}

Object Reduce_Bignum (x) Object x; {
    register i;
    register struct S_Bignum *p = BIGNUM(x);

    if (p->usize > 2 || (p->usize == 2 && p->data[1] >= 32768))
	return x;
    i = Bignum_To_Integer (x);
    if (!FIXNUM_FITS(i))
	return x;
    return Make_Fixnum (i);
}

Bignum_Mult_In_Place (x, n) register struct S_Bignum *x; {
    register i = x->usize;
    register gran_t *p = x->data;
    register j;
    register unsigned k = 0;
    
    for (j = 0; j < i; ++j) {
	k += n * *p;
        *p++ = k;
	k >>= 16;
    }
    if (k) {
	if (i >= x->size)
	    Panic ("Bignum_Mult_In_Place");
	*p++ = k;
	x->usize++;
    }
}

Bignum_Add_In_Place (x, n) register struct S_Bignum *x; {
    register i = x->usize;
    register gran_t *p = x->data;
    register j = 0;
    register unsigned k = n;

    if (i == 0) goto extend;
    k += *p;
    *p++ = k;
    while (k >>= 16) {
	if (++j >= i) {
 extend:
	    if (i >= x->size)
		Panic ("Bignum_Add_In_Place");
	    *p++ = k;
	    x->usize++;
	    return;
	}
	k += *p;
	*p++ = k;
    }
}

Bignum_Div_In_Place (x, n) register struct S_Bignum *x; {
    register i = x->usize;
    register gran_t *p = x->data + i;
    register unsigned k = 0;
    for ( ; i; --i) {
	k <<= 16;
	k += *--p;
	*p = k / n;
	k %= n;
    }
    Bignum_Normalize_In_Place (x);
    return k;
}

Bignum_Normalize_In_Place (x) register struct S_Bignum *x; {
    register i = x->usize;
    register gran_t *p = x->data + i;
    while (i && !*--p)
	--i;
    x->usize = i;
    if (!i)
	x->minusp = False;
}

Print_Bignum (port, x) Object port, x; {
    register char *buf, *p;
    register size;
    register struct S_Bignum *big;
    
    if (Bignum_Zero (x)) {
	Printf (port, "0");
	return;
    }
    
    size = BIGNUM(x)->usize * 5 + 3;
    buf = alloca (size + 1);
    p = buf + size;
    *p = 0;

    size = (sizeof (struct S_Bignum) - sizeof (gran_t))
	+ BIGNUM(x)->usize * sizeof (gran_t);
    big = (struct S_Bignum *)alloca (size);
    bcopy ((char *)POINTER(x), (char *)big, size);
    big->size = BIGNUM(x)->usize;

    while (big->usize) {
	register unsigned bigdig = Bignum_Div_In_Place (big, 10000);
	*--p = '0' + bigdig % 10;
	bigdig /= 10;
	*--p = '0' + bigdig % 10;
	bigdig /= 10;
	*--p = '0' + bigdig % 10;
	bigdig /= 10;
	*--p = '0' + bigdig;
    }
    while (*p == '0')
	++p;
    if (Truep (BIGNUM(x)->minusp))
	Printf (port, "-");
    Format (port, p, strlen (p), 0, (Object *)0);
}

Bignum_To_Integer (x) Object x; {
    unsigned n = 0;
    int s = BIGNUM(x)->usize;

    if (s) {
	n = BIGNUM(x)->data[0];
	if (s > 1) {
	    n |= BIGNUM(x)->data[1] << 16;
	    if (s > 2)
err:		
		Primitive_Error ("integer out of range: ~s", x);
	}
    }
    if (Truep (BIGNUM(x)->minusp)) {
	if (n > (~(unsigned)0 >> 1) + 1)
	    goto err;
	return -n;
    } else {
	if (n > ~(unsigned)0 >> 1)
	    goto err;
	return n;
    }
}

Object Integer_To_Bignum (i) {
    Object big = Make_Uninitialized_Bignum (2);
    unsigned n = i;

    if (i < 0) {
	BIGNUM(big)->minusp = True;
	n = -i;
    }
    BIGNUM(big)->data[0] = n;
    BIGNUM(big)->data[1] = n >> 16;
    BIGNUM(big)->usize = 2;
    Bignum_Normalize_In_Place (BIGNUM(big));
    return big;
}

Object Unsigned_To_Bignum (i) unsigned i; {
    Object big = Make_Uninitialized_Bignum (2);

    BIGNUM(big)->data[0] = i;
    BIGNUM(big)->data[1] = i >> 16;
    BIGNUM(big)->usize = 2;
    Bignum_Normalize_In_Place (BIGNUM(big));
    return big;
}

Object Double_To_Bignum (d) double d; {         /* Truncates the double */
    Object big;
    int expo, size;
    double mantissa = frexp (d, &expo);
    register gran_t *p;
    
    if (expo <= 0 || mantissa == 0.0)
	return Make_Uninitialized_Bignum (0);
    size = (expo + (16-1)) / 16;
    big = Make_Uninitialized_Bignum (size);
    BIGNUM(big)->usize = size;
    if (mantissa < 0.0) {
	BIGNUM(big)->minusp = True;
	mantissa = -mantissa;
    }
    p = BIGNUM(big)->data;
    bzero ((char *)p, size * sizeof (gran_t));
    p += size;
    if (expo &= (16-1))
	mantissa = ldexp (mantissa, expo - 16);
    while (mantissa != 0.0) {
	if (--size < 0)
	    break;		/* inexact */
	mantissa *= 65536.0;
	*--p = (int)mantissa;
	mantissa -= *p;
    }
    Bignum_Normalize_In_Place (BIGNUM(big)); /* Probably not needed */
    return Reduce_Bignum (big);
}

double Bignum_To_Double (x) Object x; {   /* error if it ain't fit */
    double rx = 0.0;
    register i = BIGNUM(x)->usize;
    register gran_t *p = BIGNUM(x)->data + i;

    for (i = BIGNUM(x)->usize; --i >= 0; ) {
	if (rx >= HUGE / 65536.0)
	    Primitive_Error ("cannot coerce to real: ~s", x);
	rx *= 65536.0;
	rx += *--p;
    }
    if (Truep (BIGNUM(x)->minusp))
	rx = -rx;
    return rx;
}

Bignum_Zero (x) Object x; {
    return BIGNUM(x)->usize == 0;
}

Bignum_Negative (x) Object x; {
    return Truep (BIGNUM(x)->minusp);
}

Bignum_Positive (x) Object x; {
    return !Truep (BIGNUM(x)->minusp) && BIGNUM(x)->usize != 0;
}

Bignum_Even (x) Object x; {
    return BIGNUM(x)->usize == 0 || (BIGNUM(x)->data[0] & 1) == 0;
}

Object Bignum_Abs (x) Object x; {
    Object big;

    big = Copy_Bignum (x);
    BIGNUM(big)->minusp = False;
    return big;
}

Bignum_Mantissa_Cmp (x, y) register struct S_Bignum *x, *y; {
    register i = x->usize;
    if (i < y->usize)
	return -1;
    else if (i > y->usize)
	return 1;
    else {
	register gran_t *xbuf = x->data + i;
	register gran_t *ybuf = y->data + i;
	for ( ; i; --i) {
	    register n;
	    if (n = (int)*--xbuf - (int)*--ybuf)
		return n;
	}
	return 0;
    }
}

Bignum_Cmp (x, y) register struct S_Bignum *x, *y; {
    register xm = Truep (x->minusp);
    register ym = Truep (y->minusp);
    if (xm) {
	if (ym)
	    return -Bignum_Mantissa_Cmp (x, y);
	else return -1;
    } else {
	if (ym)
	    return 1;
	else return Bignum_Mantissa_Cmp (x, y);
    }
}

Bignum_Equal (x, y) Object x, y; {
    return Bignum_Cmp (BIGNUM(x), BIGNUM(y)) == 0;
}

Bignum_Less (x, y) Object x, y; {
    return Bignum_Cmp (BIGNUM(x), BIGNUM(y)) < 0;
}

Bignum_Greater (x, y) Object x, y; {
    return Bignum_Cmp (BIGNUM(x), BIGNUM(y)) > 0;
}

Bignum_Eq_Less (x, y) Object x, y; {
    return Bignum_Cmp (BIGNUM(x), BIGNUM(y)) <= 0;
}

Bignum_Eq_Greater (x, y) Object x, y; {
    return Bignum_Cmp (BIGNUM(x), BIGNUM(y)) >= 0;
}

Object General_Bignum_Plus_Minus (x, y, neg) Object x, y; {
    Object big;
    int size, xsize, ysize, xminusp, yminusp;
    GC_Node2;

    GC_Link2 (x,y);
    xsize = BIGNUM(x)->usize;
    ysize = BIGNUM(y)->usize;
    xminusp = Truep (BIGNUM(x)->minusp);
    yminusp = Truep (BIGNUM(y)->minusp);
    if (neg)
	yminusp = !yminusp;
    size = xsize > ysize ? xsize : ysize;
    if (xminusp == yminusp)
	size++;
    big = Make_Uninitialized_Bignum (size);
    BIGNUM(big)->usize = size;
    GC_Unlink;

    if (xminusp == yminusp) {
	/* Add x and y */
	register unsigned k = 0;
	register i;
	register gran_t *xbuf = BIGNUM(x)->data;
	register gran_t *ybuf = BIGNUM(y)->data;
	register gran_t *zbuf = BIGNUM(big)->data;
	for (i = 0; i < size; ++i) {
	    if (i < xsize)
		k += *xbuf++;
	    if (i < ysize)
		k += *ybuf++;
	    *zbuf++ = k;
	    k >>= 16;
	}
    } else {
	if (Bignum_Mantissa_Cmp (BIGNUM(x), BIGNUM(y)) < 0) {
	    Object temp = x;
	    x = y; y = temp;
	    xsize = ysize;
	    ysize = BIGNUM(y)->usize;
	    xminusp = yminusp;
	}
	/* Subtract y from x */
    {
	register unsigned k = 1;
	register i;
	register gran_t *xbuf = BIGNUM(x)->data;
	register gran_t *ybuf = BIGNUM(y)->data;
	register gran_t *zbuf = BIGNUM(big)->data;
	for (i = 0; i < size; ++i) {
	    if (i < xsize)
		k += *xbuf++;
	    else Panic ("General_Bignum_Plus_Minus");
	    if (i < ysize)
		k += ~*ybuf++ & 0xFFFF;
	    else k += 0xFFFF;
	    *zbuf++ = k;
	    k >>= 16;
	}
    }
    }
    BIGNUM(big)->minusp = xminusp ? True : False;
    Bignum_Normalize_In_Place (BIGNUM(big));
    return Reduce_Bignum (big);
}

Object Bignum_Plus (x, y) Object x, y; {   /* bignum + bignum */
    return General_Bignum_Plus_Minus (x, y, 0);
}

Object Bignum_Minus (x, y) Object x, y; {   /* bignum - bignum */
    return General_Bignum_Plus_Minus (x, y, 1);
}

Object Bignum_Fixnum_Multiply (x, y) Object x, y; {   /* bignum * fixnum */
    Object big;
    register size, xsize, i;
    register gran_t *xbuf, *zbuf;
    int yn = FIXNUM(y);
    register unsigned yl, yh;
    GC_Node;
    
    GC_Link (x);
    xsize = BIGNUM(x)->usize;
    size = xsize + 2;
    big = Make_Uninitialized_Bignum (size);
    BIGNUM(big)->usize = size;
    if (Truep (BIGNUM(x)->minusp) != (yn < 0))
	BIGNUM(big)->minusp = True;
    bzero ((char *)BIGNUM(big)->data, size * sizeof (gran_t));
    xbuf = BIGNUM(x)->data;
    if (yn < 0)
	yn = -yn;
    yl = yn & 0xFFFF;
    yh = yn >> 16;
    zbuf = BIGNUM(big)->data;
    for (i = 0; i < xsize; ++i) {
	register unsigned xf = xbuf[i];
	register unsigned k = 0;
	register gran_t *r = zbuf + i;
	k += xf * yl + *r;
	*r++ = k;
	k >>= 16;
	k += xf * yh + *r;
	*r++ = k;
	k >>= 16;
	*r = k;
    }
    GC_Unlink;
    Bignum_Normalize_In_Place (BIGNUM(big));
    return Reduce_Bignum (big);
}

Object Bignum_Multiply (x, y) Object x, y; {   /* bignum * bignum */
    Object big;
    register size, xsize, ysize, i, j;
    register gran_t *xbuf, *ybuf, *zbuf;
    GC_Node2;
    
    GC_Link2 (x, y);
    xsize = BIGNUM(x)->usize;
    ysize = BIGNUM(y)->usize;
    size = xsize + ysize;
    big = Make_Uninitialized_Bignum (size);
    BIGNUM(big)->usize = size;
    if (!EQ(BIGNUM(x)->minusp, BIGNUM(y)->minusp))
	BIGNUM(big)->minusp = True;
    bzero ((char *)BIGNUM(big)->data, size * sizeof (gran_t));
    xbuf = BIGNUM(x)->data;
    ybuf = BIGNUM(y)->data;
    zbuf = BIGNUM(big)->data;
    for (i = 0; i < xsize; ++i) {
	register unsigned xf = xbuf[i];
	register unsigned k = 0;
	register gran_t *p = ybuf;
	register gran_t *r = zbuf + i;
	for (j = 0; j < ysize; ++j) {
	    k += xf * *p++ + *r;
	    *r++ = k;
	    k >>= 16;
	}
	*r = k;
    }
    GC_Unlink;
    Bignum_Normalize_In_Place (BIGNUM(big));
    return Reduce_Bignum (big);
}

/* Returns cons cell (quotient . remainder); cdr is a fixnum
 */
Object Bignum_Fixnum_Divide (x, y) Object x, y; {   /* bignum / fixnum */
    Object big;
    register xsize, i;
    register gran_t *xbuf, *zbuf;
    int yn = FIXNUM(y);
    int xminusp, yminusp = 0;
    register unsigned rem;
    GC_Node;

    GC_Link (x);
    if (yn < 0) {
	yn = -yn;
	yminusp = 1;
    }
    if (yn > 0xFFFF) {
	big = Integer_To_Bignum (FIXNUM(y));
	GC_Unlink;
	return Bignum_Divide (x, big);
    }
    xsize = BIGNUM(x)->usize;
    big = Make_Uninitialized_Bignum (xsize);
    BIGNUM(big)->usize = xsize;
    xminusp = Truep (BIGNUM(x)->minusp);
    if (xminusp != yminusp)
	BIGNUM(big)->minusp = True;
    xbuf = BIGNUM(x)->data;
    zbuf = BIGNUM(big)->data;
    rem = 0;
    for (i = xsize; --i >= 0; ) {
	rem <<= 16;
	rem += xbuf[i];
	zbuf[i] = rem / yn;
	rem %= yn;
    }
    GC_Unlink;
    Bignum_Normalize_In_Place (BIGNUM(big));
    if (xminusp)
	rem = -(int)rem;
    return Cons (Reduce_Bignum (big), Make_Fixnum ((int)rem));
}

/* Returns cons cell (quotient . remainder); cdr is a fixnum
 */
Object Bignum_Divide (x, y) Object x, y; {   /* bignum / bignum */
    struct S_Bignum *dend, *dor;
    int quotsize, dendsize, dorsize, scale;
    unsigned dor1, dor2;
    Object quot, rem;
    register gran_t *qp, *dendp;
    GC_Node2;

    if (BIGNUM(y)->usize < 2)
	return Bignum_Fixnum_Divide (x, Make_Fixnum (Bignum_To_Integer (y)));

    GC_Link2 (x, y);
    quotsize = BIGNUM(x)->usize - BIGNUM(y)->usize + 1;
    if (quotsize < 0)
	quotsize = 0;
    quot = Make_Uninitialized_Bignum (quotsize);
    GC_Unlink;

    dendsize = (sizeof (struct S_Bignum) - sizeof (gran_t))
	+ (BIGNUM(x)->usize + 1) * sizeof (gran_t);
    dend = (struct S_Bignum *)alloca (dendsize);
    bcopy ((char *)POINTER(x), (char *)dend, dendsize);
    dend->size = BIGNUM(x)->usize + 1;

    if (quotsize == 0 || Bignum_Mantissa_Cmp (dend, BIGNUM(y)) < 0)
	goto zero;

    dorsize = (sizeof (struct S_Bignum) - sizeof (gran_t))
	+ BIGNUM (y)->usize * sizeof (gran_t);
    dor = (struct S_Bignum *)alloca (dorsize);
    bcopy ((char *)POINTER(y), (char *)dor, dorsize);
    dor->size = dorsize = BIGNUM(y)->usize;

    scale = 65536 / (dor->data[dor->usize - 1] + 1);
    Bignum_Mult_In_Place (dend, scale);
    if (dend->usize < dend->size)
	dend->data[dend->usize++] = 0;
    Bignum_Mult_In_Place (dor, scale);

    BIGNUM(quot)->usize = BIGNUM(quot)->size;
    qp = BIGNUM(quot)->data + BIGNUM(quot)->size;
    dendp = dend->data + dend->usize;
    dor1 = dor->data[dor->usize - 1];
    dor2 = dor->data[dor->usize - 2];
        
    while (qp > BIGNUM(quot)->data) {
	unsigned msw, guess;
	int k;
	register gran_t *dep, *dop, *edop;
	
	msw = dendp[-1] << 16 | dendp[-2];
	guess = msw / dor1;
	if (guess >= 65536)	/* [65535, 0, 0] / [65535, 65535] */
	    guess = 65535;
	for (;;) {
	    unsigned d1, d2, d3;
	    d3 = dor2 * guess;
	    d2 = dor1 * guess + (d3 >> 16);
	    d3 &= 0xFFFF;
	    d1 = d2 >> 16;
	    d2 &= 0xFFFF;
	    if (d1 < dendp[-1] || (d1 == dendp[-1] &&
				   (d2 < dendp[-2] || (d2 == dendp[-2] &&
						       d3 <= dendp[-3]))))
		break;
	    --guess;
	}
	--dendp;
	k = 0;
	dep = dendp - dorsize;
	for (dop = dor->data, edop = dop + dor->usize; dop < edop; ) {
	    register unsigned prod = *dop++ * guess;
	    k += *dep;
	    k -= prod & 0xFFFF;
	    *dep++ = k;
	    k >>= 16;
	    k -= prod >> 16;
	}
	k += *dep;
	*dep = k;
	if (k < 0) {
	    k = 0;
	    dep = dendp - dorsize;
	    for (dop = dor->data, edop = dop + dor->usize; dop < edop; ) {
		k += *dep + *dop++;
		*dep++ = k;
		k >>= 16;
	    }
	    k += *dep;
	    *dep = k;
	    --guess;
	}
	*--qp = guess;
    }
    
    if (Bignum_Div_In_Place (dend, scale))
	Panic ("Bignum_Div scale");
 zero:
    if (Truep (dend->minusp = BIGNUM(x)->minusp) != Truep (BIGNUM(y)->minusp))
	BIGNUM(quot)->minusp = True;
    Bignum_Normalize_In_Place (BIGNUM(quot));
    Bignum_Normalize_In_Place (dend);
    GC_Link (quot);
    rem = Reduce_Bignum (Copy_S_Bignum (dend));
    GC_Unlink;
    return Cons (Reduce_Bignum (quot), rem);
}