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