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