From: Ken Thompson Date: Tue, 2 Dec 2008 01:22:05 +0000 (-0800) Subject: multi precision floating point X-Git-Tag: weekly.2009-11-06~2594 X-Git-Url: http://www.git.cypherpunks.su/?a=commitdiff_plain;h=3fa46106017b0eed128e83d4ce084c43efc14d5f;p=gostls13.git multi precision floating point R=r OCL=20185 CL=20185 --- diff --git a/src/cmd/gc/Makefile b/src/cmd/gc/Makefile index 502f37146d..2a764b7aed 100644 --- a/src/cmd/gc/Makefile +++ b/src/cmd/gc/Makefile @@ -9,6 +9,7 @@ LIB=\ HFILES=\ go.h\ + mparith.h\ y.tab.h\ YFILES=\ diff --git a/src/cmd/gc/const.c b/src/cmd/gc/const.c index 7ad2110af4..3a3fef191d 100644 --- a/src/cmd/gc/const.c +++ b/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; diff --git a/src/cmd/gc/go.h b/src/cmd/gc/go.h index ce1d4cee50..2654e1c564 100644 --- a/src/cmd/gc/go.h +++ b/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<> 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 diff --git a/src/cmd/gc/lex.c b/src/cmd/gc/lex.c index 7264b3cad6..f818a641fa 100644 --- a/src/cmd/gc/lex.c +++ b/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); } diff --git a/src/cmd/gc/mparith1.c b/src/cmd/gc/mparith1.c index 98fa661b46..2bf20d60ee 100644 --- a/src/cmd/gc/mparith1.c +++ b/src/cmd/gc/mparith1.c @@ -2,9 +2,7 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -#include -#include -#include "go.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,21 +268,28 @@ 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); } // // fixed point input // required syntax is [+-][0[x]]d* -// +// void mpatofix(Mpint *a, char *as) { @@ -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); +} diff --git a/src/cmd/gc/mparith2.c b/src/cmd/gc/mparith2.c index 186437602d..f7c2ea2c2b 100644 --- a/src/cmd/gc/mparith2.c +++ b/src/cmd/gc/mparith2.c @@ -2,7 +2,7 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -#include "go.h" +#include "go.h" // // return the significant @@ -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; ineg ^ 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; ineg = neg; +} + int mptestfix(Mpint *a) { diff --git a/src/cmd/gc/mparith3.c b/src/cmd/gc/mparith3.c index 2a0a1c6c2e..1bf39c9fb4 100644 --- a/src/cmd/gc/mparith3.c +++ b/src/cmd/gc/mparith3.c @@ -2,52 +2,266 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -#include "go.h" +#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<val.a[i]; + s -= Mpscale; + } + if(s > 0) + v = (v<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; }