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