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