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