The Tor Browser / file comparison

comparison: js/src/dtoa.c

js/src/dtoa.c

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

Mercurial Repository: The Tor Browser

Europalab SCM Repository Server

mercurial