1 /* $NetBSD: strtod.c,v 1.45.2.1 2005/04/19 13:35:54 tron Exp $ */
2
3 /****************************************************************
4 *
5 * The author of this software is David M. Gay.
6 *
7 * Copyright (c) 1991 by AT&T.
8 *
9 * Permission to use, copy, modify, and distribute this software for any
10 * purpose without fee is hereby granted, provided that this entire notice
11 * is included in all copies of any software which is or includes a copy
12 * or modification of this software and in all copies of the supporting
13 * documentation for such software.
14 *
15 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
16 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
17 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
18 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
19 *
20 ***************************************************************/
21
22 /* Please send bug reports to
23 David M. Gay
24 AT&T Bell Laboratories, Room 2C-463
25 600 Mountain Avenue
26 Murray Hill, NJ 07974-2070
27 U.S.A.
28 dmg@research.att.com or research!dmg
29 */
30
31 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
32 *
33 * This strtod returns a nearest machine number to the input decimal
34 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
35 * broken by the IEEE round-even rule. Otherwise ties are broken by
36 * biased rounding (add half and chop).
37 *
38 * Inspired loosely by William D. Clinger's paper "How to Read Floating
39 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
40 *
41 * Modifications:
42 *
43 * 1. We only require IEEE, IBM, or VAX double-precision
44 * arithmetic (not IEEE double-extended).
45 * 2. We get by with floating-point arithmetic in a case that
46 * Clinger missed -- when we're computing d * 10^n
47 * for a small integer d and the integer n is not too
48 * much larger than 22 (the maximum integer k for which
49 * we can represent 10^k exactly), we may be able to
50 * compute (d*10^k) * 10^(e-k) with just one roundoff.
51 * 3. Rather than a bit-at-a-time adjustment of the binary
52 * result in the hard case, we use floating-point
53 * arithmetic to determine the adjustment to within
54 * one bit; only in really hard cases do we need to
55 * compute a second residual.
56 * 4. Because of 3., we don't need a large table of powers of 10
57 * for ten-to-e (just some small tables, e.g. of 10^k
58 * for 0 <= k <= 22).
59 */
60
61 /*
62 * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least
63 * significant byte has the lowest address.
64 * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most
65 * significant byte has the lowest address.
66 * #define Long int on machines with 32-bit ints and 64-bit longs.
67 * #define Sudden_Underflow for IEEE-format machines without gradual
68 * underflow (i.e., that flush to zero on underflow).
69 * #define IBM for IBM mainframe-style floating-point arithmetic.
70 * #define VAX for VAX-style floating-point arithmetic.
71 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
72 * #define No_leftright to omit left-right logic in fast floating-point
73 * computation of dtoa.
74 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
75 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
76 * that use extended-precision instructions to compute rounded
77 * products and quotients) with IBM.
78 * #define ROUND_BIASED for IEEE-format with biased rounding.
79 * #define Inaccurate_Divide for IEEE-format with correctly rounded
80 * products but inaccurate quotients, e.g., for Intel i860.
81 * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision
82 * integer arithmetic. Whether this speeds things up or slows things
83 * down depends on the machine and the number being converted.
84 * #define KR_headers for old-style C function headers.
85 * #define Bad_float_h if your system lacks a float.h or if it does not
86 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
87 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
88 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
89 * if memory is available and otherwise does something you deem
90 * appropriate. If MALLOC is undefined, malloc will be invoked
91 * directly -- and assumed always to succeed.
92 */
93
94 #ifdef ANDROID_CHANGES
95 #include <pthread.h>
96 #define mutex_lock(x) pthread_mutex_lock(x)
97 #define mutex_unlock(x) pthread_mutex_unlock(x)
98 #endif
99
100 #include <sys/cdefs.h>
101 #if defined(LIBC_SCCS) && !defined(lint)
102 __RCSID("$NetBSD: strtod.c,v 1.45.2.1 2005/04/19 13:35:54 tron Exp $");
103 #endif /* LIBC_SCCS and not lint */
104
105 #define Unsigned_Shifts
106 #if defined(__m68k__) || defined(__sparc__) || defined(__i386__) || \
107 defined(__mips__) || defined(__ns32k__) || defined(__alpha__) || \
108 defined(__powerpc__) || defined(__sh__) || defined(__x86_64__) || \
109 defined(__hppa__) || \
110 (defined(__arm__) && defined(__VFP_FP__))
111 #include <endian.h>
112 #if BYTE_ORDER == BIG_ENDIAN
113 #define IEEE_BIG_ENDIAN
114 #else
115 #define IEEE_LITTLE_ENDIAN
116 #endif
117 #endif
118
119 #if defined(__arm__) && !defined(__VFP_FP__)
120 /*
121 * Although the CPU is little endian the FP has different
122 * byte and word endianness. The byte order is still little endian
123 * but the word order is big endian.
124 */
125 #define IEEE_BIG_ENDIAN
126 #endif
127
128 #ifdef __vax__
129 #define VAX
130 #endif
131
132 #if defined(__hppa__) || defined(__mips__) || defined(__sh__)
133 #define NAN_WORD0 0x7ff40000
134 #else
135 #define NAN_WORD0 0x7ff80000
136 #endif
137 #define NAN_WORD1 0
138
139 #define Long int32_t
140 #define ULong u_int32_t
141
142 #ifdef DEBUG
143 #include "stdio.h"
144 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
145 #endif
146
147 #ifdef __cplusplus
148 #include "malloc.h"
149 #include "memory.h"
150 #else
151 #ifndef KR_headers
152 #include "stdlib.h"
153 #include "string.h"
154 #ifndef ANDROID_CHANGES
155 #include "locale.h"
156 #endif /* ANDROID_CHANGES */
157 #else
158 #include "malloc.h"
159 #include "memory.h"
160 #endif
161 #endif
162 #ifndef ANDROID_CHANGES
163 #include "extern.h"
164 #include "reentrant.h"
165 #endif /* ANDROID_CHANGES */
166
167 #ifdef MALLOC
168 #ifdef KR_headers
169 extern char *MALLOC();
170 #else
171 extern void *MALLOC(size_t);
172 #endif
173 #else
174 #define MALLOC malloc
175 #endif
176
177 #include "ctype.h"
178 #include "errno.h"
179 #include "float.h"
180
181 #ifndef __MATH_H__
182 #include "math.h"
183 #endif
184
185 #ifdef __cplusplus
186 extern "C" {
187 #endif
188
189 #ifndef CONST
190 #ifdef KR_headers
191 #define CONST /* blank */
192 #else
193 #define CONST const
194 #endif
195 #endif
196
197 #ifdef Unsigned_Shifts
198 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
199 #else
200 #define Sign_Extend(a,b) /*no-op*/
201 #endif
202
203 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + \
204 defined(IBM) != 1
205 Exactly one of IEEE_LITTLE_ENDIAN IEEE_BIG_ENDIAN, VAX, or
206 IBM should be defined.
207 #endif
208
209 typedef union {
210 double d;
211 ULong ul[2];
212 } _double;
213 #define value(x) ((x).d)
214 #ifdef IEEE_LITTLE_ENDIAN
215 #define word0(x) ((x).ul[1])
216 #define word1(x) ((x).ul[0])
217 #else
218 #define word0(x) ((x).ul[0])
219 #define word1(x) ((x).ul[1])
220 #endif
221
222 /* The following definition of Storeinc is appropriate for MIPS processors.
223 * An alternative that might be better on some machines is
224 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
225 */
226 #if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm__)
227 #define Storeinc(a,b,c) \
228 (((u_short *)(void *)a)[1] = \
229 (u_short)b, ((u_short *)(void *)a)[0] = (u_short)c, a++)
230 #else
231 #define Storeinc(a,b,c) \
232 (((u_short *)(void *)a)[0] = \
233 (u_short)b, ((u_short *)(void *)a)[1] = (u_short)c, a++)
234 #endif
235
236 /* #define P DBL_MANT_DIG */
237 /* Ten_pmax = floor(P*log(2)/log(5)) */
238 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
239 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
240 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
241
242 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN)
243 #define Exp_shift 20
244 #define Exp_shift1 20
245 #define Exp_msk1 0x100000
246 #define Exp_msk11 0x100000
247 #define Exp_mask 0x7ff00000
248 #define P 53
249 #define Bias 1023
250 #define IEEE_Arith
251 #define Emin (-1022)
252 #define Exp_1 0x3ff00000
253 #define Exp_11 0x3ff00000
254 #define Ebits 11
255 #define Frac_mask 0xfffff
256 #define Frac_mask1 0xfffff
257 #define Ten_pmax 22
258 #define Bletch 0x10
259 #define Bndry_mask 0xfffff
260 #define Bndry_mask1 0xfffff
261 #define LSB 1
262 #define Sign_bit 0x80000000
263 #define Log2P 1
264 #define Tiny0 0
265 #define Tiny1 1
266 #define Quick_max 14
267 #define Int_max 14
268 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
269 #else
270 #undef Sudden_Underflow
271 #define Sudden_Underflow
272 #ifdef IBM
273 #define Exp_shift 24
274 #define Exp_shift1 24
275 #define Exp_msk1 0x1000000
276 #define Exp_msk11 0x1000000
277 #define Exp_mask 0x7f000000
278 #define P 14
279 #define Bias 65
280 #define Exp_1 0x41000000
281 #define Exp_11 0x41000000
282 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
283 #define Frac_mask 0xffffff
284 #define Frac_mask1 0xffffff
285 #define Bletch 4
286 #define Ten_pmax 22
287 #define Bndry_mask 0xefffff
288 #define Bndry_mask1 0xffffff
289 #define LSB 1
290 #define Sign_bit 0x80000000
291 #define Log2P 4
292 #define Tiny0 0x100000
293 #define Tiny1 0
294 #define Quick_max 14
295 #define Int_max 15
296 #else /* VAX */
297 #define Exp_shift 23
298 #define Exp_shift1 7
299 #define Exp_msk1 0x80
300 #define Exp_msk11 0x800000
301 #define Exp_mask 0x7f80
302 #define P 56
303 #define Bias 129
304 #define Exp_1 0x40800000
305 #define Exp_11 0x4080
306 #define Ebits 8
307 #define Frac_mask 0x7fffff
308 #define Frac_mask1 0xffff007f
309 #define Ten_pmax 24
310 #define Bletch 2
311 #define Bndry_mask 0xffff007f
312 #define Bndry_mask1 0xffff007f
313 #define LSB 0x10000
314 #define Sign_bit 0x8000
315 #define Log2P 1
316 #define Tiny0 0x80
317 #define Tiny1 0
318 #define Quick_max 15
319 #define Int_max 15
320 #endif
321 #endif
322
323 #ifndef IEEE_Arith
324 #define ROUND_BIASED
325 #endif
326
327 #ifdef RND_PRODQUOT
328 #define rounded_product(a,b) a = rnd_prod(a, b)
329 #define rounded_quotient(a,b) a = rnd_quot(a, b)
330 #ifdef KR_headers
331 extern double rnd_prod(), rnd_quot();
332 #else
333 extern double rnd_prod(double, double), rnd_quot(double, double);
334 #endif
335 #else
336 #define rounded_product(a,b) a *= b
337 #define rounded_quotient(a,b) a /= b
338 #endif
339
340 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
341 #define Big1 0xffffffff
342
343 #ifndef Just_16
344 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
345 * This makes some inner loops simpler and sometimes saves work
346 * during multiplications, but it often seems to make things slightly
347 * slower. Hence the default is now to store 32 bits per Long.
348 */
349 #ifndef Pack_32
350 #define Pack_32
351 #endif
352 #endif
353
354 #define Kmax 15
355
356 #ifdef __cplusplus
357 extern "C" double strtod(const char *s00, char **se);
358 extern "C" char *__dtoa(double d, int mode, int ndigits,
359 int *decpt, int *sign, char **rve);
360 #endif
361
362 struct
363 Bigint {
364 struct Bigint *next;
365 int k, maxwds, sign, wds;
366 ULong x[1];
367 };
368
369 typedef struct Bigint Bigint;
370
371 static Bigint *freelist[Kmax+1];
372
373 #ifdef ANDROID_CHANGES
374 static pthread_mutex_t freelist_mutex = PTHREAD_MUTEX_INITIALIZER;
375 #else
376 #ifdef _REENTRANT
377 static mutex_t freelist_mutex = MUTEX_INITIALIZER;
378 #endif
379 #endif
380
381 /* Special value used to indicate an invalid Bigint value,
382 * e.g. when a memory allocation fails. The idea is that we
383 * want to avoid introducing NULL checks everytime a bigint
384 * computation is performed. Also the NULL value can also be
385 * already used to indicate "value not initialized yet" and
386 * returning NULL might alter the execution code path in
387 * case of OOM.
388 */
389 #define BIGINT_INVALID ((Bigint *)&bigint_invalid_value)
390
391 static const Bigint bigint_invalid_value;
392
393
394 /* Return BIGINT_INVALID on allocation failure.
395 *
396 * Most of the code here depends on the fact that this function
397 * never returns NULL.
398 */
399 static Bigint *
Balloc(k)400 Balloc
401 #ifdef KR_headers
402 (k) int k;
403 #else
404 (int k)
405 #endif
406 {
407 int x;
408 Bigint *rv;
409
410 mutex_lock(&freelist_mutex);
411
412 if ((rv = freelist[k]) != NULL) {
413 freelist[k] = rv->next;
414 }
415 else {
416 x = 1 << k;
417 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(Long));
418 if (rv == NULL) {
419 rv = BIGINT_INVALID;
420 goto EXIT;
421 }
422 rv->k = k;
423 rv->maxwds = x;
424 }
425 rv->sign = rv->wds = 0;
426 EXIT:
427 mutex_unlock(&freelist_mutex);
428
429 return rv;
430 }
431
432 static void
Bfree(v)433 Bfree
434 #ifdef KR_headers
435 (v) Bigint *v;
436 #else
437 (Bigint *v)
438 #endif
439 {
440 if (v && v != BIGINT_INVALID) {
441 mutex_lock(&freelist_mutex);
442
443 v->next = freelist[v->k];
444 freelist[v->k] = v;
445
446 mutex_unlock(&freelist_mutex);
447 }
448 }
449
450 #define Bcopy_valid(x,y) memcpy(&(x)->sign, &(y)->sign, \
451 (y)->wds*sizeof(Long) + 2*sizeof(int))
452
453 #define Bcopy(x,y) Bcopy_ptr(&(x),(y))
454
455 static void
Bcopy_ptr(Bigint ** px,Bigint * y)456 Bcopy_ptr(Bigint **px, Bigint *y)
457 {
458 if (*px == BIGINT_INVALID)
459 return; /* no space to store copy */
460 if (y == BIGINT_INVALID) {
461 Bfree(*px); /* invalid input */
462 *px = BIGINT_INVALID;
463 } else {
464 Bcopy_valid(*px,y);
465 }
466 }
467
468 static Bigint *
multadd(b,m,a)469 multadd
470 #ifdef KR_headers
471 (b, m, a) Bigint *b; int m, a;
472 #else
473 (Bigint *b, int m, int a) /* multiply by m and add a */
474 #endif
475 {
476 int i, wds;
477 ULong *x, y;
478 #ifdef Pack_32
479 ULong xi, z;
480 #endif
481 Bigint *b1;
482
483 if (b == BIGINT_INVALID)
484 return b;
485
486 wds = b->wds;
487 x = b->x;
488 i = 0;
489 do {
490 #ifdef Pack_32
491 xi = *x;
492 y = (xi & 0xffff) * m + a;
493 z = (xi >> 16) * m + (y >> 16);
494 a = (int)(z >> 16);
495 *x++ = (z << 16) + (y & 0xffff);
496 #else
497 y = *x * m + a;
498 a = (int)(y >> 16);
499 *x++ = y & 0xffff;
500 #endif
501 }
502 while(++i < wds);
503 if (a) {
504 if (wds >= b->maxwds) {
505 b1 = Balloc(b->k+1);
506 if (b1 == BIGINT_INVALID) {
507 Bfree(b);
508 return b1;
509 }
510 Bcopy_valid(b1, b);
511 Bfree(b);
512 b = b1;
513 }
514 b->x[wds++] = a;
515 b->wds = wds;
516 }
517 return b;
518 }
519
520 static Bigint *
s2b(s,nd0,nd,y9)521 s2b
522 #ifdef KR_headers
523 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
524 #else
525 (CONST char *s, int nd0, int nd, ULong y9)
526 #endif
527 {
528 Bigint *b;
529 int i, k;
530 Long x, y;
531
532 x = (nd + 8) / 9;
533 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
534 #ifdef Pack_32
535 b = Balloc(k);
536 if (b == BIGINT_INVALID)
537 return b;
538 b->x[0] = y9;
539 b->wds = 1;
540 #else
541 b = Balloc(k+1);
542 if (b == BIGINT_INVALID)
543 return b;
544
545 b->x[0] = y9 & 0xffff;
546 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
547 #endif
548
549 i = 9;
550 if (9 < nd0) {
551 s += 9;
552 do b = multadd(b, 10, *s++ - '0');
553 while(++i < nd0);
554 s++;
555 }
556 else
557 s += 10;
558 for(; i < nd; i++)
559 b = multadd(b, 10, *s++ - '0');
560 return b;
561 }
562
563 static int
hi0bits(x)564 hi0bits
565 #ifdef KR_headers
566 (x) ULong x;
567 #else
568 (ULong x)
569 #endif
570 {
571 int k = 0;
572
573 if (!(x & 0xffff0000)) {
574 k = 16;
575 x <<= 16;
576 }
577 if (!(x & 0xff000000)) {
578 k += 8;
579 x <<= 8;
580 }
581 if (!(x & 0xf0000000)) {
582 k += 4;
583 x <<= 4;
584 }
585 if (!(x & 0xc0000000)) {
586 k += 2;
587 x <<= 2;
588 }
589 if (!(x & 0x80000000)) {
590 k++;
591 if (!(x & 0x40000000))
592 return 32;
593 }
594 return k;
595 }
596
597 static int
lo0bits(y)598 lo0bits
599 #ifdef KR_headers
600 (y) ULong *y;
601 #else
602 (ULong *y)
603 #endif
604 {
605 int k;
606 ULong x = *y;
607
608 if (x & 7) {
609 if (x & 1)
610 return 0;
611 if (x & 2) {
612 *y = x >> 1;
613 return 1;
614 }
615 *y = x >> 2;
616 return 2;
617 }
618 k = 0;
619 if (!(x & 0xffff)) {
620 k = 16;
621 x >>= 16;
622 }
623 if (!(x & 0xff)) {
624 k += 8;
625 x >>= 8;
626 }
627 if (!(x & 0xf)) {
628 k += 4;
629 x >>= 4;
630 }
631 if (!(x & 0x3)) {
632 k += 2;
633 x >>= 2;
634 }
635 if (!(x & 1)) {
636 k++;
637 x >>= 1;
638 if (!x & 1)
639 return 32;
640 }
641 *y = x;
642 return k;
643 }
644
645 static Bigint *
i2b(i)646 i2b
647 #ifdef KR_headers
648 (i) int i;
649 #else
650 (int i)
651 #endif
652 {
653 Bigint *b;
654
655 b = Balloc(1);
656 if (b != BIGINT_INVALID) {
657 b->x[0] = i;
658 b->wds = 1;
659 }
660 return b;
661 }
662
663 static Bigint *
mult(a,b)664 mult
665 #ifdef KR_headers
666 (a, b) Bigint *a, *b;
667 #else
668 (Bigint *a, Bigint *b)
669 #endif
670 {
671 Bigint *c;
672 int k, wa, wb, wc;
673 ULong carry, y, z;
674 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
675 #ifdef Pack_32
676 ULong z2;
677 #endif
678
679 if (a == BIGINT_INVALID || b == BIGINT_INVALID)
680 return BIGINT_INVALID;
681
682 if (a->wds < b->wds) {
683 c = a;
684 a = b;
685 b = c;
686 }
687 k = a->k;
688 wa = a->wds;
689 wb = b->wds;
690 wc = wa + wb;
691 if (wc > a->maxwds)
692 k++;
693 c = Balloc(k);
694 if (c == BIGINT_INVALID)
695 return c;
696 for(x = c->x, xa = x + wc; x < xa; x++)
697 *x = 0;
698 xa = a->x;
699 xae = xa + wa;
700 xb = b->x;
701 xbe = xb + wb;
702 xc0 = c->x;
703 #ifdef Pack_32
704 for(; xb < xbe; xb++, xc0++) {
705 if ((y = *xb & 0xffff) != 0) {
706 x = xa;
707 xc = xc0;
708 carry = 0;
709 do {
710 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
711 carry = z >> 16;
712 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
713 carry = z2 >> 16;
714 Storeinc(xc, z2, z);
715 }
716 while(x < xae);
717 *xc = carry;
718 }
719 if ((y = *xb >> 16) != 0) {
720 x = xa;
721 xc = xc0;
722 carry = 0;
723 z2 = *xc;
724 do {
725 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
726 carry = z >> 16;
727 Storeinc(xc, z, z2);
728 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
729 carry = z2 >> 16;
730 }
731 while(x < xae);
732 *xc = z2;
733 }
734 }
735 #else
736 for(; xb < xbe; xc0++) {
737 if (y = *xb++) {
738 x = xa;
739 xc = xc0;
740 carry = 0;
741 do {
742 z = *x++ * y + *xc + carry;
743 carry = z >> 16;
744 *xc++ = z & 0xffff;
745 }
746 while(x < xae);
747 *xc = carry;
748 }
749 }
750 #endif
751 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
752 c->wds = wc;
753 return c;
754 }
755
756 static Bigint *p5s;
757 static pthread_mutex_t p5s_mutex = PTHREAD_MUTEX_INITIALIZER;
758
759 static Bigint *
pow5mult(b,k)760 pow5mult
761 #ifdef KR_headers
762 (b, k) Bigint *b; int k;
763 #else
764 (Bigint *b, int k)
765 #endif
766 {
767 Bigint *b1, *p5, *p51;
768 int i;
769 static const int p05[3] = { 5, 25, 125 };
770
771 if (b == BIGINT_INVALID)
772 return b;
773
774 if ((i = k & 3) != 0)
775 b = multadd(b, p05[i-1], 0);
776
777 if (!(k = (unsigned int) k >> 2))
778 return b;
779 mutex_lock(&p5s_mutex);
780 if (!(p5 = p5s)) {
781 /* first time */
782 p5 = i2b(625);
783 if (p5 == BIGINT_INVALID) {
784 Bfree(b);
785 mutex_unlock(&p5s_mutex);
786 return p5;
787 }
788 p5s = p5;
789 p5->next = 0;
790 }
791 for(;;) {
792 if (k & 1) {
793 b1 = mult(b, p5);
794 Bfree(b);
795 b = b1;
796 }
797 if (!(k = (unsigned int) k >> 1))
798 break;
799 if (!(p51 = p5->next)) {
800 p51 = mult(p5,p5);
801 if (p51 == BIGINT_INVALID) {
802 Bfree(b);
803 mutex_unlock(&p5s_mutex);
804 return p51;
805 }
806 p5->next = p51;
807 p51->next = 0;
808 }
809 p5 = p51;
810 }
811 mutex_unlock(&p5s_mutex);
812 return b;
813 }
814
815 static Bigint *
lshift(b,k)816 lshift
817 #ifdef KR_headers
818 (b, k) Bigint *b; int k;
819 #else
820 (Bigint *b, int k)
821 #endif
822 {
823 int i, k1, n, n1;
824 Bigint *b1;
825 ULong *x, *x1, *xe, z;
826
827 if (b == BIGINT_INVALID)
828 return b;
829
830 #ifdef Pack_32
831 n = (unsigned int)k >> 5;
832 #else
833 n = (unsigned int)k >> 4;
834 #endif
835 k1 = b->k;
836 n1 = n + b->wds + 1;
837 for(i = b->maxwds; n1 > i; i <<= 1)
838 k1++;
839 b1 = Balloc(k1);
840 if (b1 == BIGINT_INVALID) {
841 Bfree(b);
842 return b1;
843 }
844 x1 = b1->x;
845 for(i = 0; i < n; i++)
846 *x1++ = 0;
847 x = b->x;
848 xe = x + b->wds;
849 #ifdef Pack_32
850 if (k &= 0x1f) {
851 k1 = 32 - k;
852 z = 0;
853 do {
854 *x1++ = *x << k | z;
855 z = *x++ >> k1;
856 }
857 while(x < xe);
858 if ((*x1 = z) != 0)
859 ++n1;
860 }
861 #else
862 if (k &= 0xf) {
863 k1 = 16 - k;
864 z = 0;
865 do {
866 *x1++ = *x << k & 0xffff | z;
867 z = *x++ >> k1;
868 }
869 while(x < xe);
870 if (*x1 = z)
871 ++n1;
872 }
873 #endif
874 else do
875 *x1++ = *x++;
876 while(x < xe);
877 b1->wds = n1 - 1;
878 Bfree(b);
879 return b1;
880 }
881
882 static int
cmp(a,b)883 cmp
884 #ifdef KR_headers
885 (a, b) Bigint *a, *b;
886 #else
887 (Bigint *a, Bigint *b)
888 #endif
889 {
890 ULong *xa, *xa0, *xb, *xb0;
891 int i, j;
892
893 if (a == BIGINT_INVALID || b == BIGINT_INVALID)
894 #ifdef DEBUG
895 Bug("cmp called with a or b invalid");
896 #else
897 return 0; /* equal - the best we can do right now */
898 #endif
899
900 i = a->wds;
901 j = b->wds;
902 #ifdef DEBUG
903 if (i > 1 && !a->x[i-1])
904 Bug("cmp called with a->x[a->wds-1] == 0");
905 if (j > 1 && !b->x[j-1])
906 Bug("cmp called with b->x[b->wds-1] == 0");
907 #endif
908 if (i -= j)
909 return i;
910 xa0 = a->x;
911 xa = xa0 + j;
912 xb0 = b->x;
913 xb = xb0 + j;
914 for(;;) {
915 if (*--xa != *--xb)
916 return *xa < *xb ? -1 : 1;
917 if (xa <= xa0)
918 break;
919 }
920 return 0;
921 }
922
923 static Bigint *
diff(a,b)924 diff
925 #ifdef KR_headers
926 (a, b) Bigint *a, *b;
927 #else
928 (Bigint *a, Bigint *b)
929 #endif
930 {
931 Bigint *c;
932 int i, wa, wb;
933 Long borrow, y; /* We need signed shifts here. */
934 ULong *xa, *xae, *xb, *xbe, *xc;
935 #ifdef Pack_32
936 Long z;
937 #endif
938
939 if (a == BIGINT_INVALID || b == BIGINT_INVALID)
940 return BIGINT_INVALID;
941
942 i = cmp(a,b);
943 if (!i) {
944 c = Balloc(0);
945 if (c != BIGINT_INVALID) {
946 c->wds = 1;
947 c->x[0] = 0;
948 }
949 return c;
950 }
951 if (i < 0) {
952 c = a;
953 a = b;
954 b = c;
955 i = 1;
956 }
957 else
958 i = 0;
959 c = Balloc(a->k);
960 if (c == BIGINT_INVALID)
961 return c;
962 c->sign = i;
963 wa = a->wds;
964 xa = a->x;
965 xae = xa + wa;
966 wb = b->wds;
967 xb = b->x;
968 xbe = xb + wb;
969 xc = c->x;
970 borrow = 0;
971 #ifdef Pack_32
972 do {
973 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
974 borrow = (ULong)y >> 16;
975 Sign_Extend(borrow, y);
976 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
977 borrow = (ULong)z >> 16;
978 Sign_Extend(borrow, z);
979 Storeinc(xc, z, y);
980 }
981 while(xb < xbe);
982 while(xa < xae) {
983 y = (*xa & 0xffff) + borrow;
984 borrow = (ULong)y >> 16;
985 Sign_Extend(borrow, y);
986 z = (*xa++ >> 16) + borrow;
987 borrow = (ULong)z >> 16;
988 Sign_Extend(borrow, z);
989 Storeinc(xc, z, y);
990 }
991 #else
992 do {
993 y = *xa++ - *xb++ + borrow;
994 borrow = y >> 16;
995 Sign_Extend(borrow, y);
996 *xc++ = y & 0xffff;
997 }
998 while(xb < xbe);
999 while(xa < xae) {
1000 y = *xa++ + borrow;
1001 borrow = y >> 16;
1002 Sign_Extend(borrow, y);
1003 *xc++ = y & 0xffff;
1004 }
1005 #endif
1006 while(!*--xc)
1007 wa--;
1008 c->wds = wa;
1009 return c;
1010 }
1011
1012 static double
ulp(_x)1013 ulp
1014 #ifdef KR_headers
1015 (_x) double _x;
1016 #else
1017 (double _x)
1018 #endif
1019 {
1020 _double x;
1021 Long L;
1022 _double a;
1023
1024 value(x) = _x;
1025 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1026 #ifndef Sudden_Underflow
1027 if (L > 0) {
1028 #endif
1029 #ifdef IBM
1030 L |= Exp_msk1 >> 4;
1031 #endif
1032 word0(a) = L;
1033 word1(a) = 0;
1034 #ifndef Sudden_Underflow
1035 }
1036 else {
1037 L = (ULong)-L >> Exp_shift;
1038 if (L < Exp_shift) {
1039 word0(a) = 0x80000 >> L;
1040 word1(a) = 0;
1041 }
1042 else {
1043 word0(a) = 0;
1044 L -= Exp_shift;
1045 word1(a) = L >= 31 ? 1 : 1 << (31 - L);
1046 }
1047 }
1048 #endif
1049 return value(a);
1050 }
1051
1052 static double
b2d(a,e)1053 b2d
1054 #ifdef KR_headers
1055 (a, e) Bigint *a; int *e;
1056 #else
1057 (Bigint *a, int *e)
1058 #endif
1059 {
1060 ULong *xa, *xa0, w, y, z;
1061 int k;
1062 _double d;
1063 #ifdef VAX
1064 ULong d0, d1;
1065 #else
1066 #define d0 word0(d)
1067 #define d1 word1(d)
1068 #endif
1069
1070 if (a == BIGINT_INVALID)
1071 return NAN;
1072
1073 xa0 = a->x;
1074 xa = xa0 + a->wds;
1075 y = *--xa;
1076 #ifdef DEBUG
1077 if (!y) Bug("zero y in b2d");
1078 #endif
1079 k = hi0bits(y);
1080 *e = 32 - k;
1081 #ifdef Pack_32
1082 if (k < Ebits) {
1083 d0 = Exp_1 | y >> (Ebits - k);
1084 w = xa > xa0 ? *--xa : 0;
1085 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1086 goto ret_d;
1087 }
1088 z = xa > xa0 ? *--xa : 0;
1089 if (k -= Ebits) {
1090 d0 = Exp_1 | y << k | z >> (32 - k);
1091 y = xa > xa0 ? *--xa : 0;
1092 d1 = z << k | y >> (32 - k);
1093 }
1094 else {
1095 d0 = Exp_1 | y;
1096 d1 = z;
1097 }
1098 #else
1099 if (k < Ebits + 16) {
1100 z = xa > xa0 ? *--xa : 0;
1101 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1102 w = xa > xa0 ? *--xa : 0;
1103 y = xa > xa0 ? *--xa : 0;
1104 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1105 goto ret_d;
1106 }
1107 z = xa > xa0 ? *--xa : 0;
1108 w = xa > xa0 ? *--xa : 0;
1109 k -= Ebits + 16;
1110 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1111 y = xa > xa0 ? *--xa : 0;
1112 d1 = w << k + 16 | y << k;
1113 #endif
1114 ret_d:
1115 #ifdef VAX
1116 word0(d) = d0 >> 16 | d0 << 16;
1117 word1(d) = d1 >> 16 | d1 << 16;
1118 #else
1119 #undef d0
1120 #undef d1
1121 #endif
1122 return value(d);
1123 }
1124
1125 static Bigint *
d2b(_d,e,bits)1126 d2b
1127 #ifdef KR_headers
1128 (_d, e, bits) double d; int *e, *bits;
1129 #else
1130 (double _d, int *e, int *bits)
1131 #endif
1132 {
1133 Bigint *b;
1134 int de, i, k;
1135 ULong *x, y, z;
1136 _double d;
1137 #ifdef VAX
1138 ULong d0, d1;
1139 #endif
1140
1141 value(d) = _d;
1142 #ifdef VAX
1143 d0 = word0(d) >> 16 | word0(d) << 16;
1144 d1 = word1(d) >> 16 | word1(d) << 16;
1145 #else
1146 #define d0 word0(d)
1147 #define d1 word1(d)
1148 #endif
1149
1150 #ifdef Pack_32
1151 b = Balloc(1);
1152 #else
1153 b = Balloc(2);
1154 #endif
1155 if (b == BIGINT_INVALID)
1156 return b;
1157 x = b->x;
1158
1159 z = d0 & Frac_mask;
1160 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1161 #ifdef Sudden_Underflow
1162 de = (int)(d0 >> Exp_shift);
1163 #ifndef IBM
1164 z |= Exp_msk11;
1165 #endif
1166 #else
1167 if ((de = (int)(d0 >> Exp_shift)) != 0)
1168 z |= Exp_msk1;
1169 #endif
1170 #ifdef Pack_32
1171 if ((y = d1) != 0) {
1172 if ((k = lo0bits(&y)) != 0) {
1173 x[0] = y | z << (32 - k);
1174 z >>= k;
1175 }
1176 else
1177 x[0] = y;
1178 i = b->wds = (x[1] = z) ? 2 : 1;
1179 }
1180 else {
1181 #ifdef DEBUG
1182 if (!z)
1183 Bug("Zero passed to d2b");
1184 #endif
1185 k = lo0bits(&z);
1186 x[0] = z;
1187 i = b->wds = 1;
1188 k += 32;
1189 }
1190 #else
1191 if (y = d1) {
1192 if (k = lo0bits(&y))
1193 if (k >= 16) {
1194 x[0] = y | z << 32 - k & 0xffff;
1195 x[1] = z >> k - 16 & 0xffff;
1196 x[2] = z >> k;
1197 i = 2;
1198 }
1199 else {
1200 x[0] = y & 0xffff;
1201 x[1] = y >> 16 | z << 16 - k & 0xffff;
1202 x[2] = z >> k & 0xffff;
1203 x[3] = z >> k+16;
1204 i = 3;
1205 }
1206 else {
1207 x[0] = y & 0xffff;
1208 x[1] = y >> 16;
1209 x[2] = z & 0xffff;
1210 x[3] = z >> 16;
1211 i = 3;
1212 }
1213 }
1214 else {
1215 #ifdef DEBUG
1216 if (!z)
1217 Bug("Zero passed to d2b");
1218 #endif
1219 k = lo0bits(&z);
1220 if (k >= 16) {
1221 x[0] = z;
1222 i = 0;
1223 }
1224 else {
1225 x[0] = z & 0xffff;
1226 x[1] = z >> 16;
1227 i = 1;
1228 }
1229 k += 32;
1230 }
1231 while(!x[i])
1232 --i;
1233 b->wds = i + 1;
1234 #endif
1235 #ifndef Sudden_Underflow
1236 if (de) {
1237 #endif
1238 #ifdef IBM
1239 *e = (de - Bias - (P-1) << 2) + k;
1240 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1241 #else
1242 *e = de - Bias - (P-1) + k;
1243 *bits = P - k;
1244 #endif
1245 #ifndef Sudden_Underflow
1246 }
1247 else {
1248 *e = de - Bias - (P-1) + 1 + k;
1249 #ifdef Pack_32
1250 *bits = 32*i - hi0bits(x[i-1]);
1251 #else
1252 *bits = (i+2)*16 - hi0bits(x[i]);
1253 #endif
1254 }
1255 #endif
1256 return b;
1257 }
1258 #undef d0
1259 #undef d1
1260
1261 static double
ratio(a,b)1262 ratio
1263 #ifdef KR_headers
1264 (a, b) Bigint *a, *b;
1265 #else
1266 (Bigint *a, Bigint *b)
1267 #endif
1268 {
1269 _double da, db;
1270 int k, ka, kb;
1271
1272 if (a == BIGINT_INVALID || b == BIGINT_INVALID)
1273 return NAN; /* for lack of better value ? */
1274
1275 value(da) = b2d(a, &ka);
1276 value(db) = b2d(b, &kb);
1277 #ifdef Pack_32
1278 k = ka - kb + 32*(a->wds - b->wds);
1279 #else
1280 k = ka - kb + 16*(a->wds - b->wds);
1281 #endif
1282 #ifdef IBM
1283 if (k > 0) {
1284 word0(da) += (k >> 2)*Exp_msk1;
1285 if (k &= 3)
1286 da *= 1 << k;
1287 }
1288 else {
1289 k = -k;
1290 word0(db) += (k >> 2)*Exp_msk1;
1291 if (k &= 3)
1292 db *= 1 << k;
1293 }
1294 #else
1295 if (k > 0)
1296 word0(da) += k*Exp_msk1;
1297 else {
1298 k = -k;
1299 word0(db) += k*Exp_msk1;
1300 }
1301 #endif
1302 return value(da) / value(db);
1303 }
1304
1305 static CONST double
1306 tens[] = {
1307 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1308 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1309 1e20, 1e21, 1e22
1310 #ifdef VAX
1311 , 1e23, 1e24
1312 #endif
1313 };
1314
1315 #ifdef IEEE_Arith
1316 static CONST double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1317 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1318 #define n_bigtens 5
1319 #else
1320 #ifdef IBM
1321 static CONST double bigtens[] = { 1e16, 1e32, 1e64 };
1322 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1323 #define n_bigtens 3
1324 #else
1325 static CONST double bigtens[] = { 1e16, 1e32 };
1326 static CONST double tinytens[] = { 1e-16, 1e-32 };
1327 #define n_bigtens 2
1328 #endif
1329 #endif
1330
1331 double
strtod(s00,se)1332 strtod
1333 #ifdef KR_headers
1334 (s00, se) CONST char *s00; char **se;
1335 #else
1336 (CONST char *s00, char **se)
1337 #endif
1338 {
1339 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1340 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1341 CONST char *s, *s0, *s1;
1342 double aadj, aadj1, adj;
1343 _double rv, rv0;
1344 Long L;
1345 ULong y, z;
1346 Bigint *bb1, *bd0;
1347 Bigint *bb = NULL, *bd = NULL, *bs = NULL, *delta = NULL;/* pacify gcc */
1348
1349 #ifdef ANDROID_CHANGES
1350 CONST char decimal_point = '.';
1351 #else /* ANDROID_CHANGES */
1352 #ifndef KR_headers
1353 CONST char decimal_point = localeconv()->decimal_point[0];
1354 #else
1355 CONST char decimal_point = '.';
1356 #endif
1357
1358 #endif /* ANDROID_CHANGES */
1359
1360 sign = nz0 = nz = 0;
1361 value(rv) = 0.;
1362
1363
1364 for(s = s00; isspace((unsigned char) *s); s++)
1365 ;
1366
1367 if (*s == '-') {
1368 sign = 1;
1369 s++;
1370 } else if (*s == '+') {
1371 s++;
1372 }
1373
1374 if (*s == '\0') {
1375 s = s00;
1376 goto ret;
1377 }
1378
1379 /* "INF" or "INFINITY" */
1380 if (tolower((unsigned char)*s) == 'i' && strncasecmp(s, "inf", 3) == 0) {
1381 if (strncasecmp(s + 3, "inity", 5) == 0)
1382 s += 8;
1383 else
1384 s += 3;
1385
1386 value(rv) = HUGE_VAL;
1387 goto ret;
1388 }
1389
1390 #ifdef IEEE_Arith
1391 /* "NAN" or "NAN(n-char-sequence-opt)" */
1392 if (tolower((unsigned char)*s) == 'n' && strncasecmp(s, "nan", 3) == 0) {
1393 /* Build a quiet NaN. */
1394 word0(rv) = NAN_WORD0;
1395 word1(rv) = NAN_WORD1;
1396 s+= 3;
1397
1398 /* Don't interpret (n-char-sequence-opt), for now. */
1399 if (*s == '(') {
1400 s0 = s;
1401 for (s++; *s != ')' && *s != '\0'; s++)
1402 ;
1403 if (*s == ')')
1404 s++; /* Skip over closing paren ... */
1405 else
1406 s = s0; /* ... otherwise go back. */
1407 }
1408
1409 goto ret;
1410 }
1411 #endif
1412
1413 if (*s == '0') {
1414 nz0 = 1;
1415 while(*++s == '0') ;
1416 if (!*s)
1417 goto ret;
1418 }
1419 s0 = s;
1420 y = z = 0;
1421 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1422 if (nd < 9)
1423 y = 10*y + c - '0';
1424 else if (nd < 16)
1425 z = 10*z + c - '0';
1426 nd0 = nd;
1427 if (c == decimal_point) {
1428 c = *++s;
1429 if (!nd) {
1430 for(; c == '0'; c = *++s)
1431 nz++;
1432 if (c > '0' && c <= '9') {
1433 s0 = s;
1434 nf += nz;
1435 nz = 0;
1436 goto have_dig;
1437 }
1438 goto dig_done;
1439 }
1440 for(; c >= '0' && c <= '9'; c = *++s) {
1441 have_dig:
1442 nz++;
1443 if (c -= '0') {
1444 nf += nz;
1445 for(i = 1; i < nz; i++)
1446 if (nd++ < 9)
1447 y *= 10;
1448 else if (nd <= DBL_DIG + 1)
1449 z *= 10;
1450 if (nd++ < 9)
1451 y = 10*y + c;
1452 else if (nd <= DBL_DIG + 1)
1453 z = 10*z + c;
1454 nz = 0;
1455 }
1456 }
1457 }
1458 dig_done:
1459 e = 0;
1460 if (c == 'e' || c == 'E') {
1461 if (!nd && !nz && !nz0) {
1462 s = s00;
1463 goto ret;
1464 }
1465 s00 = s;
1466 esign = 0;
1467 switch(c = *++s) {
1468 case '-':
1469 esign = 1;
1470 /* FALLTHROUGH */
1471 case '+':
1472 c = *++s;
1473 }
1474 if (c >= '0' && c <= '9') {
1475 while(c == '0')
1476 c = *++s;
1477 if (c > '0' && c <= '9') {
1478 L = c - '0';
1479 s1 = s;
1480 while((c = *++s) >= '0' && c <= '9')
1481 L = 10*L + c - '0';
1482 if (s - s1 > 8 || L > 19999)
1483 /* Avoid confusion from exponents
1484 * so large that e might overflow.
1485 */
1486 e = 19999; /* safe for 16 bit ints */
1487 else
1488 e = (int)L;
1489 if (esign)
1490 e = -e;
1491 }
1492 else
1493 e = 0;
1494 }
1495 else
1496 s = s00;
1497 }
1498 if (!nd) {
1499 if (!nz && !nz0)
1500 s = s00;
1501 goto ret;
1502 }
1503 e1 = e -= nf;
1504
1505 /* Now we have nd0 digits, starting at s0, followed by a
1506 * decimal point, followed by nd-nd0 digits. The number we're
1507 * after is the integer represented by those digits times
1508 * 10**e */
1509
1510 if (!nd0)
1511 nd0 = nd;
1512 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1513 value(rv) = y;
1514 if (k > 9)
1515 value(rv) = tens[k - 9] * value(rv) + z;
1516 bd0 = 0;
1517 if (nd <= DBL_DIG
1518 #ifndef RND_PRODQUOT
1519 && FLT_ROUNDS == 1
1520 #endif
1521 ) {
1522 if (!e)
1523 goto ret;
1524 if (e > 0) {
1525 if (e <= Ten_pmax) {
1526 #ifdef VAX
1527 goto vax_ovfl_check;
1528 #else
1529 /* value(rv) = */ rounded_product(value(rv),
1530 tens[e]);
1531 goto ret;
1532 #endif
1533 }
1534 i = DBL_DIG - nd;
1535 if (e <= Ten_pmax + i) {
1536 /* A fancier test would sometimes let us do
1537 * this for larger i values.
1538 */
1539 e -= i;
1540 value(rv) *= tens[i];
1541 #ifdef VAX
1542 /* VAX exponent range is so narrow we must
1543 * worry about overflow here...
1544 */
1545 vax_ovfl_check:
1546 word0(rv) -= P*Exp_msk1;
1547 /* value(rv) = */ rounded_product(value(rv),
1548 tens[e]);
1549 if ((word0(rv) & Exp_mask)
1550 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1551 goto ovfl;
1552 word0(rv) += P*Exp_msk1;
1553 #else
1554 /* value(rv) = */ rounded_product(value(rv),
1555 tens[e]);
1556 #endif
1557 goto ret;
1558 }
1559 }
1560 #ifndef Inaccurate_Divide
1561 else if (e >= -Ten_pmax) {
1562 /* value(rv) = */ rounded_quotient(value(rv),
1563 tens[-e]);
1564 goto ret;
1565 }
1566 #endif
1567 }
1568 e1 += nd - k;
1569
1570 /* Get starting approximation = rv * 10**e1 */
1571
1572 if (e1 > 0) {
1573 if ((i = e1 & 15) != 0)
1574 value(rv) *= tens[i];
1575 if (e1 &= ~15) {
1576 if (e1 > DBL_MAX_10_EXP) {
1577 ovfl:
1578 errno = ERANGE;
1579 value(rv) = HUGE_VAL;
1580 if (bd0)
1581 goto retfree;
1582 goto ret;
1583 }
1584 if ((e1 = (unsigned int)e1 >> 4) != 0) {
1585 for(j = 0; e1 > 1; j++,
1586 e1 = (unsigned int)e1 >> 1)
1587 if (e1 & 1)
1588 value(rv) *= bigtens[j];
1589 /* The last multiplication could overflow. */
1590 word0(rv) -= P*Exp_msk1;
1591 value(rv) *= bigtens[j];
1592 if ((z = word0(rv) & Exp_mask)
1593 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1594 goto ovfl;
1595 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1596 /* set to largest number */
1597 /* (Can't trust DBL_MAX) */
1598 word0(rv) = Big0;
1599 word1(rv) = Big1;
1600 }
1601 else
1602 word0(rv) += P*Exp_msk1;
1603 }
1604 }
1605 }
1606 else if (e1 < 0) {
1607 e1 = -e1;
1608 if ((i = e1 & 15) != 0)
1609 value(rv) /= tens[i];
1610 if (e1 &= ~15) {
1611 e1 = (unsigned int)e1 >> 4;
1612 if (e1 >= 1 << n_bigtens)
1613 goto undfl;
1614 for(j = 0; e1 > 1; j++,
1615 e1 = (unsigned int)e1 >> 1)
1616 if (e1 & 1)
1617 value(rv) *= tinytens[j];
1618 /* The last multiplication could underflow. */
1619 value(rv0) = value(rv);
1620 value(rv) *= tinytens[j];
1621 if (!value(rv)) {
1622 value(rv) = 2.*value(rv0);
1623 value(rv) *= tinytens[j];
1624 if (!value(rv)) {
1625 undfl:
1626 value(rv) = 0.;
1627 errno = ERANGE;
1628 if (bd0)
1629 goto retfree;
1630 goto ret;
1631 }
1632 word0(rv) = Tiny0;
1633 word1(rv) = Tiny1;
1634 /* The refinement below will clean
1635 * this approximation up.
1636 */
1637 }
1638 }
1639 }
1640
1641 /* Now the hard part -- adjusting rv to the correct value.*/
1642
1643 /* Put digits into bd: true value = bd * 10^e */
1644
1645 bd0 = s2b(s0, nd0, nd, y);
1646
1647 for(;;) {
1648 bd = Balloc(bd0->k);
1649 Bcopy(bd, bd0);
1650 bb = d2b(value(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
1651 bs = i2b(1);
1652
1653 if (e >= 0) {
1654 bb2 = bb5 = 0;
1655 bd2 = bd5 = e;
1656 }
1657 else {
1658 bb2 = bb5 = -e;
1659 bd2 = bd5 = 0;
1660 }
1661 if (bbe >= 0)
1662 bb2 += bbe;
1663 else
1664 bd2 -= bbe;
1665 bs2 = bb2;
1666 #ifdef Sudden_Underflow
1667 #ifdef IBM
1668 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1669 #else
1670 j = P + 1 - bbbits;
1671 #endif
1672 #else
1673 i = bbe + bbbits - 1; /* logb(rv) */
1674 if (i < Emin) /* denormal */
1675 j = bbe + (P-Emin);
1676 else
1677 j = P + 1 - bbbits;
1678 #endif
1679 bb2 += j;
1680 bd2 += j;
1681 i = bb2 < bd2 ? bb2 : bd2;
1682 if (i > bs2)
1683 i = bs2;
1684 if (i > 0) {
1685 bb2 -= i;
1686 bd2 -= i;
1687 bs2 -= i;
1688 }
1689 if (bb5 > 0) {
1690 bs = pow5mult(bs, bb5);
1691 bb1 = mult(bs, bb);
1692 Bfree(bb);
1693 bb = bb1;
1694 }
1695 if (bb2 > 0)
1696 bb = lshift(bb, bb2);
1697 if (bd5 > 0)
1698 bd = pow5mult(bd, bd5);
1699 if (bd2 > 0)
1700 bd = lshift(bd, bd2);
1701 if (bs2 > 0)
1702 bs = lshift(bs, bs2);
1703 delta = diff(bb, bd);
1704 dsign = delta->sign;
1705 delta->sign = 0;
1706 i = cmp(delta, bs);
1707 if (i < 0) {
1708 /* Error is less than half an ulp -- check for
1709 * special case of mantissa a power of two.
1710 */
1711 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1712 break;
1713 delta = lshift(delta,Log2P);
1714 if (cmp(delta, bs) > 0)
1715 goto drop_down;
1716 break;
1717 }
1718 if (i == 0) {
1719 /* exactly half-way between */
1720 if (dsign) {
1721 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1722 && word1(rv) == 0xffffffff) {
1723 /*boundary case -- increment exponent*/
1724 word0(rv) = (word0(rv) & Exp_mask)
1725 + Exp_msk1
1726 #ifdef IBM
1727 | Exp_msk1 >> 4
1728 #endif
1729 ;
1730 word1(rv) = 0;
1731 break;
1732 }
1733 }
1734 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1735 drop_down:
1736 /* boundary case -- decrement exponent */
1737 #ifdef Sudden_Underflow
1738 L = word0(rv) & Exp_mask;
1739 #ifdef IBM
1740 if (L < Exp_msk1)
1741 #else
1742 if (L <= Exp_msk1)
1743 #endif
1744 goto undfl;
1745 L -= Exp_msk1;
1746 #else
1747 L = (word0(rv) & Exp_mask) - Exp_msk1;
1748 #endif
1749 word0(rv) = L | Bndry_mask1;
1750 word1(rv) = 0xffffffff;
1751 #ifdef IBM
1752 goto cont;
1753 #else
1754 break;
1755 #endif
1756 }
1757 #ifndef ROUND_BIASED
1758 if (!(word1(rv) & LSB))
1759 break;
1760 #endif
1761 if (dsign)
1762 value(rv) += ulp(value(rv));
1763 #ifndef ROUND_BIASED
1764 else {
1765 value(rv) -= ulp(value(rv));
1766 #ifndef Sudden_Underflow
1767 if (!value(rv))
1768 goto undfl;
1769 #endif
1770 }
1771 #endif
1772 break;
1773 }
1774 if ((aadj = ratio(delta, bs)) <= 2.) {
1775 if (dsign)
1776 aadj = aadj1 = 1.;
1777 else if (word1(rv) || word0(rv) & Bndry_mask) {
1778 #ifndef Sudden_Underflow
1779 if (word1(rv) == Tiny1 && !word0(rv))
1780 goto undfl;
1781 #endif
1782 aadj = 1.;
1783 aadj1 = -1.;
1784 }
1785 else {
1786 /* special case -- power of FLT_RADIX to be */
1787 /* rounded down... */
1788
1789 if (aadj < 2./FLT_RADIX)
1790 aadj = 1./FLT_RADIX;
1791 else
1792 aadj *= 0.5;
1793 aadj1 = -aadj;
1794 }
1795 }
1796 else {
1797 aadj *= 0.5;
1798 aadj1 = dsign ? aadj : -aadj;
1799 #ifdef Check_FLT_ROUNDS
1800 switch(FLT_ROUNDS) {
1801 case 2: /* towards +infinity */
1802 aadj1 -= 0.5;
1803 break;
1804 case 0: /* towards 0 */
1805 case 3: /* towards -infinity */
1806 aadj1 += 0.5;
1807 }
1808 #else
1809 if (FLT_ROUNDS == 0)
1810 aadj1 += 0.5;
1811 #endif
1812 }
1813 y = word0(rv) & Exp_mask;
1814
1815 /* Check for overflow */
1816
1817 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1818 value(rv0) = value(rv);
1819 word0(rv) -= P*Exp_msk1;
1820 adj = aadj1 * ulp(value(rv));
1821 value(rv) += adj;
1822 if ((word0(rv) & Exp_mask) >=
1823 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1824 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1825 goto ovfl;
1826 word0(rv) = Big0;
1827 word1(rv) = Big1;
1828 goto cont;
1829 }
1830 else
1831 word0(rv) += P*Exp_msk1;
1832 }
1833 else {
1834 #ifdef Sudden_Underflow
1835 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1836 value(rv0) = value(rv);
1837 word0(rv) += P*Exp_msk1;
1838 adj = aadj1 * ulp(value(rv));
1839 value(rv) += adj;
1840 #ifdef IBM
1841 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1842 #else
1843 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1844 #endif
1845 {
1846 if (word0(rv0) == Tiny0
1847 && word1(rv0) == Tiny1)
1848 goto undfl;
1849 word0(rv) = Tiny0;
1850 word1(rv) = Tiny1;
1851 goto cont;
1852 }
1853 else
1854 word0(rv) -= P*Exp_msk1;
1855 }
1856 else {
1857 adj = aadj1 * ulp(value(rv));
1858 value(rv) += adj;
1859 }
1860 #else
1861 /* Compute adj so that the IEEE rounding rules will
1862 * correctly round rv + adj in some half-way cases.
1863 * If rv * ulp(rv) is denormalized (i.e.,
1864 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1865 * trouble from bits lost to denormalization;
1866 * example: 1.2e-307 .
1867 */
1868 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1869 aadj1 = (double)(int)(aadj + 0.5);
1870 if (!dsign)
1871 aadj1 = -aadj1;
1872 }
1873 adj = aadj1 * ulp(value(rv));
1874 value(rv) += adj;
1875 #endif
1876 }
1877 z = word0(rv) & Exp_mask;
1878 if (y == z) {
1879 /* Can we stop now? */
1880 L = aadj;
1881 aadj -= L;
1882 /* The tolerances below are conservative. */
1883 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1884 if (aadj < .4999999 || aadj > .5000001)
1885 break;
1886 }
1887 else if (aadj < .4999999/FLT_RADIX)
1888 break;
1889 }
1890 cont:
1891 Bfree(bb);
1892 Bfree(bd);
1893 Bfree(bs);
1894 Bfree(delta);
1895 }
1896 retfree:
1897 Bfree(bb);
1898 Bfree(bd);
1899 Bfree(bs);
1900 Bfree(bd0);
1901 Bfree(delta);
1902 ret:
1903 if (se)
1904 /* LINTED interface specification */
1905 *se = (char *)s;
1906 return sign ? -value(rv) : value(rv);
1907 }
1908
1909 static int
quorem(b,S)1910 quorem
1911 #ifdef KR_headers
1912 (b, S) Bigint *b, *S;
1913 #else
1914 (Bigint *b, Bigint *S)
1915 #endif
1916 {
1917 int n;
1918 Long borrow, y;
1919 ULong carry, q, ys;
1920 ULong *bx, *bxe, *sx, *sxe;
1921 #ifdef Pack_32
1922 Long z;
1923 ULong si, zs;
1924 #endif
1925
1926 if (b == BIGINT_INVALID || S == BIGINT_INVALID)
1927 return 0;
1928
1929 n = S->wds;
1930 #ifdef DEBUG
1931 /*debug*/ if (b->wds > n)
1932 /*debug*/ Bug("oversize b in quorem");
1933 #endif
1934 if (b->wds < n)
1935 return 0;
1936 sx = S->x;
1937 sxe = sx + --n;
1938 bx = b->x;
1939 bxe = bx + n;
1940 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1941 #ifdef DEBUG
1942 /*debug*/ if (q > 9)
1943 /*debug*/ Bug("oversized quotient in quorem");
1944 #endif
1945 if (q) {
1946 borrow = 0;
1947 carry = 0;
1948 do {
1949 #ifdef Pack_32
1950 si = *sx++;
1951 ys = (si & 0xffff) * q + carry;
1952 zs = (si >> 16) * q + (ys >> 16);
1953 carry = zs >> 16;
1954 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1955 borrow = (ULong)y >> 16;
1956 Sign_Extend(borrow, y);
1957 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1958 borrow = (ULong)z >> 16;
1959 Sign_Extend(borrow, z);
1960 Storeinc(bx, z, y);
1961 #else
1962 ys = *sx++ * q + carry;
1963 carry = ys >> 16;
1964 y = *bx - (ys & 0xffff) + borrow;
1965 borrow = y >> 16;
1966 Sign_Extend(borrow, y);
1967 *bx++ = y & 0xffff;
1968 #endif
1969 }
1970 while(sx <= sxe);
1971 if (!*bxe) {
1972 bx = b->x;
1973 while(--bxe > bx && !*bxe)
1974 --n;
1975 b->wds = n;
1976 }
1977 }
1978 if (cmp(b, S) >= 0) {
1979 q++;
1980 borrow = 0;
1981 carry = 0;
1982 bx = b->x;
1983 sx = S->x;
1984 do {
1985 #ifdef Pack_32
1986 si = *sx++;
1987 ys = (si & 0xffff) + carry;
1988 zs = (si >> 16) + (ys >> 16);
1989 carry = zs >> 16;
1990 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1991 borrow = (ULong)y >> 16;
1992 Sign_Extend(borrow, y);
1993 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1994 borrow = (ULong)z >> 16;
1995 Sign_Extend(borrow, z);
1996 Storeinc(bx, z, y);
1997 #else
1998 ys = *sx++ + carry;
1999 carry = ys >> 16;
2000 y = *bx - (ys & 0xffff) + borrow;
2001 borrow = y >> 16;
2002 Sign_Extend(borrow, y);
2003 *bx++ = y & 0xffff;
2004 #endif
2005 }
2006 while(sx <= sxe);
2007 bx = b->x;
2008 bxe = bx + n;
2009 if (!*bxe) {
2010 while(--bxe > bx && !*bxe)
2011 --n;
2012 b->wds = n;
2013 }
2014 }
2015 return q;
2016 }
2017
2018 /* freedtoa(s) must be used to free values s returned by dtoa
2019 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2020 * but for consistency with earlier versions of dtoa, it is optional
2021 * when MULTIPLE_THREADS is not defined.
2022 */
2023
2024 void
2025 #ifdef KR_headers
freedtoa(s)2026 freedtoa(s) char *s;
2027 #else
2028 freedtoa(char *s)
2029 #endif
2030 {
2031 free(s);
2032 }
2033
2034
2035
2036 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2037 *
2038 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2039 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
2040 *
2041 * Modifications:
2042 * 1. Rather than iterating, we use a simple numeric overestimate
2043 * to determine k = floor(log10(d)). We scale relevant
2044 * quantities using O(log2(k)) rather than O(k) multiplications.
2045 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2046 * try to generate digits strictly left to right. Instead, we
2047 * compute with fewer bits and propagate the carry if necessary
2048 * when rounding the final digit up. This is often faster.
2049 * 3. Under the assumption that input will be rounded nearest,
2050 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2051 * That is, we allow equality in stopping tests when the
2052 * round-nearest rule will give the same floating-point value
2053 * as would satisfaction of the stopping test with strict
2054 * inequality.
2055 * 4. We remove common factors of powers of 2 from relevant
2056 * quantities.
2057 * 5. When converting floating-point integers less than 1e16,
2058 * we use floating-point arithmetic rather than resorting
2059 * to multiple-precision integers.
2060 * 6. When asked to produce fewer than 15 digits, we first try
2061 * to get by with floating-point arithmetic; we resort to
2062 * multiple-precision integer arithmetic only if we cannot
2063 * guarantee that the floating-point calculation has given
2064 * the correctly rounded result. For k requested digits and
2065 * "uniformly" distributed input, the probability is
2066 * something like 10^(k-15) that we must resort to the Long
2067 * calculation.
2068 */
2069
2070 __LIBC_HIDDEN__ char *
__dtoa(_d,mode,ndigits,decpt,sign,rve)2071 __dtoa
2072 #ifdef KR_headers
2073 (_d, mode, ndigits, decpt, sign, rve)
2074 double _d; int mode, ndigits, *decpt, *sign; char **rve;
2075 #else
2076 (double _d, int mode, int ndigits, int *decpt, int *sign, char **rve)
2077 #endif
2078 {
2079 /* Arguments ndigits, decpt, sign are similar to those
2080 of ecvt and fcvt; trailing zeros are suppressed from
2081 the returned string. If not null, *rve is set to point
2082 to the end of the return value. If d is +-Infinity or NaN,
2083 then *decpt is set to 9999.
2084
2085 mode:
2086 0 ==> shortest string that yields d when read in
2087 and rounded to nearest.
2088 1 ==> like 0, but with Steele & White stopping rule;
2089 e.g. with IEEE P754 arithmetic , mode 0 gives
2090 1e23 whereas mode 1 gives 9.999999999999999e22.
2091 2 ==> max(1,ndigits) significant digits. This gives a
2092 return value similar to that of ecvt, except
2093 that trailing zeros are suppressed.
2094 3 ==> through ndigits past the decimal point. This
2095 gives a return value similar to that from fcvt,
2096 except that trailing zeros are suppressed, and
2097 ndigits can be negative.
2098 4-9 should give the same return values as 2-3, i.e.,
2099 4 <= mode <= 9 ==> same return as mode
2100 2 + (mode & 1). These modes are mainly for
2101 debugging; often they run slower but sometimes
2102 faster than modes 2-3.
2103 4,5,8,9 ==> left-to-right digit generation.
2104 6-9 ==> don't try fast floating-point estimate
2105 (if applicable).
2106
2107 Values of mode other than 0-9 are treated as mode 0.
2108
2109 Sufficient space is allocated to the return value
2110 to hold the suppressed trailing zeros.
2111 */
2112
2113 int bbits, b2, b5, be, dig, i, ieps, ilim0,
2114 j, jj1, k, k0, k_check, leftright, m2, m5, s2, s5,
2115 try_quick;
2116 int ilim = 0, ilim1 = 0, spec_case = 0; /* pacify gcc */
2117 Long L;
2118 #ifndef Sudden_Underflow
2119 int denorm;
2120 ULong x;
2121 #endif
2122 Bigint *b, *b1, *delta, *mhi, *S;
2123 Bigint *mlo = NULL; /* pacify gcc */
2124 double ds;
2125 char *s, *s0;
2126 Bigint *result = NULL;
2127 int result_k = 0;
2128 _double d, d2, eps;
2129
2130 value(d) = _d;
2131
2132 if (word0(d) & Sign_bit) {
2133 /* set sign for everything, including 0's and NaNs */
2134 *sign = 1;
2135 word0(d) &= ~Sign_bit; /* clear sign bit */
2136 }
2137 else
2138 *sign = 0;
2139
2140 #if defined(IEEE_Arith) + defined(VAX)
2141 #ifdef IEEE_Arith
2142 if ((word0(d) & Exp_mask) == Exp_mask)
2143 #else
2144 if (word0(d) == 0x8000)
2145 #endif
2146 {
2147 /* Infinity or NaN */
2148 *decpt = 9999;
2149 s =
2150 #ifdef IEEE_Arith
2151 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
2152 #endif
2153 "NaN";
2154 result = Balloc(strlen(s)+1);
2155 if (result == BIGINT_INVALID)
2156 return NULL;
2157 s0 = (char *)(void *)result;
2158 strcpy(s0, s);
2159 if (rve)
2160 *rve =
2161 #ifdef IEEE_Arith
2162 s0[3] ? s0 + 8 :
2163 #endif
2164 s0 + 3;
2165 return s0;
2166 }
2167 #endif
2168 #ifdef IBM
2169 value(d) += 0; /* normalize */
2170 #endif
2171 if (!value(d)) {
2172 *decpt = 1;
2173 result = Balloc(2);
2174 if (result == BIGINT_INVALID)
2175 return NULL;
2176 s0 = (char *)(void *)result;
2177 strcpy(s0, "0");
2178 if (rve)
2179 *rve = s0 + 1;
2180 return s0;
2181 }
2182
2183 b = d2b(value(d), &be, &bbits);
2184 #ifdef Sudden_Underflow
2185 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
2186 #else
2187 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) {
2188 #endif
2189 value(d2) = value(d);
2190 word0(d2) &= Frac_mask1;
2191 word0(d2) |= Exp_11;
2192 #ifdef IBM
2193 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
2194 value(d2) /= 1 << j;
2195 #endif
2196
2197 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2198 * log10(x) = log(x) / log(10)
2199 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2200 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2201 *
2202 * This suggests computing an approximation k to log10(d) by
2203 *
2204 * k = (i - Bias)*0.301029995663981
2205 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2206 *
2207 * We want k to be too large rather than too small.
2208 * The error in the first-order Taylor series approximation
2209 * is in our favor, so we just round up the constant enough
2210 * to compensate for any error in the multiplication of
2211 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2212 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2213 * adding 1e-13 to the constant term more than suffices.
2214 * Hence we adjust the constant term to 0.1760912590558.
2215 * (We could get a more accurate k by invoking log10,
2216 * but this is probably not worthwhile.)
2217 */
2218
2219 i -= Bias;
2220 #ifdef IBM
2221 i <<= 2;
2222 i += j;
2223 #endif
2224 #ifndef Sudden_Underflow
2225 denorm = 0;
2226 }
2227 else {
2228 /* d is denormalized */
2229
2230 i = bbits + be + (Bias + (P-1) - 1);
2231 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
2232 : word1(d) << (32 - i);
2233 value(d2) = x;
2234 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2235 i -= (Bias + (P-1) - 1) + 1;
2236 denorm = 1;
2237 }
2238 #endif
2239 ds = (value(d2)-1.5)*0.289529654602168 + 0.1760912590558 +
2240 i*0.301029995663981;
2241 k = (int)ds;
2242 if (ds < 0. && ds != k)
2243 k--; /* want k = floor(ds) */
2244 k_check = 1;
2245 if (k >= 0 && k <= Ten_pmax) {
2246 if (value(d) < tens[k])
2247 k--;
2248 k_check = 0;
2249 }
2250 j = bbits - i - 1;
2251 if (j >= 0) {
2252 b2 = 0;
2253 s2 = j;
2254 }
2255 else {
2256 b2 = -j;
2257 s2 = 0;
2258 }
2259 if (k >= 0) {
2260 b5 = 0;
2261 s5 = k;
2262 s2 += k;
2263 }
2264 else {
2265 b2 -= k;
2266 b5 = -k;
2267 s5 = 0;
2268 }
2269 if (mode < 0 || mode > 9)
2270 mode = 0;
2271 try_quick = 1;
2272 if (mode > 5) {
2273 mode -= 4;
2274 try_quick = 0;
2275 }
2276 leftright = 1;
2277 switch(mode) {
2278 case 0:
2279 case 1:
2280 ilim = ilim1 = -1;
2281 i = 18;
2282 ndigits = 0;
2283 break;
2284 case 2:
2285 leftright = 0;
2286 /* FALLTHROUGH */
2287 case 4:
2288 if (ndigits <= 0)
2289 ndigits = 1;
2290 ilim = ilim1 = i = ndigits;
2291 break;
2292 case 3:
2293 leftright = 0;
2294 /* FALLTHROUGH */
2295 case 5:
2296 i = ndigits + k + 1;
2297 ilim = i;
2298 ilim1 = i - 1;
2299 if (i <= 0)
2300 i = 1;
2301 }
2302 j = sizeof(ULong);
2303 for(result_k = 0; (int)(sizeof(Bigint) - sizeof(ULong)) + j <= i;
2304 j <<= 1) result_k++;
2305 // this is really a ugly hack, the code uses Balloc
2306 // instead of malloc, but casts the result into a char*
2307 // it seems the only reason to do that is due to the
2308 // complicated way the block size need to be computed
2309 // buuurk....
2310 result = Balloc(result_k);
2311 if (result == BIGINT_INVALID) {
2312 Bfree(b);
2313 return NULL;
2314 }
2315 s = s0 = (char *)(void *)result;
2316
2317 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2318
2319 /* Try to get by with floating-point arithmetic. */
2320
2321 i = 0;
2322 value(d2) = value(d);
2323 k0 = k;
2324 ilim0 = ilim;
2325 ieps = 2; /* conservative */
2326 if (k > 0) {
2327 ds = tens[k&0xf];
2328 j = (unsigned int)k >> 4;
2329 if (j & Bletch) {
2330 /* prevent overflows */
2331 j &= Bletch - 1;
2332 value(d) /= bigtens[n_bigtens-1];
2333 ieps++;
2334 }
2335 for(; j; j = (unsigned int)j >> 1, i++)
2336 if (j & 1) {
2337 ieps++;
2338 ds *= bigtens[i];
2339 }
2340 value(d) /= ds;
2341 }
2342 else if ((jj1 = -k) != 0) {
2343 value(d) *= tens[jj1 & 0xf];
2344 for(j = (unsigned int)jj1 >> 4; j;
2345 j = (unsigned int)j >> 1, i++)
2346 if (j & 1) {
2347 ieps++;
2348 value(d) *= bigtens[i];
2349 }
2350 }
2351 if (k_check && value(d) < 1. && ilim > 0) {
2352 if (ilim1 <= 0)
2353 goto fast_failed;
2354 ilim = ilim1;
2355 k--;
2356 value(d) *= 10.;
2357 ieps++;
2358 }
2359 value(eps) = ieps*value(d) + 7.;
2360 word0(eps) -= (P-1)*Exp_msk1;
2361 if (ilim == 0) {
2362 S = mhi = 0;
2363 value(d) -= 5.;
2364 if (value(d) > value(eps))
2365 goto one_digit;
2366 if (value(d) < -value(eps))
2367 goto no_digits;
2368 goto fast_failed;
2369 }
2370 #ifndef No_leftright
2371 if (leftright) {
2372 /* Use Steele & White method of only
2373 * generating digits needed.
2374 */
2375 value(eps) = 0.5/tens[ilim-1] - value(eps);
2376 for(i = 0;;) {
2377 L = value(d);
2378 value(d) -= L;
2379 *s++ = '0' + (int)L;
2380 if (value(d) < value(eps))
2381 goto ret1;
2382 if (1. - value(d) < value(eps))
2383 goto bump_up;
2384 if (++i >= ilim)
2385 break;
2386 value(eps) *= 10.;
2387 value(d) *= 10.;
2388 }
2389 }
2390 else {
2391 #endif
2392 /* Generate ilim digits, then fix them up. */
2393 value(eps) *= tens[ilim-1];
2394 for(i = 1;; i++, value(d) *= 10.) {
2395 L = value(d);
2396 value(d) -= L;
2397 *s++ = '0' + (int)L;
2398 if (i == ilim) {
2399 if (value(d) > 0.5 + value(eps))
2400 goto bump_up;
2401 else if (value(d) < 0.5 - value(eps)) {
2402 while(*--s == '0');
2403 s++;
2404 goto ret1;
2405 }
2406 break;
2407 }
2408 }
2409 #ifndef No_leftright
2410 }
2411 #endif
2412 fast_failed:
2413 s = s0;
2414 value(d) = value(d2);
2415 k = k0;
2416 ilim = ilim0;
2417 }
2418
2419 /* Do we have a "small" integer? */
2420
2421 if (be >= 0 && k <= Int_max) {
2422 /* Yes. */
2423 ds = tens[k];
2424 if (ndigits < 0 && ilim <= 0) {
2425 S = mhi = 0;
2426 if (ilim < 0 || value(d) <= 5*ds)
2427 goto no_digits;
2428 goto one_digit;
2429 }
2430 for(i = 1;; i++) {
2431 L = value(d) / ds;
2432 value(d) -= L*ds;
2433 #ifdef Check_FLT_ROUNDS
2434 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2435 if (value(d) < 0) {
2436 L--;
2437 value(d) += ds;
2438 }
2439 #endif
2440 *s++ = '0' + (int)L;
2441 if (i == ilim) {
2442 value(d) += value(d);
2443 if (value(d) > ds || (value(d) == ds && L & 1)) {
2444 bump_up:
2445 while(*--s == '9')
2446 if (s == s0) {
2447 k++;
2448 *s = '0';
2449 break;
2450 }
2451 ++*s++;
2452 }
2453 break;
2454 }
2455 if (!(value(d) *= 10.))
2456 break;
2457 }
2458 goto ret1;
2459 }
2460
2461 m2 = b2;
2462 m5 = b5;
2463 mhi = mlo = 0;
2464 if (leftright) {
2465 if (mode < 2) {
2466 i =
2467 #ifndef Sudden_Underflow
2468 denorm ? be + (Bias + (P-1) - 1 + 1) :
2469 #endif
2470 #ifdef IBM
2471 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2472 #else
2473 1 + P - bbits;
2474 #endif
2475 }
2476 else {
2477 j = ilim - 1;
2478 if (m5 >= j)
2479 m5 -= j;
2480 else {
2481 s5 += j -= m5;
2482 b5 += j;
2483 m5 = 0;
2484 }
2485 if ((i = ilim) < 0) {
2486 m2 -= i;
2487 i = 0;
2488 }
2489 }
2490 b2 += i;
2491 s2 += i;
2492 mhi = i2b(1);
2493 }
2494 if (m2 > 0 && s2 > 0) {
2495 i = m2 < s2 ? m2 : s2;
2496 b2 -= i;
2497 m2 -= i;
2498 s2 -= i;
2499 }
2500 if (b5 > 0) {
2501 if (leftright) {
2502 if (m5 > 0) {
2503 mhi = pow5mult(mhi, m5);
2504 b1 = mult(mhi, b);
2505 Bfree(b);
2506 b = b1;
2507 }
2508 if ((j = b5 - m5) != 0)
2509 b = pow5mult(b, j);
2510 }
2511 else
2512 b = pow5mult(b, b5);
2513 }
2514 S = i2b(1);
2515 if (s5 > 0)
2516 S = pow5mult(S, s5);
2517
2518 /* Check for special case that d is a normalized power of 2. */
2519
2520 if (mode < 2) {
2521 if (!word1(d) && !(word0(d) & Bndry_mask)
2522 #ifndef Sudden_Underflow
2523 && word0(d) & Exp_mask
2524 #endif
2525 ) {
2526 /* The special case */
2527 b2 += Log2P;
2528 s2 += Log2P;
2529 spec_case = 1;
2530 }
2531 else
2532 spec_case = 0;
2533 }
2534
2535 /* Arrange for convenient computation of quotients:
2536 * shift left if necessary so divisor has 4 leading 0 bits.
2537 *
2538 * Perhaps we should just compute leading 28 bits of S once
2539 * and for all and pass them and a shift to quorem, so it
2540 * can do shifts and ors to compute the numerator for q.
2541 */
2542 if (S == BIGINT_INVALID) {
2543 i = 0;
2544 } else {
2545 #ifdef Pack_32
2546 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0)
2547 i = 32 - i;
2548 #else
2549 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
2550 i = 16 - i;
2551 #endif
2552 }
2553
2554 if (i > 4) {
2555 i -= 4;
2556 b2 += i;
2557 m2 += i;
2558 s2 += i;
2559 }
2560 else if (i < 4) {
2561 i += 28;
2562 b2 += i;
2563 m2 += i;
2564 s2 += i;
2565 }
2566 if (b2 > 0)
2567 b = lshift(b, b2);
2568 if (s2 > 0)
2569 S = lshift(S, s2);
2570 if (k_check) {
2571 if (cmp(b,S) < 0) {
2572 k--;
2573 b = multadd(b, 10, 0); /* we botched the k estimate */
2574 if (leftright)
2575 mhi = multadd(mhi, 10, 0);
2576 ilim = ilim1;
2577 }
2578 }
2579 if (ilim <= 0 && mode > 2) {
2580 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2581 /* no digits, fcvt style */
2582 no_digits:
2583 k = -1 - ndigits;
2584 goto ret;
2585 }
2586 one_digit:
2587 *s++ = '1';
2588 k++;
2589 goto ret;
2590 }
2591 if (leftright) {
2592 if (m2 > 0)
2593 mhi = lshift(mhi, m2);
2594
2595 /* Compute mlo -- check for special case
2596 * that d is a normalized power of 2.
2597 */
2598
2599 mlo = mhi;
2600 if (spec_case) {
2601 mhi = Balloc(mhi->k);
2602 Bcopy(mhi, mlo);
2603 mhi = lshift(mhi, Log2P);
2604 }
2605
2606 for(i = 1;;i++) {
2607 dig = quorem(b,S) + '0';
2608 /* Do we yet have the shortest decimal string
2609 * that will round to d?
2610 */
2611 j = cmp(b, mlo);
2612 delta = diff(S, mhi);
2613 jj1 = delta->sign ? 1 : cmp(b, delta);
2614 Bfree(delta);
2615 #ifndef ROUND_BIASED
2616 if (jj1 == 0 && !mode && !(word1(d) & 1)) {
2617 if (dig == '9')
2618 goto round_9_up;
2619 if (j > 0)
2620 dig++;
2621 *s++ = dig;
2622 goto ret;
2623 }
2624 #endif
2625 if (j < 0 || (j == 0 && !mode
2626 #ifndef ROUND_BIASED
2627 && !(word1(d) & 1)
2628 #endif
2629 )) {
2630 if (jj1 > 0) {
2631 b = lshift(b, 1);
2632 jj1 = cmp(b, S);
2633 if ((jj1 > 0 || (jj1 == 0 && dig & 1))
2634 && dig++ == '9')
2635 goto round_9_up;
2636 }
2637 *s++ = dig;
2638 goto ret;
2639 }
2640 if (jj1 > 0) {
2641 if (dig == '9') { /* possible if i == 1 */
2642 round_9_up:
2643 *s++ = '9';
2644 goto roundoff;
2645 }
2646 *s++ = dig + 1;
2647 goto ret;
2648 }
2649 *s++ = dig;
2650 if (i == ilim)
2651 break;
2652 b = multadd(b, 10, 0);
2653 if (mlo == mhi)
2654 mlo = mhi = multadd(mhi, 10, 0);
2655 else {
2656 mlo = multadd(mlo, 10, 0);
2657 mhi = multadd(mhi, 10, 0);
2658 }
2659 }
2660 }
2661 else
2662 for(i = 1;; i++) {
2663 *s++ = dig = quorem(b,S) + '0';
2664 if (i >= ilim)
2665 break;
2666 b = multadd(b, 10, 0);
2667 }
2668
2669 /* Round off last digit */
2670
2671 b = lshift(b, 1);
2672 j = cmp(b, S);
2673 if (j > 0 || (j == 0 && dig & 1)) {
2674 roundoff:
2675 while(*--s == '9')
2676 if (s == s0) {
2677 k++;
2678 *s++ = '1';
2679 goto ret;
2680 }
2681 ++*s++;
2682 }
2683 else {
2684 while(*--s == '0');
2685 s++;
2686 }
2687 ret:
2688 Bfree(S);
2689 if (mhi) {
2690 if (mlo && mlo != mhi)
2691 Bfree(mlo);
2692 Bfree(mhi);
2693 }
2694 ret1:
2695 Bfree(b);
2696 if (s == s0) { /* don't return empty string */
2697 *s++ = '0';
2698 k = 0;
2699 }
2700 *s = 0;
2701 *decpt = k + 1;
2702 if (rve)
2703 *rve = s;
2704 return s0;
2705 }
2706 #ifdef __cplusplus
2707 }
2708 #endif
2709