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 const 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 int j, k, *r;
3617
3618 j = sizeof(ULong);
3619 for(k = 0;
3620 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + (size_t)j <= (size_t)i;
3621 j <<= 1)
3622 k++;
3623 r = (int*)Balloc(k);
3624 *r = k;
3625 return
3626 #ifndef MULTIPLE_THREADS
3627 dtoa_result =
3628 #endif
3629 (char *)(r+1);
3630 }
3631
3632 static char *
3633 #ifdef KR_headers
nrv_alloc(s,rve,n)3634 nrv_alloc(s, rve, n) char *s, **rve; int n;
3635 #else
3636 nrv_alloc(const char *s, char **rve, int n)
3637 #endif
3638 {
3639 char *rv, *t;
3640
3641 t = rv = rv_alloc(n);
3642 while((*t = *s++)) t++;
3643 if (rve)
3644 *rve = t;
3645 return rv;
3646 }
3647
3648 /* freedtoa(s) must be used to free values s returned by dtoa
3649 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
3650 * but for consistency with earlier versions of dtoa, it is optional
3651 * when MULTIPLE_THREADS is not defined.
3652 */
3653
3654 ZEND_API void
3655 #ifdef KR_headers
zend_freedtoa(s)3656 zend_freedtoa(s) char *s;
3657 #else
3658 zend_freedtoa(char *s)
3659 #endif
3660 {
3661 Bigint *b = (Bigint *)((int *)s - 1);
3662 b->maxwds = 1 << (b->k = *(int*)b);
3663 Bfree(b);
3664 #ifndef MULTIPLE_THREADS
3665 if (s == dtoa_result)
3666 dtoa_result = 0;
3667 #endif
3668 }
3669
3670 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
3671 *
3672 * Inspired by "How to Print Floating-Point Numbers Accurately" by
3673 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
3674 *
3675 * Modifications:
3676 * 1. Rather than iterating, we use a simple numeric overestimate
3677 * to determine k = floor(log10(d)). We scale relevant
3678 * quantities using O(log2(k)) rather than O(k) multiplications.
3679 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
3680 * try to generate digits strictly left to right. Instead, we
3681 * compute with fewer bits and propagate the carry if necessary
3682 * when rounding the final digit up. This is often faster.
3683 * 3. Under the assumption that input will be rounded nearest,
3684 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
3685 * That is, we allow equality in stopping tests when the
3686 * round-nearest rule will give the same floating-point value
3687 * as would satisfaction of the stopping test with strict
3688 * inequality.
3689 * 4. We remove common factors of powers of 2 from relevant
3690 * quantities.
3691 * 5. When converting floating-point integers less than 1e16,
3692 * we use floating-point arithmetic rather than resorting
3693 * to multiple-precision integers.
3694 * 6. When asked to produce fewer than 15 digits, we first try
3695 * to get by with floating-point arithmetic; we resort to
3696 * multiple-precision integer arithmetic only if we cannot
3697 * guarantee that the floating-point calculation has given
3698 * the correctly rounded result. For k requested digits and
3699 * "uniformly" distributed input, the probability is
3700 * something like 10^(k-15) that we must resort to the Long
3701 * calculation.
3702 */
3703
zend_dtoa(double dd,int mode,int ndigits,int * decpt,bool * sign,char ** rve)3704 ZEND_API char *zend_dtoa(double dd, int mode, int ndigits, int *decpt, bool *sign, char **rve)
3705 {
3706 /* Arguments ndigits, decpt, sign are similar to those
3707 of ecvt and fcvt; trailing zeros are suppressed from
3708 the returned string. If not null, *rve is set to point
3709 to the end of the return value. If d is +-Infinity or NaN,
3710 then *decpt is set to 9999.
3711
3712 mode:
3713 0 ==> shortest string that yields d when read in
3714 and rounded to nearest.
3715 1 ==> like 0, but with Steele & White stopping rule;
3716 e.g. with IEEE P754 arithmetic , mode 0 gives
3717 1e23 whereas mode 1 gives 9.999999999999999e22.
3718 2 ==> max(1,ndigits) significant digits. This gives a
3719 return value similar to that of ecvt, except
3720 that trailing zeros are suppressed.
3721 3 ==> through ndigits past the decimal point. This
3722 gives a return value similar to that from fcvt,
3723 except that trailing zeros are suppressed, and
3724 ndigits can be negative.
3725 4,5 ==> similar to 2 and 3, respectively, but (in
3726 round-nearest mode) with the tests of mode 0 to
3727 possibly return a shorter string that rounds to d.
3728 With IEEE arithmetic and compilation with
3729 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
3730 as modes 2 and 3 when FLT_ROUNDS != 1.
3731 6-9 ==> Debugging modes similar to mode - 4: don't try
3732 fast floating-point estimate (if applicable).
3733
3734 Values of mode other than 0-9 are treated as mode 0.
3735
3736 Sufficient space is allocated to the return value
3737 to hold the suppressed trailing zeros.
3738 */
3739
3740 int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1,
3741 j, j1 = 0, k, k0, k_check, leftright, m2, m5, s2, s5,
3742 spec_case = 0, try_quick;
3743 Long L;
3744 #ifndef Sudden_Underflow
3745 int denorm;
3746 ULong x;
3747 #endif
3748 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
3749 U d2, eps, u;
3750 double ds;
3751 char *s, *s0;
3752 #ifndef No_leftright
3753 #ifdef IEEE_Arith
3754 U eps1;
3755 #endif
3756 #endif
3757 #ifdef SET_INEXACT
3758 int inexact, oldinexact;
3759 #endif
3760 #ifdef Honor_FLT_ROUNDS /*{*/
3761 int Rounding;
3762 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
3763 Rounding = Flt_Rounds;
3764 #else /*}{*/
3765 Rounding = 1;
3766 switch(fegetround()) {
3767 case FE_TOWARDZERO: Rounding = 0; break;
3768 case FE_UPWARD: Rounding = 2; break;
3769 case FE_DOWNWARD: Rounding = 3;
3770 }
3771 #endif /*}}*/
3772 #endif /*}*/
3773
3774 #ifndef MULTIPLE_THREADS
3775 if (dtoa_result) {
3776 zend_freedtoa(dtoa_result);
3777 dtoa_result = 0;
3778 }
3779 #endif
3780
3781 u.d = dd;
3782 if (word0(&u) & Sign_bit) {
3783 /* set sign for everything, including 0's and NaNs */
3784 *sign = 1;
3785 word0(&u) &= ~Sign_bit; /* clear sign bit */
3786 }
3787 else
3788 *sign = 0;
3789
3790 #if defined(IEEE_Arith) + defined(VAX)
3791 #ifdef IEEE_Arith
3792 if ((word0(&u) & Exp_mask) == Exp_mask)
3793 #else
3794 if (word0(&u) == 0x8000)
3795 #endif
3796 {
3797 /* Infinity or NaN */
3798 *decpt = 9999;
3799 #ifdef IEEE_Arith
3800 if (!word1(&u) && !(word0(&u) & 0xfffff))
3801 return nrv_alloc("Infinity", rve, 8);
3802 #endif
3803 return nrv_alloc("NaN", rve, 3);
3804 }
3805 #endif
3806 #ifdef IBM
3807 dval(&u) += 0; /* normalize */
3808 #endif
3809 if (!dval(&u)) {
3810 *decpt = 1;
3811 return nrv_alloc("0", rve, 1);
3812 }
3813
3814 #ifdef SET_INEXACT
3815 try_quick = oldinexact = get_inexact();
3816 inexact = 1;
3817 #endif
3818 #ifdef Honor_FLT_ROUNDS
3819 if (Rounding >= 2) {
3820 if (*sign)
3821 Rounding = Rounding == 2 ? 0 : 2;
3822 else
3823 if (Rounding != 2)
3824 Rounding = 0;
3825 }
3826 #endif
3827
3828 b = d2b(&u, &be, &bbits);
3829 #ifdef Sudden_Underflow
3830 i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
3831 #else
3832 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
3833 #endif
3834 dval(&d2) = dval(&u);
3835 word0(&d2) &= Frac_mask1;
3836 word0(&d2) |= Exp_11;
3837 #ifdef IBM
3838 if (j = 11 - hi0bits(word0(&d2) & Frac_mask))
3839 dval(&d2) /= 1 << j;
3840 #endif
3841
3842 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
3843 * log10(x) = log(x) / log(10)
3844 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
3845 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
3846 *
3847 * This suggests computing an approximation k to log10(d) by
3848 *
3849 * k = (i - Bias)*0.301029995663981
3850 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
3851 *
3852 * We want k to be too large rather than too small.
3853 * The error in the first-order Taylor series approximation
3854 * is in our favor, so we just round up the constant enough
3855 * to compensate for any error in the multiplication of
3856 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
3857 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
3858 * adding 1e-13 to the constant term more than suffices.
3859 * Hence we adjust the constant term to 0.1760912590558.
3860 * (We could get a more accurate k by invoking log10,
3861 * but this is probably not worthwhile.)
3862 */
3863
3864 i -= Bias;
3865 #ifdef IBM
3866 i <<= 2;
3867 i += j;
3868 #endif
3869 #ifndef Sudden_Underflow
3870 denorm = 0;
3871 }
3872 else {
3873 /* d is denormalized */
3874
3875 i = bbits + be + (Bias + (P-1) - 1);
3876 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
3877 : word1(&u) << (32 - i);
3878 dval(&d2) = x;
3879 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
3880 i -= (Bias + (P-1) - 1) + 1;
3881 denorm = 1;
3882 }
3883 #endif
3884 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
3885 k = (int)ds;
3886 if (ds < 0. && ds != k)
3887 k--; /* want k = floor(ds) */
3888 k_check = 1;
3889 if (k >= 0 && k <= Ten_pmax) {
3890 if (dval(&u) < tens[k])
3891 k--;
3892 k_check = 0;
3893 }
3894 j = bbits - i - 1;
3895 if (j >= 0) {
3896 b2 = 0;
3897 s2 = j;
3898 }
3899 else {
3900 b2 = -j;
3901 s2 = 0;
3902 }
3903 if (k >= 0) {
3904 b5 = 0;
3905 s5 = k;
3906 s2 += k;
3907 }
3908 else {
3909 b2 -= k;
3910 b5 = -k;
3911 s5 = 0;
3912 }
3913 if (mode < 0 || mode > 9)
3914 mode = 0;
3915
3916 #ifndef SET_INEXACT
3917 #ifdef Check_FLT_ROUNDS
3918 try_quick = Rounding == 1;
3919 #else
3920 try_quick = 1;
3921 #endif
3922 #endif /*SET_INEXACT*/
3923
3924 if (mode > 5) {
3925 mode -= 4;
3926 try_quick = 0;
3927 }
3928 leftright = 1;
3929 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
3930 /* silence erroneous "gcc -Wall" warning. */
3931 switch(mode) {
3932 case 0:
3933 case 1:
3934 i = 18;
3935 ndigits = 0;
3936 break;
3937 case 2:
3938 leftright = 0;
3939 ZEND_FALLTHROUGH;
3940 case 4:
3941 if (ndigits <= 0)
3942 ndigits = 1;
3943 ilim = ilim1 = i = ndigits;
3944 break;
3945 case 3:
3946 leftright = 0;
3947 ZEND_FALLTHROUGH;
3948 case 5:
3949 i = ndigits + k + 1;
3950 ilim = i;
3951 ilim1 = i - 1;
3952 if (i <= 0)
3953 i = 1;
3954 }
3955 s = s0 = rv_alloc(i);
3956
3957 #ifdef Honor_FLT_ROUNDS
3958 if (mode > 1 && Rounding != 1)
3959 leftright = 0;
3960 #endif
3961
3962 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
3963
3964 /* Try to get by with floating-point arithmetic. */
3965
3966 i = 0;
3967 dval(&d2) = dval(&u);
3968 k0 = k;
3969 ilim0 = ilim;
3970 ieps = 2; /* conservative */
3971 if (k > 0) {
3972 ds = tens[k&0xf];
3973 j = k >> 4;
3974 if (j & Bletch) {
3975 /* prevent overflows */
3976 j &= Bletch - 1;
3977 dval(&u) /= bigtens[n_bigtens-1];
3978 ieps++;
3979 }
3980 for(; j; j >>= 1, i++)
3981 if (j & 1) {
3982 ieps++;
3983 ds *= bigtens[i];
3984 }
3985 dval(&u) /= ds;
3986 }
3987 else if ((j1 = -k)) {
3988 dval(&u) *= tens[j1 & 0xf];
3989 for(j = j1 >> 4; j; j >>= 1, i++)
3990 if (j & 1) {
3991 ieps++;
3992 dval(&u) *= bigtens[i];
3993 }
3994 }
3995 if (k_check && dval(&u) < 1. && ilim > 0) {
3996 if (ilim1 <= 0)
3997 goto fast_failed;
3998 ilim = ilim1;
3999 k--;
4000 dval(&u) *= 10.;
4001 ieps++;
4002 }
4003 dval(&eps) = ieps*dval(&u) + 7.;
4004 word0(&eps) -= (P-1)*Exp_msk1;
4005 if (ilim == 0) {
4006 S = mhi = 0;
4007 dval(&u) -= 5.;
4008 if (dval(&u) > dval(&eps))
4009 goto one_digit;
4010 if (dval(&u) < -dval(&eps))
4011 goto no_digits;
4012 goto fast_failed;
4013 }
4014 #ifndef No_leftright
4015 if (leftright) {
4016 /* Use Steele & White method of only
4017 * generating digits needed.
4018 */
4019 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
4020 #ifdef IEEE_Arith
4021 if (k0 < 0 && j1 >= 307) {
4022 eps1.d = 1.01e256; /* 1.01 allows roundoff in the next few lines */
4023 word0(&eps1) -= Exp_msk1 * (Bias+P-1);
4024 dval(&eps1) *= tens[j1 & 0xf];
4025 for(i = 0, j = (j1-256) >> 4; j; j >>= 1, i++)
4026 if (j & 1)
4027 dval(&eps1) *= bigtens[i];
4028 if (eps.d < eps1.d)
4029 eps.d = eps1.d;
4030 }
4031 #endif
4032 for(i = 0;;) {
4033 L = dval(&u);
4034 dval(&u) -= L;
4035 *s++ = '0' + (int)L;
4036 if (1. - dval(&u) < dval(&eps))
4037 goto bump_up;
4038 if (dval(&u) < dval(&eps))
4039 goto ret1;
4040 if (++i >= ilim)
4041 break;
4042 dval(&eps) *= 10.;
4043 dval(&u) *= 10.;
4044 }
4045 }
4046 else {
4047 #endif
4048 /* Generate ilim digits, then fix them up. */
4049 dval(&eps) *= tens[ilim-1];
4050 for(i = 1;; i++, dval(&u) *= 10.) {
4051 L = (Long)(dval(&u));
4052 if (!(dval(&u) -= L))
4053 ilim = i;
4054 *s++ = '0' + (int)L;
4055 if (i == ilim) {
4056 if (dval(&u) > 0.5 + dval(&eps))
4057 goto bump_up;
4058 else if (dval(&u) < 0.5 - dval(&eps)) {
4059 while(*--s == '0');
4060 s++;
4061 goto ret1;
4062 }
4063 break;
4064 }
4065 }
4066 #ifndef No_leftright
4067 }
4068 #endif
4069 fast_failed:
4070 s = s0;
4071 dval(&u) = dval(&d2);
4072 k = k0;
4073 ilim = ilim0;
4074 }
4075
4076 /* Do we have a "small" integer? */
4077
4078 if (be >= 0 && k <= Int_max) {
4079 /* Yes. */
4080 ds = tens[k];
4081 if (ndigits < 0 && ilim <= 0) {
4082 S = mhi = 0;
4083 if (ilim < 0 || dval(&u) <= 5*ds)
4084 goto no_digits;
4085 goto one_digit;
4086 }
4087 for(i = 1;; i++, dval(&u) *= 10.) {
4088 L = (Long)(dval(&u) / ds);
4089 dval(&u) -= L*ds;
4090 #ifdef Check_FLT_ROUNDS
4091 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
4092 if (dval(&u) < 0) {
4093 L--;
4094 dval(&u) += ds;
4095 }
4096 #endif
4097 *s++ = '0' + (int)L;
4098 if (!dval(&u)) {
4099 #ifdef SET_INEXACT
4100 inexact = 0;
4101 #endif
4102 break;
4103 }
4104 if (i == ilim) {
4105 #ifdef Honor_FLT_ROUNDS
4106 if (mode > 1)
4107 switch(Rounding) {
4108 case 0: goto ret1;
4109 case 2: goto bump_up;
4110 }
4111 #endif
4112 dval(&u) += dval(&u);
4113 #ifdef ROUND_BIASED
4114 if (dval(&u) >= ds)
4115 #else
4116 if (dval(&u) > ds || (dval(&u) == ds && L & 1))
4117 #endif
4118 {
4119 bump_up:
4120 while(*--s == '9')
4121 if (s == s0) {
4122 k++;
4123 *s = '0';
4124 break;
4125 }
4126 ++*s++;
4127 }
4128 break;
4129 }
4130 }
4131 goto ret1;
4132 }
4133
4134 m2 = b2;
4135 m5 = b5;
4136 mhi = mlo = 0;
4137 if (leftright) {
4138 i =
4139 #ifndef Sudden_Underflow
4140 denorm ? be + (Bias + (P-1) - 1 + 1) :
4141 #endif
4142 #ifdef IBM
4143 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
4144 #else
4145 1 + P - bbits;
4146 #endif
4147 b2 += i;
4148 s2 += i;
4149 mhi = i2b(1);
4150 }
4151 if (m2 > 0 && s2 > 0) {
4152 i = m2 < s2 ? m2 : s2;
4153 b2 -= i;
4154 m2 -= i;
4155 s2 -= i;
4156 }
4157 if (b5 > 0) {
4158 if (leftright) {
4159 if (m5 > 0) {
4160 mhi = pow5mult(mhi, m5);
4161 b1 = mult(mhi, b);
4162 Bfree(b);
4163 b = b1;
4164 }
4165 if ((j = b5 - m5))
4166 b = pow5mult(b, j);
4167 }
4168 else
4169 b = pow5mult(b, b5);
4170 }
4171 S = i2b(1);
4172 if (s5 > 0)
4173 S = pow5mult(S, s5);
4174
4175 /* Check for special case that d is a normalized power of 2. */
4176
4177 spec_case = 0;
4178 if ((mode < 2 || leftright)
4179 #ifdef Honor_FLT_ROUNDS
4180 && Rounding == 1
4181 #endif
4182 ) {
4183 if (!word1(&u) && !(word0(&u) & Bndry_mask)
4184 #ifndef Sudden_Underflow
4185 && word0(&u) & (Exp_mask & ~Exp_msk1)
4186 #endif
4187 ) {
4188 /* The special case */
4189 b2 += Log2P;
4190 s2 += Log2P;
4191 spec_case = 1;
4192 }
4193 }
4194
4195 /* Arrange for convenient computation of quotients:
4196 * shift left if necessary so divisor has 4 leading 0 bits.
4197 *
4198 * Perhaps we should just compute leading 28 bits of S once
4199 * and for all and pass them and a shift to quorem, so it
4200 * can do shifts and ORs to compute the numerator for q.
4201 */
4202 i = dshift(S, s2);
4203 b2 += i;
4204 m2 += i;
4205 s2 += i;
4206 if (b2 > 0)
4207 b = lshift(b, b2);
4208 if (s2 > 0)
4209 S = lshift(S, s2);
4210 if (k_check) {
4211 if (cmp(b,S) < 0) {
4212 k--;
4213 b = multadd(b, 10, 0); /* we botched the k estimate */
4214 if (leftright)
4215 mhi = multadd(mhi, 10, 0);
4216 ilim = ilim1;
4217 }
4218 }
4219 if (ilim <= 0 && (mode == 3 || mode == 5)) {
4220 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
4221 /* no digits, fcvt style */
4222 no_digits:
4223 k = -1 - ndigits;
4224 goto ret;
4225 }
4226 one_digit:
4227 *s++ = '1';
4228 k++;
4229 goto ret;
4230 }
4231 if (leftright) {
4232 if (m2 > 0)
4233 mhi = lshift(mhi, m2);
4234
4235 /* Compute mlo -- check for special case
4236 * that d is a normalized power of 2.
4237 */
4238
4239 mlo = mhi;
4240 if (spec_case) {
4241 mhi = Balloc(mhi->k);
4242 Bcopy(mhi, mlo);
4243 mhi = lshift(mhi, Log2P);
4244 }
4245
4246 for(i = 1;;i++) {
4247 dig = quorem(b,S) + '0';
4248 /* Do we yet have the shortest decimal string
4249 * that will round to d?
4250 */
4251 j = cmp(b, mlo);
4252 delta = diff(S, mhi);
4253 j1 = delta->sign ? 1 : cmp(b, delta);
4254 Bfree(delta);
4255 #ifndef ROUND_BIASED
4256 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
4257 #ifdef Honor_FLT_ROUNDS
4258 && Rounding >= 1
4259 #endif
4260 ) {
4261 if (dig == '9')
4262 goto round_9_up;
4263 if (j > 0)
4264 dig++;
4265 #ifdef SET_INEXACT
4266 else if (!b->x[0] && b->wds <= 1)
4267 inexact = 0;
4268 #endif
4269 *s++ = dig;
4270 goto ret;
4271 }
4272 #endif
4273 if (j < 0 || (j == 0 && mode != 1
4274 #ifndef ROUND_BIASED
4275 && !(word1(&u) & 1)
4276 #endif
4277 )) {
4278 if (!b->x[0] && b->wds <= 1) {
4279 #ifdef SET_INEXACT
4280 inexact = 0;
4281 #endif
4282 goto accept_dig;
4283 }
4284 #ifdef Honor_FLT_ROUNDS
4285 if (mode > 1)
4286 switch(Rounding) {
4287 case 0: goto accept_dig;
4288 case 2: goto keep_dig;
4289 }
4290 #endif /*Honor_FLT_ROUNDS*/
4291 if (j1 > 0) {
4292 b = lshift(b, 1);
4293 j1 = cmp(b, S);
4294 #ifdef ROUND_BIASED
4295 if (j1 >= 0 /*)*/
4296 #else
4297 if ((j1 > 0 || (j1 == 0 && dig & 1))
4298 #endif
4299 && dig++ == '9')
4300 goto round_9_up;
4301 }
4302 accept_dig:
4303 *s++ = dig;
4304 goto ret;
4305 }
4306 if (j1 > 0) {
4307 #ifdef Honor_FLT_ROUNDS
4308 if (!Rounding)
4309 goto accept_dig;
4310 #endif
4311 if (dig == '9') { /* possible if i == 1 */
4312 round_9_up:
4313 *s++ = '9';
4314 goto roundoff;
4315 }
4316 *s++ = dig + 1;
4317 goto ret;
4318 }
4319 #ifdef Honor_FLT_ROUNDS
4320 keep_dig:
4321 #endif
4322 *s++ = dig;
4323 if (i == ilim)
4324 break;
4325 b = multadd(b, 10, 0);
4326 if (mlo == mhi)
4327 mlo = mhi = multadd(mhi, 10, 0);
4328 else {
4329 mlo = multadd(mlo, 10, 0);
4330 mhi = multadd(mhi, 10, 0);
4331 }
4332 }
4333 }
4334 else
4335 for(i = 1;; i++) {
4336 *s++ = dig = quorem(b,S) + '0';
4337 if (!b->x[0] && b->wds <= 1) {
4338 #ifdef SET_INEXACT
4339 inexact = 0;
4340 #endif
4341 goto ret;
4342 }
4343 if (i >= ilim)
4344 break;
4345 b = multadd(b, 10, 0);
4346 }
4347
4348 /* Round off last digit */
4349
4350 #ifdef Honor_FLT_ROUNDS
4351 switch(Rounding) {
4352 case 0: goto trimzeros;
4353 case 2: goto roundoff;
4354 }
4355 #endif
4356 b = lshift(b, 1);
4357 j = cmp(b, S);
4358 #ifdef ROUND_BIASED
4359 if (j >= 0)
4360 #else
4361 if (j > 0 || (j == 0 && dig & 1))
4362 #endif
4363 {
4364 roundoff:
4365 while(*--s == '9')
4366 if (s == s0) {
4367 k++;
4368 *s++ = '1';
4369 goto ret;
4370 }
4371 ++*s++;
4372 }
4373 else {
4374 #ifdef Honor_FLT_ROUNDS
4375 trimzeros:
4376 #endif
4377 while(*--s == '0');
4378 s++;
4379 }
4380 ret:
4381 Bfree(S);
4382 if (mhi) {
4383 if (mlo && mlo != mhi)
4384 Bfree(mlo);
4385 Bfree(mhi);
4386 }
4387 ret1:
4388 #ifdef SET_INEXACT
4389 if (inexact) {
4390 if (!oldinexact) {
4391 word0(&u) = Exp_1 + (70 << Exp_shift);
4392 word1(&u) = 0;
4393 dval(&u) += 1.;
4394 }
4395 }
4396 else if (!oldinexact)
4397 clear_inexact();
4398 #endif
4399 Bfree(b);
4400 *s = 0;
4401 *decpt = k + 1;
4402 if (rve)
4403 *rve = s;
4404 return s0;
4405 }
4406
zend_hex_strtod(const char * str,const char ** endptr)4407 ZEND_API double zend_hex_strtod(const char *str, const char **endptr)
4408 {
4409 const char *s = str;
4410 char c;
4411 int any = 0;
4412 double value = 0;
4413
4414 if (*s == '0' && (s[1] == 'x' || s[1] == 'X')) {
4415 s += 2;
4416 }
4417
4418 while ((c = *s++)) {
4419 if (c >= '0' && c <= '9') {
4420 c -= '0';
4421 } else if (c >= 'A' && c <= 'F') {
4422 c -= 'A' - 10;
4423 } else if (c >= 'a' && c <= 'f') {
4424 c -= 'a' - 10;
4425 } else {
4426 break;
4427 }
4428
4429 any = 1;
4430 value = value * 16 + c;
4431 }
4432
4433 if (endptr != NULL) {
4434 *endptr = any ? s - 1 : str;
4435 }
4436
4437 return value;
4438 }
4439
zend_oct_strtod(const char * str,const char ** endptr)4440 ZEND_API double zend_oct_strtod(const char *str, const char **endptr)
4441 {
4442 const char *s = str;
4443 char c;
4444 double value = 0;
4445 int any = 0;
4446
4447 if (str[0] == '\0') {
4448 if (endptr != NULL) {
4449 *endptr = str;
4450 }
4451 return 0.0;
4452 }
4453
4454 while ((c = *s++)) {
4455 if (c < '0' || c > '7') {
4456 /* break and return the current value if the number is not well-formed
4457 * that's what Linux strtol() does
4458 */
4459 break;
4460 }
4461 value = value * 8 + c - '0';
4462 any = 1;
4463 }
4464
4465 if (endptr != NULL) {
4466 *endptr = any ? s - 1 : str;
4467 }
4468
4469 return value;
4470 }
4471
zend_bin_strtod(const char * str,const char ** endptr)4472 ZEND_API double zend_bin_strtod(const char *str, const char **endptr)
4473 {
4474 const char *s = str;
4475 char c;
4476 double value = 0;
4477 int any = 0;
4478
4479 if ('0' == *s && ('b' == s[1] || 'B' == s[1])) {
4480 s += 2;
4481 }
4482
4483 while ((c = *s++)) {
4484 /*
4485 * Verify the validity of the current character as a base-2 digit. In
4486 * the event that an invalid digit is found, halt the conversion and
4487 * return the portion which has been converted thus far.
4488 */
4489 if ('0' == c || '1' == c)
4490 value = value * 2 + c - '0';
4491 else
4492 break;
4493
4494 any = 1;
4495 }
4496
4497 /*
4498 * As with many strtoX implementations, should the subject sequence be
4499 * empty or not well-formed, no conversion is performed and the original
4500 * value of str is stored in *endptr, provided that endptr is not a null
4501 * pointer.
4502 */
4503 if (NULL != endptr) {
4504 *endptr = (char *)(any ? s - 1 : str);
4505 }
4506
4507 return value;
4508 }
4509
zend_gcvt(double value,int ndigit,char dec_point,char exponent,char * buf)4510 ZEND_API char *zend_gcvt(double value, int ndigit, char dec_point, char exponent, char *buf)
4511 {
4512 char *digits, *dst, *src;
4513 int i, decpt;
4514 bool sign;
4515 int mode = ndigit >= 0 ? 2 : 0;
4516
4517 if (mode == 0) {
4518 ndigit = 17;
4519 }
4520 digits = zend_dtoa(value, mode, ndigit, &decpt, &sign, NULL);
4521 if (decpt == 9999) {
4522 /*
4523 * Infinity or NaN, convert to inf or nan with sign.
4524 * We assume the buffer is at least ndigit long.
4525 */
4526 snprintf(buf, ndigit + 1, "%s%s", (sign && *digits == 'I') ? "-" : "", *digits == 'I' ? "INF" : "NAN");
4527 zend_freedtoa(digits);
4528 return (buf);
4529 }
4530
4531 dst = buf;
4532 if (sign) {
4533 *dst++ = '-';
4534 }
4535
4536 if ((decpt >= 0 && decpt > ndigit) || decpt < -3) { /* use E-style */
4537 /* exponential format (e.g. 1.2345e+13) */
4538 if (--decpt < 0) {
4539 sign = 1;
4540 decpt = -decpt;
4541 } else {
4542 sign = 0;
4543 }
4544 src = digits;
4545 *dst++ = *src++;
4546 *dst++ = dec_point;
4547 if (*src == '\0') {
4548 *dst++ = '0';
4549 } else {
4550 do {
4551 *dst++ = *src++;
4552 } while (*src != '\0');
4553 }
4554 *dst++ = exponent;
4555 if (sign) {
4556 *dst++ = '-';
4557 } else {
4558 *dst++ = '+';
4559 }
4560 if (decpt < 10) {
4561 *dst++ = '0' + decpt;
4562 *dst = '\0';
4563 } else {
4564 /* XXX - optimize */
4565 int n;
4566 for (n = decpt, i = 0; (n /= 10) != 0; i++);
4567 dst[i + 1] = '\0';
4568 while (decpt != 0) {
4569 dst[i--] = '0' + decpt % 10;
4570 decpt /= 10;
4571 }
4572 }
4573 } else if (decpt < 0) {
4574 /* standard format 0. */
4575 *dst++ = '0'; /* zero before decimal point */
4576 *dst++ = dec_point;
4577 do {
4578 *dst++ = '0';
4579 } while (++decpt < 0);
4580 src = digits;
4581 while (*src != '\0') {
4582 *dst++ = *src++;
4583 }
4584 *dst = '\0';
4585 } else {
4586 /* standard format */
4587 for (i = 0, src = digits; i < decpt; i++) {
4588 if (*src != '\0') {
4589 *dst++ = *src++;
4590 } else {
4591 *dst++ = '0';
4592 }
4593 }
4594 if (*src != '\0') {
4595 if (src == digits) {
4596 *dst++ = '0'; /* zero before decimal point */
4597 }
4598 *dst++ = dec_point;
4599 for (i = decpt; digits[i] != '\0'; i++) {
4600 *dst++ = digits[i];
4601 }
4602 }
4603 *dst = '\0';
4604 }
4605 zend_freedtoa(digits);
4606 return (buf);
4607 }
4608
destroy_freelist(void)4609 static void destroy_freelist(void)
4610 {
4611 int i;
4612 Bigint *tmp;
4613
4614 ACQUIRE_DTOA_LOCK(0)
4615 for (i = 0; i <= Kmax; i++) {
4616 Bigint **listp = &freelist[i];
4617 while ((tmp = *listp) != NULL) {
4618 *listp = tmp->next;
4619 free(tmp);
4620 }
4621 freelist[i] = NULL;
4622 }
4623 FREE_DTOA_LOCK(0)
4624 }
4625
free_p5s(void)4626 static void free_p5s(void)
4627 {
4628 Bigint **listp, *tmp;
4629
4630 ACQUIRE_DTOA_LOCK(1)
4631 listp = &p5s;
4632 while ((tmp = *listp) != NULL) {
4633 *listp = tmp->next;
4634 free(tmp);
4635 }
4636 FREE_DTOA_LOCK(1)
4637 }
4638
