HFILES=\
go.h\
+ mparith.h\
y.tab.h\
YFILES=\
}
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;
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;
typedef struct Mpflt Mpflt;
struct Mpflt
{
- double val;
- uchar ovf;
+ Mpint val;
+ short exp;
};
typedef struct Val Val;
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);
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
void mpdivfltflt(Mpflt *a, Mpflt *b);
void mpnegflt(Mpflt *a);
double mpgetflt(Mpflt *a);
+int Fconv(Fmt*);
+void mpnorm(Mpflt *a);
/*
* subr.c
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();
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);
}
// 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"
+#include "go.h"
/// uses arithmetic
char buf[500];
Mpflt c;
- sprint(buf, "%B", a);
+ snprint(buf, sizeof(buf), "%B", a);
mpatoflt(&c, buf);
return mpcmpfltflt(&c, b);
}
char buf[500];
Mpflt c;
- sprint(buf, "%B", b);
+ snprint(buf, sizeof(buf), "%B", b);
mpatoflt(&c, buf);
return mpcmpfltflt(a, &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
void
mpmovefixflt(Mpflt *a, Mpint *b)
{
- mpmovecflt(a, mpgetfixflt(b));
+ a->val = *b;
+ a->exp = 0;
+ mpnorm(a);
}
void
*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);
}
//
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 */
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)
{
*--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);
+}
// 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
}
}
+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
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)
{
return;
}
- while(s >= Mpscale) {
- mplshw(a);
- s -= Mpscale;
- }
- while(s > 0) {
- mplsh(a);
- s--;
- }
+ mpshiftfix(a, s);
}
void
return;
}
- while(s >= Mpscale) {
- mprshw(a);
- s -= Mpscale;
- }
- while(s > 0) {
- mprsh(a);
- s--;
- }
+ mpshiftfix(a, -s);
}
void
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)
{
}
}
+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)
{
// 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<<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;
}