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_8087 for IEEE-arithmetic machines where the least
66 * significant byte has the lowest address.
67 * #define IEEE_MC68k 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. Unless Trust_FLT_ROUNDS
76 * is also #defined, fegetround() will be queried for the rounding mode.
77 * Note that both FLT_ROUNDS and fegetround() are specified by the C99
78 * standard (and are specified to be consistent, with fesetround()
79 * affecting the value of FLT_ROUNDS), but that some (Linux) systems
80 * do not work correctly in this regard, so using fegetround() is more
81 * portable than using FLT_FOUNDS directly.
82 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
83 * and Honor_FLT_ROUNDS is not #defined.
84 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
85 * that use extended-precision instructions to compute rounded
86 * products and quotients) with IBM.
87 * #define ROUND_BIASED for IEEE-format with biased rounding.
88 * #define Inaccurate_Divide for IEEE-format with correctly rounded
89 * products but inaccurate quotients, e.g., for Intel i860.
90 * #define NO_LONG_LONG on machines that do not have a "long long"
91 * integer type (of >= 64 bits). On such machines, you can
92 * #define Just_16 to store 16 bits per 32-bit Long when doing
93 * high-precision integer arithmetic. Whether this speeds things
94 * up or slows things down depends on the machine and the number
95 * being converted. If long long is available and the name is
96 * something other than "long long", #define Llong to be the name,
97 * and if "unsigned Llong" does not work as an unsigned version of
98 * Llong, #define #ULLong to be the corresponding unsigned type.
99 * #define KR_headers for old-style C function headers.
100 * #define Bad_float_h if your system lacks a float.h or if it does not
101 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
102 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
103 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
104 * if memory is available and otherwise does something you deem
105 * appropriate. If MALLOC is undefined, malloc will be invoked
106 * directly -- and assumed always to succeed. Similarly, if you
107 * want something other than the system's free() to be called to
108 * recycle memory acquired from MALLOC, #define FREE to be the
109 * name of the alternate routine. (FREE or free is only called in
110 * pathological cases, e.g., in a dtoa call after a dtoa return in
111 * mode 3 with thousands of digits requested.)
112 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
113 * memory allocations from a private pool of memory when possible.
114 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
115 * unless #defined to be a different length. This default length
116 * suffices to get rid of MALLOC calls except for unusual cases,
117 * such as decimal-to-binary conversion of a very long string of
118 * digits. The longest string dtoa can return is about 751 bytes
119 * long. For conversions by strtod of strings of 800 digits and
120 * all dtoa conversions in single-threaded executions with 8-byte
121 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
122 * pointers, PRIVATE_MEM >= 7112 appears adequate.
123 * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
124 * #defined automatically on IEEE systems. On such systems,
125 * when INFNAN_CHECK is #defined, strtod checks
126 * for Infinity and NaN (case insensitively). On some systems
127 * (e.g., some HP systems), it may be necessary to #define NAN_WORD0
128 * appropriately -- to the most significant word of a quiet NaN.
129 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
130 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
131 * strtod also accepts (case insensitively) strings of the form
132 * NaN(x), where x is a string of hexadecimal digits and spaces;
133 * if there is only one string of hexadecimal digits, it is taken
134 * for the 52 fraction bits of the resulting NaN; if there are two
135 * or more strings of hex digits, the first is for the high 20 bits,
136 * the second and subsequent for the low 32 bits, with intervening
137 * white space ignored; but if this results in none of the 52
138 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
139 * and NAN_WORD1 are used instead.
140 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
141 * multiple threads. In this case, you must provide (or suitably
142 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
143 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
144 * in pow5mult, ensures lazy evaluation of only one copy of high
145 * powers of 5; omitting this lock would introduce a small
146 * probability of wasting memory, but would otherwise be harmless.)
147 * You must also invoke freedtoa(s) to free the value s returned by
148 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
149 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
150 * avoids underflows on inputs whose result does not underflow.
151 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
152 * floating-point numbers and flushes underflows to zero rather
153 * than implementing gradual underflow, then you must also #define
154 * Sudden_Underflow.
155 * #define USE_LOCALE to use the current locale's decimal_point value.
156 * #define SET_INEXACT if IEEE arithmetic is being used and extra
157 * computation should be done to set the inexact flag when the
158 * result is inexact and avoid setting inexact when the result
159 * is exact. In this case, dtoa.c must be compiled in
160 * an environment, perhaps provided by #include "dtoa.c" in a
161 * suitable wrapper, that defines two functions,
162 * int get_inexact(void);
163 * void clear_inexact(void);
164 * such that get_inexact() returns a nonzero value if the
165 * inexact bit is already set, and clear_inexact() sets the
166 * inexact bit to 0. When SET_INEXACT is #defined, strtod
167 * also does extra computations to set the underflow and overflow
168 * flags when appropriate (i.e., when the result is tiny and
169 * inexact or when it is a numeric value rounded to +-infinity).
170 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
171 * the result overflows to +-Infinity or underflows to 0.
172 * #define NO_HEX_FP to omit recognition of hexadecimal floating-point
173 * values by strtod.
174 * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now)
175 * to disable logic for "fast" testing of very long input strings
176 * to strtod. This testing proceeds by initially truncating the
177 * input string, then if necessary comparing the whole string with
178 * a decimal expansion to decide close cases. This logic is only
179 * used for input more than STRTOD_DIGLIM digits long (default 40).
180 */
181
182 #if defined _MSC_VER && _MSC_VER == 1800
183 // TODO(scottmg): VS2013 RC ICEs on a bunch of functions in this file.
184 // This should be removed after RTM. See http://crbug.com/288948.
185 #pragma optimize("", off)
186 #endif
187
188 #define IEEE_8087
189 #define NO_HEX_FP
190
191 #ifndef Long
192 #if __LP64__
193 #define Long int
194 #else
195 #define Long long
196 #endif
197 #endif
198 #ifndef ULong
199 typedef unsigned Long ULong;
200 #endif
201
202 #ifdef DEBUG
203 #include "stdio.h"
204 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
205 #endif
206
207 #include "stdlib.h"
208 #include "string.h"
209
210 #ifdef USE_LOCALE
211 #include "locale.h"
212 #endif
213
214 #ifdef Honor_FLT_ROUNDS
215 #ifndef Trust_FLT_ROUNDS
216 #include <fenv.h>
217 #endif
218 #endif
219
220 #ifdef MALLOC
221 #ifdef KR_headers
222 extern char *MALLOC();
223 #else
224 extern void *MALLOC(size_t);
225 #endif
226 #else
227 #define MALLOC malloc
228 #endif
229
230 #ifndef Omit_Private_Memory
231 #ifndef PRIVATE_MEM
232 #define PRIVATE_MEM 2304
233 #endif
234 #define PRIVATE_mem ((unsigned)((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)))
235 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
236 #endif
237
238 #undef IEEE_Arith
239 #undef Avoid_Underflow
240 #ifdef IEEE_MC68k
241 #define IEEE_Arith
242 #endif
243 #ifdef IEEE_8087
244 #define IEEE_Arith
245 #endif
246
247 #ifdef IEEE_Arith
248 #ifndef NO_INFNAN_CHECK
249 #undef INFNAN_CHECK
250 #define INFNAN_CHECK
251 #endif
252 #else
253 #undef INFNAN_CHECK
254 #define NO_STRTOD_BIGCOMP
255 #endif
256
257 #include "errno.h"
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 #ifndef __MATH_H__
293 #include "math.h"
294 #endif
295
296 namespace dmg_fp {
297
298 #ifndef CONST
299 #ifdef KR_headers
300 #define CONST /* blank */
301 #else
302 #define CONST const
303 #endif
304 #endif
305
306 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
307 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
308 #endif
309
310 typedef union { double d; ULong L[2]; } U;
311
312 #ifdef IEEE_8087
313 #define word0(x) (x)->L[1]
314 #define word1(x) (x)->L[0]
315 #else
316 #define word0(x) (x)->L[0]
317 #define word1(x) (x)->L[1]
318 #endif
319 #define dval(x) (x)->d
320
321 #ifndef STRTOD_DIGLIM
322 #define STRTOD_DIGLIM 40
323 #endif
324
325 #ifdef DIGLIM_DEBUG
326 extern int strtod_diglim;
327 #else
328 #define strtod_diglim STRTOD_DIGLIM
329 #endif
330
331 /* The following definition of Storeinc is appropriate for MIPS processors.
332 * An alternative that might be better on some machines is
333 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
334 */
335 #if defined(IEEE_8087) + defined(VAX)
336 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
337 ((unsigned short *)a)[0] = (unsigned short)c, a++)
338 #else
339 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
340 ((unsigned short *)a)[1] = (unsigned short)c, a++)
341 #endif
342
343 /* #define P DBL_MANT_DIG */
344 /* Ten_pmax = floor(P*log(2)/log(5)) */
345 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
346 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
347 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
348
349 #ifdef IEEE_Arith
350 #define Exp_shift 20
351 #define Exp_shift1 20
352 #define Exp_msk1 0x100000
353 #define Exp_msk11 0x100000
354 #define Exp_mask 0x7ff00000
355 #define P 53
356 #define Nbits 53
357 #define Bias 1023
358 #define Emax 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 #undef Check_FLT_ROUNDS
393 #define Check_FLT_ROUNDS
394 #else
395 #define Rounding Flt_Rounds
396 #endif
397
398 #else /* ifndef IEEE_Arith */
399 #undef Check_FLT_ROUNDS
400 #undef Honor_FLT_ROUNDS
401 #undef SET_INEXACT
402 #undef Sudden_Underflow
403 #define Sudden_Underflow
404 #ifdef IBM
405 #undef Flt_Rounds
406 #define Flt_Rounds 0
407 #define Exp_shift 24
408 #define Exp_shift1 24
409 #define Exp_msk1 0x1000000
410 #define Exp_msk11 0x1000000
411 #define Exp_mask 0x7f000000
412 #define P 14
413 #define Nbits 56
414 #define Bias 65
415 #define Emax 248
416 #define Emin (-260)
417 #define Exp_1 0x41000000
418 #define Exp_11 0x41000000
419 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
420 #define Frac_mask 0xffffff
421 #define Frac_mask1 0xffffff
422 #define Bletch 4
423 #define Ten_pmax 22
424 #define Bndry_mask 0xefffff
425 #define Bndry_mask1 0xffffff
426 #define LSB 1
427 #define Sign_bit 0x80000000
428 #define Log2P 4
429 #define Tiny0 0x100000
430 #define Tiny1 0
431 #define Quick_max 14
432 #define Int_max 15
433 #else /* VAX */
434 #undef Flt_Rounds
435 #define Flt_Rounds 1
436 #define Exp_shift 23
437 #define Exp_shift1 7
438 #define Exp_msk1 0x80
439 #define Exp_msk11 0x800000
440 #define Exp_mask 0x7f80
441 #define P 56
442 #define Nbits 56
443 #define Bias 129
444 #define Emax 126
445 #define Emin (-129)
446 #define Exp_1 0x40800000
447 #define Exp_11 0x4080
448 #define Ebits 8
449 #define Frac_mask 0x7fffff
450 #define Frac_mask1 0xffff007f
451 #define Ten_pmax 24
452 #define Bletch 2
453 #define Bndry_mask 0xffff007f
454 #define Bndry_mask1 0xffff007f
455 #define LSB 0x10000
456 #define Sign_bit 0x8000
457 #define Log2P 1
458 #define Tiny0 0x80
459 #define Tiny1 0
460 #define Quick_max 15
461 #define Int_max 15
462 #endif /* IBM, VAX */
463 #endif /* IEEE_Arith */
464
465 #ifndef IEEE_Arith
466 #define ROUND_BIASED
467 #endif
468
469 #ifdef RND_PRODQUOT
470 #define rounded_product(a,b) a = rnd_prod(a, b)
471 #define rounded_quotient(a,b) a = rnd_quot(a, b)
472 #ifdef KR_headers
473 extern double rnd_prod(), rnd_quot();
474 #else
475 extern double rnd_prod(double, double), rnd_quot(double, double);
476 #endif
477 #else
478 #define rounded_product(a,b) a *= b
479 #define rounded_quotient(a,b) a /= b
480 #endif
481
482 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
483 #define Big1 0xffffffff
484
485 #ifndef Pack_32
486 #define Pack_32
487 #endif
488
489 typedef struct BCinfo BCinfo;
490 struct
491 BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflchk; };
492
493 #ifdef KR_headers
494 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
495 #else
496 #define FFFFFFFF 0xffffffffUL
497 #endif
498
499 #ifdef NO_LONG_LONG
500 #undef ULLong
501 #ifdef Just_16
502 #undef Pack_32
503 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
504 * This makes some inner loops simpler and sometimes saves work
505 * during multiplications, but it often seems to make things slightly
506 * slower. Hence the default is now to store 32 bits per Long.
507 */
508 #endif
509 #else /* long long available */
510 #ifndef Llong
511 #define Llong long long
512 #endif
513 #ifndef ULLong
514 #define ULLong unsigned Llong
515 #endif
516 #endif /* NO_LONG_LONG */
517
518 #ifndef MULTIPLE_THREADS
519 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
520 #define FREE_DTOA_LOCK(n) /*nothing*/
521 #endif
522
523 #define Kmax 7
524
525 double strtod(const char *s00, char **se);
526 char *dtoa(double d, int mode, int ndigits,
527 int *decpt, int *sign, char **rve);
528
529 struct
530 Bigint {
531 struct Bigint *next;
532 int k, maxwds, sign, wds;
533 ULong x[1];
534 };
535
536 typedef struct Bigint Bigint;
537
538 static Bigint *freelist[Kmax+1];
539
540 static Bigint *
541 Balloc
542 #ifdef KR_headers
543 (k) int k;
544 #else
545 (int k)
546 #endif
547 {
548 int x;
549 Bigint *rv;
550 #ifndef Omit_Private_Memory
551 unsigned int len;
552 #endif
553
554 ACQUIRE_DTOA_LOCK(0);
555 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
556 /* but this case seems very unlikely. */
557 if (k <= Kmax && (rv = freelist[k]))
558 freelist[k] = rv->next;
559 else {
560 x = 1 << k;
561 #ifdef Omit_Private_Memory
562 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
563 #else
564 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
565 /sizeof(double);
566 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
567 rv = (Bigint*)pmem_next;
568 pmem_next += len;
569 }
570 else
571 rv = (Bigint*)MALLOC(len*sizeof(double));
572 #endif
573 rv->k = k;
574 rv->maxwds = x;
575 }
576 FREE_DTOA_LOCK(0);
577 rv->sign = rv->wds = 0;
578 return rv;
579 }
580
581 static void
582 Bfree
583 #ifdef KR_headers
584 (v) Bigint *v;
585 #else
586 (Bigint *v)
587 #endif
588 {
589 if (v) {
590 if (v->k > Kmax)
591 #ifdef FREE
592 FREE((void*)v);
593 #else
594 free((void*)v);
595 #endif
596 else {
597 ACQUIRE_DTOA_LOCK(0);
598 v->next = freelist[v->k];
599 freelist[v->k] = v;
600 FREE_DTOA_LOCK(0);
601 }
602 }
603 }
604
605 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
606 y->wds*sizeof(Long) + 2*sizeof(int))
607
608 static Bigint *
609 multadd
610 #ifdef KR_headers
611 (b, m, a) Bigint *b; int m, a;
612 #else
613 (Bigint *b, int m, int a) /* multiply by m and add a */
614 #endif
615 {
616 int i, wds;
617 #ifdef ULLong
618 ULong *x;
619 ULLong carry, y;
620 #else
621 ULong carry, *x, y;
622 #ifdef Pack_32
623 ULong xi, z;
624 #endif
625 #endif
626 Bigint *b1;
627
628 wds = b->wds;
629 x = b->x;
630 i = 0;
631 carry = a;
632 do {
633 #ifdef ULLong
634 y = *x * (ULLong)m + carry;
635 carry = y >> 32;
636 *x++ = y & FFFFFFFF;
637 #else
638 #ifdef Pack_32
639 xi = *x;
640 y = (xi & 0xffff) * m + carry;
641 z = (xi >> 16) * m + (y >> 16);
642 carry = z >> 16;
643 *x++ = (z << 16) + (y & 0xffff);
644 #else
645 y = *x * m + carry;
646 carry = y >> 16;
647 *x++ = y & 0xffff;
648 #endif
649 #endif
650 }
651 while(++i < wds);
652 if (carry) {
653 if (wds >= b->maxwds) {
654 b1 = Balloc(b->k+1);
655 Bcopy(b1, b);
656 Bfree(b);
657 b = b1;
658 }
659 b->x[wds++] = carry;
660 b->wds = wds;
661 }
662 return b;
663 }
664
665 static Bigint *
666 s2b
667 #ifdef KR_headers
668 (s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9;
669 #else
670 (CONST char *s, int nd0, int nd, ULong y9, int dplen)
671 #endif
672 {
673 Bigint *b;
674 int i, k;
675 Long x, y;
676
677 x = (nd + 8) / 9;
678 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
679 #ifdef Pack_32
680 b = Balloc(k);
681 b->x[0] = y9;
682 b->wds = 1;
683 #else
684 b = Balloc(k+1);
685 b->x[0] = y9 & 0xffff;
686 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
687 #endif
688
689 i = 9;
690 if (9 < nd0) {
691 s += 9;
692 do b = multadd(b, 10, *s++ - '0');
693 while(++i < nd0);
694 s += dplen;
695 }
696 else
697 s += dplen + 9;
698 for(; i < nd; i++)
699 b = multadd(b, 10, *s++ - '0');
700 return b;
701 }
702
703 static int
704 hi0bits
705 #ifdef KR_headers
706 (x) ULong x;
707 #else
708 (ULong x)
709 #endif
710 {
711 int k = 0;
712
713 if (!(x & 0xffff0000)) {
714 k = 16;
715 x <<= 16;
716 }
717 if (!(x & 0xff000000)) {
718 k += 8;
719 x <<= 8;
720 }
721 if (!(x & 0xf0000000)) {
722 k += 4;
723 x <<= 4;
724 }
725 if (!(x & 0xc0000000)) {
726 k += 2;
727 x <<= 2;
728 }
729 if (!(x & 0x80000000)) {
730 k++;
731 if (!(x & 0x40000000))
732 return 32;
733 }
734 return k;
735 }
736
737 static int
738 lo0bits
739 #ifdef KR_headers
740 (y) ULong *y;
741 #else
742 (ULong *y)
743 #endif
744 {
745 int k;
746 ULong x = *y;
747
748 if (x & 7) {
749 if (x & 1)
750 return 0;
751 if (x & 2) {
752 *y = x >> 1;
753 return 1;
754 }
755 *y = x >> 2;
756 return 2;
757 }
758 k = 0;
759 if (!(x & 0xffff)) {
760 k = 16;
761 x >>= 16;
762 }
763 if (!(x & 0xff)) {
764 k += 8;
765 x >>= 8;
766 }
767 if (!(x & 0xf)) {
768 k += 4;
769 x >>= 4;
770 }
771 if (!(x & 0x3)) {
772 k += 2;
773 x >>= 2;
774 }
775 if (!(x & 1)) {
776 k++;
777 x >>= 1;
778 if (!x)
779 return 32;
780 }
781 *y = x;
782 return k;
783 }
784
785 static Bigint *
786 i2b
787 #ifdef KR_headers
788 (i) int i;
789 #else
790 (int i)
791 #endif
792 {
793 Bigint *b;
794
795 b = Balloc(1);
796 b->x[0] = i;
797 b->wds = 1;
798 return b;
799 }
800
801 static Bigint *
802 mult
803 #ifdef KR_headers
804 (a, b) Bigint *a, *b;
805 #else
806 (Bigint *a, Bigint *b)
807 #endif
808 {
809 Bigint *c;
810 int k, wa, wb, wc;
811 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
812 ULong y;
813 #ifdef ULLong
814 ULLong carry, z;
815 #else
816 ULong carry, z;
817 #ifdef Pack_32
818 ULong z2;
819 #endif
820 #endif
821
822 if (a->wds < b->wds) {
823 c = a;
824 a = b;
825 b = c;
826 }
827 k = a->k;
828 wa = a->wds;
829 wb = b->wds;
830 wc = wa + wb;
831 if (wc > a->maxwds)
832 k++;
833 c = Balloc(k);
834 for(x = c->x, xa = x + wc; x < xa; x++)
835 *x = 0;
836 xa = a->x;
837 xae = xa + wa;
838 xb = b->x;
839 xbe = xb + wb;
840 xc0 = c->x;
841 #ifdef ULLong
842 for(; xb < xbe; xc0++) {
843 if ((y = *xb++)) {
844 x = xa;
845 xc = xc0;
846 carry = 0;
847 do {
848 z = *x++ * (ULLong)y + *xc + carry;
849 carry = z >> 32;
850 *xc++ = z & FFFFFFFF;
851 }
852 while(x < xae);
853 *xc = carry;
854 }
855 }
856 #else
857 #ifdef Pack_32
858 for(; xb < xbe; xb++, xc0++) {
859 if (y = *xb & 0xffff) {
860 x = xa;
861 xc = xc0;
862 carry = 0;
863 do {
864 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
865 carry = z >> 16;
866 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
867 carry = z2 >> 16;
868 Storeinc(xc, z2, z);
869 }
870 while(x < xae);
871 *xc = carry;
872 }
873 if (y = *xb >> 16) {
874 x = xa;
875 xc = xc0;
876 carry = 0;
877 z2 = *xc;
878 do {
879 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
880 carry = z >> 16;
881 Storeinc(xc, z, z2);
882 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
883 carry = z2 >> 16;
884 }
885 while(x < xae);
886 *xc = z2;
887 }
888 }
889 #else
890 for(; xb < xbe; xc0++) {
891 if (y = *xb++) {
892 x = xa;
893 xc = xc0;
894 carry = 0;
895 do {
896 z = *x++ * y + *xc + carry;
897 carry = z >> 16;
898 *xc++ = z & 0xffff;
899 }
900 while(x < xae);
901 *xc = carry;
902 }
903 }
904 #endif
905 #endif
906 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
907 c->wds = wc;
908 return c;
909 }
910
911 static Bigint *p5s;
912
913 static Bigint *
914 pow5mult
915 #ifdef KR_headers
916 (b, k) Bigint *b; int k;
917 #else
918 (Bigint *b, int k)
919 #endif
920 {
921 Bigint *b1, *p5, *p51;
922 int i;
923 static int p05[3] = { 5, 25, 125 };
924
925 if ((i = k & 3))
926 b = multadd(b, p05[i-1], 0);
927
928 if (!(k >>= 2))
929 return b;
930 if (!(p5 = p5s)) {
931 /* first time */
932 #ifdef MULTIPLE_THREADS
933 ACQUIRE_DTOA_LOCK(1);
934 if (!(p5 = p5s)) {
935 p5 = p5s = i2b(625);
936 p5->next = 0;
937 }
938 FREE_DTOA_LOCK(1);
939 #else
940 p5 = p5s = i2b(625);
941 p5->next = 0;
942 #endif
943 }
944 for(;;) {
945 if (k & 1) {
946 b1 = mult(b, p5);
947 Bfree(b);
948 b = b1;
949 }
950 if (!(k >>= 1))
951 break;
952 if (!(p51 = p5->next)) {
953 #ifdef MULTIPLE_THREADS
954 ACQUIRE_DTOA_LOCK(1);
955 if (!(p51 = p5->next)) {
956 p51 = p5->next = mult(p5,p5);
957 p51->next = 0;
958 }
959 FREE_DTOA_LOCK(1);
960 #else
961 p51 = p5->next = mult(p5,p5);
962 p51->next = 0;
963 #endif
964 }
965 p5 = p51;
966 }
967 return b;
968 }
969
970 static Bigint *
971 lshift
972 #ifdef KR_headers
973 (b, k) Bigint *b; int k;
974 #else
975 (Bigint *b, int k)
976 #endif
977 {
978 int i, k1, n, n1;
979 Bigint *b1;
980 ULong *x, *x1, *xe, z;
981
982 #ifdef Pack_32
983 n = k >> 5;
984 #else
985 n = k >> 4;
986 #endif
987 k1 = b->k;
988 n1 = n + b->wds + 1;
989 for(i = b->maxwds; n1 > i; i <<= 1)
990 k1++;
991 b1 = Balloc(k1);
992 x1 = b1->x;
993 for(i = 0; i < n; i++)
994 *x1++ = 0;
995 x = b->x;
996 xe = x + b->wds;
997 #ifdef Pack_32
998 if (k &= 0x1f) {
999 k1 = 32 - k;
1000 z = 0;
1001 do {
1002 *x1++ = *x << k | z;
1003 z = *x++ >> k1;
1004 }
1005 while(x < xe);
1006 if ((*x1 = z))
1007 ++n1;
1008 }
1009 #else
1010 if (k &= 0xf) {
1011 k1 = 16 - k;
1012 z = 0;
1013 do {
1014 *x1++ = *x << k & 0xffff | z;
1015 z = *x++ >> k1;
1016 }
1017 while(x < xe);
1018 if (*x1 = z)
1019 ++n1;
1020 }
1021 #endif
1022 else do
1023 *x1++ = *x++;
1024 while(x < xe);
1025 b1->wds = n1 - 1;
1026 Bfree(b);
1027 return b1;
1028 }
1029
1030 static int
1031 cmp
1032 #ifdef KR_headers
1033 (a, b) Bigint *a, *b;
1034 #else
1035 (Bigint *a, Bigint *b)
1036 #endif
1037 {
1038 ULong *xa, *xa0, *xb, *xb0;
1039 int i, j;
1040
1041 i = a->wds;
1042 j = b->wds;
1043 #ifdef DEBUG
1044 if (i > 1 && !a->x[i-1])
1045 Bug("cmp called with a->x[a->wds-1] == 0");
1046 if (j > 1 && !b->x[j-1])
1047 Bug("cmp called with b->x[b->wds-1] == 0");
1048 #endif
1049 if (i -= j)
1050 return i;
1051 xa0 = a->x;
1052 xa = xa0 + j;
1053 xb0 = b->x;
1054 xb = xb0 + j;
1055 for(;;) {
1056 if (*--xa != *--xb)
1057 return *xa < *xb ? -1 : 1;
1058 if (xa <= xa0)
1059 break;
1060 }
1061 return 0;
1062 }
1063
1064 static Bigint *
1065 diff
1066 #ifdef KR_headers
1067 (a, b) Bigint *a, *b;
1068 #else
1069 (Bigint *a, Bigint *b)
1070 #endif
1071 {
1072 Bigint *c;
1073 int i, wa, wb;
1074 ULong *xa, *xae, *xb, *xbe, *xc;
1075 #ifdef ULLong
1076 ULLong borrow, y;
1077 #else
1078 ULong borrow, y;
1079 #ifdef Pack_32
1080 ULong z;
1081 #endif
1082 #endif
1083
1084 i = cmp(a,b);
1085 if (!i) {
1086 c = Balloc(0);
1087 c->wds = 1;
1088 c->x[0] = 0;
1089 return c;
1090 }
1091 if (i < 0) {
1092 c = a;
1093 a = b;
1094 b = c;
1095 i = 1;
1096 }
1097 else
1098 i = 0;
1099 c = Balloc(a->k);
1100 c->sign = i;
1101 wa = a->wds;
1102 xa = a->x;
1103 xae = xa + wa;
1104 wb = b->wds;
1105 xb = b->x;
1106 xbe = xb + wb;
1107 xc = c->x;
1108 borrow = 0;
1109 #ifdef ULLong
1110 do {
1111 y = (ULLong)*xa++ - *xb++ - borrow;
1112 borrow = y >> 32 & (ULong)1;
1113 *xc++ = y & FFFFFFFF;
1114 }
1115 while(xb < xbe);
1116 while(xa < xae) {
1117 y = *xa++ - borrow;
1118 borrow = y >> 32 & (ULong)1;
1119 *xc++ = y & FFFFFFFF;
1120 }
1121 #else
1122 #ifdef Pack_32
1123 do {
1124 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1125 borrow = (y & 0x10000) >> 16;
1126 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1127 borrow = (z & 0x10000) >> 16;
1128 Storeinc(xc, z, y);
1129 }
1130 while(xb < xbe);
1131 while(xa < xae) {
1132 y = (*xa & 0xffff) - borrow;
1133 borrow = (y & 0x10000) >> 16;
1134 z = (*xa++ >> 16) - borrow;
1135 borrow = (z & 0x10000) >> 16;
1136 Storeinc(xc, z, y);
1137 }
1138 #else
1139 do {
1140 y = *xa++ - *xb++ - borrow;
1141 borrow = (y & 0x10000) >> 16;
1142 *xc++ = y & 0xffff;
1143 }
1144 while(xb < xbe);
1145 while(xa < xae) {
1146 y = *xa++ - borrow;
1147 borrow = (y & 0x10000) >> 16;
1148 *xc++ = y & 0xffff;
1149 }
1150 #endif
1151 #endif
1152 while(!*--xc)
1153 wa--;
1154 c->wds = wa;
1155 return c;
1156 }
1157
1158 static double
1159 ulp
1160 #ifdef KR_headers
1161 (x) U *x;
1162 #else
1163 (U *x)
1164 #endif
1165 {
1166 Long L;
1167 U u;
1168
1169 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1170 #ifndef Avoid_Underflow
1171 #ifndef Sudden_Underflow
1172 if (L > 0) {
1173 #endif
1174 #endif
1175 #ifdef IBM
1176 L |= Exp_msk1 >> 4;
1177 #endif
1178 word0(&u) = L;
1179 word1(&u) = 0;
1180 #ifndef Avoid_Underflow
1181 #ifndef Sudden_Underflow
1182 }
1183 else {
1184 L = -L >> Exp_shift;
1185 if (L < Exp_shift) {
1186 word0(&u) = 0x80000 >> L;
1187 word1(&u) = 0;
1188 }
1189 else {
1190 word0(&u) = 0;
1191 L -= Exp_shift;
1192 word1(&u) = L >= 31 ? 1 : 1 << 31 - L;
1193 }
1194 }
1195 #endif
1196 #endif
1197 return dval(&u);
1198 }
1199
1200 static double
1201 b2d
1202 #ifdef KR_headers
1203 (a, e) Bigint *a; int *e;
1204 #else
1205 (Bigint *a, int *e)
1206 #endif
1207 {
1208 ULong *xa, *xa0, w, y, z;
1209 int k;
1210 U d;
1211 #ifdef VAX
1212 ULong d0, d1;
1213 #else
1214 #define d0 word0(&d)
1215 #define d1 word1(&d)
1216 #endif
1217
1218 xa0 = a->x;
1219 xa = xa0 + a->wds;
1220 y = *--xa;
1221 #ifdef DEBUG
1222 if (!y) Bug("zero y in b2d");
1223 #endif
1224 k = hi0bits(y);
1225 *e = 32 - k;
1226 #ifdef Pack_32
1227 if (k < Ebits) {
1228 d0 = Exp_1 | y >> (Ebits - k);
1229 w = xa > xa0 ? *--xa : 0;
1230 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1231 goto ret_d;
1232 }
1233 z = xa > xa0 ? *--xa : 0;
1234 if (k -= Ebits) {
1235 d0 = Exp_1 | y << k | z >> (32 - k);
1236 y = xa > xa0 ? *--xa : 0;
1237 d1 = z << k | y >> (32 - k);
1238 }
1239 else {
1240 d0 = Exp_1 | y;
1241 d1 = z;
1242 }
1243 #else
1244 if (k < Ebits + 16) {
1245 z = xa > xa0 ? *--xa : 0;
1246 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1247 w = xa > xa0 ? *--xa : 0;
1248 y = xa > xa0 ? *--xa : 0;
1249 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1250 goto ret_d;
1251 }
1252 z = xa > xa0 ? *--xa : 0;
1253 w = xa > xa0 ? *--xa : 0;
1254 k -= Ebits + 16;
1255 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1256 y = xa > xa0 ? *--xa : 0;
1257 d1 = w << k + 16 | y << k;
1258 #endif
1259 ret_d:
1260 #ifdef VAX
1261 word0(&d) = d0 >> 16 | d0 << 16;
1262 word1(&d) = d1 >> 16 | d1 << 16;
1263 #else
1264 #undef d0
1265 #undef d1
1266 #endif
1267 return dval(&d);
1268 }
1269
1270 static Bigint *
1271 d2b
1272 #ifdef KR_headers
1273 (d, e, bits) U *d; int *e, *bits;
1274 #else
1275 (U *d, int *e, int *bits)
1276 #endif
1277 {
1278 Bigint *b;
1279 int de, k;
1280 ULong *x, y, z;
1281 #ifndef Sudden_Underflow
1282 int i;
1283 #endif
1284 #ifdef VAX
1285 ULong d0, d1;
1286 d0 = word0(d) >> 16 | word0(d) << 16;
1287 d1 = word1(d) >> 16 | word1(d) << 16;
1288 #else
1289 #define d0 word0(d)
1290 #define d1 word1(d)
1291 #endif
1292
1293 #ifdef Pack_32
1294 b = Balloc(1);
1295 #else
1296 b = Balloc(2);
1297 #endif
1298 x = b->x;
1299
1300 z = d0 & Frac_mask;
1301 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1302 #ifdef Sudden_Underflow
1303 de = (int)(d0 >> Exp_shift);
1304 #ifndef IBM
1305 z |= Exp_msk11;
1306 #endif
1307 #else
1308 if ((de = (int)(d0 >> Exp_shift)))
1309 z |= Exp_msk1;
1310 #endif
1311 #ifdef Pack_32
1312 if ((y = d1)) {
1313 if ((k = lo0bits(&y))) {
1314 x[0] = y | z << (32 - k);
1315 z >>= k;
1316 }
1317 else
1318 x[0] = y;
1319 #ifndef Sudden_Underflow
1320 i =
1321 #endif
1322 b->wds = (x[1] = z) ? 2 : 1;
1323 }
1324 else {
1325 k = lo0bits(&z);
1326 x[0] = z;
1327 #ifndef Sudden_Underflow
1328 i =
1329 #endif
1330 b->wds = 1;
1331 k += 32;
1332 }
1333 #else
1334 if (y = d1) {
1335 if (k = lo0bits(&y))
1336 if (k >= 16) {
1337 x[0] = y | z << 32 - k & 0xffff;
1338 x[1] = z >> k - 16 & 0xffff;
1339 x[2] = z >> k;
1340 i = 2;
1341 }
1342 else {
1343 x[0] = y & 0xffff;
1344 x[1] = y >> 16 | z << 16 - k & 0xffff;
1345 x[2] = z >> k & 0xffff;
1346 x[3] = z >> k+16;
1347 i = 3;
1348 }
1349 else {
1350 x[0] = y & 0xffff;
1351 x[1] = y >> 16;
1352 x[2] = z & 0xffff;
1353 x[3] = z >> 16;
1354 i = 3;
1355 }
1356 }
1357 else {
1358 #ifdef DEBUG
1359 if (!z)
1360 Bug("Zero passed to d2b");
1361 #endif
1362 k = lo0bits(&z);
1363 if (k >= 16) {
1364 x[0] = z;
1365 i = 0;
1366 }
1367 else {
1368 x[0] = z & 0xffff;
1369 x[1] = z >> 16;
1370 i = 1;
1371 }
1372 k += 32;
1373 }
1374 while(!x[i])
1375 --i;
1376 b->wds = i + 1;
1377 #endif
1378 #ifndef Sudden_Underflow
1379 if (de) {
1380 #endif
1381 #ifdef IBM
1382 *e = (de - Bias - (P-1) << 2) + k;
1383 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1384 #else
1385 *e = de - Bias - (P-1) + k;
1386 *bits = P - k;
1387 #endif
1388 #ifndef Sudden_Underflow
1389 }
1390 else {
1391 *e = de - Bias - (P-1) + 1 + k;
1392 #ifdef Pack_32
1393 *bits = 32*i - hi0bits(x[i-1]);
1394 #else
1395 *bits = (i+2)*16 - hi0bits(x[i]);
1396 #endif
1397 }
1398 #endif
1399 return b;
1400 }
1401 #undef d0
1402 #undef d1
1403
1404 static double
1405 ratio
1406 #ifdef KR_headers
1407 (a, b) Bigint *a, *b;
1408 #else
1409 (Bigint *a, Bigint *b)
1410 #endif
1411 {
1412 U da, db;
1413 int k, ka, kb;
1414
1415 dval(&da) = b2d(a, &ka);
1416 dval(&db) = b2d(b, &kb);
1417 #ifdef Pack_32
1418 k = ka - kb + 32*(a->wds - b->wds);
1419 #else
1420 k = ka - kb + 16*(a->wds - b->wds);
1421 #endif
1422 #ifdef IBM
1423 if (k > 0) {
1424 word0(&da) += (k >> 2)*Exp_msk1;
1425 if (k &= 3)
1426 dval(&da) *= 1 << k;
1427 }
1428 else {
1429 k = -k;
1430 word0(&db) += (k >> 2)*Exp_msk1;
1431 if (k &= 3)
1432 dval(&db) *= 1 << k;
1433 }
1434 #else
1435 if (k > 0)
1436 word0(&da) += k*Exp_msk1;
1437 else {
1438 k = -k;
1439 word0(&db) += k*Exp_msk1;
1440 }
1441 #endif
1442 return dval(&da) / dval(&db);
1443 }
1444
1445 static CONST double
1446 tens[] = {
1447 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1448 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1449 1e20, 1e21, 1e22
1450 #ifdef VAX
1451 , 1e23, 1e24
1452 #endif
1453 };
1454
1455 static CONST double
1456 #ifdef IEEE_Arith
1457 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1458 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1459 #ifdef Avoid_Underflow
1460 9007199254740992.*9007199254740992.e-256
1461 /* = 2^106 * 1e-256 */
1462 #else
1463 1e-256
1464 #endif
1465 };
1466 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1467 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1468 #define Scale_Bit 0x10
1469 #define n_bigtens 5
1470 #else
1471 #ifdef IBM
1472 bigtens[] = { 1e16, 1e32, 1e64 };
1473 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1474 #define n_bigtens 3
1475 #else
1476 bigtens[] = { 1e16, 1e32 };
1477 static CONST double tinytens[] = { 1e-16, 1e-32 };
1478 #define n_bigtens 2
1479 #endif
1480 #endif
1481
1482 #undef Need_Hexdig
1483 #ifdef INFNAN_CHECK
1484 #ifndef No_Hex_NaN
1485 #define Need_Hexdig
1486 #endif
1487 #endif
1488
1489 #ifndef Need_Hexdig
1490 #ifndef NO_HEX_FP
1491 #define Need_Hexdig
1492 #endif
1493 #endif
1494
1495 #ifdef Need_Hexdig /*{*/
1496 static unsigned char hexdig[256];
1497
1498 static void
1499 #ifdef KR_headers
1500 htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc;
1501 #else
1502 htinit(unsigned char *h, unsigned char *s, int inc)
1503 #endif
1504 {
1505 int i, j;
1506 for(i = 0; (j = s[i]) !=0; i++)
1507 h[j] = i + inc;
1508 }
1509
1510 static void
1511 #ifdef KR_headers
hexdig_init()1512 hexdig_init()
1513 #else
1514 hexdig_init(void)
1515 #endif
1516 {
1517 #define USC (unsigned char *)
1518 htinit(hexdig, USC "0123456789", 0x10);
1519 htinit(hexdig, USC "abcdef", 0x10 + 10);
1520 htinit(hexdig, USC "ABCDEF", 0x10 + 10);
1521 }
1522 #endif /* } Need_Hexdig */
1523
1524 #ifdef INFNAN_CHECK
1525
1526 #ifndef NAN_WORD0
1527 #define NAN_WORD0 0x7ff80000
1528 #endif
1529
1530 #ifndef NAN_WORD1
1531 #define NAN_WORD1 0
1532 #endif
1533
1534 static int
1535 match
1536 #ifdef KR_headers
1537 (sp, t) char **sp, *t;
1538 #else
1539 (CONST char **sp, CONST char *t)
1540 #endif
1541 {
1542 int c, d;
1543 CONST char *s = *sp;
1544
1545 while((d = *t++)) {
1546 if ((c = *++s) >= 'A' && c <= 'Z')
1547 c += 'a' - 'A';
1548 if (c != d)
1549 return 0;
1550 }
1551 *sp = s + 1;
1552 return 1;
1553 }
1554
1555 #ifndef No_Hex_NaN
1556 static void
1557 hexnan
1558 #ifdef KR_headers
1559 (rvp, sp) U *rvp; CONST char **sp;
1560 #else
1561 (U *rvp, CONST char **sp)
1562 #endif
1563 {
1564 ULong c, x[2];
1565 CONST char *s;
1566 int c1, havedig, udx0, xshift;
1567
1568 if (!hexdig['0'])
1569 hexdig_init();
1570 x[0] = x[1] = 0;
1571 havedig = xshift = 0;
1572 udx0 = 1;
1573 s = *sp;
1574 /* allow optional initial 0x or 0X */
1575 while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
1576 ++s;
1577 if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
1578 s += 2;
1579 while((c = *(CONST unsigned char*)++s)) {
1580 if ((c1 = hexdig[c]))
1581 c = c1 & 0xf;
1582 else if (c <= ' ') {
1583 if (udx0 && havedig) {
1584 udx0 = 0;
1585 xshift = 1;
1586 }
1587 continue;
1588 }
1589 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
1590 else if (/*(*/ c == ')' && havedig) {
1591 *sp = s + 1;
1592 break;
1593 }
1594 else
1595 return; /* invalid form: don't change *sp */
1596 #else
1597 else {
1598 do {
1599 if (/*(*/ c == ')') {
1600 *sp = s + 1;
1601 break;
1602 }
1603 } while((c = *++s));
1604 break;
1605 }
1606 #endif
1607 havedig = 1;
1608 if (xshift) {
1609 xshift = 0;
1610 x[0] = x[1];
1611 x[1] = 0;
1612 }
1613 if (udx0)
1614 x[0] = (x[0] << 4) | (x[1] >> 28);
1615 x[1] = (x[1] << 4) | c;
1616 }
1617 if ((x[0] &= 0xfffff) || x[1]) {
1618 word0(rvp) = Exp_mask | x[0];
1619 word1(rvp) = x[1];
1620 }
1621 }
1622 #endif /*No_Hex_NaN*/
1623 #endif /* INFNAN_CHECK */
1624
1625 #ifdef Pack_32
1626 #define ULbits 32
1627 #define kshift 5
1628 #define kmask 31
1629 #else
1630 #define ULbits 16
1631 #define kshift 4
1632 #define kmask 15
1633 #endif
1634 #ifndef NO_HEX_FP /*{*/
1635
1636 static void
1637 #ifdef KR_headers
1638 rshift(b, k) Bigint *b; int k;
1639 #else
1640 rshift(Bigint *b, int k)
1641 #endif
1642 {
1643 ULong *x, *x1, *xe, y;
1644 int n;
1645
1646 x = x1 = b->x;
1647 n = k >> kshift;
1648 if (n < b->wds) {
1649 xe = x + b->wds;
1650 x += n;
1651 if (k &= kmask) {
1652 n = 32 - k;
1653 y = *x++ >> k;
1654 while(x < xe) {
1655 *x1++ = (y | (*x << n)) & 0xffffffff;
1656 y = *x++ >> k;
1657 }
1658 if ((*x1 = y) !=0)
1659 x1++;
1660 }
1661 else
1662 while(x < xe)
1663 *x1++ = *x++;
1664 }
1665 if ((b->wds = x1 - b->x) == 0)
1666 b->x[0] = 0;
1667 }
1668
1669 static ULong
1670 #ifdef KR_headers
1671 any_on(b, k) Bigint *b; int k;
1672 #else
1673 any_on(Bigint *b, int k)
1674 #endif
1675 {
1676 int n, nwds;
1677 ULong *x, *x0, x1, x2;
1678
1679 x = b->x;
1680 nwds = b->wds;
1681 n = k >> kshift;
1682 if (n > nwds)
1683 n = nwds;
1684 else if (n < nwds && (k &= kmask)) {
1685 x1 = x2 = x[n];
1686 x1 >>= k;
1687 x1 <<= k;
1688 if (x1 != x2)
1689 return 1;
1690 }
1691 x0 = x;
1692 x += n;
1693 while(x > x0)
1694 if (*--x)
1695 return 1;
1696 return 0;
1697 }
1698
1699 enum { /* rounding values: same as FLT_ROUNDS */
1700 Round_zero = 0,
1701 Round_near = 1,
1702 Round_up = 2,
1703 Round_down = 3
1704 };
1705
1706 static Bigint *
1707 #ifdef KR_headers
1708 increment(b) Bigint *b;
1709 #else
1710 increment(Bigint *b)
1711 #endif
1712 {
1713 ULong *x, *xe;
1714 Bigint *b1;
1715
1716 x = b->x;
1717 xe = x + b->wds;
1718 do {
1719 if (*x < (ULong)0xffffffffL) {
1720 ++*x;
1721 return b;
1722 }
1723 *x++ = 0;
1724 } while(x < xe);
1725 {
1726 if (b->wds >= b->maxwds) {
1727 b1 = Balloc(b->k+1);
1728 Bcopy(b1,b);
1729 Bfree(b);
1730 b = b1;
1731 }
1732 b->x[b->wds++] = 1;
1733 }
1734 return b;
1735 }
1736
1737 void
1738 #ifdef KR_headers
1739 gethex(sp, rvp, rounding, sign)
1740 CONST char **sp; U *rvp; int rounding, sign;
1741 #else
1742 gethex( CONST char **sp, U *rvp, int rounding, int sign)
1743 #endif
1744 {
1745 Bigint *b;
1746 CONST unsigned char *decpt, *s0, *s, *s1;
1747 Long e, e1;
1748 ULong L, lostbits, *x;
1749 int big, denorm, esign, havedig, k, n, nbits, up, zret;
1750 #ifdef IBM
1751 int j;
1752 #endif
1753 enum {
1754 #ifdef IEEE_Arith /*{{*/
1755 emax = 0x7fe - Bias - P + 1,
1756 emin = Emin - P + 1
1757 #else /*}{*/
1758 emin = Emin - P,
1759 #ifdef VAX
1760 emax = 0x7ff - Bias - P + 1
1761 #endif
1762 #ifdef IBM
1763 emax = 0x7f - Bias - P
1764 #endif
1765 #endif /*}}*/
1766 };
1767 #ifdef USE_LOCALE
1768 int i;
1769 #ifdef NO_LOCALE_CACHE
1770 const unsigned char *decimalpoint = (unsigned char*)
1771 localeconv()->decimal_point;
1772 #else
1773 const unsigned char *decimalpoint;
1774 static unsigned char *decimalpoint_cache;
1775 if (!(s0 = decimalpoint_cache)) {
1776 s0 = (unsigned char*)localeconv()->decimal_point;
1777 if ((decimalpoint_cache = (unsigned char*)
1778 MALLOC(strlen((CONST char*)s0) + 1))) {
1779 strcpy((char*)decimalpoint_cache, (CONST char*)s0);
1780 s0 = decimalpoint_cache;
1781 }
1782 }
1783 decimalpoint = s0;
1784 #endif
1785 #endif
1786
1787 if (!hexdig['0'])
1788 hexdig_init();
1789 havedig = 0;
1790 s0 = *(CONST unsigned char **)sp + 2;
1791 while(s0[havedig] == '0')
1792 havedig++;
1793 s0 += havedig;
1794 s = s0;
1795 decpt = 0;
1796 zret = 0;
1797 e = 0;
1798 if (hexdig[*s])
1799 havedig++;
1800 else {
1801 zret = 1;
1802 #ifdef USE_LOCALE
1803 for(i = 0; decimalpoint[i]; ++i) {
1804 if (s[i] != decimalpoint[i])
1805 goto pcheck;
1806 }
1807 decpt = s += i;
1808 #else
1809 if (*s != '.')
1810 goto pcheck;
1811 decpt = ++s;
1812 #endif
1813 if (!hexdig[*s])
1814 goto pcheck;
1815 while(*s == '0')
1816 s++;
1817 if (hexdig[*s])
1818 zret = 0;
1819 havedig = 1;
1820 s0 = s;
1821 }
1822 while(hexdig[*s])
1823 s++;
1824 #ifdef USE_LOCALE
1825 if (*s == *decimalpoint && !decpt) {
1826 for(i = 1; decimalpoint[i]; ++i) {
1827 if (s[i] != decimalpoint[i])
1828 goto pcheck;
1829 }
1830 decpt = s += i;
1831 #else
1832 if (*s == '.' && !decpt) {
1833 decpt = ++s;
1834 #endif
1835 while(hexdig[*s])
1836 s++;
1837 }/*}*/
1838 if (decpt)
1839 e = -(((Long)(s-decpt)) << 2);
1840 pcheck:
1841 s1 = s;
1842 big = esign = 0;
1843 switch(*s) {
1844 case 'p':
1845 case 'P':
1846 switch(*++s) {
1847 case '-':
1848 esign = 1;
1849 /* no break */
1850 case '+':
1851 s++;
1852 }
1853 if ((n = hexdig[*s]) == 0 || n > 0x19) {
1854 s = s1;
1855 break;
1856 }
1857 e1 = n - 0x10;
1858 while((n = hexdig[*++s]) !=0 && n <= 0x19) {
1859 if (e1 & 0xf8000000)
1860 big = 1;
1861 e1 = 10*e1 + n - 0x10;
1862 }
1863 if (esign)
1864 e1 = -e1;
1865 e += e1;
1866 }
1867 *sp = (char*)s;
1868 if (!havedig)
1869 *sp = (char*)s0 - 1;
1870 if (zret)
1871 goto retz1;
1872 if (big) {
1873 if (esign) {
1874 #ifdef IEEE_Arith
1875 switch(rounding) {
1876 case Round_up:
1877 if (sign)
1878 break;
1879 goto ret_tiny;
1880 case Round_down:
1881 if (!sign)
1882 break;
1883 goto ret_tiny;
1884 }
1885 #endif
1886 goto retz;
1887 #ifdef IEEE_Arith
1888 ret_tiny:
1889 #ifndef NO_ERRNO
1890 errno = ERANGE;
1891 #endif
1892 word0(rvp) = 0;
1893 word1(rvp) = 1;
1894 return;
1895 #endif /* IEEE_Arith */
1896 }
1897 switch(rounding) {
1898 case Round_near:
1899 goto ovfl1;
1900 case Round_up:
1901 if (!sign)
1902 goto ovfl1;
1903 goto ret_big;
1904 case Round_down:
1905 if (sign)
1906 goto ovfl1;
1907 goto ret_big;
1908 }
1909 ret_big:
1910 word0(rvp) = Big0;
1911 word1(rvp) = Big1;
1912 return;
1913 }
1914 n = s1 - s0 - 1;
1915 for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
1916 k++;
1917 b = Balloc(k);
1918 x = b->x;
1919 n = 0;
1920 L = 0;
1921 #ifdef USE_LOCALE
1922 for(i = 0; decimalpoint[i+1]; ++i);
1923 #endif
1924 while(s1 > s0) {
1925 #ifdef USE_LOCALE
1926 if (*--s1 == decimalpoint[i]) {
1927 s1 -= i;
1928 continue;
1929 }
1930 #else
1931 if (*--s1 == '.')
1932 continue;
1933 #endif
1934 if (n == ULbits) {
1935 *x++ = L;
1936 L = 0;
1937 n = 0;
1938 }
1939 L |= (hexdig[*s1] & 0x0f) << n;
1940 n += 4;
1941 }
1942 *x++ = L;
1943 b->wds = n = x - b->x;
1944 n = ULbits*n - hi0bits(L);
1945 nbits = Nbits;
1946 lostbits = 0;
1947 x = b->x;
1948 if (n > nbits) {
1949 n -= nbits;
1950 if (any_on(b,n)) {
1951 lostbits = 1;
1952 k = n - 1;
1953 if (x[k>>kshift] & 1 << (k & kmask)) {
1954 lostbits = 2;
1955 if (k > 0 && any_on(b,k))
1956 lostbits = 3;
1957 }
1958 }
1959 rshift(b, n);
1960 e += n;
1961 }
1962 else if (n < nbits) {
1963 n = nbits - n;
1964 b = lshift(b, n);
1965 e -= n;
1966 x = b->x;
1967 }
1968 if (e > Emax) {
1969 ovfl:
1970 Bfree(b);
1971 ovfl1:
1972 #ifndef NO_ERRNO
1973 errno = ERANGE;
1974 #endif
1975 word0(rvp) = Exp_mask;
1976 word1(rvp) = 0;
1977 return;
1978 }
1979 denorm = 0;
1980 if (e < emin) {
1981 denorm = 1;
1982 n = emin - e;
1983 if (n >= nbits) {
1984 #ifdef IEEE_Arith /*{*/
1985 switch (rounding) {
1986 case Round_near:
1987 if (n == nbits && (n < 2 || any_on(b,n-1)))
1988 goto ret_tiny;
1989 break;
1990 case Round_up:
1991 if (!sign)
1992 goto ret_tiny;
1993 break;
1994 case Round_down:
1995 if (sign)
1996 goto ret_tiny;
1997 }
1998 #endif /* } IEEE_Arith */
1999 Bfree(b);
2000 retz:
2001 #ifndef NO_ERRNO
2002 errno = ERANGE;
2003 #endif
2004 retz1:
2005 rvp->d = 0.;
2006 return;
2007 }
2008 k = n - 1;
2009 if (lostbits)
2010 lostbits = 1;
2011 else if (k > 0)
2012 lostbits = any_on(b,k);
2013 if (x[k>>kshift] & 1 << (k & kmask))
2014 lostbits |= 2;
2015 nbits -= n;
2016 rshift(b,n);
2017 e = emin;
2018 }
2019 if (lostbits) {
2020 up = 0;
2021 switch(rounding) {
2022 case Round_zero:
2023 break;
2024 case Round_near:
2025 if (lostbits & 2
2026 && (lostbits & 1) | (x[0] & 1))
2027 up = 1;
2028 break;
2029 case Round_up:
2030 up = 1 - sign;
2031 break;
2032 case Round_down:
2033 up = sign;
2034 }
2035 if (up) {
2036 k = b->wds;
2037 b = increment(b);
2038 x = b->x;
2039 if (denorm) {
2040 #if 0
2041 if (nbits == Nbits - 1
2042 && x[nbits >> kshift] & 1 << (nbits & kmask))
2043 denorm = 0; /* not currently used */
2044 #endif
2045 }
2046 else if (b->wds > k
2047 || ((n = nbits & kmask) !=0
2048 && hi0bits(x[k-1]) < 32-n)) {
2049 rshift(b,1);
2050 if (++e > Emax)
2051 goto ovfl;
2052 }
2053 }
2054 }
2055 #ifdef IEEE_Arith
2056 if (denorm)
2057 word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0;
2058 else
2059 word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20);
2060 word1(rvp) = b->x[0];
2061 #endif
2062 #ifdef IBM
2063 if ((j = e & 3)) {
2064 k = b->x[0] & ((1 << j) - 1);
2065 rshift(b,j);
2066 if (k) {
2067 switch(rounding) {
2068 case Round_up:
2069 if (!sign)
2070 increment(b);
2071 break;
2072 case Round_down:
2073 if (sign)
2074 increment(b);
2075 break;
2076 case Round_near:
2077 j = 1 << (j-1);
2078 if (k & j && ((k & (j-1)) | lostbits))
2079 increment(b);
2080 }
2081 }
2082 }
2083 e >>= 2;
2084 word0(rvp) = b->x[1] | ((e + 65 + 13) << 24);
2085 word1(rvp) = b->x[0];
2086 #endif
2087 #ifdef VAX
2088 /* The next two lines ignore swap of low- and high-order 2 bytes. */
2089 /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */
2090 /* word1(rvp) = b->x[0]; */
2091 word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16);
2092 word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16);
2093 #endif
2094 Bfree(b);
2095 }
2096 #endif /*}!NO_HEX_FP*/
2097
2098 static int
2099 #ifdef KR_headers
2100 dshift(b, p2) Bigint *b; int p2;
2101 #else
2102 dshift(Bigint *b, int p2)
2103 #endif
2104 {
2105 int rv = hi0bits(b->x[b->wds-1]) - 4;
2106 if (p2 > 0)
2107 rv -= p2;
2108 return rv & kmask;
2109 }
2110
2111 static int
2112 quorem
2113 #ifdef KR_headers
2114 (b, S) Bigint *b, *S;
2115 #else
2116 (Bigint *b, Bigint *S)
2117 #endif
2118 {
2119 int n;
2120 ULong *bx, *bxe, q, *sx, *sxe;
2121 #ifdef ULLong
2122 ULLong borrow, carry, y, ys;
2123 #else
2124 ULong borrow, carry, y, ys;
2125 #ifdef Pack_32
2126 ULong si, z, zs;
2127 #endif
2128 #endif
2129
2130 n = S->wds;
2131 #ifdef DEBUG
2132 /*debug*/ if (b->wds > n)
2133 /*debug*/ Bug("oversize b in quorem");
2134 #endif
2135 if (b->wds < n)
2136 return 0;
2137 sx = S->x;
2138 sxe = sx + --n;
2139 bx = b->x;
2140 bxe = bx + n;
2141 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
2142 #ifdef DEBUG
2143 /*debug*/ if (q > 9)
2144 /*debug*/ Bug("oversized quotient in quorem");
2145 #endif
2146 if (q) {
2147 borrow = 0;
2148 carry = 0;
2149 do {
2150 #ifdef ULLong
2151 ys = *sx++ * (ULLong)q + carry;
2152 carry = ys >> 32;
2153 y = *bx - (ys & FFFFFFFF) - borrow;
2154 borrow = y >> 32 & (ULong)1;
2155 *bx++ = y & FFFFFFFF;
2156 #else
2157 #ifdef Pack_32
2158 si = *sx++;
2159 ys = (si & 0xffff) * q + carry;
2160 zs = (si >> 16) * q + (ys >> 16);
2161 carry = zs >> 16;
2162 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2163 borrow = (y & 0x10000) >> 16;
2164 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2165 borrow = (z & 0x10000) >> 16;
2166 Storeinc(bx, z, y);
2167 #else
2168 ys = *sx++ * q + carry;
2169 carry = ys >> 16;
2170 y = *bx - (ys & 0xffff) - borrow;
2171 borrow = (y & 0x10000) >> 16;
2172 *bx++ = y & 0xffff;
2173 #endif
2174 #endif
2175 }
2176 while(sx <= sxe);
2177 if (!*bxe) {
2178 bx = b->x;
2179 while(--bxe > bx && !*bxe)
2180 --n;
2181 b->wds = n;
2182 }
2183 }
2184 if (cmp(b, S) >= 0) {
2185 q++;
2186 borrow = 0;
2187 carry = 0;
2188 bx = b->x;
2189 sx = S->x;
2190 do {
2191 #ifdef ULLong
2192 ys = *sx++ + carry;
2193 carry = ys >> 32;
2194 y = *bx - (ys & FFFFFFFF) - borrow;
2195 borrow = y >> 32 & (ULong)1;
2196 *bx++ = y & FFFFFFFF;
2197 #else
2198 #ifdef Pack_32
2199 si = *sx++;
2200 ys = (si & 0xffff) + carry;
2201 zs = (si >> 16) + (ys >> 16);
2202 carry = zs >> 16;
2203 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2204 borrow = (y & 0x10000) >> 16;
2205 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2206 borrow = (z & 0x10000) >> 16;
2207 Storeinc(bx, z, y);
2208 #else
2209 ys = *sx++ + carry;
2210 carry = ys >> 16;
2211 y = *bx - (ys & 0xffff) - borrow;
2212 borrow = (y & 0x10000) >> 16;
2213 *bx++ = y & 0xffff;
2214 #endif
2215 #endif
2216 }
2217 while(sx <= sxe);
2218 bx = b->x;
2219 bxe = bx + n;
2220 if (!*bxe) {
2221 while(--bxe > bx && !*bxe)
2222 --n;
2223 b->wds = n;
2224 }
2225 }
2226 return q;
2227 }
2228
2229 #ifndef NO_STRTOD_BIGCOMP
2230
2231 static void
2232 bigcomp
2233 #ifdef KR_headers
2234 (rv, s0, bc)
2235 U *rv; CONST char *s0; BCinfo *bc;
2236 #else
2237 (U *rv, CONST char *s0, BCinfo *bc)
2238 #endif
2239 {
2240 Bigint *b, *d;
2241 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;
2242
2243 dsign = bc->dsign;
2244 nd = bc->nd;
2245 nd0 = bc->nd0;
2246 p5 = nd + bc->e0 - 1;
2247 dd = speccase = 0;
2248 #ifndef Sudden_Underflow
2249 if (rv->d == 0.) { /* special case: value near underflow-to-zero */
2250 /* threshold was rounded to zero */
2251 b = i2b(1);
2252 p2 = Emin - P + 1;
2253 bbits = 1;
2254 #ifdef Avoid_Underflow
2255 word0(rv) = (P+2) << Exp_shift;
2256 #else
2257 word1(rv) = 1;
2258 #endif
2259 i = 0;
2260 #ifdef Honor_FLT_ROUNDS
2261 if (bc->rounding == 1)
2262 #endif
2263 {
2264 speccase = 1;
2265 --p2;
2266 dsign = 0;
2267 goto have_i;
2268 }
2269 }
2270 else
2271 #endif
2272 b = d2b(rv, &p2, &bbits);
2273 #ifdef Avoid_Underflow
2274 p2 -= bc->scale;
2275 #endif
2276 /* floor(log2(rv)) == bbits - 1 + p2 */
2277 /* Check for denormal case. */
2278 i = P - bbits;
2279 if (i > (j = P - Emin - 1 + p2)) {
2280 #ifdef Sudden_Underflow
2281 Bfree(b);
2282 b = i2b(1);
2283 p2 = Emin;
2284 i = P - 1;
2285 #ifdef Avoid_Underflow
2286 word0(rv) = (1 + bc->scale) << Exp_shift;
2287 #else
2288 word0(rv) = Exp_msk1;
2289 #endif
2290 word1(rv) = 0;
2291 #else
2292 i = j;
2293 #endif
2294 }
2295 #ifdef Honor_FLT_ROUNDS
2296 if (bc->rounding != 1) {
2297 if (i > 0)
2298 b = lshift(b, i);
2299 if (dsign)
2300 b = increment(b);
2301 }
2302 else
2303 #endif
2304 {
2305 b = lshift(b, ++i);
2306 b->x[0] |= 1;
2307 }
2308 #ifndef Sudden_Underflow
2309 have_i:
2310 #endif
2311 p2 -= p5 + i;
2312 d = i2b(1);
2313 /* Arrange for convenient computation of quotients:
2314 * shift left if necessary so divisor has 4 leading 0 bits.
2315 */
2316 if (p5 > 0)
2317 d = pow5mult(d, p5);
2318 else if (p5 < 0)
2319 b = pow5mult(b, -p5);
2320 if (p2 > 0) {
2321 b2 = p2;
2322 d2 = 0;
2323 }
2324 else {
2325 b2 = 0;
2326 d2 = -p2;
2327 }
2328 i = dshift(d, d2);
2329 if ((b2 += i) > 0)
2330 b = lshift(b, b2);
2331 if ((d2 += i) > 0)
2332 d = lshift(d, d2);
2333
2334 /* Now b/d = exactly half-way between the two floating-point values */
2335 /* on either side of the input string. Compute first digit of b/d. */
2336
2337 if (!(dig = quorem(b,d))) {
2338 b = multadd(b, 10, 0); /* very unlikely */
2339 dig = quorem(b,d);
2340 }
2341
2342 /* Compare b/d with s0 */
2343
2344 for(i = 0; i < nd0; ) {
2345 if ((dd = s0[i++] - '0' - dig))
2346 goto ret;
2347 if (!b->x[0] && b->wds == 1) {
2348 if (i < nd)
2349 dd = 1;
2350 goto ret;
2351 }
2352 b = multadd(b, 10, 0);
2353 dig = quorem(b,d);
2354 }
2355 for(j = bc->dp1; i++ < nd;) {
2356 if ((dd = s0[j++] - '0' - dig))
2357 goto ret;
2358 if (!b->x[0] && b->wds == 1) {
2359 if (i < nd)
2360 dd = 1;
2361 goto ret;
2362 }
2363 b = multadd(b, 10, 0);
2364 dig = quorem(b,d);
2365 }
2366 if (b->x[0] || b->wds > 1)
2367 dd = -1;
2368 ret:
2369 Bfree(b);
2370 Bfree(d);
2371 #ifdef Honor_FLT_ROUNDS
2372 if (bc->rounding != 1) {
2373 if (dd < 0) {
2374 if (bc->rounding == 0) {
2375 if (!dsign)
2376 goto retlow1;
2377 }
2378 else if (dsign)
2379 goto rethi1;
2380 }
2381 else if (dd > 0) {
2382 if (bc->rounding == 0) {
2383 if (dsign)
2384 goto rethi1;
2385 goto ret1;
2386 }
2387 if (!dsign)
2388 goto rethi1;
2389 dval(rv) += 2.*ulp(rv);
2390 }
2391 else {
2392 bc->inexact = 0;
2393 if (dsign)
2394 goto rethi1;
2395 }
2396 }
2397 else
2398 #endif
2399 if (speccase) {
2400 if (dd <= 0)
2401 rv->d = 0.;
2402 }
2403 else if (dd < 0) {
2404 if (!dsign) /* does not happen for round-near */
2405 retlow1:
2406 dval(rv) -= ulp(rv);
2407 }
2408 else if (dd > 0) {
2409 if (dsign) {
2410 rethi1:
2411 dval(rv) += ulp(rv);
2412 }
2413 }
2414 else {
2415 /* Exact half-way case: apply round-even rule. */
2416 if (word1(rv) & 1) {
2417 if (dsign)
2418 goto rethi1;
2419 goto retlow1;
2420 }
2421 }
2422
2423 #ifdef Honor_FLT_ROUNDS
2424 ret1:
2425 #endif
2426 return;
2427 }
2428 #endif /* NO_STRTOD_BIGCOMP */
2429
2430 double
2431 strtod
2432 #ifdef KR_headers
2433 (s00, se) CONST char *s00; char **se;
2434 #else
2435 (CONST char *s00, char **se)
2436 #endif
2437 {
2438 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1;
2439 int esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
2440 CONST char *s, *s0, *s1;
2441 double aadj, aadj1;
2442 Long L;
2443 U aadj2, adj, rv, rv0;
2444 ULong y, z;
2445 BCinfo bc;
2446 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
2447 #ifdef SET_INEXACT
2448 int oldinexact;
2449 #endif
2450 #ifdef Honor_FLT_ROUNDS /*{*/
2451 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
2452 bc.rounding = Flt_Rounds;
2453 #else /*}{*/
2454 bc.rounding = 1;
2455 switch(fegetround()) {
2456 case FE_TOWARDZERO: bc.rounding = 0; break;
2457 case FE_UPWARD: bc.rounding = 2; break;
2458 case FE_DOWNWARD: bc.rounding = 3;
2459 }
2460 #endif /*}}*/
2461 #endif /*}*/
2462 #ifdef USE_LOCALE
2463 CONST char *s2;
2464 #endif
2465
2466 sign = nz0 = nz = bc.dplen = bc.uflchk = 0;
2467 dval(&rv) = 0.;
2468 for(s = s00;;s++) switch(*s) {
2469 case '-':
2470 sign = 1;
2471 /* no break */
2472 case '+':
2473 if (*++s)
2474 goto break2;
2475 /* no break */
2476 case 0:
2477 goto ret0;
2478 case '\t':
2479 case '\n':
2480 case '\v':
2481 case '\f':
2482 case '\r':
2483 case ' ':
2484 continue;
2485 default:
2486 goto break2;
2487 }
2488 break2:
2489 if (*s == '0') {
2490 #ifndef NO_HEX_FP /*{*/
2491 switch(s[1]) {
2492 case 'x':
2493 case 'X':
2494 #ifdef Honor_FLT_ROUNDS
2495 gethex(&s, &rv, bc.rounding, sign);
2496 #else
2497 gethex(&s, &rv, 1, sign);
2498 #endif
2499 goto ret;
2500 }
2501 #endif /*}*/
2502 nz0 = 1;
2503 while(*++s == '0') ;
2504 if (!*s)
2505 goto ret;
2506 }
2507 s0 = s;
2508 y = z = 0;
2509 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
2510 if (nd < 9)
2511 y = 10*y + c - '0';
2512 else if (nd < 16)
2513 z = 10*z + c - '0';
2514 nd0 = nd;
2515 bc.dp0 = bc.dp1 = s - s0;
2516 #ifdef USE_LOCALE
2517 s1 = localeconv()->decimal_point;
2518 if (c == *s1) {
2519 c = '.';
2520 if (*++s1) {
2521 s2 = s;
2522 for(;;) {
2523 if (*++s2 != *s1) {
2524 c = 0;
2525 break;
2526 }
2527 if (!*++s1) {
2528 s = s2;
2529 break;
2530 }
2531 }
2532 }
2533 }
2534 #endif
2535 if (c == '.') {
2536 c = *++s;
2537 bc.dp1 = s - s0;
2538 bc.dplen = bc.dp1 - bc.dp0;
2539 if (!nd) {
2540 for(; c == '0'; c = *++s)
2541 nz++;
2542 if (c > '0' && c <= '9') {
2543 s0 = s;
2544 nf += nz;
2545 nz = 0;
2546 goto have_dig;
2547 }
2548 goto dig_done;
2549 }
2550 for(; c >= '0' && c <= '9'; c = *++s) {
2551 have_dig:
2552 nz++;
2553 if (c -= '0') {
2554 nf += nz;
2555 for(i = 1; i < nz; i++)
2556 if (nd++ < 9)
2557 y *= 10;
2558 else if (nd <= DBL_DIG + 1)
2559 z *= 10;
2560 if (nd++ < 9)
2561 y = 10*y + c;
2562 else if (nd <= DBL_DIG + 1)
2563 z = 10*z + c;
2564 nz = 0;
2565 }
2566 }
2567 }
2568 dig_done:
2569 e = 0;
2570 if (c == 'e' || c == 'E') {
2571 if (!nd && !nz && !nz0) {
2572 goto ret0;
2573 }
2574 s00 = s;
2575 esign = 0;
2576 switch(c = *++s) {
2577 case '-':
2578 esign = 1;
2579 case '+':
2580 c = *++s;
2581 }
2582 if (c >= '0' && c <= '9') {
2583 while(c == '0')
2584 c = *++s;
2585 if (c > '0' && c <= '9') {
2586 L = c - '0';
2587 s1 = s;
2588 while((c = *++s) >= '0' && c <= '9')
2589 L = 10*L + c - '0';
2590 if (s - s1 > 8 || L > 19999)
2591 /* Avoid confusion from exponents
2592 * so large that e might overflow.
2593 */
2594 e = 19999; /* safe for 16 bit ints */
2595 else
2596 e = (int)L;
2597 if (esign)
2598 e = -e;
2599 }
2600 else
2601 e = 0;
2602 }
2603 else
2604 s = s00;
2605 }
2606 if (!nd) {
2607 if (!nz && !nz0) {
2608 #ifdef INFNAN_CHECK
2609 /* Check for Nan and Infinity */
2610 if (!bc.dplen)
2611 switch(c) {
2612 case 'i':
2613 case 'I':
2614 if (match(&s,"nf")) {
2615 --s;
2616 if (!match(&s,"inity"))
2617 ++s;
2618 word0(&rv) = 0x7ff00000;
2619 word1(&rv) = 0;
2620 goto ret;
2621 }
2622 break;
2623 case 'n':
2624 case 'N':
2625 if (match(&s, "an")) {
2626 word0(&rv) = NAN_WORD0;
2627 word1(&rv) = NAN_WORD1;
2628 #ifndef No_Hex_NaN
2629 if (*s == '(') /*)*/
2630 hexnan(&rv, &s);
2631 #endif
2632 goto ret;
2633 }
2634 }
2635 #endif /* INFNAN_CHECK */
2636 ret0:
2637 s = s00;
2638 sign = 0;
2639 }
2640 goto ret;
2641 }
2642 bc.e0 = e1 = e -= nf;
2643
2644 /* Now we have nd0 digits, starting at s0, followed by a
2645 * decimal point, followed by nd-nd0 digits. The number we're
2646 * after is the integer represented by those digits times
2647 * 10**e */
2648
2649 if (!nd0)
2650 nd0 = nd;
2651 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
2652 dval(&rv) = y;
2653 if (k > 9) {
2654 #ifdef SET_INEXACT
2655 if (k > DBL_DIG)
2656 oldinexact = get_inexact();
2657 #endif
2658 dval(&rv) = tens[k - 9] * dval(&rv) + z;
2659 }
2660 bd0 = 0;
2661 if (nd <= DBL_DIG
2662 #ifndef RND_PRODQUOT
2663 #ifndef Honor_FLT_ROUNDS
2664 && Flt_Rounds == 1
2665 #endif
2666 #endif
2667 ) {
2668 if (!e)
2669 goto ret;
2670 if (e > 0) {
2671 if (e <= Ten_pmax) {
2672 #ifdef VAX
2673 goto vax_ovfl_check;
2674 #else
2675 #ifdef Honor_FLT_ROUNDS
2676 /* round correctly FLT_ROUNDS = 2 or 3 */
2677 if (sign) {
2678 rv.d = -rv.d;
2679 sign = 0;
2680 }
2681 #endif
2682 /* rv = */ rounded_product(dval(&rv), tens[e]);
2683 goto ret;
2684 #endif
2685 }
2686 i = DBL_DIG - nd;
2687 if (e <= Ten_pmax + i) {
2688 /* A fancier test would sometimes let us do
2689 * this for larger i values.
2690 */
2691 #ifdef Honor_FLT_ROUNDS
2692 /* round correctly FLT_ROUNDS = 2 or 3 */
2693 if (sign) {
2694 rv.d = -rv.d;
2695 sign = 0;
2696 }
2697 #endif
2698 e -= i;
2699 dval(&rv) *= tens[i];
2700 #ifdef VAX
2701 /* VAX exponent range is so narrow we must
2702 * worry about overflow here...
2703 */
2704 vax_ovfl_check:
2705 word0(&rv) -= P*Exp_msk1;
2706 /* rv = */ rounded_product(dval(&rv), tens[e]);
2707 if ((word0(&rv) & Exp_mask)
2708 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
2709 goto ovfl;
2710 word0(&rv) += P*Exp_msk1;
2711 #else
2712 /* rv = */ rounded_product(dval(&rv), tens[e]);
2713 #endif
2714 goto ret;
2715 }
2716 }
2717 #ifndef Inaccurate_Divide
2718 else if (e >= -Ten_pmax) {
2719 #ifdef Honor_FLT_ROUNDS
2720 /* round correctly FLT_ROUNDS = 2 or 3 */
2721 if (sign) {
2722 rv.d = -rv.d;
2723 sign = 0;
2724 }
2725 #endif
2726 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
2727 goto ret;
2728 }
2729 #endif
2730 }
2731 e1 += nd - k;
2732
2733 #ifdef IEEE_Arith
2734 #ifdef SET_INEXACT
2735 bc.inexact = 1;
2736 if (k <= DBL_DIG)
2737 oldinexact = get_inexact();
2738 #endif
2739 #ifdef Avoid_Underflow
2740 bc.scale = 0;
2741 #endif
2742 #ifdef Honor_FLT_ROUNDS
2743 if (bc.rounding >= 2) {
2744 if (sign)
2745 bc.rounding = bc.rounding == 2 ? 0 : 2;
2746 else
2747 if (bc.rounding != 2)
2748 bc.rounding = 0;
2749 }
2750 #endif
2751 #endif /*IEEE_Arith*/
2752
2753 /* Get starting approximation = rv * 10**e1 */
2754
2755 if (e1 > 0) {
2756 if ((i = e1 & 15))
2757 dval(&rv) *= tens[i];
2758 if (e1 &= ~15) {
2759 if (e1 > DBL_MAX_10_EXP) {
2760 ovfl:
2761 #ifndef NO_ERRNO
2762 errno = ERANGE;
2763 #endif
2764 /* Can't trust HUGE_VAL */
2765 #ifdef IEEE_Arith
2766 #ifdef Honor_FLT_ROUNDS
2767 switch(bc.rounding) {
2768 case 0: /* toward 0 */
2769 case 3: /* toward -infinity */
2770 word0(&rv) = Big0;
2771 word1(&rv) = Big1;
2772 break;
2773 default:
2774 word0(&rv) = Exp_mask;
2775 word1(&rv) = 0;
2776 }
2777 #else /*Honor_FLT_ROUNDS*/
2778 word0(&rv) = Exp_mask;
2779 word1(&rv) = 0;
2780 #endif /*Honor_FLT_ROUNDS*/
2781 #ifdef SET_INEXACT
2782 /* set overflow bit */
2783 dval(&rv0) = 1e300;
2784 dval(&rv0) *= dval(&rv0);
2785 #endif
2786 #else /*IEEE_Arith*/
2787 word0(&rv) = Big0;
2788 word1(&rv) = Big1;
2789 #endif /*IEEE_Arith*/
2790 goto ret;
2791 }
2792 e1 >>= 4;
2793 for(j = 0; e1 > 1; j++, e1 >>= 1)
2794 if (e1 & 1)
2795 dval(&rv) *= bigtens[j];
2796 /* The last multiplication could overflow. */
2797 word0(&rv) -= P*Exp_msk1;
2798 dval(&rv) *= bigtens[j];
2799 if ((z = word0(&rv) & Exp_mask)
2800 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
2801 goto ovfl;
2802 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
2803 /* set to largest number */
2804 /* (Can't trust DBL_MAX) */
2805 word0(&rv) = Big0;
2806 word1(&rv) = Big1;
2807 }
2808 else
2809 word0(&rv) += P*Exp_msk1;
2810 }
2811 }
2812 else if (e1 < 0) {
2813 e1 = -e1;
2814 if ((i = e1 & 15))
2815 dval(&rv) /= tens[i];
2816 if (e1 >>= 4) {
2817 if (e1 >= 1 << n_bigtens)
2818 goto undfl;
2819 #ifdef Avoid_Underflow
2820 if (e1 & Scale_Bit)
2821 bc.scale = 2*P;
2822 for(j = 0; e1 > 0; j++, e1 >>= 1)
2823 if (e1 & 1)
2824 dval(&rv) *= tinytens[j];
2825 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
2826 >> Exp_shift)) > 0) {
2827 /* scaled rv is denormal; clear j low bits */
2828 if (j >= 32) {
2829 word1(&rv) = 0;
2830 if (j >= 53)
2831 word0(&rv) = (P+2)*Exp_msk1;
2832 else
2833 word0(&rv) &= 0xffffffff << (j-32);
2834 }
2835 else
2836 word1(&rv) &= 0xffffffff << j;
2837 }
2838 #else
2839 for(j = 0; e1 > 1; j++, e1 >>= 1)
2840 if (e1 & 1)
2841 dval(&rv) *= tinytens[j];
2842 /* The last multiplication could underflow. */
2843 dval(&rv0) = dval(&rv);
2844 dval(&rv) *= tinytens[j];
2845 if (!dval(&rv)) {
2846 dval(&rv) = 2.*dval(&rv0);
2847 dval(&rv) *= tinytens[j];
2848 #endif
2849 if (!dval(&rv)) {
2850 undfl:
2851 dval(&rv) = 0.;
2852 #ifndef NO_ERRNO
2853 errno = ERANGE;
2854 #endif
2855 goto ret;
2856 }
2857 #ifndef Avoid_Underflow
2858 word0(&rv) = Tiny0;
2859 word1(&rv) = Tiny1;
2860 /* The refinement below will clean
2861 * this approximation up.
2862 */
2863 }
2864 #endif
2865 }
2866 }
2867
2868 /* Now the hard part -- adjusting rv to the correct value.*/
2869
2870 /* Put digits into bd: true value = bd * 10^e */
2871
2872 bc.nd = nd;
2873 #ifndef NO_STRTOD_BIGCOMP
2874 bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */
2875 /* to silence an erroneous warning about bc.nd0 */
2876 /* possibly not being initialized. */
2877 if (nd > strtod_diglim) {
2878 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */
2879 /* minimum number of decimal digits to distinguish double values */
2880 /* in IEEE arithmetic. */
2881 i = j = 18;
2882 if (i > nd0)
2883 j += bc.dplen;
2884 for(;;) {
2885 if (--j <= bc.dp1 && j >= bc.dp0)
2886 j = bc.dp0 - 1;
2887 if (s0[j] != '0')
2888 break;
2889 --i;
2890 }
2891 e += nd - i;
2892 nd = i;
2893 if (nd0 > nd)
2894 nd0 = nd;
2895 if (nd < 9) { /* must recompute y */
2896 y = 0;
2897 for(i = 0; i < nd0; ++i)
2898 y = 10*y + s0[i] - '0';
2899 for(j = bc.dp1; i < nd; ++i)
2900 y = 10*y + s0[j++] - '0';
2901 }
2902 }
2903 #endif
2904 bd0 = s2b(s0, nd0, nd, y, bc.dplen);
2905
2906 for(;;) {
2907 bd = Balloc(bd0->k);
2908 Bcopy(bd, bd0);
2909 bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
2910 bs = i2b(1);
2911
2912 if (e >= 0) {
2913 bb2 = bb5 = 0;
2914 bd2 = bd5 = e;
2915 }
2916 else {
2917 bb2 = bb5 = -e;
2918 bd2 = bd5 = 0;
2919 }
2920 if (bbe >= 0)
2921 bb2 += bbe;
2922 else
2923 bd2 -= bbe;
2924 bs2 = bb2;
2925 #ifdef Honor_FLT_ROUNDS
2926 if (bc.rounding != 1)
2927 bs2++;
2928 #endif
2929 #ifdef Avoid_Underflow
2930 j = bbe - bc.scale;
2931 i = j + bbbits - 1; /* logb(rv) */
2932 if (i < Emin) /* denormal */
2933 j += P - Emin;
2934 else
2935 j = P + 1 - bbbits;
2936 #else /*Avoid_Underflow*/
2937 #ifdef Sudden_Underflow
2938 #ifdef IBM
2939 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
2940 #else
2941 j = P + 1 - bbbits;
2942 #endif
2943 #else /*Sudden_Underflow*/
2944 j = bbe;
2945 i = j + bbbits - 1; /* logb(rv) */
2946 if (i < Emin) /* denormal */
2947 j += P - Emin;
2948 else
2949 j = P + 1 - bbbits;
2950 #endif /*Sudden_Underflow*/
2951 #endif /*Avoid_Underflow*/
2952 bb2 += j;
2953 bd2 += j;
2954 #ifdef Avoid_Underflow
2955 bd2 += bc.scale;
2956 #endif
2957 i = bb2 < bd2 ? bb2 : bd2;
2958 if (i > bs2)
2959 i = bs2;
2960 if (i > 0) {
2961 bb2 -= i;
2962 bd2 -= i;
2963 bs2 -= i;
2964 }
2965 if (bb5 > 0) {
2966 bs = pow5mult(bs, bb5);
2967 bb1 = mult(bs, bb);
2968 Bfree(bb);
2969 bb = bb1;
2970 }
2971 if (bb2 > 0)
2972 bb = lshift(bb, bb2);
2973 if (bd5 > 0)
2974 bd = pow5mult(bd, bd5);
2975 if (bd2 > 0)
2976 bd = lshift(bd, bd2);
2977 if (bs2 > 0)
2978 bs = lshift(bs, bs2);
2979 delta = diff(bb, bd);
2980 bc.dsign = delta->sign;
2981 delta->sign = 0;
2982 i = cmp(delta, bs);
2983 #ifndef NO_STRTOD_BIGCOMP
2984 if (bc.nd > nd && i <= 0) {
2985 if (bc.dsign)
2986 break; /* Must use bigcomp(). */
2987 #ifdef Honor_FLT_ROUNDS
2988 if (bc.rounding != 1) {
2989 if (i < 0)
2990 break;
2991 }
2992 else
2993 #endif
2994 {
2995 bc.nd = nd;
2996 i = -1; /* Discarded digits make delta smaller. */
2997 }
2998 }
2999 #endif
3000 #ifdef Honor_FLT_ROUNDS
3001 if (bc.rounding != 1) {
3002 if (i < 0) {
3003 /* Error is less than an ulp */
3004 if (!delta->x[0] && delta->wds <= 1) {
3005 /* exact */
3006 #ifdef SET_INEXACT
3007 bc.inexact = 0;
3008 #endif
3009 break;
3010 }
3011 if (bc.rounding) {
3012 if (bc.dsign) {
3013 adj.d = 1.;
3014 goto apply_adj;
3015 }
3016 }
3017 else if (!bc.dsign) {
3018 adj.d = -1.;
3019 if (!word1(&rv)
3020 && !(word0(&rv) & Frac_mask)) {
3021 y = word0(&rv) & Exp_mask;
3022 #ifdef Avoid_Underflow
3023 if (!bc.scale || y > 2*P*Exp_msk1)
3024 #else
3025 if (y)
3026 #endif
3027 {
3028 delta = lshift(delta,Log2P);
3029 if (cmp(delta, bs) <= 0)
3030 adj.d = -0.5;
3031 }
3032 }
3033 apply_adj:
3034 #ifdef Avoid_Underflow
3035 if (bc.scale && (y = word0(&rv) & Exp_mask)
3036 <= 2*P*Exp_msk1)
3037 word0(&adj) += (2*P+1)*Exp_msk1 - y;
3038 #else
3039 #ifdef Sudden_Underflow
3040 if ((word0(&rv) & Exp_mask) <=
3041 P*Exp_msk1) {
3042 word0(&rv) += P*Exp_msk1;
3043 dval(&rv) += adj.d*ulp(dval(&rv));
3044 word0(&rv) -= P*Exp_msk1;
3045 }
3046 else
3047 #endif /*Sudden_Underflow*/
3048 #endif /*Avoid_Underflow*/
3049 dval(&rv) += adj.d*ulp(&rv);
3050 }
3051 break;
3052 }
3053 adj.d = ratio(delta, bs);
3054 if (adj.d < 1.)
3055 adj.d = 1.;
3056 if (adj.d <= 0x7ffffffe) {
3057 /* adj = rounding ? ceil(adj) : floor(adj); */
3058 y = adj.d;
3059 if (y != adj.d) {
3060 if (!((bc.rounding>>1) ^ bc.dsign))
3061 y++;
3062 adj.d = y;
3063 }
3064 }
3065 #ifdef Avoid_Underflow
3066 if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3067 word0(&adj) += (2*P+1)*Exp_msk1 - y;
3068 #else
3069 #ifdef Sudden_Underflow
3070 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3071 word0(&rv) += P*Exp_msk1;
3072 adj.d *= ulp(dval(&rv));
3073 if (bc.dsign)
3074 dval(&rv) += adj.d;
3075 else
3076 dval(&rv) -= adj.d;
3077 word0(&rv) -= P*Exp_msk1;
3078 goto cont;
3079 }
3080 #endif /*Sudden_Underflow*/
3081 #endif /*Avoid_Underflow*/
3082 adj.d *= ulp(&rv);
3083 if (bc.dsign) {
3084 if (word0(&rv) == Big0 && word1(&rv) == Big1)
3085 goto ovfl;
3086 dval(&rv) += adj.d;
3087 }
3088 else
3089 dval(&rv) -= adj.d;
3090 goto cont;
3091 }
3092 #endif /*Honor_FLT_ROUNDS*/
3093
3094 if (i < 0) {
3095 /* Error is less than half an ulp -- check for
3096 * special case of mantissa a power of two.
3097 */
3098 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask
3099 #ifdef IEEE_Arith
3100 #ifdef Avoid_Underflow
3101 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
3102 #else
3103 || (word0(&rv) & Exp_mask) <= Exp_msk1
3104 #endif
3105 #endif
3106 ) {
3107 #ifdef SET_INEXACT
3108 if (!delta->x[0] && delta->wds <= 1)
3109 bc.inexact = 0;
3110 #endif
3111 break;
3112 }
3113 if (!delta->x[0] && delta->wds <= 1) {
3114 /* exact result */
3115 #ifdef SET_INEXACT
3116 bc.inexact = 0;
3117 #endif
3118 break;
3119 }
3120 delta = lshift(delta,Log2P);
3121 if (cmp(delta, bs) > 0)
3122 goto drop_down;
3123 break;
3124 }
3125 if (i == 0) {
3126 /* exactly half-way between */
3127 if (bc.dsign) {
3128 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
3129 && word1(&rv) == (
3130 #ifdef Avoid_Underflow
3131 (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3132 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
3133 #endif
3134 0xffffffff)) {
3135 /*boundary case -- increment exponent*/
3136 word0(&rv) = (word0(&rv) & Exp_mask)
3137 + Exp_msk1
3138 #ifdef IBM
3139 | Exp_msk1 >> 4
3140 #endif
3141 ;
3142 word1(&rv) = 0;
3143 #ifdef Avoid_Underflow
3144 bc.dsign = 0;
3145 #endif
3146 break;
3147 }
3148 }
3149 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
3150 drop_down:
3151 /* boundary case -- decrement exponent */
3152 #ifdef Sudden_Underflow /*{{*/
3153 L = word0(&rv) & Exp_mask;
3154 #ifdef IBM
3155 if (L < Exp_msk1)
3156 #else
3157 #ifdef Avoid_Underflow
3158 if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
3159 #else
3160 if (L <= Exp_msk1)
3161 #endif /*Avoid_Underflow*/
3162 #endif /*IBM*/
3163 {
3164 if (bc.nd >nd) {
3165 bc.uflchk = 1;
3166 break;
3167 }
3168 goto undfl;
3169 }
3170 L -= Exp_msk1;
3171 #else /*Sudden_Underflow}{*/
3172 #ifdef Avoid_Underflow
3173 if (bc.scale) {
3174 L = word0(&rv) & Exp_mask;
3175 if (L <= (2*P+1)*Exp_msk1) {
3176 if (L > (P+2)*Exp_msk1)
3177 /* round even ==> */
3178 /* accept rv */
3179 break;
3180 /* rv = smallest denormal */
3181 if (bc.nd >nd) {
3182 bc.uflchk = 1;
3183 break;
3184 }
3185 goto undfl;
3186 }
3187 }
3188 #endif /*Avoid_Underflow*/
3189 L = (word0(&rv) & Exp_mask) - Exp_msk1;
3190 #endif /*Sudden_Underflow}}*/
3191 word0(&rv) = L | Bndry_mask1;
3192 word1(&rv) = 0xffffffff;
3193 #ifdef IBM
3194 goto cont;
3195 #else
3196 break;
3197 #endif
3198 }
3199 #ifndef ROUND_BIASED
3200 if (!(word1(&rv) & LSB))
3201 break;
3202 #endif
3203 if (bc.dsign)
3204 dval(&rv) += ulp(&rv);
3205 #ifndef ROUND_BIASED
3206 else {
3207 dval(&rv) -= ulp(&rv);
3208 #ifndef Sudden_Underflow
3209 if (!dval(&rv)) {
3210 if (bc.nd >nd) {
3211 bc.uflchk = 1;
3212 break;
3213 }
3214 goto undfl;
3215 }
3216 #endif
3217 }
3218 #ifdef Avoid_Underflow
3219 bc.dsign = 1 - bc.dsign;
3220 #endif
3221 #endif
3222 break;
3223 }
3224 if ((aadj = ratio(delta, bs)) <= 2.) {
3225 if (bc.dsign)
3226 aadj = aadj1 = 1.;
3227 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
3228 #ifndef Sudden_Underflow
3229 if (word1(&rv) == Tiny1 && !word0(&rv)) {
3230 if (bc.nd >nd) {
3231 bc.uflchk = 1;
3232 break;
3233 }
3234 goto undfl;
3235 }
3236 #endif
3237 aadj = 1.;
3238 aadj1 = -1.;
3239 }
3240 else {
3241 /* special case -- power of FLT_RADIX to be */
3242 /* rounded down... */
3243
3244 if (aadj < 2./FLT_RADIX)
3245 aadj = 1./FLT_RADIX;
3246 else
3247 aadj *= 0.5;
3248 aadj1 = -aadj;
3249 }
3250 }
3251 else {
3252 aadj *= 0.5;
3253 aadj1 = bc.dsign ? aadj : -aadj;
3254 #ifdef Check_FLT_ROUNDS
3255 switch(bc.rounding) {
3256 case 2: /* towards +infinity */
3257 aadj1 -= 0.5;
3258 break;
3259 case 0: /* towards 0 */
3260 case 3: /* towards -infinity */
3261 aadj1 += 0.5;
3262 }
3263 #else
3264 if (Flt_Rounds == 0)
3265 aadj1 += 0.5;
3266 #endif /*Check_FLT_ROUNDS*/
3267 }
3268 y = word0(&rv) & Exp_mask;
3269
3270 /* Check for overflow */
3271
3272 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
3273 dval(&rv0) = dval(&rv);
3274 word0(&rv) -= P*Exp_msk1;
3275 adj.d = aadj1 * ulp(&rv);
3276 dval(&rv) += adj.d;
3277 if ((word0(&rv) & Exp_mask) >=
3278 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
3279 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
3280 goto ovfl;
3281 word0(&rv) = Big0;
3282 word1(&rv) = Big1;
3283 goto cont;
3284 }
3285 else
3286 word0(&rv) += P*Exp_msk1;
3287 }
3288 else {
3289 #ifdef Avoid_Underflow
3290 if (bc.scale && y <= 2*P*Exp_msk1) {
3291 if (aadj <= 0x7fffffff) {
3292 if ((z = aadj) <= 0)
3293 z = 1;
3294 aadj = z;
3295 aadj1 = bc.dsign ? aadj : -aadj;
3296 }
3297 dval(&aadj2) = aadj1;
3298 word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
3299 aadj1 = dval(&aadj2);
3300 }
3301 adj.d = aadj1 * ulp(&rv);
3302 dval(&rv) += adj.d;
3303 #else
3304 #ifdef Sudden_Underflow
3305 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3306 dval(&rv0) = dval(&rv);
3307 word0(&rv) += P*Exp_msk1;
3308 adj.d = aadj1 * ulp(&rv);
3309 dval(&rv) += adj.d;
3310 #ifdef IBM
3311 if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
3312 #else
3313 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
3314 #endif
3315 {
3316 if (word0(&rv0) == Tiny0
3317 && word1(&rv0) == Tiny1) {
3318 if (bc.nd >nd) {
3319 bc.uflchk = 1;
3320 break;
3321 }
3322 goto undfl;
3323 }
3324 word0(&rv) = Tiny0;
3325 word1(&rv) = Tiny1;
3326 goto cont;
3327 }
3328 else
3329 word0(&rv) -= P*Exp_msk1;
3330 }
3331 else {
3332 adj.d = aadj1 * ulp(&rv);
3333 dval(&rv) += adj.d;
3334 }
3335 #else /*Sudden_Underflow*/
3336 /* Compute adj so that the IEEE rounding rules will
3337 * correctly round rv + adj in some half-way cases.
3338 * If rv * ulp(rv) is denormalized (i.e.,
3339 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
3340 * trouble from bits lost to denormalization;
3341 * example: 1.2e-307 .
3342 */
3343 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
3344 aadj1 = (double)(int)(aadj + 0.5);
3345 if (!bc.dsign)
3346 aadj1 = -aadj1;
3347 }
3348 adj.d = aadj1 * ulp(&rv);
3349 dval(&rv) += adj.d;
3350 #endif /*Sudden_Underflow*/
3351 #endif /*Avoid_Underflow*/
3352 }
3353 z = word0(&rv) & Exp_mask;
3354 #ifndef SET_INEXACT
3355 if (bc.nd == nd) {
3356 #ifdef Avoid_Underflow
3357 if (!bc.scale)
3358 #endif
3359 if (y == z) {
3360 /* Can we stop now? */
3361 L = (Long)aadj;
3362 aadj -= L;
3363 /* The tolerances below are conservative. */
3364 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
3365 if (aadj < .4999999 || aadj > .5000001)
3366 break;
3367 }
3368 else if (aadj < .4999999/FLT_RADIX)
3369 break;
3370 }
3371 }
3372 #endif
3373 cont:
3374 Bfree(bb);
3375 Bfree(bd);
3376 Bfree(bs);
3377 Bfree(delta);
3378 }
3379 Bfree(bb);
3380 Bfree(bd);
3381 Bfree(bs);
3382 Bfree(bd0);
3383 Bfree(delta);
3384 #ifndef NO_STRTOD_BIGCOMP
3385 if (bc.nd > nd)
3386 bigcomp(&rv, s0, &bc);
3387 #endif
3388 #ifdef SET_INEXACT
3389 if (bc.inexact) {
3390 if (!oldinexact) {
3391 word0(&rv0) = Exp_1 + (70 << Exp_shift);
3392 word1(&rv0) = 0;
3393 dval(&rv0) += 1.;
3394 }
3395 }
3396 else if (!oldinexact)
3397 clear_inexact();
3398 #endif
3399 #ifdef Avoid_Underflow
3400 if (bc.scale) {
3401 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
3402 word1(&rv0) = 0;
3403 dval(&rv) *= dval(&rv0);
3404 #ifndef NO_ERRNO
3405 /* try to avoid the bug of testing an 8087 register value */
3406 #ifdef IEEE_Arith
3407 if (!(word0(&rv) & Exp_mask))
3408 #else
3409 if (word0(&rv) == 0 && word1(&rv) == 0)
3410 #endif
3411 errno = ERANGE;
3412 #endif
3413 }
3414 #endif /* Avoid_Underflow */
3415 #ifdef SET_INEXACT
3416 if (bc.inexact && !(word0(&rv) & Exp_mask)) {
3417 /* set underflow bit */
3418 dval(&rv0) = 1e-300;
3419 dval(&rv0) *= dval(&rv0);
3420 }
3421 #endif
3422 ret:
3423 if (se)
3424 *se = (char *)s;
3425 return sign ? -dval(&rv) : dval(&rv);
3426 }
3427
3428 #ifndef MULTIPLE_THREADS
3429 static char *dtoa_result;
3430 #endif
3431
3432 static char *
3433 #ifdef KR_headers
3434 rv_alloc(i) int i;
3435 #else
3436 rv_alloc(int i)
3437 #endif
3438 {
3439 int j, k, *r;
3440
3441 j = sizeof(ULong);
3442 for(k = 0;
3443 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (size_t)i;
3444 j <<= 1)
3445 k++;
3446 r = (int*)Balloc(k);
3447 *r = k;
3448 return
3449 #ifndef MULTIPLE_THREADS
3450 dtoa_result =
3451 #endif
3452 (char *)(r+1);
3453 }
3454
3455 static char *
3456 #ifdef KR_headers
3457 nrv_alloc(s, rve, n) char *s, **rve; int n;
3458 #else
3459 nrv_alloc(CONST char *s, char **rve, int n)
3460 #endif
3461 {
3462 char *rv, *t;
3463
3464 t = rv = rv_alloc(n);
3465 while((*t = *s++)) t++;
3466 if (rve)
3467 *rve = t;
3468 return rv;
3469 }
3470
3471 /* freedtoa(s) must be used to free values s returned by dtoa
3472 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
3473 * but for consistency with earlier versions of dtoa, it is optional
3474 * when MULTIPLE_THREADS is not defined.
3475 */
3476
3477 void
3478 #ifdef KR_headers
3479 freedtoa(s) char *s;
3480 #else
3481 freedtoa(char *s)
3482 #endif
3483 {
3484 Bigint *b = (Bigint *)((int *)s - 1);
3485 b->maxwds = 1 << (b->k = *(int*)b);
3486 Bfree(b);
3487 #ifndef MULTIPLE_THREADS
3488 if (s == dtoa_result)
3489 dtoa_result = 0;
3490 #endif
3491 }
3492
3493 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
3494 *
3495 * Inspired by "How to Print Floating-Point Numbers Accurately" by
3496 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
3497 *
3498 * Modifications:
3499 * 1. Rather than iterating, we use a simple numeric overestimate
3500 * to determine k = floor(log10(d)). We scale relevant
3501 * quantities using O(log2(k)) rather than O(k) multiplications.
3502 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
3503 * try to generate digits strictly left to right. Instead, we
3504 * compute with fewer bits and propagate the carry if necessary
3505 * when rounding the final digit up. This is often faster.
3506 * 3. Under the assumption that input will be rounded nearest,
3507 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
3508 * That is, we allow equality in stopping tests when the
3509 * round-nearest rule will give the same floating-point value
3510 * as would satisfaction of the stopping test with strict
3511 * inequality.
3512 * 4. We remove common factors of powers of 2 from relevant
3513 * quantities.
3514 * 5. When converting floating-point integers less than 1e16,
3515 * we use floating-point arithmetic rather than resorting
3516 * to multiple-precision integers.
3517 * 6. When asked to produce fewer than 15 digits, we first try
3518 * to get by with floating-point arithmetic; we resort to
3519 * multiple-precision integer arithmetic only if we cannot
3520 * guarantee that the floating-point calculation has given
3521 * the correctly rounded result. For k requested digits and
3522 * "uniformly" distributed input, the probability is
3523 * something like 10^(k-15) that we must resort to the Long
3524 * calculation.
3525 */
3526
3527 char *
3528 dtoa
3529 #ifdef KR_headers
3530 (dd, mode, ndigits, decpt, sign, rve)
3531 double dd; int mode, ndigits, *decpt, *sign; char **rve;
3532 #else
3533 (double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
3534 #endif
3535 {
3536 /* Arguments ndigits, decpt, sign are similar to those
3537 of ecvt and fcvt; trailing zeros are suppressed from
3538 the returned string. If not null, *rve is set to point
3539 to the end of the return value. If d is +-Infinity or NaN,
3540 then *decpt is set to 9999.
3541
3542 mode:
3543 0 ==> shortest string that yields d when read in
3544 and rounded to nearest.
3545 1 ==> like 0, but with Steele & White stopping rule;
3546 e.g. with IEEE P754 arithmetic , mode 0 gives
3547 1e23 whereas mode 1 gives 9.999999999999999e22.
3548 2 ==> max(1,ndigits) significant digits. This gives a
3549 return value similar to that of ecvt, except
3550 that trailing zeros are suppressed.
3551 3 ==> through ndigits past the decimal point. This
3552 gives a return value similar to that from fcvt,
3553 except that trailing zeros are suppressed, and
3554 ndigits can be negative.
3555 4,5 ==> similar to 2 and 3, respectively, but (in
3556 round-nearest mode) with the tests of mode 0 to
3557 possibly return a shorter string that rounds to d.
3558 With IEEE arithmetic and compilation with
3559 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
3560 as modes 2 and 3 when FLT_ROUNDS != 1.
3561 6-9 ==> Debugging modes similar to mode - 4: don't try
3562 fast floating-point estimate (if applicable).
3563
3564 Values of mode other than 0-9 are treated as mode 0.
3565
3566 Sufficient space is allocated to the return value
3567 to hold the suppressed trailing zeros.
3568 */
3569
3570 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
3571 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
3572 spec_case, try_quick;
3573 Long L;
3574 #ifndef Sudden_Underflow
3575 int denorm;
3576 ULong x;
3577 #endif
3578 Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
3579 U d2, eps, u;
3580 double ds;
3581 char *s, *s0;
3582 #ifdef SET_INEXACT
3583 int inexact, oldinexact;
3584 #endif
3585 #ifdef Honor_FLT_ROUNDS /*{*/
3586 int Rounding;
3587 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
3588 Rounding = Flt_Rounds;
3589 #else /*}{*/
3590 Rounding = 1;
3591 switch(fegetround()) {
3592 case FE_TOWARDZERO: Rounding = 0; break;
3593 case FE_UPWARD: Rounding = 2; break;
3594 case FE_DOWNWARD: Rounding = 3;
3595 }
3596 #endif /*}}*/
3597 #endif /*}*/
3598
3599 #ifndef MULTIPLE_THREADS
3600 if (dtoa_result) {
3601 freedtoa(dtoa_result);
3602 dtoa_result = 0;
3603 }
3604 #endif
3605
3606 u.d = dd;
3607 if (word0(&u) & Sign_bit) {
3608 /* set sign for everything, including 0's and NaNs */
3609 *sign = 1;
3610 word0(&u) &= ~Sign_bit; /* clear sign bit */
3611 }
3612 else
3613 *sign = 0;
3614
3615 #if defined(IEEE_Arith) + defined(VAX)
3616 #ifdef IEEE_Arith
3617 if ((word0(&u) & Exp_mask) == Exp_mask)
3618 #else
3619 if (word0(&u) == 0x8000)
3620 #endif
3621 {
3622 /* Infinity or NaN */
3623 *decpt = 9999;
3624 #ifdef IEEE_Arith
3625 if (!word1(&u) && !(word0(&u) & 0xfffff))
3626 return nrv_alloc("Infinity", rve, 8);
3627 #endif
3628 return nrv_alloc("NaN", rve, 3);
3629 }
3630 #endif
3631 #ifdef IBM
3632 dval(&u) += 0; /* normalize */
3633 #endif
3634 if (!dval(&u)) {
3635 *decpt = 1;
3636 return nrv_alloc("0", rve, 1);
3637 }
3638
3639 #ifdef SET_INEXACT
3640 try_quick = oldinexact = get_inexact();
3641 inexact = 1;
3642 #endif
3643 #ifdef Honor_FLT_ROUNDS
3644 if (Rounding >= 2) {
3645 if (*sign)
3646 Rounding = Rounding == 2 ? 0 : 2;
3647 else
3648 if (Rounding != 2)
3649 Rounding = 0;
3650 }
3651 #endif
3652
3653 b = d2b(&u, &be, &bbits);
3654 #ifdef Sudden_Underflow
3655 i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
3656 #else
3657 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
3658 #endif
3659 dval(&d2) = dval(&u);
3660 word0(&d2) &= Frac_mask1;
3661 word0(&d2) |= Exp_11;
3662 #ifdef IBM
3663 if (j = 11 - hi0bits(word0(&d2) & Frac_mask))
3664 dval(&d2) /= 1 << j;
3665 #endif
3666
3667 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
3668 * log10(x) = log(x) / log(10)
3669 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
3670 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
3671 *
3672 * This suggests computing an approximation k to log10(d) by
3673 *
3674 * k = (i - Bias)*0.301029995663981
3675 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
3676 *
3677 * We want k to be too large rather than too small.
3678 * The error in the first-order Taylor series approximation
3679 * is in our favor, so we just round up the constant enough
3680 * to compensate for any error in the multiplication of
3681 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
3682 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
3683 * adding 1e-13 to the constant term more than suffices.
3684 * Hence we adjust the constant term to 0.1760912590558.
3685 * (We could get a more accurate k by invoking log10,
3686 * but this is probably not worthwhile.)
3687 */
3688
3689 i -= Bias;
3690 #ifdef IBM
3691 i <<= 2;
3692 i += j;
3693 #endif
3694 #ifndef Sudden_Underflow
3695 denorm = 0;
3696 }
3697 else {
3698 /* d is denormalized */
3699
3700 i = bbits + be + (Bias + (P-1) - 1);
3701 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
3702 : word1(&u) << (32 - i);
3703 dval(&d2) = x;
3704 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
3705 i -= (Bias + (P-1) - 1) + 1;
3706 denorm = 1;
3707 }
3708 #endif
3709 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
3710 k = (int)ds;
3711 if (ds < 0. && ds != k)
3712 k--; /* want k = floor(ds) */
3713 k_check = 1;
3714 if (k >= 0 && k <= Ten_pmax) {
3715 if (dval(&u) < tens[k])
3716 k--;
3717 k_check = 0;
3718 }
3719 j = bbits - i - 1;
3720 if (j >= 0) {
3721 b2 = 0;
3722 s2 = j;
3723 }
3724 else {
3725 b2 = -j;
3726 s2 = 0;
3727 }
3728 if (k >= 0) {
3729 b5 = 0;
3730 s5 = k;
3731 s2 += k;
3732 }
3733 else {
3734 b2 -= k;
3735 b5 = -k;
3736 s5 = 0;
3737 }
3738 if (mode < 0 || mode > 9)
3739 mode = 0;
3740
3741 #ifndef SET_INEXACT
3742 #ifdef Check_FLT_ROUNDS
3743 try_quick = Rounding == 1;
3744 #else
3745 try_quick = 1;
3746 #endif
3747 #endif /*SET_INEXACT*/
3748
3749 if (mode > 5) {
3750 mode -= 4;
3751 try_quick = 0;
3752 }
3753 leftright = 1;
3754 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
3755 /* silence erroneous "gcc -Wall" warning. */
3756 switch(mode) {
3757 case 0:
3758 case 1:
3759 i = 18;
3760 ndigits = 0;
3761 break;
3762 case 2:
3763 leftright = 0;
3764 /* no break */
3765 case 4:
3766 if (ndigits <= 0)
3767 ndigits = 1;
3768 ilim = ilim1 = i = ndigits;
3769 break;
3770 case 3:
3771 leftright = 0;
3772 /* no break */
3773 case 5:
3774 i = ndigits + k + 1;
3775 ilim = i;
3776 ilim1 = i - 1;
3777 if (i <= 0)
3778 i = 1;
3779 }
3780 s = s0 = rv_alloc(i);
3781
3782 #ifdef Honor_FLT_ROUNDS
3783 if (mode > 1 && Rounding != 1)
3784 leftright = 0;
3785 #endif
3786
3787 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
3788
3789 /* Try to get by with floating-point arithmetic. */
3790
3791 i = 0;
3792 dval(&d2) = dval(&u);
3793 k0 = k;
3794 ilim0 = ilim;
3795 ieps = 2; /* conservative */
3796 if (k > 0) {
3797 ds = tens[k&0xf];
3798 j = k >> 4;
3799 if (j & Bletch) {
3800 /* prevent overflows */
3801 j &= Bletch - 1;
3802 dval(&u) /= bigtens[n_bigtens-1];
3803 ieps++;
3804 }
3805 for(; j; j >>= 1, i++)
3806 if (j & 1) {
3807 ieps++;
3808 ds *= bigtens[i];
3809 }
3810 dval(&u) /= ds;
3811 }
3812 else if ((j1 = -k)) {
3813 dval(&u) *= tens[j1 & 0xf];
3814 for(j = j1 >> 4; j; j >>= 1, i++)
3815 if (j & 1) {
3816 ieps++;
3817 dval(&u) *= bigtens[i];
3818 }
3819 }
3820 if (k_check && dval(&u) < 1. && ilim > 0) {
3821 if (ilim1 <= 0)
3822 goto fast_failed;
3823 ilim = ilim1;
3824 k--;
3825 dval(&u) *= 10.;
3826 ieps++;
3827 }
3828 dval(&eps) = ieps*dval(&u) + 7.;
3829 word0(&eps) -= (P-1)*Exp_msk1;
3830 if (ilim == 0) {
3831 S = mhi = 0;
3832 dval(&u) -= 5.;
3833 if (dval(&u) > dval(&eps))
3834 goto one_digit;
3835 if (dval(&u) < -dval(&eps))
3836 goto no_digits;
3837 goto fast_failed;
3838 }
3839 #ifndef No_leftright
3840 if (leftright) {
3841 /* Use Steele & White method of only
3842 * generating digits needed.
3843 */
3844 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
3845 for(i = 0;;) {
3846 L = dval(&u);
3847 dval(&u) -= L;
3848 *s++ = '0' + (int)L;
3849 if (dval(&u) < dval(&eps))
3850 goto ret1;
3851 if (1. - dval(&u) < dval(&eps))
3852 goto bump_up;
3853 if (++i >= ilim)
3854 break;
3855 dval(&eps) *= 10.;
3856 dval(&u) *= 10.;
3857 }
3858 }
3859 else {
3860 #endif
3861 /* Generate ilim digits, then fix them up. */
3862 dval(&eps) *= tens[ilim-1];
3863 for(i = 1;; i++, dval(&u) *= 10.) {
3864 L = (Long)(dval(&u));
3865 if (!(dval(&u) -= L))
3866 ilim = i;
3867 *s++ = '0' + (int)L;
3868 if (i == ilim) {
3869 if (dval(&u) > 0.5 + dval(&eps))
3870 goto bump_up;
3871 else if (dval(&u) < 0.5 - dval(&eps)) {
3872 while(*--s == '0') {}
3873 s++;
3874 goto ret1;
3875 }
3876 break;
3877 }
3878 }
3879 #ifndef No_leftright
3880 }
3881 #endif
3882 fast_failed:
3883 s = s0;
3884 dval(&u) = dval(&d2);
3885 k = k0;
3886 ilim = ilim0;
3887 }
3888
3889 /* Do we have a "small" integer? */
3890
3891 if (be >= 0 && k <= Int_max) {
3892 /* Yes. */
3893 ds = tens[k];
3894 if (ndigits < 0 && ilim <= 0) {
3895 S = mhi = 0;
3896 if (ilim < 0 || dval(&u) <= 5*ds)
3897 goto no_digits;
3898 goto one_digit;
3899 }
3900 for(i = 1; i <= k + 1; i++, dval(&u) *= 10.) {
3901 L = (Long)(dval(&u) / ds);
3902 dval(&u) -= L*ds;
3903 #ifdef Check_FLT_ROUNDS
3904 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
3905 if (dval(&u) < 0) {
3906 L--;
3907 dval(&u) += ds;
3908 }
3909 #endif
3910 *s++ = '0' + (int)L;
3911 if (!dval(&u)) {
3912 #ifdef SET_INEXACT
3913 inexact = 0;
3914 #endif
3915 break;
3916 }
3917 if (i == ilim) {
3918 #ifdef Honor_FLT_ROUNDS
3919 if (mode > 1)
3920 switch(Rounding) {
3921 case 0: goto ret1;
3922 case 2: goto bump_up;
3923 }
3924 #endif
3925 dval(&u) += dval(&u);
3926 if (dval(&u) > ds || (dval(&u) == ds && L & 1)) {
3927 bump_up:
3928 while(*--s == '9')
3929 if (s == s0) {
3930 k++;
3931 *s = '0';
3932 break;
3933 }
3934 ++*s++;
3935 }
3936 break;
3937 }
3938 }
3939 goto ret1;
3940 }
3941
3942 m2 = b2;
3943 m5 = b5;
3944 mhi = mlo = 0;
3945 if (leftright) {
3946 i =
3947 #ifndef Sudden_Underflow
3948 denorm ? be + (Bias + (P-1) - 1 + 1) :
3949 #endif
3950 #ifdef IBM
3951 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
3952 #else
3953 1 + P - bbits;
3954 #endif
3955 b2 += i;
3956 s2 += i;
3957 mhi = i2b(1);
3958 }
3959 if (m2 > 0 && s2 > 0) {
3960 i = m2 < s2 ? m2 : s2;
3961 b2 -= i;
3962 m2 -= i;
3963 s2 -= i;
3964 }
3965 if (b5 > 0) {
3966 if (leftright) {
3967 if (m5 > 0) {
3968 mhi = pow5mult(mhi, m5);
3969 b1 = mult(mhi, b);
3970 Bfree(b);
3971 b = b1;
3972 }
3973 if ((j = b5 - m5))
3974 b = pow5mult(b, j);
3975 }
3976 else
3977 b = pow5mult(b, b5);
3978 }
3979 S = i2b(1);
3980 if (s5 > 0)
3981 S = pow5mult(S, s5);
3982
3983 /* Check for special case that d is a normalized power of 2. */
3984
3985 spec_case = 0;
3986 if ((mode < 2 || leftright)
3987 #ifdef Honor_FLT_ROUNDS
3988 && Rounding == 1
3989 #endif
3990 ) {
3991 if (!word1(&u) && !(word0(&u) & Bndry_mask)
3992 #ifndef Sudden_Underflow
3993 && word0(&u) & (Exp_mask & ~Exp_msk1)
3994 #endif
3995 ) {
3996 /* The special case */
3997 b2 += Log2P;
3998 s2 += Log2P;
3999 spec_case = 1;
4000 }
4001 }
4002
4003 /* Arrange for convenient computation of quotients:
4004 * shift left if necessary so divisor has 4 leading 0 bits.
4005 *
4006 * Perhaps we should just compute leading 28 bits of S once
4007 * and for all and pass them and a shift to quorem, so it
4008 * can do shifts and ors to compute the numerator for q.
4009 */
4010 #ifdef Pack_32
4011 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
4012 i = 32 - i;
4013 #define iInc 28
4014 #else
4015 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
4016 i = 16 - i;
4017 #define iInc 12
4018 #endif
4019 i = dshift(S, s2);
4020 b2 += i;
4021 m2 += i;
4022 s2 += i;
4023 if (b2 > 0)
4024 b = lshift(b, b2);
4025 if (s2 > 0)
4026 S = lshift(S, s2);
4027 if (k_check) {
4028 if (cmp(b,S) < 0) {
4029 k--;
4030 b = multadd(b, 10, 0); /* we botched the k estimate */
4031 if (leftright)
4032 mhi = multadd(mhi, 10, 0);
4033 ilim = ilim1;
4034 }
4035 }
4036 if (ilim <= 0 && (mode == 3 || mode == 5)) {
4037 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
4038 /* no digits, fcvt style */
4039 no_digits:
4040 k = -1 - ndigits;
4041 goto ret;
4042 }
4043 one_digit:
4044 *s++ = '1';
4045 k++;
4046 goto ret;
4047 }
4048 if (leftright) {
4049 if (m2 > 0)
4050 mhi = lshift(mhi, m2);
4051
4052 /* Compute mlo -- check for special case
4053 * that d is a normalized power of 2.
4054 */
4055
4056 mlo = mhi;
4057 if (spec_case) {
4058 mhi = Balloc(mhi->k);
4059 Bcopy(mhi, mlo);
4060 mhi = lshift(mhi, Log2P);
4061 }
4062
4063 for(i = 1;;i++) {
4064 dig = quorem(b,S) + '0';
4065 /* Do we yet have the shortest decimal string
4066 * that will round to d?
4067 */
4068 j = cmp(b, mlo);
4069 delta = diff(S, mhi);
4070 j1 = delta->sign ? 1 : cmp(b, delta);
4071 Bfree(delta);
4072 #ifndef ROUND_BIASED
4073 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
4074 #ifdef Honor_FLT_ROUNDS
4075 && Rounding >= 1
4076 #endif
4077 ) {
4078 if (dig == '9')
4079 goto round_9_up;
4080 if (j > 0)
4081 dig++;
4082 #ifdef SET_INEXACT
4083 else if (!b->x[0] && b->wds <= 1)
4084 inexact = 0;
4085 #endif
4086 *s++ = dig;
4087 goto ret;
4088 }
4089 #endif
4090 if (j < 0 || (j == 0 && mode != 1
4091 #ifndef ROUND_BIASED
4092 && !(word1(&u) & 1)
4093 #endif
4094 )) {
4095 if (!b->x[0] && b->wds <= 1) {
4096 #ifdef SET_INEXACT
4097 inexact = 0;
4098 #endif
4099 goto accept_dig;
4100 }
4101 #ifdef Honor_FLT_ROUNDS
4102 if (mode > 1)
4103 switch(Rounding) {
4104 case 0: goto accept_dig;
4105 case 2: goto keep_dig;
4106 }
4107 #endif /*Honor_FLT_ROUNDS*/
4108 if (j1 > 0) {
4109 b = lshift(b, 1);
4110 j1 = cmp(b, S);
4111 if ((j1 > 0 || (j1 == 0 && dig & 1))
4112 && dig++ == '9')
4113 goto round_9_up;
4114 }
4115 accept_dig:
4116 *s++ = dig;
4117 goto ret;
4118 }
4119 if (j1 > 0) {
4120 #ifdef Honor_FLT_ROUNDS
4121 if (!Rounding)
4122 goto accept_dig;
4123 #endif
4124 if (dig == '9') { /* possible if i == 1 */
4125 round_9_up:
4126 *s++ = '9';
4127 goto roundoff;
4128 }
4129 *s++ = dig + 1;
4130 goto ret;
4131 }
4132 #ifdef Honor_FLT_ROUNDS
4133 keep_dig:
4134 #endif
4135 *s++ = dig;
4136 if (i == ilim)
4137 break;
4138 b = multadd(b, 10, 0);
4139 if (mlo == mhi)
4140 mlo = mhi = multadd(mhi, 10, 0);
4141 else {
4142 mlo = multadd(mlo, 10, 0);
4143 mhi = multadd(mhi, 10, 0);
4144 }
4145 }
4146 }
4147 else
4148 for(i = 1;; i++) {
4149 *s++ = dig = quorem(b,S) + '0';
4150 if (!b->x[0] && b->wds <= 1) {
4151 #ifdef SET_INEXACT
4152 inexact = 0;
4153 #endif
4154 goto ret;
4155 }
4156 if (i >= ilim)
4157 break;
4158 b = multadd(b, 10, 0);
4159 }
4160
4161 /* Round off last digit */
4162
4163 #ifdef Honor_FLT_ROUNDS
4164 switch(Rounding) {
4165 case 0: goto trimzeros;
4166 case 2: goto roundoff;
4167 }
4168 #endif
4169 b = lshift(b, 1);
4170 j = cmp(b, S);
4171 if (j > 0 || (j == 0 && dig & 1)) {
4172 roundoff:
4173 while(*--s == '9')
4174 if (s == s0) {
4175 k++;
4176 *s++ = '1';
4177 goto ret;
4178 }
4179 ++*s++;
4180 }
4181 else {
4182 #ifdef Honor_FLT_ROUNDS
4183 trimzeros:
4184 #endif
4185 while(*--s == '0') {}
4186 s++;
4187 }
4188 ret:
4189 Bfree(S);
4190 if (mhi) {
4191 if (mlo && mlo != mhi)
4192 Bfree(mlo);
4193 Bfree(mhi);
4194 }
4195 ret1:
4196 #ifdef SET_INEXACT
4197 if (inexact) {
4198 if (!oldinexact) {
4199 word0(&u) = Exp_1 + (70 << Exp_shift);
4200 word1(&u) = 0;
4201 dval(&u) += 1.;
4202 }
4203 }
4204 else if (!oldinexact)
4205 clear_inexact();
4206 #endif
4207 Bfree(b);
4208 *s = 0;
4209 *decpt = k + 1;
4210 if (rve)
4211 *rve = s;
4212 return s0;
4213 }
4214
4215 } // namespace dmg_fp
4216