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