Ruby 3.5.0dev (2025-04-25 revision c772d2691dabd15f2b6fcb8bddce64c4385b3b23)
dtoa.c
1/****************************************************************
2 *
3 * The author of this software is David M. Gay.
4 *
5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
6 *
7 * Permission to use, copy, modify, and distribute this software for any
8 * purpose without fee is hereby granted, provided that this entire notice
9 * is included in all copies of any software which is or includes a copy
10 * or modification of this software and in all copies of the supporting
11 * documentation for such software.
12 *
13 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
15 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
17 *
18 ***************************************************************/
19
20/* Please send bug reports to David M. Gay (dmg at acm dot org,
21 * with " at " changed at "@" and " dot " changed to "."). */
22
23/* On a machine with IEEE extended-precision registers, it is
24 * necessary to specify double-precision (53-bit) rounding precision
25 * before invoking strtod or dtoa. If the machine uses (the equivalent
26 * of) Intel 80x87 arithmetic, the call
27 * _control87(PC_53, MCW_PC);
28 * does this with many compilers. Whether this or another call is
29 * appropriate depends on the compiler; for this to work, it may be
30 * necessary to #include "float.h" or another system-dependent header
31 * file.
32 */
33
34/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
35 *
36 * This strtod returns a nearest machine number to the input decimal
37 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
38 * broken by the IEEE round-even rule. Otherwise ties are broken by
39 * biased rounding (add half and chop).
40 *
41 * Inspired loosely by William D. Clinger's paper "How to Read Floating
42 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
43 *
44 * Modifications:
45 *
46 * 1. We only require IEEE, IBM, or VAX double-precision
47 * arithmetic (not IEEE double-extended).
48 * 2. We get by with floating-point arithmetic in a case that
49 * Clinger missed -- when we're computing d * 10^n
50 * for a small integer d and the integer n is not too
51 * much larger than 22 (the maximum integer k for which
52 * we can represent 10^k exactly), we may be able to
53 * compute (d*10^k) * 10^(e-k) with just one roundoff.
54 * 3. Rather than a bit-at-a-time adjustment of the binary
55 * result in the hard case, we use floating-point
56 * arithmetic to determine the adjustment to within
57 * one bit; only in really hard cases do we need to
58 * compute a second residual.
59 * 4. Because of 3., we don't need a large table of powers of 10
60 * for ten-to-e (just some small tables, e.g. of 10^k
61 * for 0 <= k <= 22).
62 */
63
64/*
65 * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least
66 * significant byte has the lowest address.
67 * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most
68 * significant byte has the lowest address.
69 * #define Long int on machines with 32-bit ints and 64-bit longs.
70 * #define IBM for IBM mainframe-style floating-point arithmetic.
71 * #define VAX for VAX-style floating-point arithmetic (D_floating).
72 * #define No_leftright to omit left-right logic in fast floating-point
73 * computation of dtoa.
74 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
75 * and strtod and dtoa should round accordingly.
76 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
77 * and Honor_FLT_ROUNDS is not #defined.
78 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
79 * that use extended-precision instructions to compute rounded
80 * products and quotients) with IBM.
81 * #define ROUND_BIASED for IEEE-format with biased rounding.
82 * #define Inaccurate_Divide for IEEE-format with correctly rounded
83 * products but inaccurate quotients, e.g., for Intel i860.
84 * #define NO_LONG_LONG on machines that do not have a "long long"
85 * integer type (of >= 64 bits). On such machines, you can
86 * #define Just_16 to store 16 bits per 32-bit Long when doing
87 * high-precision integer arithmetic. Whether this speeds things
88 * up or slows things down depends on the machine and the number
89 * being converted. If long long is available and the name is
90 * something other than "long long", #define Llong to be the name,
91 * and if "unsigned Llong" does not work as an unsigned version of
92 * Llong, #define #ULLong to be the corresponding unsigned type.
93 * #define KR_headers for old-style C function headers.
94 * #define Bad_float_h if your system lacks a float.h or if it does not
95 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
96 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
97 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
98 * if memory is available and otherwise does something you deem
99 * appropriate. If MALLOC is undefined, malloc will be invoked
100 * directly -- and assumed always to succeed.
101 * #define INFNAN_CHECK on IEEE systems to cause strtod to check for
102 * Infinity and NaN (case insensitively). On some systems (e.g.,
103 * some HP systems), it may be necessary to #define NAN_WORD0
104 * appropriately -- to the most significant word of a quiet NaN.
105 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
106 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
107 * strtod also accepts (case insensitively) strings of the form
108 * NaN(x), where x is a string of hexadecimal digits and spaces;
109 * if there is only one string of hexadecimal digits, it is taken
110 * for the 52 fraction bits of the resulting NaN; if there are two
111 * or more strings of hex digits, the first is for the high 20 bits,
112 * the second and subsequent for the low 32 bits, with intervening
113 * white space ignored; but if this results in none of the 52
114 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
115 * and NAN_WORD1 are used instead.
116 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
117 * multiple threads. In this case, you must provide (or suitably
118 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
119 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
120 * in pow5mult, ensures lazy evaluation of only one copy of high
121 * powers of 5; omitting this lock would introduce a small
122 * probability of wasting memory, but would otherwise be harmless.)
123 * You must also invoke freedtoa(s) to free the value s returned by
124 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
125 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
126 * avoids underflows on inputs whose result does not underflow.
127 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
128 * floating-point numbers and flushes underflows to zero rather
129 * than implementing gradual underflow, then you must also #define
130 * Sudden_Underflow.
131 * #define YES_ALIAS to permit aliasing certain double values with
132 * arrays of ULongs. This leads to slightly better code with
133 * some compilers and was always used prior to 19990916, but it
134 * is not strictly legal and can cause trouble with aggressively
135 * optimizing compilers (e.g., gcc 2.95.1 under -O2).
136 * #define USE_LOCALE to use the current locale's decimal_point value.
137 * #define SET_INEXACT if IEEE arithmetic is being used and extra
138 * computation should be done to set the inexact flag when the
139 * result is inexact and avoid setting inexact when the result
140 * is exact. In this case, dtoa.c must be compiled in
141 * an environment, perhaps provided by #include "dtoa.c" in a
142 * suitable wrapper, that defines two functions,
143 * int get_inexact(void);
144 * void clear_inexact(void);
145 * such that get_inexact() returns a nonzero value if the
146 * inexact bit is already set, and clear_inexact() sets the
147 * inexact bit to 0. When SET_INEXACT is #defined, strtod
148 * also does extra computations to set the underflow and overflow
149 * flags when appropriate (i.e., when the result is tiny and
150 * inexact or when it is a numeric value rounded to +-infinity).
151 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
152 * the result overflows to +-Infinity or underflows to 0.
153 */
154
155#ifdef WORDS_BIGENDIAN
156#define IEEE_BIG_ENDIAN
157#else
158#define IEEE_LITTLE_ENDIAN
159#endif
160
161#ifdef __vax__
162#define VAX
163#undef IEEE_BIG_ENDIAN
164#undef IEEE_LITTLE_ENDIAN
165#endif
166
167#if defined(__arm__) && !defined(__VFP_FP__)
168#define IEEE_BIG_ENDIAN
169#undef IEEE_LITTLE_ENDIAN
170#endif
171
172#undef Long
173#undef ULong
174
175#include <assert.h>
176#include <limits.h>
177#include <stddef.h>
178#include <stdint.h>
179
180#if (INT_MAX >> 30) && !(INT_MAX >> 31)
181#define Long int
182#define ULong unsigned int
183#elif (LONG_MAX >> 30) && !(LONG_MAX >> 31)
184#define Long long int
185#define ULong unsigned long int
186#else
187#error No 32bit integer
188#endif
189
190#if defined(HAVE_LONG_LONG) && (HAVE_LONG_LONG)
191#define Llong LONG_LONG
192#else
193#define NO_LONG_LONG
194#endif
195
196#ifdef DEBUG
197#include <stdio.h>
198#define Bug(x) {fprintf(stderr, "%s\n", (x)); exit(EXIT_FAILURE);}
199#endif
200
201#ifndef ISDIGIT
202#include <ctype.h>
203#define ISDIGIT(c) isdigit(c)
204#endif
205#include <errno.h>
206#include <stdlib.h>
207#include <string.h>
208
209#ifdef USE_LOCALE
210#include <locale.h>
211#endif
212
213#ifdef MALLOC
214extern void *MALLOC(size_t);
215#else
216#define MALLOC malloc
217#endif
218#ifdef FREE
219extern void FREE(void*);
220#else
221#define FREE free
222#endif
223#ifndef NO_SANITIZE
224#define NO_SANITIZE(x, y) y
225#endif
226
227#undef IEEE_Arith
228#undef Avoid_Underflow
229#ifdef IEEE_BIG_ENDIAN
230#define IEEE_Arith
231#endif
232#ifdef IEEE_LITTLE_ENDIAN
233#define IEEE_Arith
234#endif
235
236#ifdef Bad_float_h
237
238#ifdef IEEE_Arith
239#define DBL_DIG 15
240#define DBL_MAX_10_EXP 308
241#define DBL_MAX_EXP 1024
242#define FLT_RADIX 2
243#endif /*IEEE_Arith*/
244
245#ifdef IBM
246#define DBL_DIG 16
247#define DBL_MAX_10_EXP 75
248#define DBL_MAX_EXP 63
249#define FLT_RADIX 16
250#define DBL_MAX 7.2370055773322621e+75
251#endif
252
253#ifdef VAX
254#define DBL_DIG 16
255#define DBL_MAX_10_EXP 38
256#define DBL_MAX_EXP 127
257#define FLT_RADIX 2
258#define DBL_MAX 1.7014118346046923e+38
259#endif
260
261#ifndef LONG_MAX
262#define LONG_MAX 2147483647
263#endif
264
265#else /* ifndef Bad_float_h */
266#include <float.h>
267#endif /* Bad_float_h */
268
269#include <math.h>
270
271#ifdef __cplusplus
272extern "C" {
273#if 0
274} /* satisfy cc-mode */
275#endif
276#endif
277
278#ifndef hexdigit
279static const char hexdigit[] = "0123456789abcdef0123456789ABCDEF";
280#endif
281
282#if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + defined(IBM) != 1
283Exactly one of IEEE_LITTLE_ENDIAN, IEEE_BIG_ENDIAN, VAX, or IBM should be defined.
284#endif
285
286typedef union { double d; ULong L[2]; } U;
287
288#ifdef YES_ALIAS
289typedef double double_u;
290# define dval(x) (x)
291# ifdef IEEE_LITTLE_ENDIAN
292# define word0(x) (((ULong *)&(x))[1])
293# define word1(x) (((ULong *)&(x))[0])
294# else
295# define word0(x) (((ULong *)&(x))[0])
296# define word1(x) (((ULong *)&(x))[1])
297# endif
298#else
299typedef U double_u;
300# ifdef IEEE_LITTLE_ENDIAN
301# define word0(x) ((x).L[1])
302# define word1(x) ((x).L[0])
303# else
304# define word0(x) ((x).L[0])
305# define word1(x) ((x).L[1])
306# endif
307# define dval(x) ((x).d)
308#endif
309
310/* The following definition of Storeinc is appropriate for MIPS processors.
311 * An alternative that might be better on some machines is
312 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
313 */
314#if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm__)
315#define Storeinc(a,b,c) (((unsigned short *)(a))[1] = (unsigned short)(b), \
316((unsigned short *)(a))[0] = (unsigned short)(c), (a)++)
317#else
318#define Storeinc(a,b,c) (((unsigned short *)(a))[0] = (unsigned short)(b), \
319((unsigned short *)(a))[1] = (unsigned short)(c), (a)++)
320#endif
321
322/* #define P DBL_MANT_DIG */
323/* Ten_pmax = floor(P*log(2)/log(5)) */
324/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
325/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
326/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
327
328#ifdef IEEE_Arith
329#define Exp_shift 20
330#define Exp_shift1 20
331#define Exp_msk1 0x100000
332#define Exp_msk11 0x100000
333#define Exp_mask 0x7ff00000
334#define P 53
335#define Bias 1023
336#define Emin (-1022)
337#define Exp_1 0x3ff00000
338#define Exp_11 0x3ff00000
339#define Ebits 11
340#define Frac_mask 0xfffff
341#define Frac_mask1 0xfffff
342#define Ten_pmax 22
343#define Bletch 0x10
344#define Bndry_mask 0xfffff
345#define Bndry_mask1 0xfffff
346#define LSB 1
347#define Sign_bit 0x80000000
348#define Log2P 1
349#define Tiny0 0
350#define Tiny1 1
351#define Quick_max 14
352#define Int_max 14
353#ifndef NO_IEEE_Scale
354#define Avoid_Underflow
355#ifdef Flush_Denorm /* debugging option */
356#undef Sudden_Underflow
357#endif
358#endif
359
360#ifndef Flt_Rounds
361#ifdef FLT_ROUNDS
362#define Flt_Rounds FLT_ROUNDS
363#else
364#define Flt_Rounds 1
365#endif
366#endif /*Flt_Rounds*/
367
368#ifdef Honor_FLT_ROUNDS
369#define Rounding rounding
370#undef Check_FLT_ROUNDS
371#define Check_FLT_ROUNDS
372#else
373#define Rounding Flt_Rounds
374#endif
375
376#else /* ifndef IEEE_Arith */
377#undef Check_FLT_ROUNDS
378#undef Honor_FLT_ROUNDS
379#undef SET_INEXACT
380#undef Sudden_Underflow
381#define Sudden_Underflow
382#ifdef IBM
383#undef Flt_Rounds
384#define Flt_Rounds 0
385#define Exp_shift 24
386#define Exp_shift1 24
387#define Exp_msk1 0x1000000
388#define Exp_msk11 0x1000000
389#define Exp_mask 0x7f000000
390#define P 14
391#define Bias 65
392#define Exp_1 0x41000000
393#define Exp_11 0x41000000
394#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
395#define Frac_mask 0xffffff
396#define Frac_mask1 0xffffff
397#define Bletch 4
398#define Ten_pmax 22
399#define Bndry_mask 0xefffff
400#define Bndry_mask1 0xffffff
401#define LSB 1
402#define Sign_bit 0x80000000
403#define Log2P 4
404#define Tiny0 0x100000
405#define Tiny1 0
406#define Quick_max 14
407#define Int_max 15
408#else /* VAX */
409#undef Flt_Rounds
410#define Flt_Rounds 1
411#define Exp_shift 23
412#define Exp_shift1 7
413#define Exp_msk1 0x80
414#define Exp_msk11 0x800000
415#define Exp_mask 0x7f80
416#define P 56
417#define Bias 129
418#define Exp_1 0x40800000
419#define Exp_11 0x4080
420#define Ebits 8
421#define Frac_mask 0x7fffff
422#define Frac_mask1 0xffff007f
423#define Ten_pmax 24
424#define Bletch 2
425#define Bndry_mask 0xffff007f
426#define Bndry_mask1 0xffff007f
427#define LSB 0x10000
428#define Sign_bit 0x8000
429#define Log2P 1
430#define Tiny0 0x80
431#define Tiny1 0
432#define Quick_max 15
433#define Int_max 15
434#endif /* IBM, VAX */
435#endif /* IEEE_Arith */
436
437#ifndef IEEE_Arith
438#define ROUND_BIASED
439#endif
440
441#ifdef RND_PRODQUOT
442#define rounded_product(a,b) ((a) = rnd_prod((a), (b)))
443#define rounded_quotient(a,b) ((a) = rnd_quot((a), (b)))
444extern double rnd_prod(double, double), rnd_quot(double, double);
445#else
446#define rounded_product(a,b) ((a) *= (b))
447#define rounded_quotient(a,b) ((a) /= (b))
448#endif
449
450#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
451#define Big1 0xffffffff
452
453#ifndef Pack_32
454#define Pack_32
455#endif
456
457#define FFFFFFFF 0xffffffffUL
458
459#ifdef NO_LONG_LONG
460#undef ULLong
461#ifdef Just_16
462#undef Pack_32
463/* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
464 * This makes some inner loops simpler and sometimes saves work
465 * during multiplications, but it often seems to make things slightly
466 * slower. Hence the default is now to store 32 bits per Long.
467 */
468#endif
469#else /* long long available */
470#ifndef Llong
471#define Llong long long
472#endif
473#ifndef ULLong
474#define ULLong unsigned Llong
475#endif
476#endif /* NO_LONG_LONG */
477
478#define MULTIPLE_THREADS 1
479
480#ifndef MULTIPLE_THREADS
481#define ACQUIRE_DTOA_LOCK(n) /*nothing*/
482#define FREE_DTOA_LOCK(n) /*nothing*/
483#else
484#define ACQUIRE_DTOA_LOCK(n) /*unused right now*/
485#define FREE_DTOA_LOCK(n) /*unused right now*/
486#endif
487
488#ifndef ATOMIC_PTR_CAS
489#define ATOMIC_PTR_CAS(var, old, new) ((var) = (new), (void *)(old))
490#endif
491#ifndef LIKELY
492#define LIKELY(x) (x)
493#endif
494#ifndef UNLIKELY
495#define UNLIKELY(x) (x)
496#endif
497#ifndef ASSUME
498#define ASSUME(x) (void)(x)
499#endif
500
501#define Kmax 15
502
503struct Bigint {
504 struct Bigint *next;
505 int k, maxwds, sign, wds;
506 ULong x[1];
507};
508
509typedef struct Bigint Bigint;
510
511static Bigint *
512Balloc(int k)
513{
514 int x;
515 Bigint *rv;
516
517 x = 1 << k;
518 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
519 rv->k = k;
520 rv->maxwds = x;
521 rv->sign = rv->wds = 0;
522 return rv;
523}
524
525static void
526Bfree(Bigint *v)
527{
528 FREE(v);
529}
530
531#define Bcopy(x,y) memcpy((char *)&(x)->sign, (char *)&(y)->sign, \
532(y)->wds*sizeof(Long) + 2*sizeof(int))
533
534static Bigint *
535multadd(Bigint *b, int m, int a) /* multiply by m and add a */
536{
537 int i, wds;
538 ULong *x;
539#ifdef ULLong
540 ULLong carry, y;
541#else
542 ULong carry, y;
543#ifdef Pack_32
544 ULong xi, z;
545#endif
546#endif
547 Bigint *b1;
548
549 wds = b->wds;
550 x = b->x;
551 i = 0;
552 carry = a;
553 do {
554#ifdef ULLong
555 y = *x * (ULLong)m + carry;
556 carry = y >> 32;
557 *x++ = (ULong)(y & FFFFFFFF);
558#else
559#ifdef Pack_32
560 xi = *x;
561 y = (xi & 0xffff) * m + carry;
562 z = (xi >> 16) * m + (y >> 16);
563 carry = z >> 16;
564 *x++ = (z << 16) + (y & 0xffff);
565#else
566 y = *x * m + carry;
567 carry = y >> 16;
568 *x++ = y & 0xffff;
569#endif
570#endif
571 } while (++i < wds);
572 if (carry) {
573 if (wds >= b->maxwds) {
574 b1 = Balloc(b->k+1);
575 Bcopy(b1, b);
576 Bfree(b);
577 b = b1;
578 }
579 b->x[wds++] = (ULong)carry;
580 b->wds = wds;
581 }
582 return b;
583}
584
585static Bigint *
586s2b(const char *s, int nd0, int nd, ULong y9)
587{
588 Bigint *b;
589 int i, k;
590 Long x, y;
591
592 x = (nd + 8) / 9;
593 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
594#ifdef Pack_32
595 b = Balloc(k);
596 b->x[0] = y9;
597 b->wds = 1;
598#else
599 b = Balloc(k+1);
600 b->x[0] = y9 & 0xffff;
601 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
602#endif
603
604 i = 9;
605 if (9 < nd0) {
606 s += 9;
607 do {
608 b = multadd(b, 10, *s++ - '0');
609 } while (++i < nd0);
610 s++;
611 }
612 else
613 s += 10;
614 for (; i < nd; i++)
615 b = multadd(b, 10, *s++ - '0');
616 return b;
617}
618
619static int
620hi0bits(register ULong x)
621{
622 register int k = 0;
623
624 if (!(x & 0xffff0000)) {
625 k = 16;
626 x <<= 16;
627 }
628 if (!(x & 0xff000000)) {
629 k += 8;
630 x <<= 8;
631 }
632 if (!(x & 0xf0000000)) {
633 k += 4;
634 x <<= 4;
635 }
636 if (!(x & 0xc0000000)) {
637 k += 2;
638 x <<= 2;
639 }
640 if (!(x & 0x80000000)) {
641 k++;
642 if (!(x & 0x40000000))
643 return 32;
644 }
645 return k;
646}
647
648static int
649lo0bits(ULong *y)
650{
651 register int k;
652 register ULong x = *y;
653
654 if (x & 7) {
655 if (x & 1)
656 return 0;
657 if (x & 2) {
658 *y = x >> 1;
659 return 1;
660 }
661 *y = x >> 2;
662 return 2;
663 }
664 k = 0;
665 if (!(x & 0xffff)) {
666 k = 16;
667 x >>= 16;
668 }
669 if (!(x & 0xff)) {
670 k += 8;
671 x >>= 8;
672 }
673 if (!(x & 0xf)) {
674 k += 4;
675 x >>= 4;
676 }
677 if (!(x & 0x3)) {
678 k += 2;
679 x >>= 2;
680 }
681 if (!(x & 1)) {
682 k++;
683 x >>= 1;
684 if (!x)
685 return 32;
686 }
687 *y = x;
688 return k;
689}
690
691static Bigint *
692i2b(int i)
693{
694 Bigint *b;
695
696 b = Balloc(1);
697 b->x[0] = i;
698 b->wds = 1;
699 return b;
700}
701
702static Bigint *
703mult(Bigint *a, Bigint *b)
704{
705 Bigint *c;
706 int k, wa, wb, wc;
707 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
708 ULong y;
709#ifdef ULLong
710 ULLong carry, z;
711#else
712 ULong carry, z;
713#ifdef Pack_32
714 ULong z2;
715#endif
716#endif
717
718 if (a->wds < b->wds) {
719 c = a;
720 a = b;
721 b = c;
722 }
723 k = a->k;
724 wa = a->wds;
725 wb = b->wds;
726 wc = wa + wb;
727 if (wc > a->maxwds)
728 k++;
729 c = Balloc(k);
730 for (x = c->x, xa = x + wc; x < xa; x++)
731 *x = 0;
732 xa = a->x;
733 xae = xa + wa;
734 xb = b->x;
735 xbe = xb + wb;
736 xc0 = c->x;
737#ifdef ULLong
738 for (; xb < xbe; xc0++) {
739 if ((y = *xb++) != 0) {
740 x = xa;
741 xc = xc0;
742 carry = 0;
743 do {
744 z = *x++ * (ULLong)y + *xc + carry;
745 carry = z >> 32;
746 *xc++ = (ULong)(z & FFFFFFFF);
747 } while (x < xae);
748 *xc = (ULong)carry;
749 }
750 }
751#else
752#ifdef Pack_32
753 for (; xb < xbe; xb++, xc0++) {
754 if ((y = *xb & 0xffff) != 0) {
755 x = xa;
756 xc = xc0;
757 carry = 0;
758 do {
759 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
760 carry = z >> 16;
761 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
762 carry = z2 >> 16;
763 Storeinc(xc, z2, z);
764 } while (x < xae);
765 *xc = (ULong)carry;
766 }
767 if ((y = *xb >> 16) != 0) {
768 x = xa;
769 xc = xc0;
770 carry = 0;
771 z2 = *xc;
772 do {
773 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
774 carry = z >> 16;
775 Storeinc(xc, z, z2);
776 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
777 carry = z2 >> 16;
778 } while (x < xae);
779 *xc = z2;
780 }
781 }
782#else
783 for (; xb < xbe; xc0++) {
784 if (y = *xb++) {
785 x = xa;
786 xc = xc0;
787 carry = 0;
788 do {
789 z = *x++ * y + *xc + carry;
790 carry = z >> 16;
791 *xc++ = z & 0xffff;
792 } while (x < xae);
793 *xc = (ULong)carry;
794 }
795 }
796#endif
797#endif
798 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
799 c->wds = wc;
800 return c;
801}
802
803static Bigint *p5s;
804
805static Bigint *
806pow5mult(Bigint *b, int k)
807{
808 Bigint *b1, *p5, *p51;
809 Bigint *p5tmp;
810 int i;
811 static const int p05[3] = { 5, 25, 125 };
812
813 if ((i = k & 3) != 0)
814 b = multadd(b, p05[i-1], 0);
815
816 if (!(k >>= 2))
817 return b;
818 if (!(p5 = p5s)) {
819 /* first time */
820 ACQUIRE_DTOA_LOCK(1);
821 if (!(p5 = p5s)) {
822 p5 = i2b(625);
823 p5->next = 0;
824 p5tmp = ATOMIC_PTR_CAS(p5s, NULL, p5);
825 if (UNLIKELY(p5tmp)) {
826 Bfree(p5);
827 p5 = p5tmp;
828 }
829 }
830 FREE_DTOA_LOCK(1);
831 }
832 for (;;) {
833 if (k & 1) {
834 b1 = mult(b, p5);
835 Bfree(b);
836 b = b1;
837 }
838 if (!(k >>= 1))
839 break;
840 if (!(p51 = p5->next)) {
841 ACQUIRE_DTOA_LOCK(1);
842 if (!(p51 = p5->next)) {
843 p51 = mult(p5,p5);
844 p51->next = 0;
845 p5tmp = ATOMIC_PTR_CAS(p5->next, NULL, p51);
846 if (UNLIKELY(p5tmp)) {
847 Bfree(p51);
848 p51 = p5tmp;
849 }
850 }
851 FREE_DTOA_LOCK(1);
852 }
853 p5 = p51;
854 }
855 return b;
856}
857
858static Bigint *
859lshift(Bigint *b, int k)
860{
861 int i, k1, n, n1;
862 Bigint *b1;
863 ULong *x, *x1, *xe, z;
864
865#ifdef Pack_32
866 n = k >> 5;
867#else
868 n = k >> 4;
869#endif
870 k1 = b->k;
871 n1 = n + b->wds + 1;
872 for (i = b->maxwds; n1 > i; i <<= 1)
873 k1++;
874 b1 = Balloc(k1);
875 x1 = b1->x;
876 for (i = 0; i < n; i++)
877 *x1++ = 0;
878 x = b->x;
879 xe = x + b->wds;
880#ifdef Pack_32
881 if (k &= 0x1f) {
882 k1 = 32 - k;
883 z = 0;
884 do {
885 *x1++ = *x << k | z;
886 z = *x++ >> k1;
887 } while (x < xe);
888 if ((*x1 = z) != 0)
889 ++n1;
890 }
891#else
892 if (k &= 0xf) {
893 k1 = 16 - k;
894 z = 0;
895 do {
896 *x1++ = *x << k & 0xffff | z;
897 z = *x++ >> k1;
898 } while (x < xe);
899 if (*x1 = z)
900 ++n1;
901 }
902#endif
903 else
904 do {
905 *x1++ = *x++;
906 } while (x < xe);
907 b1->wds = n1 - 1;
908 Bfree(b);
909 return b1;
910}
911
912static int
913cmp(Bigint *a, Bigint *b)
914{
915 ULong *xa, *xa0, *xb, *xb0;
916 int i, j;
917
918 i = a->wds;
919 j = b->wds;
920#ifdef DEBUG
921 if (i > 1 && !a->x[i-1])
922 Bug("cmp called with a->x[a->wds-1] == 0");
923 if (j > 1 && !b->x[j-1])
924 Bug("cmp called with b->x[b->wds-1] == 0");
925#endif
926 if (i -= j)
927 return i;
928 xa0 = a->x;
929 xa = xa0 + j;
930 xb0 = b->x;
931 xb = xb0 + j;
932 for (;;) {
933 if (*--xa != *--xb)
934 return *xa < *xb ? -1 : 1;
935 if (xa <= xa0)
936 break;
937 }
938 return 0;
939}
940
941NO_SANITIZE("unsigned-integer-overflow", static Bigint * diff(Bigint *a, Bigint *b));
942static Bigint *
943diff(Bigint *a, Bigint *b)
944{
945 Bigint *c;
946 int i, wa, wb;
947 ULong *xa, *xae, *xb, *xbe, *xc;
948#ifdef ULLong
949 ULLong borrow, y;
950#else
951 ULong borrow, y;
952#ifdef Pack_32
953 ULong z;
954#endif
955#endif
956
957 i = cmp(a,b);
958 if (!i) {
959 c = Balloc(0);
960 c->wds = 1;
961 c->x[0] = 0;
962 return c;
963 }
964 if (i < 0) {
965 c = a;
966 a = b;
967 b = c;
968 i = 1;
969 }
970 else
971 i = 0;
972 c = Balloc(a->k);
973 c->sign = i;
974 wa = a->wds;
975 xa = a->x;
976 xae = xa + wa;
977 wb = b->wds;
978 xb = b->x;
979 xbe = xb + wb;
980 xc = c->x;
981 borrow = 0;
982#ifdef ULLong
983 do {
984 y = (ULLong)*xa++ - *xb++ - borrow;
985 borrow = y >> 32 & (ULong)1;
986 *xc++ = (ULong)(y & FFFFFFFF);
987 } while (xb < xbe);
988 while (xa < xae) {
989 y = *xa++ - borrow;
990 borrow = y >> 32 & (ULong)1;
991 *xc++ = (ULong)(y & FFFFFFFF);
992 }
993#else
994#ifdef Pack_32
995 do {
996 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
997 borrow = (y & 0x10000) >> 16;
998 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
999 borrow = (z & 0x10000) >> 16;
1000 Storeinc(xc, z, y);
1001 } while (xb < xbe);
1002 while (xa < xae) {
1003 y = (*xa & 0xffff) - borrow;
1004 borrow = (y & 0x10000) >> 16;
1005 z = (*xa++ >> 16) - borrow;
1006 borrow = (z & 0x10000) >> 16;
1007 Storeinc(xc, z, y);
1008 }
1009#else
1010 do {
1011 y = *xa++ - *xb++ - borrow;
1012 borrow = (y & 0x10000) >> 16;
1013 *xc++ = y & 0xffff;
1014 } while (xb < xbe);
1015 while (xa < xae) {
1016 y = *xa++ - borrow;
1017 borrow = (y & 0x10000) >> 16;
1018 *xc++ = y & 0xffff;
1019 }
1020#endif
1021#endif
1022 while (!*--xc)
1023 wa--;
1024 c->wds = wa;
1025 return c;
1026}
1027
1028static double
1029ulp(double x_)
1030{
1031 register Long L;
1032 double_u x, a;
1033 dval(x) = x_;
1034
1035 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1036#ifndef Avoid_Underflow
1037#ifndef Sudden_Underflow
1038 if (L > 0) {
1039#endif
1040#endif
1041#ifdef IBM
1042 L |= Exp_msk1 >> 4;
1043#endif
1044 word0(a) = L;
1045 word1(a) = 0;
1046#ifndef Avoid_Underflow
1047#ifndef Sudden_Underflow
1048 }
1049 else {
1050 L = -L >> Exp_shift;
1051 if (L < Exp_shift) {
1052 word0(a) = 0x80000 >> L;
1053 word1(a) = 0;
1054 }
1055 else {
1056 word0(a) = 0;
1057 L -= Exp_shift;
1058 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
1059 }
1060 }
1061#endif
1062#endif
1063 return dval(a);
1064}
1065
1066static double
1067b2d(Bigint *a, int *e)
1068{
1069 ULong *xa, *xa0, w, y, z;
1070 int k;
1071 double_u d;
1072#ifdef VAX
1073 ULong d0, d1;
1074#else
1075#define d0 word0(d)
1076#define d1 word1(d)
1077#endif
1078
1079 xa0 = a->x;
1080 xa = xa0 + a->wds;
1081 y = *--xa;
1082#ifdef DEBUG
1083 if (!y) Bug("zero y in b2d");
1084#endif
1085 k = hi0bits(y);
1086 *e = 32 - k;
1087#ifdef Pack_32
1088 if (k < Ebits) {
1089 d0 = Exp_1 | y >> (Ebits - k);
1090 w = xa > xa0 ? *--xa : 0;
1091 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1092 goto ret_d;
1093 }
1094 z = xa > xa0 ? *--xa : 0;
1095 if (k -= Ebits) {
1096 d0 = Exp_1 | y << k | z >> (32 - k);
1097 y = xa > xa0 ? *--xa : 0;
1098 d1 = z << k | y >> (32 - k);
1099 }
1100 else {
1101 d0 = Exp_1 | y;
1102 d1 = z;
1103 }
1104#else
1105 if (k < Ebits + 16) {
1106 z = xa > xa0 ? *--xa : 0;
1107 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1108 w = xa > xa0 ? *--xa : 0;
1109 y = xa > xa0 ? *--xa : 0;
1110 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1111 goto ret_d;
1112 }
1113 z = xa > xa0 ? *--xa : 0;
1114 w = xa > xa0 ? *--xa : 0;
1115 k -= Ebits + 16;
1116 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1117 y = xa > xa0 ? *--xa : 0;
1118 d1 = w << k + 16 | y << k;
1119#endif
1120ret_d:
1121#ifdef VAX
1122 word0(d) = d0 >> 16 | d0 << 16;
1123 word1(d) = d1 >> 16 | d1 << 16;
1124#else
1125#undef d0
1126#undef d1
1127#endif
1128 return dval(d);
1129}
1130
1131static Bigint *
1132d2b(double d_, int *e, int *bits)
1133{
1134 double_u d;
1135 Bigint *b;
1136 int de, k;
1137 ULong *x, y, z;
1138#ifndef Sudden_Underflow
1139 int i;
1140#endif
1141#ifdef VAX
1142 ULong d0, d1;
1143#endif
1144 dval(d) = d_;
1145#ifdef VAX
1146 d0 = word0(d) >> 16 | word0(d) << 16;
1147 d1 = word1(d) >> 16 | word1(d) << 16;
1148#else
1149#define d0 word0(d)
1150#define d1 word1(d)
1151#endif
1152
1153#ifdef Pack_32
1154 b = Balloc(1);
1155#else
1156 b = Balloc(2);
1157#endif
1158 x = b->x;
1159
1160 z = d0 & Frac_mask;
1161 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1162#ifdef Sudden_Underflow
1163 de = (int)(d0 >> Exp_shift);
1164#ifndef IBM
1165 z |= Exp_msk11;
1166#endif
1167#else
1168 if ((de = (int)(d0 >> Exp_shift)) != 0)
1169 z |= Exp_msk1;
1170#endif
1171#ifdef Pack_32
1172 if ((y = d1) != 0) {
1173 if ((k = lo0bits(&y)) != 0) {
1174 x[0] = y | z << (32 - k);
1175 z >>= k;
1176 }
1177 else
1178 x[0] = y;
1179#ifndef Sudden_Underflow
1180 i =
1181#endif
1182 b->wds = (x[1] = z) ? 2 : 1;
1183 }
1184 else {
1185#ifdef DEBUG
1186 if (!z)
1187 Bug("Zero passed to d2b");
1188#endif
1189 k = lo0bits(&z);
1190 x[0] = z;
1191#ifndef Sudden_Underflow
1192 i =
1193#endif
1194 b->wds = 1;
1195 k += 32;
1196 }
1197#else
1198 if (y = d1) {
1199 if (k = lo0bits(&y))
1200 if (k >= 16) {
1201 x[0] = y | z << 32 - k & 0xffff;
1202 x[1] = z >> k - 16 & 0xffff;
1203 x[2] = z >> k;
1204 i = 2;
1205 }
1206 else {
1207 x[0] = y & 0xffff;
1208 x[1] = y >> 16 | z << 16 - k & 0xffff;
1209 x[2] = z >> k & 0xffff;
1210 x[3] = z >> k+16;
1211 i = 3;
1212 }
1213 else {
1214 x[0] = y & 0xffff;
1215 x[1] = y >> 16;
1216 x[2] = z & 0xffff;
1217 x[3] = z >> 16;
1218 i = 3;
1219 }
1220 }
1221 else {
1222#ifdef DEBUG
1223 if (!z)
1224 Bug("Zero passed to d2b");
1225#endif
1226 k = lo0bits(&z);
1227 if (k >= 16) {
1228 x[0] = z;
1229 i = 0;
1230 }
1231 else {
1232 x[0] = z & 0xffff;
1233 x[1] = z >> 16;
1234 i = 1;
1235 }
1236 k += 32;
1237 }
1238 while (!x[i])
1239 --i;
1240 b->wds = i + 1;
1241#endif
1242#ifndef Sudden_Underflow
1243 if (de) {
1244#endif
1245#ifdef IBM
1246 *e = (de - Bias - (P-1) << 2) + k;
1247 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1248#else
1249 *e = de - Bias - (P-1) + k;
1250 *bits = P - k;
1251#endif
1252#ifndef Sudden_Underflow
1253 }
1254 else {
1255 *e = de - Bias - (P-1) + 1 + k;
1256#ifdef Pack_32
1257 *bits = 32*i - hi0bits(x[i-1]);
1258#else
1259 *bits = (i+2)*16 - hi0bits(x[i]);
1260#endif
1261 }
1262#endif
1263 return b;
1264}
1265#undef d0
1266#undef d1
1267
1268static double
1269ratio(Bigint *a, Bigint *b)
1270{
1271 double_u da, db;
1272 int k, ka, kb;
1273
1274 dval(da) = b2d(a, &ka);
1275 dval(db) = b2d(b, &kb);
1276#ifdef Pack_32
1277 k = ka - kb + 32*(a->wds - b->wds);
1278#else
1279 k = ka - kb + 16*(a->wds - b->wds);
1280#endif
1281#ifdef IBM
1282 if (k > 0) {
1283 word0(da) += (k >> 2)*Exp_msk1;
1284 if (k &= 3)
1285 dval(da) *= 1 << k;
1286 }
1287 else {
1288 k = -k;
1289 word0(db) += (k >> 2)*Exp_msk1;
1290 if (k &= 3)
1291 dval(db) *= 1 << k;
1292 }
1293#else
1294 if (k > 0)
1295 word0(da) += k*Exp_msk1;
1296 else {
1297 k = -k;
1298 word0(db) += k*Exp_msk1;
1299 }
1300#endif
1301 return dval(da) / dval(db);
1302}
1303
1304static const double
1305tens[] = {
1306 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1307 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1308 1e20, 1e21, 1e22
1309#ifdef VAX
1310 , 1e23, 1e24
1311#endif
1312};
1313
1314static const double
1315#ifdef IEEE_Arith
1316bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1317static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1318#ifdef Avoid_Underflow
1319 9007199254740992.*9007199254740992.e-256
1320 /* = 2^106 * 1e-53 */
1321#else
1322 1e-256
1323#endif
1324};
1325/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1326/* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1327#define Scale_Bit 0x10
1328#define n_bigtens 5
1329#else
1330#ifdef IBM
1331bigtens[] = { 1e16, 1e32, 1e64 };
1332static const double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1333#define n_bigtens 3
1334#else
1335bigtens[] = { 1e16, 1e32 };
1336static const double tinytens[] = { 1e-16, 1e-32 };
1337#define n_bigtens 2
1338#endif
1339#endif
1340
1341#ifndef IEEE_Arith
1342#undef INFNAN_CHECK
1343#endif
1344
1345#ifdef INFNAN_CHECK
1346
1347#ifndef NAN_WORD0
1348#define NAN_WORD0 0x7ff80000
1349#endif
1350
1351#ifndef NAN_WORD1
1352#define NAN_WORD1 0
1353#endif
1354
1355static int
1356match(const char **sp, char *t)
1357{
1358 int c, d;
1359 const char *s = *sp;
1360
1361 while (d = *t++) {
1362 if ((c = *++s) >= 'A' && c <= 'Z')
1363 c += 'a' - 'A';
1364 if (c != d)
1365 return 0;
1366 }
1367 *sp = s + 1;
1368 return 1;
1369}
1370
1371#ifndef No_Hex_NaN
1372static void
1373hexnan(double *rvp, const char **sp)
1374{
1375 ULong c, x[2];
1376 const char *s;
1377 int havedig, udx0, xshift;
1378
1379 x[0] = x[1] = 0;
1380 havedig = xshift = 0;
1381 udx0 = 1;
1382 s = *sp;
1383 while (c = *(const unsigned char*)++s) {
1384 if (c >= '0' && c <= '9')
1385 c -= '0';
1386 else if (c >= 'a' && c <= 'f')
1387 c += 10 - 'a';
1388 else if (c >= 'A' && c <= 'F')
1389 c += 10 - 'A';
1390 else if (c <= ' ') {
1391 if (udx0 && havedig) {
1392 udx0 = 0;
1393 xshift = 1;
1394 }
1395 continue;
1396 }
1397 else if (/*(*/ c == ')' && havedig) {
1398 *sp = s + 1;
1399 break;
1400 }
1401 else
1402 return; /* invalid form: don't change *sp */
1403 havedig = 1;
1404 if (xshift) {
1405 xshift = 0;
1406 x[0] = x[1];
1407 x[1] = 0;
1408 }
1409 if (udx0)
1410 x[0] = (x[0] << 4) | (x[1] >> 28);
1411 x[1] = (x[1] << 4) | c;
1412 }
1413 if ((x[0] &= 0xfffff) || x[1]) {
1414 word0(*rvp) = Exp_mask | x[0];
1415 word1(*rvp) = x[1];
1416 }
1417}
1418#endif /*No_Hex_NaN*/
1419#endif /* INFNAN_CHECK */
1420
1421NO_SANITIZE("unsigned-integer-overflow", double strtod(const char *s00, char **se));
1422double
1423strtod(const char *s00, char **se)
1424{
1425#ifdef Avoid_Underflow
1426 int scale;
1427#endif
1428 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1429 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1430 const char *s, *s0, *s1;
1431 double aadj, adj;
1432 double_u aadj1, rv, rv0;
1433 Long L;
1434 ULong y, z;
1435 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1436#ifdef SET_INEXACT
1437 int inexact, oldinexact;
1438#endif
1439#ifdef Honor_FLT_ROUNDS
1440 int rounding;
1441#endif
1442#ifdef USE_LOCALE
1443 const char *s2;
1444#endif
1445
1446 errno = 0;
1447 sign = nz0 = nz = 0;
1448 dval(rv) = 0.;
1449 for (s = s00;;s++)
1450 switch (*s) {
1451 case '-':
1452 sign = 1;
1453 /* no break */
1454 case '+':
1455 if (*++s)
1456 goto break2;
1457 /* no break */
1458 case 0:
1459 goto ret0;
1460 case '\t':
1461 case '\n':
1462 case '\v':
1463 case '\f':
1464 case '\r':
1465 case ' ':
1466 continue;
1467 default:
1468 goto break2;
1469 }
1470break2:
1471 if (*s == '0') {
1472 if (s[1] == 'x' || s[1] == 'X') {
1473 s0 = ++s;
1474 adj = 0;
1475 aadj = 1.0;
1476 nd0 = -4;
1477
1478 if (!*++s || (!(s1 = strchr(hexdigit, *s)) && *s != '.')) goto ret0;
1479 if (*s == '0') {
1480 while (*++s == '0');
1481 if (!*s) goto ret;
1482 s1 = strchr(hexdigit, *s);
1483 }
1484 if (s1 != NULL) {
1485 do {
1486 adj += aadj * ((s1 - hexdigit) & 15);
1487 nd0 += 4;
1488 aadj /= 16;
1489 } while (*++s && (s1 = strchr(hexdigit, *s)));
1490 }
1491
1492 if ((*s == '.') && *++s && (s1 = strchr(hexdigit, *s))) {
1493 if (nd0 < 0) {
1494 while (*s == '0') {
1495 s++;
1496 nd0 -= 4;
1497 }
1498 }
1499 for (; *s && (s1 = strchr(hexdigit, *s)); ++s) {
1500 adj += aadj * ((s1 - hexdigit) & 15);
1501 if ((aadj /= 16) == 0.0) {
1502 while (*++s && strchr(hexdigit, *s));
1503 break;
1504 }
1505 }
1506 }
1507
1508 if (*s == 'P' || *s == 'p') {
1509 dsign = 0x2C - *++s; /* +: 2B, -: 2D */
1510 if (abs(dsign) == 1) s++;
1511 else dsign = 1;
1512
1513 nd = 0;
1514 c = *s;
1515 if (c < '0' || '9' < c) goto ret0;
1516 do {
1517 nd *= 10;
1518 nd += c;
1519 nd -= '0';
1520 c = *++s;
1521 /* Float("0x0."+("0"*267)+"1fp2095") */
1522 if (nd + dsign * nd0 > 2095) {
1523 while ('0' <= c && c <= '9') c = *++s;
1524 break;
1525 }
1526 } while ('0' <= c && c <= '9');
1527 nd0 += nd * dsign;
1528 }
1529 dval(rv) = ldexp(adj, nd0);
1530 goto ret;
1531 }
1532 nz0 = 1;
1533 while (*++s == '0') ;
1534 if (!*s)
1535 goto ret;
1536 }
1537 s0 = s;
1538 y = z = 0;
1539 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1540 if (nd < 9)
1541 y = 10*y + c - '0';
1542 else if (nd < DBL_DIG + 2)
1543 z = 10*z + c - '0';
1544 nd0 = nd;
1545#ifdef USE_LOCALE
1546 s1 = localeconv()->decimal_point;
1547 if (c == *s1) {
1548 c = '.';
1549 if (*++s1) {
1550 s2 = s;
1551 for (;;) {
1552 if (*++s2 != *s1) {
1553 c = 0;
1554 break;
1555 }
1556 if (!*++s1) {
1557 s = s2;
1558 break;
1559 }
1560 }
1561 }
1562 }
1563#endif
1564 if (c == '.') {
1565 c = *++s;
1566 if (!ISDIGIT(c))
1567 goto dig_done;
1568 if (!nd) {
1569 for (; c == '0'; c = *++s)
1570 nz++;
1571 if (c > '0' && c <= '9') {
1572 s0 = s;
1573 nf += nz;
1574 nz = 0;
1575 goto have_dig;
1576 }
1577 goto dig_done;
1578 }
1579 for (; c >= '0' && c <= '9'; c = *++s) {
1580have_dig:
1581 nz++;
1582 if (nd > DBL_DIG * 4) {
1583 continue;
1584 }
1585 if (c -= '0') {
1586 nf += nz;
1587 for (i = 1; i < nz; i++)
1588 if (nd++ < 9)
1589 y *= 10;
1590 else if (nd <= DBL_DIG + 2)
1591 z *= 10;
1592 if (nd++ < 9)
1593 y = 10*y + c;
1594 else if (nd <= DBL_DIG + 2)
1595 z = 10*z + c;
1596 nz = 0;
1597 }
1598 }
1599 }
1600dig_done:
1601 e = 0;
1602 if (c == 'e' || c == 'E') {
1603 if (!nd && !nz && !nz0) {
1604 goto ret0;
1605 }
1606 s00 = s;
1607 esign = 0;
1608 switch (c = *++s) {
1609 case '-':
1610 esign = 1;
1611 case '+':
1612 c = *++s;
1613 }
1614 if (c >= '0' && c <= '9') {
1615 while (c == '0')
1616 c = *++s;
1617 if (c > '0' && c <= '9') {
1618 L = c - '0';
1619 s1 = s;
1620 while ((c = *++s) >= '0' && c <= '9')
1621 L = 10*L + c - '0';
1622 if (s - s1 > 8 || L > 19999)
1623 /* Avoid confusion from exponents
1624 * so large that e might overflow.
1625 */
1626 e = 19999; /* safe for 16 bit ints */
1627 else
1628 e = (int)L;
1629 if (esign)
1630 e = -e;
1631 }
1632 else
1633 e = 0;
1634 }
1635 else
1636 s = s00;
1637 }
1638 if (!nd) {
1639 if (!nz && !nz0) {
1640#ifdef INFNAN_CHECK
1641 /* Check for Nan and Infinity */
1642 switch (c) {
1643 case 'i':
1644 case 'I':
1645 if (match(&s,"nf")) {
1646 --s;
1647 if (!match(&s,"inity"))
1648 ++s;
1649 word0(rv) = 0x7ff00000;
1650 word1(rv) = 0;
1651 goto ret;
1652 }
1653 break;
1654 case 'n':
1655 case 'N':
1656 if (match(&s, "an")) {
1657 word0(rv) = NAN_WORD0;
1658 word1(rv) = NAN_WORD1;
1659#ifndef No_Hex_NaN
1660 if (*s == '(') /*)*/
1661 hexnan(&rv, &s);
1662#endif
1663 goto ret;
1664 }
1665 }
1666#endif /* INFNAN_CHECK */
1667ret0:
1668 s = s00;
1669 sign = 0;
1670 }
1671 goto ret;
1672 }
1673 e1 = e -= nf;
1674
1675 /* Now we have nd0 digits, starting at s0, followed by a
1676 * decimal point, followed by nd-nd0 digits. The number we're
1677 * after is the integer represented by those digits times
1678 * 10**e */
1679
1680 if (!nd0)
1681 nd0 = nd;
1682 k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;
1683 dval(rv) = y;
1684 if (k > 9) {
1685#ifdef SET_INEXACT
1686 if (k > DBL_DIG)
1687 oldinexact = get_inexact();
1688#endif
1689 dval(rv) = tens[k - 9] * dval(rv) + z;
1690 }
1691 bd0 = bb = bd = bs = delta = 0;
1692 if (nd <= DBL_DIG
1693#ifndef RND_PRODQUOT
1694#ifndef Honor_FLT_ROUNDS
1695 && Flt_Rounds == 1
1696#endif
1697#endif
1698 ) {
1699 if (!e)
1700 goto ret;
1701 if (e > 0) {
1702 if (e <= Ten_pmax) {
1703#ifdef VAX
1704 goto vax_ovfl_check;
1705#else
1706#ifdef Honor_FLT_ROUNDS
1707 /* round correctly FLT_ROUNDS = 2 or 3 */
1708 if (sign) {
1709 dval(rv) = -dval(rv);
1710 sign = 0;
1711 }
1712#endif
1713 /* rv = */ rounded_product(dval(rv), tens[e]);
1714 goto ret;
1715#endif
1716 }
1717 i = DBL_DIG - nd;
1718 if (e <= Ten_pmax + i) {
1719 /* A fancier test would sometimes let us do
1720 * this for larger i values.
1721 */
1722#ifdef Honor_FLT_ROUNDS
1723 /* round correctly FLT_ROUNDS = 2 or 3 */
1724 if (sign) {
1725 dval(rv) = -dval(rv);
1726 sign = 0;
1727 }
1728#endif
1729 e -= i;
1730 dval(rv) *= tens[i];
1731#ifdef VAX
1732 /* VAX exponent range is so narrow we must
1733 * worry about overflow here...
1734 */
1735vax_ovfl_check:
1736 word0(rv) -= P*Exp_msk1;
1737 /* rv = */ rounded_product(dval(rv), tens[e]);
1738 if ((word0(rv) & Exp_mask)
1739 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1740 goto ovfl;
1741 word0(rv) += P*Exp_msk1;
1742#else
1743 /* rv = */ rounded_product(dval(rv), tens[e]);
1744#endif
1745 goto ret;
1746 }
1747 }
1748#ifndef Inaccurate_Divide
1749 else if (e >= -Ten_pmax) {
1750#ifdef Honor_FLT_ROUNDS
1751 /* round correctly FLT_ROUNDS = 2 or 3 */
1752 if (sign) {
1753 dval(rv) = -dval(rv);
1754 sign = 0;
1755 }
1756#endif
1757 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1758 goto ret;
1759 }
1760#endif
1761 }
1762 e1 += nd - k;
1763
1764#ifdef IEEE_Arith
1765#ifdef SET_INEXACT
1766 inexact = 1;
1767 if (k <= DBL_DIG)
1768 oldinexact = get_inexact();
1769#endif
1770#ifdef Avoid_Underflow
1771 scale = 0;
1772#endif
1773#ifdef Honor_FLT_ROUNDS
1774 if ((rounding = Flt_Rounds) >= 2) {
1775 if (sign)
1776 rounding = rounding == 2 ? 0 : 2;
1777 else
1778 if (rounding != 2)
1779 rounding = 0;
1780 }
1781#endif
1782#endif /*IEEE_Arith*/
1783
1784 /* Get starting approximation = rv * 10**e1 */
1785
1786 if (e1 > 0) {
1787 if ((i = e1 & 15) != 0)
1788 dval(rv) *= tens[i];
1789 if (e1 &= ~15) {
1790 if (e1 > DBL_MAX_10_EXP) {
1791ovfl:
1792#ifndef NO_ERRNO
1793 errno = ERANGE;
1794#endif
1795 /* Can't trust HUGE_VAL */
1796#ifdef IEEE_Arith
1797#ifdef Honor_FLT_ROUNDS
1798 switch (rounding) {
1799 case 0: /* toward 0 */
1800 case 3: /* toward -infinity */
1801 word0(rv) = Big0;
1802 word1(rv) = Big1;
1803 break;
1804 default:
1805 word0(rv) = Exp_mask;
1806 word1(rv) = 0;
1807 }
1808#else /*Honor_FLT_ROUNDS*/
1809 word0(rv) = Exp_mask;
1810 word1(rv) = 0;
1811#endif /*Honor_FLT_ROUNDS*/
1812#ifdef SET_INEXACT
1813 /* set overflow bit */
1814 dval(rv0) = 1e300;
1815 dval(rv0) *= dval(rv0);
1816#endif
1817#else /*IEEE_Arith*/
1818 word0(rv) = Big0;
1819 word1(rv) = Big1;
1820#endif /*IEEE_Arith*/
1821 if (bd0)
1822 goto retfree;
1823 goto ret;
1824 }
1825 e1 >>= 4;
1826 for (j = 0; e1 > 1; j++, e1 >>= 1)
1827 if (e1 & 1)
1828 dval(rv) *= bigtens[j];
1829 /* The last multiplication could overflow. */
1830 word0(rv) -= P*Exp_msk1;
1831 dval(rv) *= bigtens[j];
1832 if ((z = word0(rv) & Exp_mask)
1833 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1834 goto ovfl;
1835 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1836 /* set to largest number */
1837 /* (Can't trust DBL_MAX) */
1838 word0(rv) = Big0;
1839 word1(rv) = Big1;
1840 }
1841 else
1842 word0(rv) += P*Exp_msk1;
1843 }
1844 }
1845 else if (e1 < 0) {
1846 e1 = -e1;
1847 if ((i = e1 & 15) != 0)
1848 dval(rv) /= tens[i];
1849 if (e1 >>= 4) {
1850 if (e1 >= 1 << n_bigtens)
1851 goto undfl;
1852#ifdef Avoid_Underflow
1853 if (e1 & Scale_Bit)
1854 scale = 2*P;
1855 for (j = 0; e1 > 0; j++, e1 >>= 1)
1856 if (e1 & 1)
1857 dval(rv) *= tinytens[j];
1858 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
1859 >> Exp_shift)) > 0) {
1860 /* scaled rv is denormal; zap j low bits */
1861 if (j >= 32) {
1862 word1(rv) = 0;
1863 if (j >= 53)
1864 word0(rv) = (P+2)*Exp_msk1;
1865 else
1866 word0(rv) &= 0xffffffff << (j-32);
1867 }
1868 else
1869 word1(rv) &= 0xffffffff << j;
1870 }
1871#else
1872 for (j = 0; e1 > 1; j++, e1 >>= 1)
1873 if (e1 & 1)
1874 dval(rv) *= tinytens[j];
1875 /* The last multiplication could underflow. */
1876 dval(rv0) = dval(rv);
1877 dval(rv) *= tinytens[j];
1878 if (!dval(rv)) {
1879 dval(rv) = 2.*dval(rv0);
1880 dval(rv) *= tinytens[j];
1881#endif
1882 if (!dval(rv)) {
1883undfl:
1884 dval(rv) = 0.;
1885#ifndef NO_ERRNO
1886 errno = ERANGE;
1887#endif
1888 if (bd0)
1889 goto retfree;
1890 goto ret;
1891 }
1892#ifndef Avoid_Underflow
1893 word0(rv) = Tiny0;
1894 word1(rv) = Tiny1;
1895 /* The refinement below will clean
1896 * this approximation up.
1897 */
1898 }
1899#endif
1900 }
1901 }
1902
1903 /* Now the hard part -- adjusting rv to the correct value.*/
1904
1905 /* Put digits into bd: true value = bd * 10^e */
1906
1907 bd0 = s2b(s0, nd0, nd, y);
1908
1909 for (;;) {
1910 bd = Balloc(bd0->k);
1911 Bcopy(bd, bd0);
1912 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
1913 bs = i2b(1);
1914
1915 if (e >= 0) {
1916 bb2 = bb5 = 0;
1917 bd2 = bd5 = e;
1918 }
1919 else {
1920 bb2 = bb5 = -e;
1921 bd2 = bd5 = 0;
1922 }
1923 if (bbe >= 0)
1924 bb2 += bbe;
1925 else
1926 bd2 -= bbe;
1927 bs2 = bb2;
1928#ifdef Honor_FLT_ROUNDS
1929 if (rounding != 1)
1930 bs2++;
1931#endif
1932#ifdef Avoid_Underflow
1933 j = bbe - scale;
1934 i = j + bbbits - 1; /* logb(rv) */
1935 if (i < Emin) /* denormal */
1936 j += P - Emin;
1937 else
1938 j = P + 1 - bbbits;
1939#else /*Avoid_Underflow*/
1940#ifdef Sudden_Underflow
1941#ifdef IBM
1942 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1943#else
1944 j = P + 1 - bbbits;
1945#endif
1946#else /*Sudden_Underflow*/
1947 j = bbe;
1948 i = j + bbbits - 1; /* logb(rv) */
1949 if (i < Emin) /* denormal */
1950 j += P - Emin;
1951 else
1952 j = P + 1 - bbbits;
1953#endif /*Sudden_Underflow*/
1954#endif /*Avoid_Underflow*/
1955 bb2 += j;
1956 bd2 += j;
1957#ifdef Avoid_Underflow
1958 bd2 += scale;
1959#endif
1960 i = bb2 < bd2 ? bb2 : bd2;
1961 if (i > bs2)
1962 i = bs2;
1963 if (i > 0) {
1964 bb2 -= i;
1965 bd2 -= i;
1966 bs2 -= i;
1967 }
1968 if (bb5 > 0) {
1969 bs = pow5mult(bs, bb5);
1970 bb1 = mult(bs, bb);
1971 Bfree(bb);
1972 bb = bb1;
1973 }
1974 if (bb2 > 0)
1975 bb = lshift(bb, bb2);
1976 if (bd5 > 0)
1977 bd = pow5mult(bd, bd5);
1978 if (bd2 > 0)
1979 bd = lshift(bd, bd2);
1980 if (bs2 > 0)
1981 bs = lshift(bs, bs2);
1982 delta = diff(bb, bd);
1983 dsign = delta->sign;
1984 delta->sign = 0;
1985 i = cmp(delta, bs);
1986#ifdef Honor_FLT_ROUNDS
1987 if (rounding != 1) {
1988 if (i < 0) {
1989 /* Error is less than an ulp */
1990 if (!delta->x[0] && delta->wds <= 1) {
1991 /* exact */
1992#ifdef SET_INEXACT
1993 inexact = 0;
1994#endif
1995 break;
1996 }
1997 if (rounding) {
1998 if (dsign) {
1999 adj = 1.;
2000 goto apply_adj;
2001 }
2002 }
2003 else if (!dsign) {
2004 adj = -1.;
2005 if (!word1(rv)
2006 && !(word0(rv) & Frac_mask)) {
2007 y = word0(rv) & Exp_mask;
2008#ifdef Avoid_Underflow
2009 if (!scale || y > 2*P*Exp_msk1)
2010#else
2011 if (y)
2012#endif
2013 {
2014 delta = lshift(delta,Log2P);
2015 if (cmp(delta, bs) <= 0)
2016 adj = -0.5;
2017 }
2018 }
2019apply_adj:
2020#ifdef Avoid_Underflow
2021 if (scale && (y = word0(rv) & Exp_mask)
2022 <= 2*P*Exp_msk1)
2023 word0(adj) += (2*P+1)*Exp_msk1 - y;
2024#else
2025#ifdef Sudden_Underflow
2026 if ((word0(rv) & Exp_mask) <=
2027 P*Exp_msk1) {
2028 word0(rv) += P*Exp_msk1;
2029 dval(rv) += adj*ulp(dval(rv));
2030 word0(rv) -= P*Exp_msk1;
2031 }
2032 else
2033#endif /*Sudden_Underflow*/
2034#endif /*Avoid_Underflow*/
2035 dval(rv) += adj*ulp(dval(rv));
2036 }
2037 break;
2038 }
2039 adj = ratio(delta, bs);
2040 if (adj < 1.)
2041 adj = 1.;
2042 if (adj <= 0x7ffffffe) {
2043 /* adj = rounding ? ceil(adj) : floor(adj); */
2044 y = adj;
2045 if (y != adj) {
2046 if (!((rounding>>1) ^ dsign))
2047 y++;
2048 adj = y;
2049 }
2050 }
2051#ifdef Avoid_Underflow
2052 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2053 word0(adj) += (2*P+1)*Exp_msk1 - y;
2054#else
2055#ifdef Sudden_Underflow
2056 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2057 word0(rv) += P*Exp_msk1;
2058 adj *= ulp(dval(rv));
2059 if (dsign)
2060 dval(rv) += adj;
2061 else
2062 dval(rv) -= adj;
2063 word0(rv) -= P*Exp_msk1;
2064 goto cont;
2065 }
2066#endif /*Sudden_Underflow*/
2067#endif /*Avoid_Underflow*/
2068 adj *= ulp(dval(rv));
2069 if (dsign)
2070 dval(rv) += adj;
2071 else
2072 dval(rv) -= adj;
2073 goto cont;
2074 }
2075#endif /*Honor_FLT_ROUNDS*/
2076
2077 if (i < 0) {
2078 /* Error is less than half an ulp -- check for
2079 * special case of mantissa a power of two.
2080 */
2081 if (dsign || word1(rv) || word0(rv) & Bndry_mask
2082#ifdef IEEE_Arith
2083#ifdef Avoid_Underflow
2084 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2085#else
2086 || (word0(rv) & Exp_mask) <= Exp_msk1
2087#endif
2088#endif
2089 ) {
2090#ifdef SET_INEXACT
2091 if (!delta->x[0] && delta->wds <= 1)
2092 inexact = 0;
2093#endif
2094 break;
2095 }
2096 if (!delta->x[0] && delta->wds <= 1) {
2097 /* exact result */
2098#ifdef SET_INEXACT
2099 inexact = 0;
2100#endif
2101 break;
2102 }
2103 delta = lshift(delta,Log2P);
2104 if (cmp(delta, bs) > 0)
2105 goto drop_down;
2106 break;
2107 }
2108 if (i == 0) {
2109 /* exactly half-way between */
2110 if (dsign) {
2111 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2112 && word1(rv) == (
2113#ifdef Avoid_Underflow
2114 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2115 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2116#endif
2117 0xffffffff)) {
2118 /*boundary case -- increment exponent*/
2119 word0(rv) = (word0(rv) & Exp_mask)
2120 + Exp_msk1
2121#ifdef IBM
2122 | Exp_msk1 >> 4
2123#endif
2124 ;
2125 word1(rv) = 0;
2126#ifdef Avoid_Underflow
2127 dsign = 0;
2128#endif
2129 break;
2130 }
2131 }
2132 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2133drop_down:
2134 /* boundary case -- decrement exponent */
2135#ifdef Sudden_Underflow /*{{*/
2136 L = word0(rv) & Exp_mask;
2137#ifdef IBM
2138 if (L < Exp_msk1)
2139#else
2140#ifdef Avoid_Underflow
2141 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
2142#else
2143 if (L <= Exp_msk1)
2144#endif /*Avoid_Underflow*/
2145#endif /*IBM*/
2146 goto undfl;
2147 L -= Exp_msk1;
2148#else /*Sudden_Underflow}{*/
2149#ifdef Avoid_Underflow
2150 if (scale) {
2151 L = word0(rv) & Exp_mask;
2152 if (L <= (2*P+1)*Exp_msk1) {
2153 if (L > (P+2)*Exp_msk1)
2154 /* round even ==> */
2155 /* accept rv */
2156 break;
2157 /* rv = smallest denormal */
2158 goto undfl;
2159 }
2160 }
2161#endif /*Avoid_Underflow*/
2162 L = (word0(rv) & Exp_mask) - Exp_msk1;
2163#endif /*Sudden_Underflow}}*/
2164 word0(rv) = L | Bndry_mask1;
2165 word1(rv) = 0xffffffff;
2166#ifdef IBM
2167 goto cont;
2168#else
2169 break;
2170#endif
2171 }
2172#ifndef ROUND_BIASED
2173 if (!(word1(rv) & LSB))
2174 break;
2175#endif
2176 if (dsign)
2177 dval(rv) += ulp(dval(rv));
2178#ifndef ROUND_BIASED
2179 else {
2180 dval(rv) -= ulp(dval(rv));
2181#ifndef Sudden_Underflow
2182 if (!dval(rv))
2183 goto undfl;
2184#endif
2185 }
2186#ifdef Avoid_Underflow
2187 dsign = 1 - dsign;
2188#endif
2189#endif
2190 break;
2191 }
2192 if ((aadj = ratio(delta, bs)) <= 2.) {
2193 if (dsign)
2194 aadj = dval(aadj1) = 1.;
2195 else if (word1(rv) || word0(rv) & Bndry_mask) {
2196#ifndef Sudden_Underflow
2197 if (word1(rv) == Tiny1 && !word0(rv))
2198 goto undfl;
2199#endif
2200 aadj = 1.;
2201 dval(aadj1) = -1.;
2202 }
2203 else {
2204 /* special case -- power of FLT_RADIX to be */
2205 /* rounded down... */
2206
2207 if (aadj < 2./FLT_RADIX)
2208 aadj = 1./FLT_RADIX;
2209 else
2210 aadj *= 0.5;
2211 dval(aadj1) = -aadj;
2212 }
2213 }
2214 else {
2215 aadj *= 0.5;
2216 dval(aadj1) = dsign ? aadj : -aadj;
2217#ifdef Check_FLT_ROUNDS
2218 switch (Rounding) {
2219 case 2: /* towards +infinity */
2220 dval(aadj1) -= 0.5;
2221 break;
2222 case 0: /* towards 0 */
2223 case 3: /* towards -infinity */
2224 dval(aadj1) += 0.5;
2225 }
2226#else
2227 if (Flt_Rounds == 0)
2228 dval(aadj1) += 0.5;
2229#endif /*Check_FLT_ROUNDS*/
2230 }
2231 y = word0(rv) & Exp_mask;
2232
2233 /* Check for overflow */
2234
2235 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2236 dval(rv0) = dval(rv);
2237 word0(rv) -= P*Exp_msk1;
2238 adj = dval(aadj1) * ulp(dval(rv));
2239 dval(rv) += adj;
2240 if ((word0(rv) & Exp_mask) >=
2241 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2242 if (word0(rv0) == Big0 && word1(rv0) == Big1)
2243 goto ovfl;
2244 word0(rv) = Big0;
2245 word1(rv) = Big1;
2246 goto cont;
2247 }
2248 else
2249 word0(rv) += P*Exp_msk1;
2250 }
2251 else {
2252#ifdef Avoid_Underflow
2253 if (scale && y <= 2*P*Exp_msk1) {
2254 if (aadj <= 0x7fffffff) {
2255 if ((z = (int)aadj) <= 0)
2256 z = 1;
2257 aadj = z;
2258 dval(aadj1) = dsign ? aadj : -aadj;
2259 }
2260 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
2261 }
2262 adj = dval(aadj1) * ulp(dval(rv));
2263 dval(rv) += adj;
2264#else
2265#ifdef Sudden_Underflow
2266 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2267 dval(rv0) = dval(rv);
2268 word0(rv) += P*Exp_msk1;
2269 adj = dval(aadj1) * ulp(dval(rv));
2270 dval(rv) += adj;
2271#ifdef IBM
2272 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
2273#else
2274 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2275#endif
2276 {
2277 if (word0(rv0) == Tiny0 && word1(rv0) == Tiny1)
2278 goto undfl;
2279 word0(rv) = Tiny0;
2280 word1(rv) = Tiny1;
2281 goto cont;
2282 }
2283 else
2284 word0(rv) -= P*Exp_msk1;
2285 }
2286 else {
2287 adj = dval(aadj1) * ulp(dval(rv));
2288 dval(rv) += adj;
2289 }
2290#else /*Sudden_Underflow*/
2291 /* Compute adj so that the IEEE rounding rules will
2292 * correctly round rv + adj in some half-way cases.
2293 * If rv * ulp(rv) is denormalized (i.e.,
2294 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2295 * trouble from bits lost to denormalization;
2296 * example: 1.2e-307 .
2297 */
2298 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
2299 dval(aadj1) = (double)(int)(aadj + 0.5);
2300 if (!dsign)
2301 dval(aadj1) = -dval(aadj1);
2302 }
2303 adj = dval(aadj1) * ulp(dval(rv));
2304 dval(rv) += adj;
2305#endif /*Sudden_Underflow*/
2306#endif /*Avoid_Underflow*/
2307 }
2308 z = word0(rv) & Exp_mask;
2309#ifndef SET_INEXACT
2310#ifdef Avoid_Underflow
2311 if (!scale)
2312#endif
2313 if (y == z) {
2314 /* Can we stop now? */
2315 L = (Long)aadj;
2316 aadj -= L;
2317 /* The tolerances below are conservative. */
2318 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2319 if (aadj < .4999999 || aadj > .5000001)
2320 break;
2321 }
2322 else if (aadj < .4999999/FLT_RADIX)
2323 break;
2324 }
2325#endif
2326cont:
2327 Bfree(bb);
2328 Bfree(bd);
2329 Bfree(bs);
2330 Bfree(delta);
2331 }
2332#ifdef SET_INEXACT
2333 if (inexact) {
2334 if (!oldinexact) {
2335 word0(rv0) = Exp_1 + (70 << Exp_shift);
2336 word1(rv0) = 0;
2337 dval(rv0) += 1.;
2338 }
2339 }
2340 else if (!oldinexact)
2341 clear_inexact();
2342#endif
2343#ifdef Avoid_Underflow
2344 if (scale) {
2345 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
2346 word1(rv0) = 0;
2347 dval(rv) *= dval(rv0);
2348#ifndef NO_ERRNO
2349 /* try to avoid the bug of testing an 8087 register value */
2350 if (word0(rv) == 0 && word1(rv) == 0)
2351 errno = ERANGE;
2352#endif
2353 }
2354#endif /* Avoid_Underflow */
2355#ifdef SET_INEXACT
2356 if (inexact && !(word0(rv) & Exp_mask)) {
2357 /* set underflow bit */
2358 dval(rv0) = 1e-300;
2359 dval(rv0) *= dval(rv0);
2360 }
2361#endif
2362retfree:
2363 Bfree(bb);
2364 Bfree(bd);
2365 Bfree(bs);
2366 Bfree(bd0);
2367 Bfree(delta);
2368ret:
2369 if (se)
2370 *se = (char *)s;
2371 return sign ? -dval(rv) : dval(rv);
2372}
2373
2374NO_SANITIZE("unsigned-integer-overflow", static int quorem(Bigint *b, Bigint *S));
2375static int
2376quorem(Bigint *b, Bigint *S)
2377{
2378 int n;
2379 ULong *bx, *bxe, q, *sx, *sxe;
2380#ifdef ULLong
2381 ULLong borrow, carry, y, ys;
2382#else
2383 ULong borrow, carry, y, ys;
2384#ifdef Pack_32
2385 ULong si, z, zs;
2386#endif
2387#endif
2388
2389 n = S->wds;
2390#ifdef DEBUG
2391 /*debug*/ if (b->wds > n)
2392 /*debug*/ Bug("oversize b in quorem");
2393#endif
2394 if (b->wds < n)
2395 return 0;
2396 sx = S->x;
2397 sxe = sx + --n;
2398 bx = b->x;
2399 bxe = bx + n;
2400 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
2401#ifdef DEBUG
2402 /*debug*/ if (q > 9)
2403 /*debug*/ Bug("oversized quotient in quorem");
2404#endif
2405 if (q) {
2406 borrow = 0;
2407 carry = 0;
2408 do {
2409#ifdef ULLong
2410 ys = *sx++ * (ULLong)q + carry;
2411 carry = ys >> 32;
2412 y = *bx - (ys & FFFFFFFF) - borrow;
2413 borrow = y >> 32 & (ULong)1;
2414 *bx++ = (ULong)(y & FFFFFFFF);
2415#else
2416#ifdef Pack_32
2417 si = *sx++;
2418 ys = (si & 0xffff) * q + carry;
2419 zs = (si >> 16) * q + (ys >> 16);
2420 carry = zs >> 16;
2421 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2422 borrow = (y & 0x10000) >> 16;
2423 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2424 borrow = (z & 0x10000) >> 16;
2425 Storeinc(bx, z, y);
2426#else
2427 ys = *sx++ * q + carry;
2428 carry = ys >> 16;
2429 y = *bx - (ys & 0xffff) - borrow;
2430 borrow = (y & 0x10000) >> 16;
2431 *bx++ = y & 0xffff;
2432#endif
2433#endif
2434 } while (sx <= sxe);
2435 if (!*bxe) {
2436 bx = b->x;
2437 while (--bxe > bx && !*bxe)
2438 --n;
2439 b->wds = n;
2440 }
2441 }
2442 if (cmp(b, S) >= 0) {
2443 q++;
2444 borrow = 0;
2445 carry = 0;
2446 bx = b->x;
2447 sx = S->x;
2448 do {
2449#ifdef ULLong
2450 ys = *sx++ + carry;
2451 carry = ys >> 32;
2452 y = *bx - (ys & FFFFFFFF) - borrow;
2453 borrow = y >> 32 & (ULong)1;
2454 *bx++ = (ULong)(y & FFFFFFFF);
2455#else
2456#ifdef Pack_32
2457 si = *sx++;
2458 ys = (si & 0xffff) + carry;
2459 zs = (si >> 16) + (ys >> 16);
2460 carry = zs >> 16;
2461 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2462 borrow = (y & 0x10000) >> 16;
2463 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2464 borrow = (z & 0x10000) >> 16;
2465 Storeinc(bx, z, y);
2466#else
2467 ys = *sx++ + carry;
2468 carry = ys >> 16;
2469 y = *bx - (ys & 0xffff) - borrow;
2470 borrow = (y & 0x10000) >> 16;
2471 *bx++ = y & 0xffff;
2472#endif
2473#endif
2474 } while (sx <= sxe);
2475 bx = b->x;
2476 bxe = bx + n;
2477 if (!*bxe) {
2478 while (--bxe > bx && !*bxe)
2479 --n;
2480 b->wds = n;
2481 }
2482 }
2483 return q;
2484}
2485
2486#ifndef MULTIPLE_THREADS
2487static char *dtoa_result;
2488#endif
2489
2490#ifndef MULTIPLE_THREADS
2491static char *
2492rv_alloc(int i)
2493{
2494 return dtoa_result = MALLOC(i);
2495}
2496#else
2497#define rv_alloc(i) MALLOC(i)
2498#endif
2499
2500static char *
2501nrv_alloc(const char *s, char **rve, size_t n)
2502{
2503 char *rv, *t;
2504
2505 t = rv = rv_alloc(n);
2506 while ((*t = *s++) != 0) t++;
2507 if (rve)
2508 *rve = t;
2509 return rv;
2510}
2511
2512#define rv_strdup(s, rve) nrv_alloc((s), (rve), strlen(s)+1)
2513
2514#ifndef MULTIPLE_THREADS
2515/* freedtoa(s) must be used to free values s returned by dtoa
2516 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
2517 * but for consistency with earlier versions of dtoa, it is optional
2518 * when MULTIPLE_THREADS is not defined.
2519 */
2520
2521static void
2522freedtoa(char *s)
2523{
2524 FREE(s);
2525}
2526#endif
2527
2528static const char INFSTR[] = "Infinity";
2529static const char NANSTR[] = "NaN";
2530static const char ZEROSTR[] = "0";
2531
2532/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2533 *
2534 * Inspired by "How to Print Floating-Point Numbers Accurately" by
2535 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
2536 *
2537 * Modifications:
2538 * 1. Rather than iterating, we use a simple numeric overestimate
2539 * to determine k = floor(log10(d)). We scale relevant
2540 * quantities using O(log2(k)) rather than O(k) multiplications.
2541 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2542 * try to generate digits strictly left to right. Instead, we
2543 * compute with fewer bits and propagate the carry if necessary
2544 * when rounding the final digit up. This is often faster.
2545 * 3. Under the assumption that input will be rounded nearest,
2546 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2547 * That is, we allow equality in stopping tests when the
2548 * round-nearest rule will give the same floating-point value
2549 * as would satisfaction of the stopping test with strict
2550 * inequality.
2551 * 4. We remove common factors of powers of 2 from relevant
2552 * quantities.
2553 * 5. When converting floating-point integers less than 1e16,
2554 * we use floating-point arithmetic rather than resorting
2555 * to multiple-precision integers.
2556 * 6. When asked to produce fewer than 15 digits, we first try
2557 * to get by with floating-point arithmetic; we resort to
2558 * multiple-precision integer arithmetic only if we cannot
2559 * guarantee that the floating-point calculation has given
2560 * the correctly rounded result. For k requested digits and
2561 * "uniformly" distributed input, the probability is
2562 * something like 10^(k-15) that we must resort to the Long
2563 * calculation.
2564 */
2565
2566char *
2567dtoa(double d_, int mode, int ndigits, int *decpt, int *sign, char **rve)
2568{
2569 /* Arguments ndigits, decpt, sign are similar to those
2570 of ecvt and fcvt; trailing zeros are suppressed from
2571 the returned string. If not null, *rve is set to point
2572 to the end of the return value. If d is +-Infinity or NaN,
2573 then *decpt is set to 9999.
2574
2575 mode:
2576 0 ==> shortest string that yields d when read in
2577 and rounded to nearest.
2578 1 ==> like 0, but with Steele & White stopping rule;
2579 e.g. with IEEE P754 arithmetic , mode 0 gives
2580 1e23 whereas mode 1 gives 9.999999999999999e22.
2581 2 ==> max(1,ndigits) significant digits. This gives a
2582 return value similar to that of ecvt, except
2583 that trailing zeros are suppressed.
2584 3 ==> through ndigits past the decimal point. This
2585 gives a return value similar to that from fcvt,
2586 except that trailing zeros are suppressed, and
2587 ndigits can be negative.
2588 4,5 ==> similar to 2 and 3, respectively, but (in
2589 round-nearest mode) with the tests of mode 0 to
2590 possibly return a shorter string that rounds to d.
2591 With IEEE arithmetic and compilation with
2592 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2593 as modes 2 and 3 when FLT_ROUNDS != 1.
2594 6-9 ==> Debugging modes similar to mode - 4: don't try
2595 fast floating-point estimate (if applicable).
2596
2597 Values of mode other than 0-9 are treated as mode 0.
2598
2599 Sufficient space is allocated to the return value
2600 to hold the suppressed trailing zeros.
2601 */
2602
2603 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2604 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2605 spec_case, try_quick, half = 0;
2606 Long L;
2607#ifndef Sudden_Underflow
2608 int denorm;
2609 ULong x;
2610#endif
2611 Bigint *b, *b1, *delta, *mlo = 0, *mhi = 0, *S;
2612 double ds;
2613 double_u d, d2, eps;
2614 char *s, *s0;
2615#ifdef Honor_FLT_ROUNDS
2616 int rounding;
2617#endif
2618#ifdef SET_INEXACT
2619 int inexact, oldinexact;
2620#endif
2621
2622 dval(d) = d_;
2623
2624#ifndef MULTIPLE_THREADS
2625 if (dtoa_result) {
2626 freedtoa(dtoa_result);
2627 dtoa_result = 0;
2628 }
2629#endif
2630
2631 if (word0(d) & Sign_bit) {
2632 /* set sign for everything, including 0's and NaNs */
2633 *sign = 1;
2634 word0(d) &= ~Sign_bit; /* clear sign bit */
2635 }
2636 else
2637 *sign = 0;
2638
2639#if defined(IEEE_Arith) + defined(VAX)
2640#ifdef IEEE_Arith
2641 if ((word0(d) & Exp_mask) == Exp_mask)
2642#else
2643 if (word0(d) == 0x8000)
2644#endif
2645 {
2646 /* Infinity or NaN */
2647 *decpt = 9999;
2648#ifdef IEEE_Arith
2649 if (!word1(d) && !(word0(d) & 0xfffff))
2650 return rv_strdup(INFSTR, rve);
2651#endif
2652 return rv_strdup(NANSTR, rve);
2653 }
2654#endif
2655#ifdef IBM
2656 dval(d) += 0; /* normalize */
2657#endif
2658 if (!dval(d)) {
2659 *decpt = 1;
2660 return rv_strdup(ZEROSTR, rve);
2661 }
2662
2663#ifdef SET_INEXACT
2664 try_quick = oldinexact = get_inexact();
2665 inexact = 1;
2666#endif
2667#ifdef Honor_FLT_ROUNDS
2668 if ((rounding = Flt_Rounds) >= 2) {
2669 if (*sign)
2670 rounding = rounding == 2 ? 0 : 2;
2671 else
2672 if (rounding != 2)
2673 rounding = 0;
2674 }
2675#endif
2676
2677 b = d2b(dval(d), &be, &bbits);
2678#ifdef Sudden_Underflow
2679 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
2680#else
2681 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) {
2682#endif
2683 dval(d2) = dval(d);
2684 word0(d2) &= Frac_mask1;
2685 word0(d2) |= Exp_11;
2686#ifdef IBM
2687 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
2688 dval(d2) /= 1 << j;
2689#endif
2690
2691 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2692 * log10(x) = log(x) / log(10)
2693 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2694 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2695 *
2696 * This suggests computing an approximation k to log10(d) by
2697 *
2698 * k = (i - Bias)*0.301029995663981
2699 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2700 *
2701 * We want k to be too large rather than too small.
2702 * The error in the first-order Taylor series approximation
2703 * is in our favor, so we just round up the constant enough
2704 * to compensate for any error in the multiplication of
2705 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2706 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2707 * adding 1e-13 to the constant term more than suffices.
2708 * Hence we adjust the constant term to 0.1760912590558.
2709 * (We could get a more accurate k by invoking log10,
2710 * but this is probably not worthwhile.)
2711 */
2712
2713 i -= Bias;
2714#ifdef IBM
2715 i <<= 2;
2716 i += j;
2717#endif
2718#ifndef Sudden_Underflow
2719 denorm = 0;
2720 }
2721 else {
2722 /* d is denormalized */
2723
2724 i = bbits + be + (Bias + (P-1) - 1);
2725 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
2726 : word1(d) << (32 - i);
2727 dval(d2) = x;
2728 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2729 i -= (Bias + (P-1) - 1) + 1;
2730 denorm = 1;
2731 }
2732#endif
2733 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2734 k = (int)ds;
2735 if (ds < 0. && ds != k)
2736 k--; /* want k = floor(ds) */
2737 k_check = 1;
2738 if (k >= 0 && k <= Ten_pmax) {
2739 if (dval(d) < tens[k])
2740 k--;
2741 k_check = 0;
2742 }
2743 j = bbits - i - 1;
2744 if (j >= 0) {
2745 b2 = 0;
2746 s2 = j;
2747 }
2748 else {
2749 b2 = -j;
2750 s2 = 0;
2751 }
2752 if (k >= 0) {
2753 b5 = 0;
2754 s5 = k;
2755 s2 += k;
2756 }
2757 else {
2758 b2 -= k;
2759 b5 = -k;
2760 s5 = 0;
2761 }
2762 if (mode < 0 || mode > 9)
2763 mode = 0;
2764
2765#ifndef SET_INEXACT
2766#ifdef Check_FLT_ROUNDS
2767 try_quick = Rounding == 1;
2768#else
2769 try_quick = 1;
2770#endif
2771#endif /*SET_INEXACT*/
2772
2773 if (mode > 5) {
2774 mode -= 4;
2775 try_quick = 0;
2776 }
2777 leftright = 1;
2778 ilim = ilim1 = -1;
2779 switch (mode) {
2780 case 0:
2781 case 1:
2782 i = 18;
2783 ndigits = 0;
2784 break;
2785 case 2:
2786 leftright = 0;
2787 /* no break */
2788 case 4:
2789 if (ndigits <= 0)
2790 ndigits = 1;
2791 ilim = ilim1 = i = ndigits;
2792 break;
2793 case 3:
2794 leftright = 0;
2795 /* no break */
2796 case 5:
2797 i = ndigits + k + 1;
2798 ilim = i;
2799 ilim1 = i - 1;
2800 if (i <= 0)
2801 i = 1;
2802 }
2803 s = s0 = rv_alloc(i+1);
2804
2805#ifdef Honor_FLT_ROUNDS
2806 if (mode > 1 && rounding != 1)
2807 leftright = 0;
2808#endif
2809
2810 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2811
2812 /* Try to get by with floating-point arithmetic. */
2813
2814 i = 0;
2815 dval(d2) = dval(d);
2816 k0 = k;
2817 ilim0 = ilim;
2818 ieps = 2; /* conservative */
2819 if (k > 0) {
2820 ds = tens[k&0xf];
2821 j = k >> 4;
2822 if (j & Bletch) {
2823 /* prevent overflows */
2824 j &= Bletch - 1;
2825 dval(d) /= bigtens[n_bigtens-1];
2826 ieps++;
2827 }
2828 for (; j; j >>= 1, i++)
2829 if (j & 1) {
2830 ieps++;
2831 ds *= bigtens[i];
2832 }
2833 dval(d) /= ds;
2834 }
2835 else if ((j1 = -k) != 0) {
2836 dval(d) *= tens[j1 & 0xf];
2837 for (j = j1 >> 4; j; j >>= 1, i++)
2838 if (j & 1) {
2839 ieps++;
2840 dval(d) *= bigtens[i];
2841 }
2842 }
2843 if (k_check && dval(d) < 1. && ilim > 0) {
2844 if (ilim1 <= 0)
2845 goto fast_failed;
2846 ilim = ilim1;
2847 k--;
2848 dval(d) *= 10.;
2849 ieps++;
2850 }
2851 dval(eps) = ieps*dval(d) + 7.;
2852 word0(eps) -= (P-1)*Exp_msk1;
2853 if (ilim == 0) {
2854 S = mhi = 0;
2855 dval(d) -= 5.;
2856 if (dval(d) > dval(eps))
2857 goto one_digit;
2858 if (dval(d) < -dval(eps))
2859 goto no_digits;
2860 goto fast_failed;
2861 }
2862#ifndef No_leftright
2863 if (leftright) {
2864 /* Use Steele & White method of only
2865 * generating digits needed.
2866 */
2867 dval(eps) = 0.5/tens[ilim-1] - dval(eps);
2868 for (i = 0;;) {
2869 L = (int)dval(d);
2870 dval(d) -= L;
2871 *s++ = '0' + (int)L;
2872 if (dval(d) < dval(eps))
2873 goto ret1;
2874 if (1. - dval(d) < dval(eps))
2875 goto bump_up;
2876 if (++i >= ilim)
2877 break;
2878 dval(eps) *= 10.;
2879 dval(d) *= 10.;
2880 }
2881 }
2882 else {
2883#endif
2884 /* Generate ilim digits, then fix them up. */
2885 dval(eps) *= tens[ilim-1];
2886 for (i = 1;; i++, dval(d) *= 10.) {
2887 L = (Long)(dval(d));
2888 if (!(dval(d) -= L))
2889 ilim = i;
2890 *s++ = '0' + (int)L;
2891 if (i == ilim) {
2892 if (dval(d) > 0.5 + dval(eps))
2893 goto bump_up;
2894 else if (dval(d) < 0.5 - dval(eps)) {
2895 while (*--s == '0') ;
2896 s++;
2897 goto ret1;
2898 }
2899 half = 1;
2900 if ((*(s-1) - '0') & 1) {
2901 goto bump_up;
2902 }
2903 break;
2904 }
2905 }
2906#ifndef No_leftright
2907 }
2908#endif
2909fast_failed:
2910 s = s0;
2911 dval(d) = dval(d2);
2912 k = k0;
2913 ilim = ilim0;
2914 }
2915
2916 /* Do we have a "small" integer? */
2917
2918 if (be >= 0 && k <= Int_max) {
2919 /* Yes. */
2920 ds = tens[k];
2921 if (ndigits < 0 && ilim <= 0) {
2922 S = mhi = 0;
2923 if (ilim < 0 || dval(d) <= 5*ds)
2924 goto no_digits;
2925 goto one_digit;
2926 }
2927 for (i = 1;; i++, dval(d) *= 10.) {
2928 L = (Long)(dval(d) / ds);
2929 dval(d) -= L*ds;
2930#ifdef Check_FLT_ROUNDS
2931 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2932 if (dval(d) < 0) {
2933 L--;
2934 dval(d) += ds;
2935 }
2936#endif
2937 *s++ = '0' + (int)L;
2938 if (!dval(d)) {
2939#ifdef SET_INEXACT
2940 inexact = 0;
2941#endif
2942 break;
2943 }
2944 if (i == ilim) {
2945#ifdef Honor_FLT_ROUNDS
2946 if (mode > 1)
2947 switch (rounding) {
2948 case 0: goto ret1;
2949 case 2: goto bump_up;
2950 }
2951#endif
2952 dval(d) += dval(d);
2953 if (dval(d) > ds || (dval(d) == ds && (L & 1))) {
2954bump_up:
2955 while (*--s == '9')
2956 if (s == s0) {
2957 k++;
2958 *s = '0';
2959 break;
2960 }
2961 ++*s++;
2962 }
2963 break;
2964 }
2965 }
2966 goto ret1;
2967 }
2968
2969 m2 = b2;
2970 m5 = b5;
2971 if (leftright) {
2972 i =
2973#ifndef Sudden_Underflow
2974 denorm ? be + (Bias + (P-1) - 1 + 1) :
2975#endif
2976#ifdef IBM
2977 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2978#else
2979 1 + P - bbits;
2980#endif
2981 b2 += i;
2982 s2 += i;
2983 mhi = i2b(1);
2984 }
2985 if (m2 > 0 && s2 > 0) {
2986 i = m2 < s2 ? m2 : s2;
2987 b2 -= i;
2988 m2 -= i;
2989 s2 -= i;
2990 }
2991 if (b5 > 0) {
2992 if (leftright) {
2993 if (m5 > 0) {
2994 mhi = pow5mult(mhi, m5);
2995 b1 = mult(mhi, b);
2996 Bfree(b);
2997 b = b1;
2998 }
2999 if ((j = b5 - m5) != 0)
3000 b = pow5mult(b, j);
3001 }
3002 else
3003 b = pow5mult(b, b5);
3004 }
3005 S = i2b(1);
3006 if (s5 > 0)
3007 S = pow5mult(S, s5);
3008
3009 /* Check for special case that d is a normalized power of 2. */
3010
3011 spec_case = 0;
3012 if ((mode < 2 || leftright)
3013#ifdef Honor_FLT_ROUNDS
3014 && rounding == 1
3015#endif
3016 ) {
3017 if (!word1(d) && !(word0(d) & Bndry_mask)
3018#ifndef Sudden_Underflow
3019 && word0(d) & (Exp_mask & ~Exp_msk1)
3020#endif
3021 ) {
3022 /* The special case */
3023 b2 += Log2P;
3024 s2 += Log2P;
3025 spec_case = 1;
3026 }
3027 }
3028
3029 /* Arrange for convenient computation of quotients:
3030 * shift left if necessary so divisor has 4 leading 0 bits.
3031 *
3032 * Perhaps we should just compute leading 28 bits of S once
3033 * and for all and pass them and a shift to quorem, so it
3034 * can do shifts and ors to compute the numerator for q.
3035 */
3036#ifdef Pack_32
3037 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0)
3038 i = 32 - i;
3039#else
3040 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) != 0)
3041 i = 16 - i;
3042#endif
3043 if (i > 4) {
3044 i -= 4;
3045 b2 += i;
3046 m2 += i;
3047 s2 += i;
3048 }
3049 else if (i < 4) {
3050 i += 28;
3051 b2 += i;
3052 m2 += i;
3053 s2 += i;
3054 }
3055 if (b2 > 0)
3056 b = lshift(b, b2);
3057 if (s2 > 0)
3058 S = lshift(S, s2);
3059 if (k_check) {
3060 if (cmp(b,S) < 0) {
3061 k--;
3062 b = multadd(b, 10, 0); /* we botched the k estimate */
3063 if (leftright)
3064 mhi = multadd(mhi, 10, 0);
3065 ilim = ilim1;
3066 }
3067 }
3068 if (ilim <= 0 && (mode == 3 || mode == 5)) {
3069 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
3070 /* no digits, fcvt style */
3071no_digits:
3072 k = -1 - ndigits;
3073 goto ret;
3074 }
3075one_digit:
3076 *s++ = '1';
3077 k++;
3078 goto ret;
3079 }
3080 if (leftright) {
3081 if (m2 > 0)
3082 mhi = lshift(mhi, m2);
3083
3084 /* Compute mlo -- check for special case
3085 * that d is a normalized power of 2.
3086 */
3087
3088 mlo = mhi;
3089 if (spec_case) {
3090 mhi = Balloc(mhi->k);
3091 Bcopy(mhi, mlo);
3092 mhi = lshift(mhi, Log2P);
3093 }
3094
3095 for (i = 1;;i++) {
3096 dig = quorem(b,S) + '0';
3097 /* Do we yet have the shortest decimal string
3098 * that will round to d?
3099 */
3100 j = cmp(b, mlo);
3101 delta = diff(S, mhi);
3102 j1 = delta->sign ? 1 : cmp(b, delta);
3103 Bfree(delta);
3104#ifndef ROUND_BIASED
3105 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
3106#ifdef Honor_FLT_ROUNDS
3107 && rounding >= 1
3108#endif
3109 ) {
3110 if (dig == '9')
3111 goto round_9_up;
3112 if (j > 0)
3113 dig++;
3114#ifdef SET_INEXACT
3115 else if (!b->x[0] && b->wds <= 1)
3116 inexact = 0;
3117#endif
3118 *s++ = dig;
3119 goto ret;
3120 }
3121#endif
3122 if (j < 0 || (j == 0 && mode != 1
3123#ifndef ROUND_BIASED
3124 && !(word1(d) & 1)
3125#endif
3126 )) {
3127 if (!b->x[0] && b->wds <= 1) {
3128#ifdef SET_INEXACT
3129 inexact = 0;
3130#endif
3131 goto accept_dig;
3132 }
3133#ifdef Honor_FLT_ROUNDS
3134 if (mode > 1)
3135 switch (rounding) {
3136 case 0: goto accept_dig;
3137 case 2: goto keep_dig;
3138 }
3139#endif /*Honor_FLT_ROUNDS*/
3140 if (j1 > 0) {
3141 b = lshift(b, 1);
3142 j1 = cmp(b, S);
3143 if ((j1 > 0 || (j1 == 0 && (dig & 1))) && dig++ == '9')
3144 goto round_9_up;
3145 }
3146accept_dig:
3147 *s++ = dig;
3148 goto ret;
3149 }
3150 if (j1 > 0) {
3151#ifdef Honor_FLT_ROUNDS
3152 if (!rounding)
3153 goto accept_dig;
3154#endif
3155 if (dig == '9') { /* possible if i == 1 */
3156round_9_up:
3157 *s++ = '9';
3158 goto roundoff;
3159 }
3160 *s++ = dig + 1;
3161 goto ret;
3162 }
3163#ifdef Honor_FLT_ROUNDS
3164keep_dig:
3165#endif
3166 *s++ = dig;
3167 if (i == ilim)
3168 break;
3169 b = multadd(b, 10, 0);
3170 if (mlo == mhi)
3171 mlo = mhi = multadd(mhi, 10, 0);
3172 else {
3173 mlo = multadd(mlo, 10, 0);
3174 mhi = multadd(mhi, 10, 0);
3175 }
3176 }
3177 }
3178 else
3179 for (i = 1;; i++) {
3180 *s++ = dig = quorem(b,S) + '0';
3181 if (!b->x[0] && b->wds <= 1) {
3182#ifdef SET_INEXACT
3183 inexact = 0;
3184#endif
3185 goto ret;
3186 }
3187 if (i >= ilim)
3188 break;
3189 b = multadd(b, 10, 0);
3190 }
3191
3192 /* Round off last digit */
3193
3194#ifdef Honor_FLT_ROUNDS
3195 switch (rounding) {
3196 case 0: goto trimzeros;
3197 case 2: goto roundoff;
3198 }
3199#endif
3200 b = lshift(b, 1);
3201 j = cmp(b, S);
3202 if (j > 0 || (j == 0 && (dig & 1))) {
3203 roundoff:
3204 while (*--s == '9')
3205 if (s == s0) {
3206 k++;
3207 *s++ = '1';
3208 goto ret;
3209 }
3210 if (!half || (*s - '0') & 1)
3211 ++*s;
3212 }
3213 else {
3214 while (*--s == '0') ;
3215 }
3216 s++;
3217ret:
3218 Bfree(S);
3219 if (mhi) {
3220 if (mlo && mlo != mhi)
3221 Bfree(mlo);
3222 Bfree(mhi);
3223 }
3224ret1:
3225#ifdef SET_INEXACT
3226 if (inexact) {
3227 if (!oldinexact) {
3228 word0(d) = Exp_1 + (70 << Exp_shift);
3229 word1(d) = 0;
3230 dval(d) += 1.;
3231 }
3232 }
3233 else if (!oldinexact)
3234 clear_inexact();
3235#endif
3236 Bfree(b);
3237 *s = 0;
3238 *decpt = k + 1;
3239 if (rve)
3240 *rve = s;
3241 return s0;
3242}
3243
3244/*-
3245 * Copyright (c) 2004-2008 David Schultz <das@FreeBSD.ORG>
3246 * All rights reserved.
3247 *
3248 * Redistribution and use in source and binary forms, with or without
3249 * modification, are permitted provided that the following conditions
3250 * are met:
3251 * 1. Redistributions of source code must retain the above copyright
3252 * notice, this list of conditions and the following disclaimer.
3253 * 2. Redistributions in binary form must reproduce the above copyright
3254 * notice, this list of conditions and the following disclaimer in the
3255 * documentation and/or other materials provided with the distribution.
3256 *
3257 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
3258 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
3259 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
3260 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
3261 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
3262 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
3263 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
3264 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
3265 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
3266 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
3267 * SUCH DAMAGE.
3268 */
3269
3270#define DBL_MANH_SIZE 20
3271#define DBL_MANL_SIZE 32
3272#define DBL_ADJ (DBL_MAX_EXP - 2)
3273#define SIGFIGS ((DBL_MANT_DIG + 3) / 4 + 1)
3274#define dexp_get(u) ((int)(word0(u) >> Exp_shift) & ~Exp_msk1)
3275#define dexp_set(u,v) (word0(u) = (((int)(word0(u)) & ~Exp_mask) | ((v) << Exp_shift)))
3276#define dmanh_get(u) ((uint32_t)(word0(u) & Frac_mask))
3277#define dmanl_get(u) ((uint32_t)word1(u))
3278
3279
3280/*
3281 * This procedure converts a double-precision number in IEEE format
3282 * into a string of hexadecimal digits and an exponent of 2. Its
3283 * behavior is bug-for-bug compatible with dtoa() in mode 2, with the
3284 * following exceptions:
3285 *
3286 * - An ndigits < 0 causes it to use as many digits as necessary to
3287 * represent the number exactly.
3288 * - The additional xdigs argument should point to either the string
3289 * "0123456789ABCDEF" or the string "0123456789abcdef", depending on
3290 * which case is desired.
3291 * - This routine does not repeat dtoa's mistake of setting decpt
3292 * to 9999 in the case of an infinity or NaN. INT_MAX is used
3293 * for this purpose instead.
3294 *
3295 * Note that the C99 standard does not specify what the leading digit
3296 * should be for non-zero numbers. For instance, 0x1.3p3 is the same
3297 * as 0x2.6p2 is the same as 0x4.cp3. This implementation always makes
3298 * the leading digit a 1. This ensures that the exponent printed is the
3299 * actual base-2 exponent, i.e., ilogb(d).
3300 *
3301 * Inputs: d, xdigs, ndigits
3302 * Outputs: decpt, sign, rve
3303 */
3304char *
3305hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign, char **rve)
3306{
3307 U u;
3308 char *s, *s0;
3309 int bufsize;
3310 uint32_t manh, manl;
3311
3312 u.d = d;
3313 if (word0(u) & Sign_bit) {
3314 /* set sign for everything, including 0's and NaNs */
3315 *sign = 1;
3316 word0(u) &= ~Sign_bit; /* clear sign bit */
3317 }
3318 else
3319 *sign = 0;
3320
3321 if (isinf(d)) { /* FP_INFINITE */
3322 *decpt = INT_MAX;
3323 return rv_strdup(INFSTR, rve);
3324 }
3325 else if (isnan(d)) { /* FP_NAN */
3326 *decpt = INT_MAX;
3327 return rv_strdup(NANSTR, rve);
3328 }
3329 else if (d == 0.0) { /* FP_ZERO */
3330 *decpt = 1;
3331 return rv_strdup(ZEROSTR, rve);
3332 }
3333 else if (dexp_get(u)) { /* FP_NORMAL */
3334 *decpt = dexp_get(u) - DBL_ADJ;
3335 }
3336 else { /* FP_SUBNORMAL */
3337 u.d *= 5.363123171977039e+154 /* 0x1p514 */;
3338 *decpt = dexp_get(u) - (514 + DBL_ADJ);
3339 }
3340
3341 if (ndigits == 0) /* dtoa() compatibility */
3342 ndigits = 1;
3343
3344 /*
3345 * If ndigits < 0, we are expected to auto-size, so we allocate
3346 * enough space for all the digits.
3347 */
3348 bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
3349 s0 = rv_alloc(bufsize+1);
3350
3351 /* Round to the desired number of digits. */
3352 if (SIGFIGS > ndigits && ndigits > 0) {
3353 float redux = 1.0f;
3354 int offset = 4 * ndigits + DBL_MAX_EXP - 4 - DBL_MANT_DIG;
3355 dexp_set(u, offset);
3356 u.d += redux;
3357 u.d -= redux;
3358 *decpt += dexp_get(u) - offset;
3359 }
3360
3361 manh = dmanh_get(u);
3362 manl = dmanl_get(u);
3363 *s0 = '1';
3364 for (s = s0 + 1; s < s0 + bufsize; s++) {
3365 *s = xdigs[(manh >> (DBL_MANH_SIZE - 4)) & 0xf];
3366 manh = (manh << 4) | (manl >> (DBL_MANL_SIZE - 4));
3367 manl <<= 4;
3368 }
3369
3370 /* If ndigits < 0, we are expected to auto-size the precision. */
3371 if (ndigits < 0) {
3372 for (ndigits = SIGFIGS; s0[ndigits - 1] == '0'; ndigits--)
3373 ;
3374 }
3375
3376 s = s0 + ndigits;
3377 *s = '\0';
3378 if (rve != NULL)
3379 *rve = s;
3380 return (s0);
3381}
3382
3383#ifdef __cplusplus
3384#if 0
3385{ /* satisfy cc-mode */
3386#endif
3387}
3388#endif
#define ISDIGIT
Old name of rb_isdigit.
Definition ctype.h:93
#define strtod(s, e)
Just another name of ruby_strtod.
Definition util.h:223
#define errno
Ractor-aware version of errno.
Definition ruby.h:388
Definition dtoa.c:503
Definition dtoa.c:286