xref: /plan9/sys/src/ape/lib/ap/stdio/_fconv.c (revision 3e12c5d1bb89fc02707907988834ef147769ddaf)
1 /* Common routines for _dtoa and strtod */
2 
3 #include "fconv.h"
4 
5 #ifdef DEBUG
6 #include <stdio.h>
7 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
8 #endif
9 
10  double
11 _tens[] = {
12 		1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
13 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
14 		1e20, 1e21, 1e22
15 #ifdef VAX
16 		, 1e23, 1e24
17 #endif
18 		};
19 
20 
21 #ifdef IEEE_Arith
22 double _bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
23 double _tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
24 #else
25 #ifdef IBM
26 double _bigtens[] = { 1e16, 1e32, 1e64 };
27 double _tinytens[] = { 1e-16, 1e-32, 1e-64 };
28 #else
29 double _bigtens[] = { 1e16, 1e32 };
30 double _tinytens[] = { 1e-16, 1e-32 };
31 #endif
32 #endif
33 
34  static Bigint *freelist[Kmax+1];
35 
36  Bigint *
_Balloc(int k)37 _Balloc(int k)
38 {
39 	int x;
40 	Bigint *rv;
41 
42 	if (rv = freelist[k]) {
43 		freelist[k] = rv->next;
44 		}
45 	else {
46 		x = 1 << k;
47 		rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long));
48 		rv->k = k;
49 		rv->maxwds = x;
50 		}
51 	rv->sign = rv->wds = 0;
52 	return rv;
53 	}
54 
55  void
_Bfree(Bigint * v)56 _Bfree(Bigint *v)
57 {
58 	if (v) {
59 		v->next = freelist[v->k];
60 		freelist[v->k] = v;
61 		}
62 	}
63 
64 
65  Bigint *
_multadd(Bigint * b,int m,int a)66 _multadd(Bigint *b, int m, int a)	/* multiply by m and add a */
67 {
68 	int i, wds;
69 	unsigned long *x, y;
70 #ifdef Pack_32
71 	unsigned long xi, z;
72 #endif
73 	Bigint *b1;
74 
75 	wds = b->wds;
76 	x = b->x;
77 	i = 0;
78 	do {
79 #ifdef Pack_32
80 		xi = *x;
81 		y = (xi & 0xffff) * m + a;
82 		z = (xi >> 16) * m + (y >> 16);
83 		a = (int)(z >> 16);
84 		*x++ = (z << 16) + (y & 0xffff);
85 #else
86 		y = *x * m + a;
87 		a = (int)(y >> 16);
88 		*x++ = y & 0xffff;
89 #endif
90 		}
91 		while(++i < wds);
92 	if (a) {
93 		if (wds >= b->maxwds) {
94 			b1 = Balloc(b->k+1);
95 			Bcopy(b1, b);
96 			Bfree(b);
97 			b = b1;
98 			}
99 		b->x[wds++] = a;
100 		b->wds = wds;
101 		}
102 	return b;
103 	}
104 
105  int
_hi0bits(register unsigned long x)106 _hi0bits(register unsigned long x)
107 {
108 	register int k = 0;
109 
110 	if (!(x & 0xffff0000)) {
111 		k = 16;
112 		x <<= 16;
113 		}
114 	if (!(x & 0xff000000)) {
115 		k += 8;
116 		x <<= 8;
117 		}
118 	if (!(x & 0xf0000000)) {
119 		k += 4;
120 		x <<= 4;
121 		}
122 	if (!(x & 0xc0000000)) {
123 		k += 2;
124 		x <<= 2;
125 		}
126 	if (!(x & 0x80000000)) {
127 		k++;
128 		if (!(x & 0x40000000))
129 			return 32;
130 		}
131 	return k;
132 	}
133 
134  static int
lo0bits(unsigned long * y)135 lo0bits(unsigned long *y)
136 {
137 	register int k;
138 	register unsigned long x = *y;
139 
140 	if (x & 7) {
141 		if (x & 1)
142 			return 0;
143 		if (x & 2) {
144 			*y = x >> 1;
145 			return 1;
146 			}
147 		*y = x >> 2;
148 		return 2;
149 		}
150 	k = 0;
151 	if (!(x & 0xffff)) {
152 		k = 16;
153 		x >>= 16;
154 		}
155 	if (!(x & 0xff)) {
156 		k += 8;
157 		x >>= 8;
158 		}
159 	if (!(x & 0xf)) {
160 		k += 4;
161 		x >>= 4;
162 		}
163 	if (!(x & 0x3)) {
164 		k += 2;
165 		x >>= 2;
166 		}
167 	if (!(x & 1)) {
168 		k++;
169 		x >>= 1;
170 		if (!x & 1)
171 			return 32;
172 		}
173 	*y = x;
174 	return k;
175 	}
176 
177  Bigint *
_i2b(int i)178 _i2b(int i)
179 {
180 	Bigint *b;
181 
182 	b = Balloc(1);
183 	b->x[0] = i;
184 	b->wds = 1;
185 	return b;
186 	}
187 
188  Bigint *
_mult(Bigint * a,Bigint * b)189 _mult(Bigint *a, Bigint *b)
190 {
191 	Bigint *c;
192 	int k, wa, wb, wc;
193 	unsigned long carry, y, z;
194 	unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
195 #ifdef Pack_32
196 	unsigned long z2;
197 #endif
198 
199 	if (a->wds < b->wds) {
200 		c = a;
201 		a = b;
202 		b = c;
203 		}
204 	k = a->k;
205 	wa = a->wds;
206 	wb = b->wds;
207 	wc = wa + wb;
208 	if (wc > a->maxwds)
209 		k++;
210 	c = Balloc(k);
211 	for(x = c->x, xa = x + wc; x < xa; x++)
212 		*x = 0;
213 	xa = a->x;
214 	xae = xa + wa;
215 	xb = b->x;
216 	xbe = xb + wb;
217 	xc0 = c->x;
218 #ifdef Pack_32
219 	for(; xb < xbe; xb++, xc0++) {
220 		if (y = *xb & 0xffff) {
221 			x = xa;
222 			xc = xc0;
223 			carry = 0;
224 			do {
225 				z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
226 				carry = z >> 16;
227 				z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
228 				carry = z2 >> 16;
229 				Storeinc(xc, z2, z);
230 				}
231 				while(x < xae);
232 			*xc = carry;
233 			}
234 		if (y = *xb >> 16) {
235 			x = xa;
236 			xc = xc0;
237 			carry = 0;
238 			z2 = *xc;
239 			do {
240 				z = (*x & 0xffff) * y + (*xc >> 16) + carry;
241 				carry = z >> 16;
242 				Storeinc(xc, z, z2);
243 				z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
244 				carry = z2 >> 16;
245 				}
246 				while(x < xae);
247 			*xc = z2;
248 			}
249 		}
250 #else
251 	for(; xb < xbe; xc0++) {
252 		if (y = *xb++) {
253 			x = xa;
254 			xc = xc0;
255 			carry = 0;
256 			do {
257 				z = *x++ * y + *xc + carry;
258 				carry = z >> 16;
259 				*xc++ = z & 0xffff;
260 				}
261 				while(x < xae);
262 			*xc = carry;
263 			}
264 		}
265 #endif
266 	for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
267 	c->wds = wc;
268 	return c;
269 	}
270 
271  static Bigint *p5s;
272 
273  Bigint *
_pow5mult(Bigint * b,int k)274 _pow5mult(Bigint *b, int k)
275 {
276 	Bigint *b1, *p5, *p51;
277 	int i;
278 	static int p05[3] = { 5, 25, 125 };
279 
280 	if (i = k & 3)
281 		b = multadd(b, p05[i-1], 0);
282 
283 	if (!(k >>= 2))
284 		return b;
285 	if (!(p5 = p5s)) {
286 		/* first time */
287 		p5 = p5s = i2b(625);
288 		p5->next = 0;
289 		}
290 	for(;;) {
291 		if (k & 1) {
292 			b1 = mult(b, p5);
293 			Bfree(b);
294 			b = b1;
295 			}
296 		if (!(k >>= 1))
297 			break;
298 		if (!(p51 = p5->next)) {
299 			p51 = p5->next = mult(p5,p5);
300 			p51->next = 0;
301 			}
302 		p5 = p51;
303 		}
304 	return b;
305 	}
306 
307  Bigint *
_lshift(Bigint * b,int k)308 _lshift(Bigint *b, int k)
309 {
310 	int i, k1, n, n1;
311 	Bigint *b1;
312 	unsigned long *x, *x1, *xe, z;
313 
314 #ifdef Pack_32
315 	n = k >> 5;
316 #else
317 	n = k >> 4;
318 #endif
319 	k1 = b->k;
320 	n1 = n + b->wds + 1;
321 	for(i = b->maxwds; n1 > i; i <<= 1)
322 		k1++;
323 	b1 = Balloc(k1);
324 	x1 = b1->x;
325 	for(i = 0; i < n; i++)
326 		*x1++ = 0;
327 	x = b->x;
328 	xe = x + b->wds;
329 #ifdef Pack_32
330 	if (k &= 0x1f) {
331 		k1 = 32 - k;
332 		z = 0;
333 		do {
334 			*x1++ = *x << k | z;
335 			z = *x++ >> k1;
336 			}
337 			while(x < xe);
338 		if (*x1 = z)
339 			++n1;
340 		}
341 #else
342 	if (k &= 0xf) {
343 		k1 = 16 - k;
344 		z = 0;
345 		do {
346 			*x1++ = *x << k  & 0xffff | z;
347 			z = *x++ >> k1;
348 			}
349 			while(x < xe);
350 		if (*x1 = z)
351 			++n1;
352 		}
353 #endif
354 	else do
355 		*x1++ = *x++;
356 		while(x < xe);
357 	b1->wds = n1 - 1;
358 	Bfree(b);
359 	return b1;
360 	}
361 
362  int
_cmp(Bigint * a,Bigint * b)363 _cmp(Bigint *a, Bigint *b)
364 {
365 	unsigned long *xa, *xa0, *xb, *xb0;
366 	int i, j;
367 
368 	i = a->wds;
369 	j = b->wds;
370 #ifdef DEBUG
371 	if (i > 1 && !a->x[i-1])
372 		Bug("cmp called with a->x[a->wds-1] == 0");
373 	if (j > 1 && !b->x[j-1])
374 		Bug("cmp called with b->x[b->wds-1] == 0");
375 #endif
376 	if (i -= j)
377 		return i;
378 	xa0 = a->x;
379 	xa = xa0 + j;
380 	xb0 = b->x;
381 	xb = xb0 + j;
382 	for(;;) {
383 		if (*--xa != *--xb)
384 			return *xa < *xb ? -1 : 1;
385 		if (xa <= xa0)
386 			break;
387 		}
388 	return 0;
389 	}
390 
391  Bigint *
_diff(Bigint * a,Bigint * b)392 _diff(Bigint *a, Bigint *b)
393 {
394 	Bigint *c;
395 	int i, wa, wb;
396 	long borrow, y;	/* We need signed shifts here. */
397 	unsigned long *xa, *xae, *xb, *xbe, *xc;
398 #ifdef Pack_32
399 	long z;
400 #endif
401 
402 	i = cmp(a,b);
403 	if (!i) {
404 		c = Balloc(0);
405 		c->wds = 1;
406 		c->x[0] = 0;
407 		return c;
408 		}
409 	if (i < 0) {
410 		c = a;
411 		a = b;
412 		b = c;
413 		i = 1;
414 		}
415 	else
416 		i = 0;
417 	c = Balloc(a->k);
418 	c->sign = i;
419 	wa = a->wds;
420 	xa = a->x;
421 	xae = xa + wa;
422 	wb = b->wds;
423 	xb = b->x;
424 	xbe = xb + wb;
425 	xc = c->x;
426 	borrow = 0;
427 #ifdef Pack_32
428 	do {
429 		y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
430 		borrow = y >> 16;
431 		Sign_Extend(borrow, y);
432 		z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
433 		borrow = z >> 16;
434 		Sign_Extend(borrow, z);
435 		Storeinc(xc, z, y);
436 		}
437 		while(xb < xbe);
438 	while(xa < xae) {
439 		y = (*xa & 0xffff) + borrow;
440 		borrow = y >> 16;
441 		Sign_Extend(borrow, y);
442 		z = (*xa++ >> 16) + borrow;
443 		borrow = z >> 16;
444 		Sign_Extend(borrow, z);
445 		Storeinc(xc, z, y);
446 		}
447 #else
448 	do {
449 		y = *xa++ - *xb++ + borrow;
450 		borrow = y >> 16;
451 		Sign_Extend(borrow, y);
452 		*xc++ = y & 0xffff;
453 		}
454 		while(xb < xbe);
455 	while(xa < xae) {
456 		y = *xa++ + borrow;
457 		borrow = y >> 16;
458 		Sign_Extend(borrow, y);
459 		*xc++ = y & 0xffff;
460 		}
461 #endif
462 	while(!*--xc)
463 		wa--;
464 	c->wds = wa;
465 	return c;
466 	}
467 
468  Bigint *
_d2b(double darg,int * e,int * bits)469 _d2b(double darg, int *e, int *bits)
470 {
471 	Bigint *b;
472 	int de, i, k;
473 	unsigned long *x, y, z;
474 	Dul d;
475 #ifdef VAX
476 	unsigned long d0, d1;
477 	d.d = darg;
478 	d0 = word0(d) >> 16 | word0(d) << 16;
479 	d1 = word1(d) >> 16 | word1(d) << 16;
480 #else
481 	d.d = darg;
482 #define d0 word0(d)
483 #define d1 word1(d)
484 #endif
485 
486 #ifdef Pack_32
487 	b = Balloc(1);
488 #else
489 	b = Balloc(2);
490 #endif
491 	x = b->x;
492 
493 	z = d0 & Frac_mask;
494 	d0 &= 0x7fffffff;	/* clear sign bit, which we ignore */
495 #ifdef Sudden_Underflow
496 	de = (int)(d0 >> Exp_shift);
497 #ifndef IBM
498 	z |= Exp_msk11;
499 #endif
500 #else
501 	if (de = (int)(d0 >> Exp_shift))
502 		z |= Exp_msk1;
503 #endif
504 #ifdef Pack_32
505 	if (y = d1) {
506 		if (k = lo0bits(&y)) {
507 			x[0] = y | z << 32 - k;
508 			z >>= k;
509 			}
510 		else
511 			x[0] = y;
512 		i = b->wds = (x[1] = z) ? 2 : 1;
513 		}
514 	else {
515 #ifdef DEBUG
516 		if (!z)
517 			Bug("Zero passed to d2b");
518 #endif
519 		k = lo0bits(&z);
520 		x[0] = z;
521 		i = b->wds = 1;
522 		k += 32;
523 		}
524 #else
525 	if (y = d1) {
526 		if (k = lo0bits(&y))
527 			if (k >= 16) {
528 				x[0] = y | z << 32 - k & 0xffff;
529 				x[1] = z >> k - 16 & 0xffff;
530 				x[2] = z >> k;
531 				i = 2;
532 				}
533 			else {
534 				x[0] = y & 0xffff;
535 				x[1] = y >> 16 | z << 16 - k & 0xffff;
536 				x[2] = z >> k & 0xffff;
537 				x[3] = z >> k+16;
538 				i = 3;
539 				}
540 		else {
541 			x[0] = y & 0xffff;
542 			x[1] = y >> 16;
543 			x[2] = z & 0xffff;
544 			x[3] = z >> 16;
545 			i = 3;
546 			}
547 		}
548 	else {
549 #ifdef DEBUG
550 		if (!z)
551 			Bug("Zero passed to d2b");
552 #endif
553 		k = lo0bits(&z);
554 		if (k >= 16) {
555 			x[0] = z;
556 			i = 0;
557 			}
558 		else {
559 			x[0] = z & 0xffff;
560 			x[1] = z >> 16;
561 			i = 1;
562 			}
563 		k += 32;
564 		}
565 	while(!x[i])
566 		--i;
567 	b->wds = i + 1;
568 #endif
569 #ifndef Sudden_Underflow
570 	if (de) {
571 #endif
572 #ifdef IBM
573 		*e = (de - Bias - (P-1) << 2) + k;
574 		*bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
575 #else
576 		*e = de - Bias - (P-1) + k;
577 		*bits = P - k;
578 #endif
579 #ifndef Sudden_Underflow
580 		}
581 	else {
582 		*e = de - Bias - (P-1) + 1 + k;
583 #ifdef Pack_32
584 		*bits = 32*i - hi0bits(x[i-1]);
585 #else
586 		*bits = (i+2)*16 - hi0bits(x[i]);
587 #endif
588 		}
589 #endif
590 	return b;
591 	}
592 #undef d0
593 #undef d1
594