// Copyright 2009 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

#include	<u.h>
#include	<libc.h>
#include	"go.h"

//
// return the significant
// words of the argument
//
static int
mplen(Mpint *a)
{
	int i, n;
	long *a1;

	n = -1;
	a1 = &a->a[0];
	for(i=0; i<Mpprec; i++) {
		if(*a1++ != 0)
			n = i;
	}
	return n+1;
}

//
// left shift mpint by one
// ignores sign
//
static void
mplsh(Mpint *a, int quiet)
{
	long *a1, x;
	int i, c;

	c = 0;
	a1 = &a->a[0];
	for(i=0; i<Mpprec; i++) {
		x = (*a1 << 1) + c;
		c = 0;
		if(x >= Mpbase) {
			x -= Mpbase;
			c = 1;
		}
		*a1++ = x;
	}
	a->ovf = c;
	if(a->ovf && !quiet)
		yyerror("constant shift overflow");
}

//
// left shift mpint by Mpscale
// ignores sign
//
static void
mplshw(Mpint *a, int quiet)
{
	long *a1;
	int i;

	a1 = &a->a[Mpprec-1];
	if(*a1) {
		a->ovf = 1;
		if(!quiet)
			yyerror("constant shift overflow");
	}
	for(i=1; i<Mpprec; i++) {
		a1[0] = a1[-1];
		a1--;
	}
	a1[0] = 0;
}

//
// right shift mpint by one
// ignores sign and overflow
//
static void
mprsh(Mpint *a)
{
	long *a1, x, lo;
	int i, c;

	c = 0;
	lo = a->a[0] & 1;
	a1 = &a->a[Mpprec];
	for(i=0; i<Mpprec; i++) {
		x = *--a1;
		*a1 = (x + c) >> 1;
		c = 0;
		if(x & 1)
			c = Mpbase;
	}
	if(a->neg && lo != 0)
		mpaddcfix(a, -1);
}

//
// right shift mpint by Mpscale
// ignores sign and overflow
//
static void
mprshw(Mpint *a)
{
	long *a1, lo;
	int i;

	lo = a->a[0];
	a1 = &a->a[0];
	for(i=1; i<Mpprec; i++) {
		a1[0] = a1[1];
		a1++;
	}
	a1[0] = 0;
	if(a->neg && lo != 0)
		mpaddcfix(a, -1);
}

//
// return the sign of (abs(a)-abs(b))
//
static int
mpcmp(Mpint *a, Mpint *b)
{
	long x, *a1, *b1;
	int i;

	if(a->ovf || b->ovf) {
		if(nsavederrors+nerrors == 0)
			yyerror("ovf in cmp");
		return 0;
	}

	a1 = &a->a[0] + Mpprec;
	b1 = &b->a[0] + Mpprec;

	for(i=0; i<Mpprec; i++) {
		x = *--a1 - *--b1;
		if(x > 0)
			return +1;
		if(x < 0)
			return -1;
	}
	return 0;
}

//
// negate a
// ignore sign and ovf
//
static void
mpneg(Mpint *a)
{
	long x, *a1;
	int i, c;

	a1 = &a->a[0];
	c = 0;
	for(i=0; i<Mpprec; i++) {
		x = -*a1 -c;
		c = 0;
		if(x < 0) {
			x += Mpbase;
			c = 1;
		}
		*a1++ = x;
	}
}

// shift left by s (or right by -s)
void
mpshiftfix(Mpint *a, int s)
{
	if(s >= 0) {
		while(s >= Mpscale) {
			mplshw(a, 0);
			s -= Mpscale;
		}
		while(s > 0) {
			mplsh(a, 0);
			s--;
		}
	} else {
		s = -s;
		while(s >= Mpscale) {
			mprshw(a);
			s -= Mpscale;
		}
		while(s > 0) {
			mprsh(a);
			s--;
		}
	}
}

/// implements fix arihmetic

void
mpaddfixfix(Mpint *a, Mpint *b, int quiet)
{
	int i, c;
	long x, *a1, *b1;

	if(a->ovf || b->ovf) {
		if(nsavederrors+nerrors == 0)
			yyerror("ovf in mpaddxx");
		a->ovf = 1;
		return;
	}

	c = 0;
	a1 = &a->a[0];
	b1 = &b->a[0];
	if(a->neg != b->neg)
		goto sub;

	// perform a+b
	for(i=0; i<Mpprec; i++) {
		x = *a1 + *b1++ + c;
		c = 0;
		if(x >= Mpbase) {
			x -= Mpbase;
			c = 1;
		}
		*a1++ = x;
	}
	a->ovf = c;
	if(a->ovf && !quiet)
		yyerror("constant addition overflow");

	return;

sub:
	// perform a-b
	switch(mpcmp(a, b)) {
	case 0:
		mpmovecfix(a, 0);
		break;

	case 1:
		for(i=0; i<Mpprec; i++) {
			x = *a1 - *b1++ - c;
			c = 0;
			if(x < 0) {
				x += Mpbase;
				c = 1;
			}
			*a1++ = x;
		}
		break;

	case -1:
		a->neg ^= 1;
		for(i=0; i<Mpprec; i++) {
			x = *b1++ - *a1 - c;
			c = 0;
			if(x < 0) {
				x += Mpbase;
				c = 1;
			}
			*a1++ = x;
		}
		break;
	}
}

void
mpmulfixfix(Mpint *a, Mpint *b)
{

	int i, j, na, nb;
	long *a1, x;
	Mpint s, q;

	if(a->ovf || b->ovf) {
		if(nsavederrors+nerrors == 0)
			yyerror("ovf in mpmulfixfix");
		a->ovf = 1;
		return;
	}

	// pick the smaller
	// to test for bits
	na = mplen(a);
	nb = mplen(b);
	if(na > nb) {
		mpmovefixfix(&s, a);
		a1 = &b->a[0];
		na = nb;
	} else {
		mpmovefixfix(&s, b);
		a1 = &a->a[0];
	}
	s.neg = 0;

	mpmovecfix(&q, 0);
	for(i=0; i<na; i++) {
		x = *a1++;
		for(j=0; j<Mpscale; j++) {
			if(x & 1)
				mpaddfixfix(&q, &s, 1);
			mplsh(&s, 1);
			x >>= 1;
		}
	}

	q.neg = a->neg ^ b->neg;
	mpmovefixfix(a, &q);
	if(a->ovf)
		yyerror("constant multiplication overflow");
}

void
mpmulfract(Mpint *a, Mpint *b)
{

	int i, j;
	long *a1, x;
	Mpint s, q;

	if(a->ovf || b->ovf) {
		if(nsavederrors+nerrors == 0)
			yyerror("ovf in mpmulflt");
		a->ovf = 1;
		return;
	}

	mpmovefixfix(&s, b);
	a1 = &a->a[Mpprec];
	s.neg = 0;
	mpmovecfix(&q, 0);

	x = *--a1;
	if(x != 0)
		yyerror("mpmulfract not normal");

	for(i=0; i<Mpprec-1; i++) {
		x = *--a1;
		if(x == 0) {
			mprshw(&s);
			continue;
		}
		for(j=0; j<Mpscale; j++) {
			x <<= 1;
			if(x & Mpbase)
				mpaddfixfix(&q, &s, 1);
			mprsh(&s);
		}
	}

	q.neg = a->neg ^ b->neg;
	mpmovefixfix(a, &q);
	if(a->ovf)
		yyerror("constant multiplication overflow");
}

void
mporfixfix(Mpint *a, Mpint *b)
{
	int i;
	long x, *a1, *b1;

	x = 0;
	if(a->ovf || b->ovf) {
		if(nsavederrors+nerrors == 0)
			yyerror("ovf in mporfixfix");
		mpmovecfix(a, 0);
		a->ovf = 1;
		return;
	}
	if(a->neg) {
		a->neg = 0;
		mpneg(a);
	}
	if(b->neg)
		mpneg(b);

	a1 = &a->a[0];
	b1 = &b->a[0];
	for(i=0; i<Mpprec; i++) {
		x = *a1 | *b1++;
		*a1++ = x;
	}

	if(b->neg)
		mpneg(b);
	if(x & Mpsign) {
		a->neg = 1;
		mpneg(a);
	}
}

void
mpandfixfix(Mpint *a, Mpint *b)
{
	int i;
	long x, *a1, *b1;

	x = 0;
	if(a->ovf || b->ovf) {
		if(nsavederrors+nerrors == 0)
			yyerror("ovf in mpandfixfix");
		mpmovecfix(a, 0);
		a->ovf = 1;
		return;
	}
	if(a->neg) {
		a->neg = 0;
		mpneg(a);
	}
	if(b->neg)
		mpneg(b);

	a1 = &a->a[0];
	b1 = &b->a[0];
	for(i=0; i<Mpprec; i++) {
		x = *a1 & *b1++;
		*a1++ = x;
	}

	if(b->neg)
		mpneg(b);
	if(x & Mpsign) {
		a->neg = 1;
		mpneg(a);
	}
}

void
mpandnotfixfix(Mpint *a, Mpint *b)
{
	int i;
	long x, *a1, *b1;

	x = 0;
	if(a->ovf || b->ovf) {
		if(nsavederrors+nerrors == 0)
			yyerror("ovf in mpandnotfixfix");
		mpmovecfix(a, 0);
		a->ovf = 1;
		return;
	}
	if(a->neg) {
		a->neg = 0;
		mpneg(a);
	}
	if(b->neg)
		mpneg(b);

	a1 = &a->a[0];
	b1 = &b->a[0];
	for(i=0; i<Mpprec; i++) {
		x = *a1 & ~*b1++;
		*a1++ = x;
	}

	if(b->neg)
		mpneg(b);
	if(x & Mpsign) {
		a->neg = 1;
		mpneg(a);
	}
}

void
mpxorfixfix(Mpint *a, Mpint *b)
{
	int i;
	long x, *a1, *b1;

	x = 0;
	if(a->ovf || b->ovf) {
		if(nsavederrors+nerrors == 0)
			yyerror("ovf in mporfixfix");
		mpmovecfix(a, 0);
		a->ovf = 1;
		return;
	}
	if(a->neg) {
		a->neg = 0;
		mpneg(a);
	}
	if(b->neg)
		mpneg(b);

	a1 = &a->a[0];
	b1 = &b->a[0];
	for(i=0; i<Mpprec; i++) {
		x = *a1 ^ *b1++;
		*a1++ = x;
	}

	if(b->neg)
		mpneg(b);
	if(x & Mpsign) {
		a->neg = 1;
		mpneg(a);
	}
}

void
mplshfixfix(Mpint *a, Mpint *b)
{
	vlong s;

	if(a->ovf || b->ovf) {
		if(nsavederrors+nerrors == 0)
			yyerror("ovf in mporfixfix");
		mpmovecfix(a, 0);
		a->ovf = 1;
		return;
	}
	s = mpgetfix(b);
	if(s < 0 || s >= Mpprec*Mpscale) {
		yyerror("stupid shift: %lld", s);
		mpmovecfix(a, 0);
		return;
	}

	mpshiftfix(a, s);
}

void
mprshfixfix(Mpint *a, Mpint *b)
{
	vlong s;

	if(a->ovf || b->ovf) {
		if(nsavederrors+nerrors == 0)
			yyerror("ovf in mprshfixfix");
		mpmovecfix(a, 0);
		a->ovf = 1;
		return;
	}
	s = mpgetfix(b);
	if(s < 0 || s >= Mpprec*Mpscale) {
		yyerror("stupid shift: %lld", s);
		if(a->neg)
			mpmovecfix(a, -1);
		else
			mpmovecfix(a, 0);
		return;
	}

	mpshiftfix(a, -s);
}

void
mpnegfix(Mpint *a)
{
	a->neg ^= 1;
}

vlong
mpgetfix(Mpint *a)
{
	vlong v;

	if(a->ovf) {
		if(nsavederrors+nerrors == 0)
			yyerror("constant overflow");
		return 0;
	}

	v = (vlong)a->a[0];
	v |= (vlong)a->a[1] << Mpscale;
	v |= (vlong)a->a[2] << (Mpscale+Mpscale);
	if(a->neg)
		v = -v;
	return v;
}

void
mpmovecfix(Mpint *a, vlong c)
{
	int i;
	long *a1;
	vlong x;

	a->neg = 0;
	a->ovf = 0;

	x = c;
	if(x < 0) {
		a->neg = 1;
		x = -x;
	}

	a1 = &a->a[0];
	for(i=0; i<Mpprec; i++) {
		*a1++ = x&Mpmask;
		x >>= Mpscale;
	}
}

void
mpdivmodfixfix(Mpint *q, Mpint *r, Mpint *n, Mpint *d)
{
	int i, ns, ds;

	ns = n->neg;
	ds = d->neg;
	n->neg = 0;
	d->neg = 0;

	mpmovefixfix(r, n);
	mpmovecfix(q, 0);

	// shift denominator until it
	// is larger than numerator
	for(i=0; i<Mpprec*Mpscale; i++) {
		if(mpcmp(d, r) > 0)
			break;
		mplsh(d, 1);
	}

	// if it never happens
	// denominator is probably zero
	if(i >= Mpprec*Mpscale) {
		q->ovf = 1;
		r->ovf = 1;
		n->neg = ns;
		d->neg = ds;
		yyerror("constant division overflow");
		return;
	}

	// shift denominator back creating
	// quotient a bit at a time
	// when done the remaining numerator
	// will be the remainder
	for(; i>0; i--) {
		mplsh(q, 1);
		mprsh(d);
		if(mpcmp(d, r) <= 0) {
			mpaddcfix(q, 1);
			mpsubfixfix(r, d);
		}
	}

	n->neg = ns;
	d->neg = ds;
	r->neg = ns;
	q->neg = ns^ds;
}

static int
iszero(Mpint *a)
{
	long *a1;
	int i;
	a1 = &a->a[0] + Mpprec;
	for(i=0; i<Mpprec; i++) {
		if(*--a1 != 0)
			return 0;
	}
	return 1;
}

void
mpdivfract(Mpint *a, Mpint *b)
{
	Mpint n, d;
	int i, j, neg;
	long *a1, x;

	mpmovefixfix(&n, a);	// numerator
	mpmovefixfix(&d, b);	// denominator
	a1 = &a->a[Mpprec];	// quotient

	neg = n.neg ^ d.neg;
	n.neg = 0;
	d.neg = 0;
	for(i=0; i<Mpprec; i++) {
		x = 0;
		for(j=0; j<Mpscale; j++) {
			x <<= 1;
			if(mpcmp(&d, &n) <= 0) {
				if(!iszero(&d))
					x |= 1;
				mpsubfixfix(&n, &d);
			}
			mprsh(&d);
		}
		*--a1 = x;
	}
	a->neg = neg;
}

int
mptestfix(Mpint *a)
{
	Mpint b;
	int r;

	mpmovecfix(&b, 0);
	r = mpcmp(a, &b);
	if(a->neg) {
		if(r > 0)
			return -1;
		if(r < 0)
			return +1;
	}
	return r;
}