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