multi precision floating point

R=r
OCL=20185
CL=20185

src/cmd/gc/Makefile

@@ -9,6 +9,7 @@ LIB=\
 
 HFILES=\
 	go.h\
+	mparith.h\
 	y.tab.h\
 
 YFILES=\

src/cmd/gc/const.c

@@ -126,10 +126,28 @@ convlit(Node *n, Type *t)
 		}
 		if(isfloat[et]) {
 			// float to float
-			if(mpcmpfltflt(n->val.u.fval, minfltval[et]) < 0)
+			double d;
+			float f;
+			Mpflt *fv;
+
+			fv = n->val.u.fval;
+			if(mpcmpfltflt(fv, minfltval[et]) < 0)
 				goto bad2;
-			if(mpcmpfltflt(n->val.u.fval, maxfltval[et]) > 0)
+			if(mpcmpfltflt(fv, maxfltval[et]) > 0)
 				goto bad2;
+//			switch(et) {
+//			case TFLOAT64:
+//				d = mpgetflt(fv);
+//				mpmovecflt(fv, d);
+//				break;
+//
+//			case TFLOAT32:
+//				d = mpgetflt(fv);
+//				f = d;
+//				d = f;
+//				mpmovecflt(fv, d);
+//				break;
+//			}
 			break;
 		}
 		goto bad1;

src/cmd/gc/go.h

@@ -76,18 +76,18 @@ struct	Array
 
 enum
 {
-	Mpscale	= 29,		/* safely smaller than bits in a long */
-	Mpprec	= 10,		/* Mpscale*Mpprec is max number of bits */
-	Mpbase	= 1L<<Mpscale,
+	Mpscale	= 29,		// safely smaller than bits in a long
+	Mpprec	= 16,		// Mpscale*Mpprec is max number of bits
+	Mpnorm	= Mpprec - 1,	// significant words in a normalized float
+	Mpbase	= 1L << Mpscale,
 	Mpsign	= Mpbase >> 1,
-	Mpmask	= Mpbase -1,
-	Debug	= 1,
+	Mpmask	= Mpbase - 1,
+	Mpdebug	= 0,
 };
 
 typedef	struct	Mpint	Mpint;
 struct	Mpint
 {
-	vlong	val;
 	long	a[Mpprec];
 	uchar	neg;
 	uchar	ovf;
@@ -96,8 +96,8 @@ struct	Mpint
 typedef	struct	Mpflt	Mpflt;
 struct	Mpflt
 {
-	double	val;
-	uchar	ovf;
+	Mpint	val;
+	short	exp;
 };
 
 typedef	struct	Val	Val;
@@ -551,7 +551,9 @@ void	mpmovecfix(Mpint *a, vlong v);
 int	mptestfix(Mpint *a);
 void	mpaddfixfix(Mpint *a, Mpint *b);
 void	mpmulfixfix(Mpint *a, Mpint *b);
+void	mpmulfract(Mpint *a, Mpint *b);
 void	mpdivmodfixfix(Mpint *q, Mpint *r, Mpint *n, Mpint *d);
+void	mpdivfract(Mpint *a, Mpint *b);
 void	mpnegfix(Mpint *a);
 void	mpandfixfix(Mpint *a, Mpint *b);
 void	mplshfixfix(Mpint *a, Mpint *b);
@@ -560,7 +562,7 @@ void	mprshfixfix(Mpint *a, Mpint *b);
 void	mpxorfixfix(Mpint *a, Mpint *b);
 void	mpcomfix(Mpint *a);
 vlong	mpgetfix(Mpint *a);
-double	mpgetfixflt(Mpint *a);
+void	mpshiftfix(Mpint *a, int s);
 
 /*
  *	mparith3.c
@@ -573,6 +575,8 @@ void	mpmulfltflt(Mpflt *a, Mpflt *b);
 void	mpdivfltflt(Mpflt *a, Mpflt *b);
 void	mpnegflt(Mpflt *a);
 double	mpgetflt(Mpflt *a);
+int	Fconv(Fmt*);
+void	mpnorm(Mpflt *a);
 
 /*
  *	subr.c

src/cmd/gc/lex.c

@@ -53,6 +53,7 @@ mainlex(int argc, char *argv[])
 	fmtinstall('Z', Zconv);		// escaped string
 	fmtinstall('L', Lconv);		// line number
 	fmtinstall('B', Bconv);		// big numbers
+	fmtinstall('F', Fconv);		// big float numbers
 	fmtinstall('W', Wconv);		// whatis numbers (Wlitint)
 
 	lexinit();
@@ -786,7 +787,7 @@ caseout:
 
 	yylval.val.u.fval = mal(sizeof(*yylval.val.u.fval));
 	mpatoflt(yylval.val.u.fval, namebuf);
-	if(yylval.val.u.fval->ovf) {
+	if(yylval.val.u.fval->val.ovf) {
 		yyerror("overflow in float constant");
 		mpmovecflt(yylval.val.u.fval, 0.0);
 	}

src/cmd/gc/mparith1.c

@@ -2,8 +2,6 @@
 // Use of this source code is governed by a BSD-style
 // license that can be found in the LICENSE file.
 
-#include <u.h>
-#include <errno.h>
 #include	"go.h"
 
 /// uses arithmetic
@@ -14,7 +12,7 @@ mpcmpfixflt(Mpint *a, Mpflt *b)
 	char buf[500];
 	Mpflt c;
 
-	sprint(buf, "%B", a);
+	snprint(buf, sizeof(buf), "%B", a);
 	mpatoflt(&c, buf);
 	return mpcmpfltflt(&c, b);
 }
@@ -25,7 +23,7 @@ mpcmpfltfix(Mpflt *a, Mpint *b)
 	char buf[500];
 	Mpflt c;
 
-	sprint(buf, "%B", b);
+	snprint(buf, sizeof(buf), "%B", b);
 	mpatoflt(&c, buf);
 	return mpcmpfltflt(a, &c);
 }
@@ -71,17 +69,17 @@ mpcmpfltc(Mpflt *b, double c)
 void
 mpsubfixfix(Mpint *a, Mpint *b)
 {
-	mpnegfix(b);
+	mpnegfix(a);
 	mpaddfixfix(a, b);
-	mpnegfix(b);
+	mpnegfix(a);
 }
 
 void
 mpsubfltflt(Mpflt *a, Mpflt *b)
 {
-	mpnegflt(b);
+	mpnegflt(a);
 	mpaddfltflt(a, b);
-	mpnegflt(b);
+	mpnegflt(a);
 }
 
 void
@@ -151,7 +149,9 @@ mpcomfix(Mpint *a)
 void
 mpmovefixflt(Mpflt *a, Mpint *b)
 {
-	mpmovecflt(a, mpgetfixflt(b));
+	a->val = *b;
+	a->exp = 0;
+	mpnorm(a);
 }
 
 void
@@ -172,25 +172,18 @@ mpmovefltflt(Mpflt *a, Mpflt *b)
 	*a = *b;
 }
 
-//
-// power of ten
-//
-static	double
-tentab[] = { 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9 };
-
-static double
-dppow10(int n)
+static	double	tab[] = { 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7 };
+static void
+mppow10flt(Mpflt *a, int p)
 {
-	int i;
-
-	if(n < 0)
-		return 1.0/dppow10(-n);
-
-	if(n < nelem(tentab))
-		return tentab[n];
-
-	i = n/2;
-	return dppow10(i) * dppow10(n-i);
+	if(p < nelem(tab)) {
+		mpmovecflt(a, tab[p]);
+		return;
+	}
+	mppow10flt(a, p>>1);
+	mpmulfltflt(a, a);
+	if(p & 1)
+		mpmulcflt(a, 10);
 }
 
 //
@@ -200,17 +193,9 @@ dppow10(int n)
 void
 mpatoflt(Mpflt *a, char *as)
 {
+	Mpflt b;
 	int dp, c, f, ef, ex, zer;
 	char *s;
-	double f64;
-
-	/* until Mpflt is really mp, use strtod to get rounding right */
-	errno = 0;
-	f64 = strtod(as, &s);
-	mpmovecflt(a, f64);
-	if(errno != 0)
-		a->ovf = 1;
-	return;
 
 	s = as;
 	dp = 0;		/* digits after decimal point */
@@ -283,14 +268,21 @@ mpatoflt(Mpflt *a, char *as)
 
 	if(dp)
 		dp--;
-	if(mpcmpfltc(a, 0.0) != 0)
-		mpmulcflt(a, dppow10(ex-dp));
+	if(mpcmpfltc(a, 0.0) != 0) {
+		if(ex >= dp) {
+			mppow10flt(&b, ex-dp);
+			mpmulfltflt(a, &b);
+		} else {
+			mppow10flt(&b, dp-ex);
+			mpdivfltflt(a, &b);
+		}
+	}
 	if(f)
 		mpnegflt(a);
 	return;
 
 bad:
-	warn("set ovf in mpatof: %s", as);
+	warn("set ovf in mpatof");
 	mpmovecflt(a, 0.0);
 }
 
@@ -410,3 +402,17 @@ Bconv(Fmt *fp)
 		*--p = '-';
 	return fmtstrcpy(fp, p);
 }
+
+int
+Fconv(Fmt *fp)
+{
+	char buf[500];
+	Mpflt *fval;
+
+	fval = va_arg(fp->args, Mpflt*);
+	if(fval->exp >= 0)
+		snprint(buf, sizeof(buf), "(%B*2^%d)", &fval->val, fval->exp);
+	else
+		snprint(buf, sizeof(buf), "(%B/2^%d)", &fval->val, -fval->exp);
+	return fmtstrcpy(fp, buf);
+}

src/cmd/gc/mparith2.c

@@ -159,6 +159,31 @@ mpneg(Mpint *a)
 	}
 }
 
+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
@@ -274,6 +299,45 @@ mpmulfixfix(Mpint *a, Mpint *b)
 		warn("set ovf in mpmulfixfix");
 }
 
+void
+mpmulfract(Mpint *a, Mpint *b)
+{
+
+	int i, j;
+	long *a1, x;
+	Mpint s, q;
+
+	if(a->ovf || b->ovf) {
+		warn("ovf in mpmulflt");
+		a->ovf = 1;
+		return;
+	}
+
+	mpmovefixfix(&s, b);
+	a1 = &a->a[Mpprec];
+	s.neg = 0;
+	mpmovecfix(&q, 0);
+
+	for(i=0; i<Mpprec; 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)
+		warn("set ovf in mpmulflt");
+}
+
 void
 mporfixfix(Mpint *a, Mpint *b)
 {
@@ -394,14 +458,7 @@ mplshfixfix(Mpint *a, Mpint *b)
 		return;
 	}
 
-	while(s >= Mpscale) {
-		mplshw(a);
-		s -= Mpscale;
-	}
-	while(s > 0) {
-		mplsh(a);
-		s--;
-	}
+	mpshiftfix(a, s);
 }
 
 void
@@ -425,14 +482,7 @@ mprshfixfix(Mpint *a, Mpint *b)
 		return;
 	}
 
-	while(s >= Mpscale) {
-		mprshw(a);
-		s -= Mpscale;
-	}
-	while(s > 0) {
-		mprsh(a);
-		s--;
-	}
+	mpshiftfix(a, -s);
 }
 
 void
@@ -459,17 +509,6 @@ mpgetfix(Mpint *a)
 	return v;
 }
 
-double
-mpgetfixflt(Mpint *a)
-{
-	// answer might not fit in intermediate vlong, so format
-	// to string and then let the string routine convert.
-	char buf[1000];
-
-	snprint(buf, sizeof buf, "%B", a);
-	return strtod(buf, nil);
-}
-
 void
 mpmovecfix(Mpint *a, vlong c)
 {
@@ -532,6 +571,36 @@ mpdivmodfixfix(Mpint *q, Mpint *r, Mpint *n, Mpint *d)
 	}
 }
 
+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) {
+				x |= 1;
+				mpsubfixfix(&n, &d);
+			}
+			mprsh(&d);
+		}
+		*--a1 = x;
+	}
+	a->neg = neg;
+}
+
 int
 mptestfix(Mpint *a)
 {

src/cmd/gc/mparith3.c

@@ -4,50 +4,264 @@
 
 #include	"go.h"
 
+/*
+ * returns the leading non-zero
+ * word of the number
+ */
+int
+sigfig(Mpflt *a)
+{
+	int i;
+
+	for(i=Mpprec-1; i>=0; i--)
+		if(a->val.a[i] != 0)
+			break;
+//print("sigfig %d %d\n", i-z+1, z);
+	return i+1;
+}
+
+/*
+ * shifts the leading non-zero
+ * word of the number to Mpnorm
+ */
+void
+mpnorm(Mpflt *a)
+{
+	int s;
+
+	s = sigfig(a);
+	if(s == 0) {
+		// zero
+		a->exp = 0;
+		a->val.neg = 0;
+		return;
+	}
+	s = (Mpnorm-s) * Mpscale;
+	mpshiftfix(&a->val, s);
+	a->exp -= s;
+}
+
 /// implements float arihmetic
 
 void
 mpaddfltflt(Mpflt *a, Mpflt *b)
 {
-	a->val += b->val;
+	int sa, sb, s;
+	Mpflt c;
+
+	if(Mpdebug)
+		print("\n%F + %F", a, b);
+
+	sa = sigfig(a);
+	sb = sigfig(b);
+
+	if(sa == 0) {
+		if(sb == 0) {
+			// zero
+			a->exp = 0;
+			a->val.neg = 0;
+			return;
+		}
+		mpmovefltflt(a, b);
+		goto out;
+	}
+	if(sb == 0)
+		goto out;
+
+	s = a->exp - b->exp;
+	if(s > 0) {
+		// a is larger, shift b right
+		mpmovefltflt(&c, b);
+		mpshiftfix(&c.val, -s);
+		mpaddfixfix(&a->val, &c.val);
+		goto out;
+	}
+	if(s < 0) {
+		// b is larger, shift a right
+		mpshiftfix(&a->val, s);
+		a->exp -= s;
+		mpaddfixfix(&a->val, &b->val);
+		goto out;
+	}
+	mpaddfixfix(&a->val, &b->val);
+
+out:
+	mpnorm(a);
+	if(Mpdebug)
+		print(" = %F\n\n", a);
 }
 
 void
 mpmulfltflt(Mpflt *a, Mpflt *b)
 {
-	a->val *= b->val;
+	int sa, sb;
+
+	if(Mpdebug)
+		print("%F\n * %F\n", a, b);
+
+	sa = sigfig(a);
+	sb = sigfig(b);
+
+	if(sa == 0 || sb == 0) {
+		// zero
+		a->exp = 0;
+		a->val.neg = 0;
+		return;
+	}
+
+	mpmulfract(&a->val, &b->val);
+	a->exp = (a->exp + b->exp) + Mpscale*Mpprec - 1;
+
+	mpnorm(a);
+	if(Mpdebug)
+		print(" = %F\n\n", a);
 }
 
 void
 mpdivfltflt(Mpflt *a, Mpflt *b)
 {
-	a->val /= b->val;
+	int sa, sb;
+	Mpflt c;
+
+	if(Mpdebug)
+		print("%F\n / %F\n", a, b);
+
+	sa = sigfig(a);
+	sb = sigfig(b);
+
+	if(sb == 0) {
+		// zero and ovfl
+		a->exp = 0;
+		a->val.neg = 0;
+		a->val.ovf = 1;
+		warn("mpdivfltflt divide by zero");
+		return;
+	}
+	if(sa == 0) {
+		// zero
+		a->exp = 0;
+		a->val.neg = 0;
+		return;
+	}
+
+	// adjust b to top
+	mpmovefltflt(&c, b);
+	mpshiftfix(&c.val, Mpscale);
+
+	// divide
+	mpdivfract(&a->val, &c.val);
+	a->exp = (a->exp-c.exp) - Mpscale*(Mpprec-1) + 1;
+
+	mpnorm(a);
+	if(Mpdebug)
+		print(" = %F\n\n", a);
 }
 
 double
 mpgetflt(Mpflt *a)
 {
-	return a->val;
+	int s, i;
+	uvlong v;
+	double f;
+
+	if(a->val.ovf)
+		warn("mpgetflt ovf");
+
+	s = sigfig(a);
+	if(s == 0)
+		return 0;
+
+	if(s != Mpnorm) {
+		warn("mpgetflt norm");
+		mpnorm(a);
+	}
+
+	while((a->val.a[Mpnorm-1] & (1L<<(Mpscale-1))) == 0) {
+		mpshiftfix(&a->val, 1);
+		a->exp -= 1;
+	}
+
+	// pick up the mantissa in a uvlong
+	s = 63;
+	v = 0;
+	for(i=Mpnorm-1; s>=Mpscale; i--) {
+		v = (v<<Mpscale) | a->val.a[i];
+		s -= Mpscale;
+	}
+	if(s > 0)
+		v = (v<<s) | (a->val.a[i]>>(Mpscale-s));
+
+	// should do this in multi precision
+	// 63 bits of mantissa being rounded to 53
+	if((v&0x3ffULL) != 0x200ULL || (v&0x400) != 0)
+		v += 0x200ULL;		// round
+	v &= ~0x3ffULL;
+	f = (double)(v);
+	f = ldexp(f, Mpnorm*Mpscale + a->exp - 63);
+
+	if(a->val.neg)
+		f = -f;
+	return f;
 }
 
 void
 mpmovecflt(Mpflt *a, double c)
 {
-	a->val = c;
+	int i;
+	double f;
+	long l;
+
+	if(Mpdebug)
+		print("\nconst %g", c);
+	mpmovecfix(&a->val, 0);
+	a->exp = 0;
+	if(c == 0)
+		goto out;
+	if(c < 0) {
+		a->val.neg = 1;
+		c = -c;
+	}
+
+	f = frexp(c, &i);
+	a->exp = i;
+
+	for(i=0; i<10; i++) {
+		f = f*Mpbase;
+		l = floor(f);
+		f = f - l;
+		a->exp -= Mpscale;
+		a->val.a[0] = l;
+		if(f == 0)
+			break;
+		mpshiftfix(&a->val, Mpscale);
+	}
+
+out:
+	mpnorm(a);
+	if(Mpdebug)
+		print(" = %F\n", a);
 }
 
 void
 mpnegflt(Mpflt *a)
 {
-	a->val = -a->val;
+	a->val.neg ^= 1;
 }
 
 int
 mptestflt(Mpflt *a)
 {
-	if(a->val < 0)
-		return -1;
-	if(a->val > 0)
-		return +1;
-	return 0;
+	int s;
+
+	if(Mpdebug)
+		print("\n%F?", a);
+	s = sigfig(a);
+	if(s != 0) {
+		s = +1;
+		if(a->val.neg)
+			s = -1;
+	}
+	if(Mpdebug)
+		print(" = %d\n", s);
+	return s;
 }