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