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