1*74a4d8c2SCharles.Forsyth #include "cc.h"
2*74a4d8c2SCharles.Forsyth
3*74a4d8c2SCharles.Forsyth enum
4*74a4d8c2SCharles.Forsyth {
5*74a4d8c2SCharles.Forsyth Mpscale = 29, /* safely smaller than bits in a long */
6*74a4d8c2SCharles.Forsyth Mpprec = 36, /* Mpscale*Mpprec sb > largest fp exp */
7*74a4d8c2SCharles.Forsyth Mpbase = 1L<<Mpscale,
8*74a4d8c2SCharles.Forsyth };
9*74a4d8c2SCharles.Forsyth
10*74a4d8c2SCharles.Forsyth typedef
11*74a4d8c2SCharles.Forsyth struct
12*74a4d8c2SCharles.Forsyth {
13*74a4d8c2SCharles.Forsyth long a[Mpprec];
14*74a4d8c2SCharles.Forsyth char ovf;
15*74a4d8c2SCharles.Forsyth } Mp;
16*74a4d8c2SCharles.Forsyth
17*74a4d8c2SCharles.Forsyth int mpatof(char*, double*);
18*74a4d8c2SCharles.Forsyth int mpatov(char *s, vlong *v);
19*74a4d8c2SCharles.Forsyth void mpint(Mp*, int);
20*74a4d8c2SCharles.Forsyth void mppow(Mp*, int, int);
21*74a4d8c2SCharles.Forsyth void mpmul(Mp*, int);
22*74a4d8c2SCharles.Forsyth void mpadd(Mp*, Mp*);
23*74a4d8c2SCharles.Forsyth int mptof(Mp*, double*);
24*74a4d8c2SCharles.Forsyth
25*74a4d8c2SCharles.Forsyth /*
26*74a4d8c2SCharles.Forsyth * convert a string, s, to floating in *d
27*74a4d8c2SCharles.Forsyth * return conversion overflow.
28*74a4d8c2SCharles.Forsyth * required syntax is [+-]d*[.]d*[e[+-]d*]
29*74a4d8c2SCharles.Forsyth */
30*74a4d8c2SCharles.Forsyth int
mpatof(char * s,double * d)31*74a4d8c2SCharles.Forsyth mpatof(char *s, double *d)
32*74a4d8c2SCharles.Forsyth {
33*74a4d8c2SCharles.Forsyth Mp a, b;
34*74a4d8c2SCharles.Forsyth int dp, c, f, ef, ex, zer;
35*74a4d8c2SCharles.Forsyth double d1, d2;
36*74a4d8c2SCharles.Forsyth
37*74a4d8c2SCharles.Forsyth dp = 0; /* digits after decimal point */
38*74a4d8c2SCharles.Forsyth f = 0; /* sign */
39*74a4d8c2SCharles.Forsyth ex = 0; /* exponent */
40*74a4d8c2SCharles.Forsyth zer = 1; /* zero */
41*74a4d8c2SCharles.Forsyth memset(&a, 0, sizeof(a));
42*74a4d8c2SCharles.Forsyth for(;;) {
43*74a4d8c2SCharles.Forsyth switch(c = *s++) {
44*74a4d8c2SCharles.Forsyth default:
45*74a4d8c2SCharles.Forsyth goto bad;
46*74a4d8c2SCharles.Forsyth case '-':
47*74a4d8c2SCharles.Forsyth f = 1;
48*74a4d8c2SCharles.Forsyth case ' ':
49*74a4d8c2SCharles.Forsyth case '\t':
50*74a4d8c2SCharles.Forsyth case '+':
51*74a4d8c2SCharles.Forsyth continue;
52*74a4d8c2SCharles.Forsyth case '.':
53*74a4d8c2SCharles.Forsyth dp = 1;
54*74a4d8c2SCharles.Forsyth continue;
55*74a4d8c2SCharles.Forsyth case '1':
56*74a4d8c2SCharles.Forsyth case '2':
57*74a4d8c2SCharles.Forsyth case '3':
58*74a4d8c2SCharles.Forsyth case '4':
59*74a4d8c2SCharles.Forsyth case '5':
60*74a4d8c2SCharles.Forsyth case '6':
61*74a4d8c2SCharles.Forsyth case '7':
62*74a4d8c2SCharles.Forsyth case '8':
63*74a4d8c2SCharles.Forsyth case '9':
64*74a4d8c2SCharles.Forsyth zer = 0;
65*74a4d8c2SCharles.Forsyth case '0':
66*74a4d8c2SCharles.Forsyth mpint(&b, c-'0');
67*74a4d8c2SCharles.Forsyth mpmul(&a, 10);
68*74a4d8c2SCharles.Forsyth mpadd(&a, &b);
69*74a4d8c2SCharles.Forsyth if(dp)
70*74a4d8c2SCharles.Forsyth dp++;
71*74a4d8c2SCharles.Forsyth continue;
72*74a4d8c2SCharles.Forsyth case 'E':
73*74a4d8c2SCharles.Forsyth case 'e':
74*74a4d8c2SCharles.Forsyth ex = 0;
75*74a4d8c2SCharles.Forsyth ef = 0;
76*74a4d8c2SCharles.Forsyth for(;;) {
77*74a4d8c2SCharles.Forsyth c = *s++;
78*74a4d8c2SCharles.Forsyth if(c == '+' || c == ' ' || c == '\t')
79*74a4d8c2SCharles.Forsyth continue;
80*74a4d8c2SCharles.Forsyth if(c == '-') {
81*74a4d8c2SCharles.Forsyth ef = 1;
82*74a4d8c2SCharles.Forsyth continue;
83*74a4d8c2SCharles.Forsyth }
84*74a4d8c2SCharles.Forsyth if(c >= '0' && c <= '9') {
85*74a4d8c2SCharles.Forsyth ex = ex*10 + (c-'0');
86*74a4d8c2SCharles.Forsyth continue;
87*74a4d8c2SCharles.Forsyth }
88*74a4d8c2SCharles.Forsyth break;
89*74a4d8c2SCharles.Forsyth }
90*74a4d8c2SCharles.Forsyth if(ef)
91*74a4d8c2SCharles.Forsyth ex = -ex;
92*74a4d8c2SCharles.Forsyth case 0:
93*74a4d8c2SCharles.Forsyth break;
94*74a4d8c2SCharles.Forsyth }
95*74a4d8c2SCharles.Forsyth break;
96*74a4d8c2SCharles.Forsyth }
97*74a4d8c2SCharles.Forsyth if(a.ovf)
98*74a4d8c2SCharles.Forsyth goto bad;
99*74a4d8c2SCharles.Forsyth if(zer) {
100*74a4d8c2SCharles.Forsyth *d = 0;
101*74a4d8c2SCharles.Forsyth return 0;
102*74a4d8c2SCharles.Forsyth }
103*74a4d8c2SCharles.Forsyth if(dp)
104*74a4d8c2SCharles.Forsyth dp--;
105*74a4d8c2SCharles.Forsyth dp -= ex;
106*74a4d8c2SCharles.Forsyth if(dp > 0) {
107*74a4d8c2SCharles.Forsyth /*
108*74a4d8c2SCharles.Forsyth * must divide by 10**dp
109*74a4d8c2SCharles.Forsyth */
110*74a4d8c2SCharles.Forsyth if(mptof(&a, &d1))
111*74a4d8c2SCharles.Forsyth goto bad;
112*74a4d8c2SCharles.Forsyth
113*74a4d8c2SCharles.Forsyth /*
114*74a4d8c2SCharles.Forsyth * trial exponent of 8**dp
115*74a4d8c2SCharles.Forsyth * 8 (being between 5 and 10)
116*74a4d8c2SCharles.Forsyth * should pick up all underflows
117*74a4d8c2SCharles.Forsyth * in the division of 5**dp.
118*74a4d8c2SCharles.Forsyth */
119*74a4d8c2SCharles.Forsyth d2 = frexp(d1, &ex);
120*74a4d8c2SCharles.Forsyth d2 = ldexp(d2, ex-3*dp);
121*74a4d8c2SCharles.Forsyth if(d2 == 0)
122*74a4d8c2SCharles.Forsyth goto bad;
123*74a4d8c2SCharles.Forsyth
124*74a4d8c2SCharles.Forsyth /*
125*74a4d8c2SCharles.Forsyth * decompose each 10 into 5*2.
126*74a4d8c2SCharles.Forsyth * create 5**dp in fixed point
127*74a4d8c2SCharles.Forsyth * and then play with the exponent
128*74a4d8c2SCharles.Forsyth * for the remaining 2**dp.
129*74a4d8c2SCharles.Forsyth * note that 5**dp will overflow
130*74a4d8c2SCharles.Forsyth * with as few as 134 input digits.
131*74a4d8c2SCharles.Forsyth */
132*74a4d8c2SCharles.Forsyth mpint(&a, 1);
133*74a4d8c2SCharles.Forsyth mppow(&a, 5, dp);
134*74a4d8c2SCharles.Forsyth if(mptof(&a, &d2))
135*74a4d8c2SCharles.Forsyth goto bad;
136*74a4d8c2SCharles.Forsyth d1 = frexp(d1/d2, &ex);
137*74a4d8c2SCharles.Forsyth d1 = ldexp(d1, ex-dp);
138*74a4d8c2SCharles.Forsyth if(d1 == 0)
139*74a4d8c2SCharles.Forsyth goto bad;
140*74a4d8c2SCharles.Forsyth } else {
141*74a4d8c2SCharles.Forsyth /*
142*74a4d8c2SCharles.Forsyth * must multiply by 10**|dp| --
143*74a4d8c2SCharles.Forsyth * just do it in fixed point.
144*74a4d8c2SCharles.Forsyth */
145*74a4d8c2SCharles.Forsyth mppow(&a, 10, -dp);
146*74a4d8c2SCharles.Forsyth if(mptof(&a, &d1))
147*74a4d8c2SCharles.Forsyth goto bad;
148*74a4d8c2SCharles.Forsyth }
149*74a4d8c2SCharles.Forsyth if(f)
150*74a4d8c2SCharles.Forsyth d1 = -d1;
151*74a4d8c2SCharles.Forsyth *d = d1;
152*74a4d8c2SCharles.Forsyth return 0;
153*74a4d8c2SCharles.Forsyth
154*74a4d8c2SCharles.Forsyth bad:
155*74a4d8c2SCharles.Forsyth return 1;
156*74a4d8c2SCharles.Forsyth }
157*74a4d8c2SCharles.Forsyth
158*74a4d8c2SCharles.Forsyth /*
159*74a4d8c2SCharles.Forsyth * convert a to floating in *d
160*74a4d8c2SCharles.Forsyth * return conversion overflow
161*74a4d8c2SCharles.Forsyth */
162*74a4d8c2SCharles.Forsyth int
mptof(Mp * a,double * d)163*74a4d8c2SCharles.Forsyth mptof(Mp *a, double *d)
164*74a4d8c2SCharles.Forsyth {
165*74a4d8c2SCharles.Forsyth double f, g;
166*74a4d8c2SCharles.Forsyth long x, *a1;
167*74a4d8c2SCharles.Forsyth int i;
168*74a4d8c2SCharles.Forsyth
169*74a4d8c2SCharles.Forsyth if(a->ovf)
170*74a4d8c2SCharles.Forsyth return 1;
171*74a4d8c2SCharles.Forsyth a1 = a->a;
172*74a4d8c2SCharles.Forsyth f = ldexp(*a1++, 0);
173*74a4d8c2SCharles.Forsyth for(i=Mpscale; i<Mpprec*Mpscale; i+=Mpscale)
174*74a4d8c2SCharles.Forsyth if(x = *a1++) {
175*74a4d8c2SCharles.Forsyth g = ldexp(x, i);
176*74a4d8c2SCharles.Forsyth /*
177*74a4d8c2SCharles.Forsyth * NOTE: the test (g==0) is plan9
178*74a4d8c2SCharles.Forsyth * specific. ansi compliant overflow
179*74a4d8c2SCharles.Forsyth * is signaled by HUGE and errno==ERANGE.
180*74a4d8c2SCharles.Forsyth * change this for your particular ldexp.
181*74a4d8c2SCharles.Forsyth */
182*74a4d8c2SCharles.Forsyth if(g == 0)
183*74a4d8c2SCharles.Forsyth return 1;
184*74a4d8c2SCharles.Forsyth f += g; /* this could bomb! */
185*74a4d8c2SCharles.Forsyth }
186*74a4d8c2SCharles.Forsyth *d = f;
187*74a4d8c2SCharles.Forsyth return 0;
188*74a4d8c2SCharles.Forsyth }
189*74a4d8c2SCharles.Forsyth
190*74a4d8c2SCharles.Forsyth /*
191*74a4d8c2SCharles.Forsyth * return a += b
192*74a4d8c2SCharles.Forsyth */
193*74a4d8c2SCharles.Forsyth void
mpadd(Mp * a,Mp * b)194*74a4d8c2SCharles.Forsyth mpadd(Mp *a, Mp *b)
195*74a4d8c2SCharles.Forsyth {
196*74a4d8c2SCharles.Forsyth int i, c;
197*74a4d8c2SCharles.Forsyth long x, *a1, *b1;
198*74a4d8c2SCharles.Forsyth
199*74a4d8c2SCharles.Forsyth if(b->ovf)
200*74a4d8c2SCharles.Forsyth a->ovf = 1;
201*74a4d8c2SCharles.Forsyth if(a->ovf)
202*74a4d8c2SCharles.Forsyth return;
203*74a4d8c2SCharles.Forsyth c = 0;
204*74a4d8c2SCharles.Forsyth a1 = a->a;
205*74a4d8c2SCharles.Forsyth b1 = b->a;
206*74a4d8c2SCharles.Forsyth for(i=0; i<Mpprec; i++) {
207*74a4d8c2SCharles.Forsyth x = *a1 + *b1++ + c;
208*74a4d8c2SCharles.Forsyth c = 0;
209*74a4d8c2SCharles.Forsyth if(x >= Mpbase) {
210*74a4d8c2SCharles.Forsyth x -= Mpbase;
211*74a4d8c2SCharles.Forsyth c = 1;
212*74a4d8c2SCharles.Forsyth }
213*74a4d8c2SCharles.Forsyth *a1++ = x;
214*74a4d8c2SCharles.Forsyth }
215*74a4d8c2SCharles.Forsyth a->ovf = c;
216*74a4d8c2SCharles.Forsyth }
217*74a4d8c2SCharles.Forsyth
218*74a4d8c2SCharles.Forsyth /*
219*74a4d8c2SCharles.Forsyth * return a = c
220*74a4d8c2SCharles.Forsyth */
221*74a4d8c2SCharles.Forsyth void
mpint(Mp * a,int c)222*74a4d8c2SCharles.Forsyth mpint(Mp *a, int c)
223*74a4d8c2SCharles.Forsyth {
224*74a4d8c2SCharles.Forsyth
225*74a4d8c2SCharles.Forsyth memset(a, 0, sizeof(*a));
226*74a4d8c2SCharles.Forsyth a->a[0] = c;
227*74a4d8c2SCharles.Forsyth }
228*74a4d8c2SCharles.Forsyth
229*74a4d8c2SCharles.Forsyth /*
230*74a4d8c2SCharles.Forsyth * return a *= c
231*74a4d8c2SCharles.Forsyth */
232*74a4d8c2SCharles.Forsyth void
mpmul(Mp * a,int c)233*74a4d8c2SCharles.Forsyth mpmul(Mp *a, int c)
234*74a4d8c2SCharles.Forsyth {
235*74a4d8c2SCharles.Forsyth Mp p;
236*74a4d8c2SCharles.Forsyth int b;
237*74a4d8c2SCharles.Forsyth
238*74a4d8c2SCharles.Forsyth memmove(&p, a, sizeof(p));
239*74a4d8c2SCharles.Forsyth if(!(c & 1))
240*74a4d8c2SCharles.Forsyth memset(a, 0, sizeof(*a));
241*74a4d8c2SCharles.Forsyth c &= ~1;
242*74a4d8c2SCharles.Forsyth for(b=2; c; b<<=1) {
243*74a4d8c2SCharles.Forsyth mpadd(&p, &p);
244*74a4d8c2SCharles.Forsyth if(c & b) {
245*74a4d8c2SCharles.Forsyth mpadd(a, &p);
246*74a4d8c2SCharles.Forsyth c &= ~b;
247*74a4d8c2SCharles.Forsyth }
248*74a4d8c2SCharles.Forsyth }
249*74a4d8c2SCharles.Forsyth }
250*74a4d8c2SCharles.Forsyth
251*74a4d8c2SCharles.Forsyth /*
252*74a4d8c2SCharles.Forsyth * return a *= b**e
253*74a4d8c2SCharles.Forsyth */
254*74a4d8c2SCharles.Forsyth void
mppow(Mp * a,int b,int e)255*74a4d8c2SCharles.Forsyth mppow(Mp *a, int b, int e)
256*74a4d8c2SCharles.Forsyth {
257*74a4d8c2SCharles.Forsyth int b1;
258*74a4d8c2SCharles.Forsyth
259*74a4d8c2SCharles.Forsyth b1 = b*b;
260*74a4d8c2SCharles.Forsyth b1 = b1*b1;
261*74a4d8c2SCharles.Forsyth while(e >= 4) {
262*74a4d8c2SCharles.Forsyth mpmul(a, b1);
263*74a4d8c2SCharles.Forsyth e -= 4;
264*74a4d8c2SCharles.Forsyth if(a->ovf)
265*74a4d8c2SCharles.Forsyth return;
266*74a4d8c2SCharles.Forsyth }
267*74a4d8c2SCharles.Forsyth while(e > 0) {
268*74a4d8c2SCharles.Forsyth mpmul(a, b);
269*74a4d8c2SCharles.Forsyth e--;
270*74a4d8c2SCharles.Forsyth }
271*74a4d8c2SCharles.Forsyth }
272*74a4d8c2SCharles.Forsyth
273*74a4d8c2SCharles.Forsyth /*
274*74a4d8c2SCharles.Forsyth * convert a string, s, to vlong in *v
275*74a4d8c2SCharles.Forsyth * return conversion overflow.
276*74a4d8c2SCharles.Forsyth * required syntax is [0[x]]d*
277*74a4d8c2SCharles.Forsyth */
278*74a4d8c2SCharles.Forsyth int
mpatov(char * s,vlong * v)279*74a4d8c2SCharles.Forsyth mpatov(char *s, vlong *v)
280*74a4d8c2SCharles.Forsyth {
281*74a4d8c2SCharles.Forsyth vlong n, nn;
282*74a4d8c2SCharles.Forsyth int c;
283*74a4d8c2SCharles.Forsyth n = 0;
284*74a4d8c2SCharles.Forsyth c = *s;
285*74a4d8c2SCharles.Forsyth if(c == '0')
286*74a4d8c2SCharles.Forsyth goto oct;
287*74a4d8c2SCharles.Forsyth while(c = *s++) {
288*74a4d8c2SCharles.Forsyth if(c >= '0' && c <= '9')
289*74a4d8c2SCharles.Forsyth nn = n*10 + c-'0';
290*74a4d8c2SCharles.Forsyth else
291*74a4d8c2SCharles.Forsyth goto bad;
292*74a4d8c2SCharles.Forsyth if(n < 0 && nn >= 0)
293*74a4d8c2SCharles.Forsyth goto bad;
294*74a4d8c2SCharles.Forsyth n = nn;
295*74a4d8c2SCharles.Forsyth }
296*74a4d8c2SCharles.Forsyth goto out;
297*74a4d8c2SCharles.Forsyth oct:
298*74a4d8c2SCharles.Forsyth s++;
299*74a4d8c2SCharles.Forsyth c = *s;
300*74a4d8c2SCharles.Forsyth if(c == 'x' || c == 'X')
301*74a4d8c2SCharles.Forsyth goto hex;
302*74a4d8c2SCharles.Forsyth while(c = *s++) {
303*74a4d8c2SCharles.Forsyth if(c >= '0' || c <= '7')
304*74a4d8c2SCharles.Forsyth nn = n*8 + c-'0';
305*74a4d8c2SCharles.Forsyth else
306*74a4d8c2SCharles.Forsyth goto bad;
307*74a4d8c2SCharles.Forsyth if(n < 0 && nn >= 0)
308*74a4d8c2SCharles.Forsyth goto bad;
309*74a4d8c2SCharles.Forsyth n = nn;
310*74a4d8c2SCharles.Forsyth }
311*74a4d8c2SCharles.Forsyth goto out;
312*74a4d8c2SCharles.Forsyth hex:
313*74a4d8c2SCharles.Forsyth s++;
314*74a4d8c2SCharles.Forsyth while(c = *s++) {
315*74a4d8c2SCharles.Forsyth if(c >= '0' && c <= '9')
316*74a4d8c2SCharles.Forsyth c += 0-'0';
317*74a4d8c2SCharles.Forsyth else
318*74a4d8c2SCharles.Forsyth if(c >= 'a' && c <= 'f')
319*74a4d8c2SCharles.Forsyth c += 10-'a';
320*74a4d8c2SCharles.Forsyth else
321*74a4d8c2SCharles.Forsyth if(c >= 'A' && c <= 'F')
322*74a4d8c2SCharles.Forsyth c += 10-'A';
323*74a4d8c2SCharles.Forsyth else
324*74a4d8c2SCharles.Forsyth goto bad;
325*74a4d8c2SCharles.Forsyth nn = n*16 + c;
326*74a4d8c2SCharles.Forsyth if(n < 0 && nn >= 0)
327*74a4d8c2SCharles.Forsyth goto bad;
328*74a4d8c2SCharles.Forsyth n = nn;
329*74a4d8c2SCharles.Forsyth }
330*74a4d8c2SCharles.Forsyth out:
331*74a4d8c2SCharles.Forsyth *v = n;
332*74a4d8c2SCharles.Forsyth return 0;
333*74a4d8c2SCharles.Forsyth
334*74a4d8c2SCharles.Forsyth bad:
335*74a4d8c2SCharles.Forsyth *v = ~0;
336*74a4d8c2SCharles.Forsyth return 1;
337*74a4d8c2SCharles.Forsyth }
338