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