/[jscoverage]/trunk/js/dtoa.c
ViewVC logotype

Contents of /trunk/js/dtoa.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 585 - (show annotations)
Sun Sep 12 15:13:23 2010 UTC (9 years, 1 month ago) by siliconforks
File MIME type: text/plain
File size: 68725 byte(s)
Update to SpiderMonkey from Firefox 3.6.9.

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

  ViewVC Help
Powered by ViewVC 1.1.24