summaryrefslogtreecommitdiff
path: root/src/cmd/gc/mparith2.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/cmd/gc/mparith2.c')
-rw-r--r--src/cmd/gc/mparith2.c683
1 files changed, 683 insertions, 0 deletions
diff --git a/src/cmd/gc/mparith2.c b/src/cmd/gc/mparith2.c
new file mode 100644
index 000000000..403255005
--- /dev/null
+++ b/src/cmd/gc/mparith2.c
@@ -0,0 +1,683 @@
+// 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 "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 and overflow
+//
+static void
+mplsh(Mpint *a)
+{
+ 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;
+ }
+}
+
+//
+// left shift mpint by Mpscale
+// ignores sign and overflow
+//
+static void
+mplshw(Mpint *a)
+{
+ long *a1;
+ int i;
+
+ a1 = &a->a[Mpprec-1];
+ 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) {
+ 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);
+ s -= Mpscale;
+ }
+ while(s > 0) {
+ mplsh(a);
+ 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 i, c;
+ long x, *a1, *b1;
+
+ if(a->ovf || b->ovf) {
+ 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)
+ yyerror("set ovf in mpaddxx");
+
+ 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) {
+ 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);
+ mplsh(&s);
+ x >>= 1;
+ }
+ }
+
+ q.neg = a->neg ^ b->neg;
+ mpmovefixfix(a, &q);
+ if(a->ovf)
+ yyerror("set ovf in mpmulfixfix");
+}
+
+void
+mpmulfract(Mpint *a, Mpint *b)
+{
+
+ int i, j;
+ long *a1, x;
+ Mpint s, q;
+
+ if(a->ovf || b->ovf) {
+ 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);
+ mprsh(&s);
+ }
+ }
+
+ q.neg = a->neg ^ b->neg;
+ mpmovefixfix(a, &q);
+ if(a->ovf)
+ yyerror("set ovf in mpmulflt");
+}
+
+void
+mporfixfix(Mpint *a, Mpint *b)
+{
+ int i;
+ long x, *a1, *b1;
+
+ if(a->ovf || b->ovf) {
+ 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;
+
+ if(a->ovf || b->ovf) {
+ 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;
+
+ if(a->ovf || b->ovf) {
+ 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;
+
+ if(a->ovf || b->ovf) {
+ 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) {
+ 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) {
+ 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) {
+ 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);
+ }
+
+ // 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("set ovf in mpdivmodfixfix");
+ 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);
+ 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;
+}