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