1 : /****************************************************************
2 : *
3 : * The author of this software is David M. Gay.
4 : *
5 : * Copyright (c) 1991 by AT&T.
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 AT&T 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
21 : David M. Gay
22 : AT&T Bell Laboratories, Room 2C-463
23 : 600 Mountain Avenue
24 : Murray Hill, NJ 07974-2070
25 : U.S.A.
26 : dmg@research.att.com or research!dmg
27 : */
28 :
29 : /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
30 : *
31 : * This strtod returns a nearest machine number to the input decimal
32 : * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
33 : * broken by the IEEE round-even rule. Otherwise ties are broken by
34 : * biased rounding (add half and chop).
35 : *
36 : * Inspired loosely by William D. Clinger's paper "How to Read Floating
37 : * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
38 : *
39 : * Modifications:
40 : *
41 : * 1. We only require IEEE, IBM, or VAX double-precision
42 : * arithmetic (not IEEE double-extended).
43 : * 2. We get by with floating-point arithmetic in a case that
44 : * Clinger missed -- when we're computing d * 10^n
45 : * for a small integer d and the integer n is not too
46 : * much larger than 22 (the maximum integer k for which
47 : * we can represent 10^k exactly), we may be able to
48 : * compute (d*10^k) * 10^(e-k) with just one roundoff.
49 : * 3. Rather than a bit-at-a-time adjustment of the binary
50 : * result in the hard case, we use floating-point
51 : * arithmetic to determine the adjustment to within
52 : * one bit; only in really hard cases do we need to
53 : * compute a second residual.
54 : * 4. Because of 3., we don't need a large table of powers of 10
55 : * for ten-to-e (just some small tables, e.g. of 10^k
56 : * for 0 <= k <= 22).
57 : */
58 :
59 : /*
60 : * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least
61 : * significant byte has the lowest address.
62 : * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most
63 : * significant byte has the lowest address.
64 : * #define Long int on machines with 32-bit ints and 64-bit longs.
65 : * #define Sudden_Underflow for IEEE-format machines without gradual
66 : * underflow (i.e., that flush to zero on underflow).
67 : * #define IBM for IBM mainframe-style floating-point arithmetic.
68 : * #define VAX for VAX-style floating-point arithmetic.
69 : * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
70 : * #define No_leftright to omit left-right logic in fast floating-point
71 : * computation of dtoa.
72 : * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
73 : * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
74 : * that use extended-precision instructions to compute rounded
75 : * products and quotients) with IBM.
76 : * #define ROUND_BIASED for IEEE-format with biased rounding.
77 : * #define Inaccurate_Divide for IEEE-format with correctly rounded
78 : * products but inaccurate quotients, e.g., for Intel i860.
79 : * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision
80 : * integer arithmetic. Whether this speeds things up or slows things
81 : * down depends on the machine and the number being converted.
82 : * #define KR_headers for old-style C function headers.
83 : * #define Bad_float_h if your system lacks a float.h or if it does not
84 : * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
85 : * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
86 : * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
87 : * if memory is available and otherwise does something you deem
88 : * appropriate. If MALLOC is undefined, malloc will be invoked
89 : * directly -- and assumed always to succeed.
90 : */
91 :
92 : #include <zend_strtod.h>
93 :
94 : #ifdef ZTS
95 : #include <TSRM.h>
96 : #endif
97 :
98 : #include <stddef.h>
99 : #include <stdio.h>
100 : #include <ctype.h>
101 : #include <stdarg.h>
102 : #include <string.h>
103 : #include <stdlib.h>
104 : #include <math.h>
105 :
106 : #ifdef HAVE_LOCALE_H
107 : #include <locale.h>
108 : #endif
109 :
110 : #ifdef HAVE_SYS_TYPES_H
111 : #include <sys/types.h>
112 : #endif
113 :
114 : #if defined(HAVE_INTTYPES_H)
115 : #include <inttypes.h>
116 : #elif defined(HAVE_STDINT_H)
117 : #include <stdint.h>
118 : #endif
119 :
120 : #ifndef HAVE_INT32_T
121 : # if SIZEOF_INT == 4
122 : typedef int int32_t;
123 : # elif SIZEOF_LONG == 4
124 : typedef long int int32_t;
125 : # endif
126 : #endif
127 :
128 : #ifndef HAVE_UINT32_T
129 : # if SIZEOF_INT == 4
130 : typedef unsigned int uint32_t;
131 : # elif SIZEOF_LONG == 4
132 : typedef unsigned long int uint32_t;
133 : # endif
134 : #endif
135 :
136 : #ifdef WORDS_BIGENDIAN
137 : #define IEEE_BIG_ENDIAN
138 : #else
139 : #define IEEE_LITTLE_ENDIAN
140 : #endif
141 :
142 : #if defined(__arm__) && !defined(__VFP_FP__)
143 : /*
144 : * * Although the CPU is little endian the FP has different
145 : * * byte and word endianness. The byte order is still little endian
146 : * * but the word order is big endian.
147 : * */
148 : #define IEEE_BIG_ENDIAN
149 : #undef IEEE_LITTLE_ENDIAN
150 : #endif
151 :
152 : #ifdef __vax__
153 : #define VAX
154 : #endif
155 :
156 : #if defined(_MSC_VER)
157 : #define int32_t __int32
158 : #define uint32_t unsigned __int32
159 : #define IEEE_LITTLE_ENDIAN
160 : #endif
161 :
162 : #define Long int32_t
163 : #define ULong uint32_t
164 :
165 : #ifdef __cplusplus
166 : #include "malloc.h"
167 : #include "memory.h"
168 : #else
169 : #ifndef KR_headers
170 : #include "stdlib.h"
171 : #include "string.h"
172 : #include "locale.h"
173 : #else
174 : #include "malloc.h"
175 : #include "memory.h"
176 : #endif
177 : #endif
178 :
179 : #ifdef MALLOC
180 : #ifdef KR_headers
181 : extern char *MALLOC();
182 : #else
183 : extern void *MALLOC(size_t);
184 : #endif
185 : #else
186 : #define MALLOC malloc
187 : #endif
188 :
189 : #include "ctype.h"
190 : #include "errno.h"
191 :
192 : #ifdef Bad_float_h
193 : #ifdef IEEE_BIG_ENDIAN
194 : #define IEEE_ARITHMETIC
195 : #endif
196 : #ifdef IEEE_LITTLE_ENDIAN
197 : #define IEEE_ARITHMETIC
198 : #endif
199 :
200 : #ifdef IEEE_ARITHMETIC
201 : #define DBL_DIG 15
202 : #define DBL_MAX_10_EXP 308
203 : #define DBL_MAX_EXP 1024
204 : #define FLT_RADIX 2
205 : #define FLT_ROUNDS 1
206 : #define DBL_MAX 1.7976931348623157e+308
207 : #endif
208 :
209 : #ifdef IBM
210 : #define DBL_DIG 16
211 : #define DBL_MAX_10_EXP 75
212 : #define DBL_MAX_EXP 63
213 : #define FLT_RADIX 16
214 : #define FLT_ROUNDS 0
215 : #define DBL_MAX 7.2370055773322621e+75
216 : #endif
217 :
218 : #ifdef VAX
219 : #define DBL_DIG 16
220 : #define DBL_MAX_10_EXP 38
221 : #define DBL_MAX_EXP 127
222 : #define FLT_RADIX 2
223 : #define FLT_ROUNDS 1
224 : #define DBL_MAX 1.7014118346046923e+38
225 : #endif
226 :
227 :
228 : #ifndef LONG_MAX
229 : #define LONG_MAX 2147483647
230 : #endif
231 : #else
232 : #include "float.h"
233 : #endif
234 : #ifndef __MATH_H__
235 : #include "math.h"
236 : #endif
237 :
238 : BEGIN_EXTERN_C()
239 :
240 : #ifndef CONST
241 : #ifdef KR_headers
242 : #define CONST /* blank */
243 : #else
244 : #define CONST const
245 : #endif
246 : #endif
247 :
248 : #ifdef Unsigned_Shifts
249 : #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
250 : #else
251 : #define Sign_Extend(a,b) /*no-op*/
252 : #endif
253 :
254 : #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + \
255 : defined(IBM) != 1
256 : Exactly one of IEEE_LITTLE_ENDIAN IEEE_BIG_ENDIAN, VAX, or
257 : IBM should be defined.
258 : #endif
259 :
260 : typedef union {
261 : double d;
262 : ULong ul[2];
263 : } _double;
264 : #define value(x) ((x).d)
265 : #ifdef IEEE_LITTLE_ENDIAN
266 : #define word0(x) ((x).ul[1])
267 : #define word1(x) ((x).ul[0])
268 : #else
269 : #define word0(x) ((x).ul[0])
270 : #define word1(x) ((x).ul[1])
271 : #endif
272 :
273 : /* The following definition of Storeinc is appropriate for MIPS processors.
274 : * * An alternative that might be better on some machines is
275 : * * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
276 : * */
277 : #if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm__)
278 : #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
279 : ((unsigned short *)a)[0] = (unsigned short)c, a++)
280 : #else
281 : #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
282 : ((unsigned short *)a)[1] = (unsigned short)c, a++)
283 : #endif
284 :
285 : /* #define P DBL_MANT_DIG */
286 : /* Ten_pmax = floor(P*log(2)/log(5)) */
287 : /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
288 : /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
289 : /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
290 :
291 : #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN)
292 : #define Exp_shift 20
293 : #define Exp_shift1 20
294 : #define Exp_msk1 0x100000
295 : #define Exp_msk11 0x100000
296 : #define Exp_mask 0x7ff00000
297 : #define P 53
298 : #define Bias 1023
299 : #define IEEE_Arith
300 : #define Emin (-1022)
301 : #define Exp_1 0x3ff00000
302 : #define Exp_11 0x3ff00000
303 : #define Ebits 11
304 : #define Frac_mask 0xfffff
305 : #define Frac_mask1 0xfffff
306 : #define Ten_pmax 22
307 : #define Bletch 0x10
308 : #define Bndry_mask 0xfffff
309 : #define Bndry_mask1 0xfffff
310 : #define LSB 1
311 : #define Sign_bit 0x80000000
312 : #define Log2P 1
313 : #define Tiny0 0
314 : #define Tiny1 1
315 : #define Quick_max 14
316 : #define Int_max 14
317 : #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
318 : #else
319 : #undef Sudden_Underflow
320 : #define Sudden_Underflow
321 : #ifdef IBM
322 : #define Exp_shift 24
323 : #define Exp_shift1 24
324 : #define Exp_msk1 0x1000000
325 : #define Exp_msk11 0x1000000
326 : #define Exp_mask 0x7f000000
327 : #define P 14
328 : #define Bias 65
329 : #define Exp_1 0x41000000
330 : #define Exp_11 0x41000000
331 : #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
332 : #define Frac_mask 0xffffff
333 : #define Frac_mask1 0xffffff
334 : #define Bletch 4
335 : #define Ten_pmax 22
336 : #define Bndry_mask 0xefffff
337 : #define Bndry_mask1 0xffffff
338 : #define LSB 1
339 : #define Sign_bit 0x80000000
340 : #define Log2P 4
341 : #define Tiny0 0x100000
342 : #define Tiny1 0
343 : #define Quick_max 14
344 : #define Int_max 15
345 : #else /* VAX */
346 : #define Exp_shift 23
347 : #define Exp_shift1 7
348 : #define Exp_msk1 0x80
349 : #define Exp_msk11 0x800000
350 : #define Exp_mask 0x7f80
351 : #define P 56
352 : #define Bias 129
353 : #define Exp_1 0x40800000
354 : #define Exp_11 0x4080
355 : #define Ebits 8
356 : #define Frac_mask 0x7fffff
357 : #define Frac_mask1 0xffff007f
358 : #define Ten_pmax 24
359 : #define Bletch 2
360 : #define Bndry_mask 0xffff007f
361 : #define Bndry_mask1 0xffff007f
362 : #define LSB 0x10000
363 : #define Sign_bit 0x8000
364 : #define Log2P 1
365 : #define Tiny0 0x80
366 : #define Tiny1 0
367 : #define Quick_max 15
368 : #define Int_max 15
369 : #endif
370 : #endif
371 :
372 : #ifndef IEEE_Arith
373 : #define ROUND_BIASED
374 : #endif
375 :
376 : #ifdef RND_PRODQUOT
377 : #define rounded_product(a,b) a = rnd_prod(a, b)
378 : #define rounded_quotient(a,b) a = rnd_quot(a, b)
379 : #ifdef KR_headers
380 : extern double rnd_prod(), rnd_quot();
381 : #else
382 : extern double rnd_prod(double, double), rnd_quot(double, double);
383 : #endif
384 : #else
385 : #define rounded_product(a,b) a *= b
386 : #define rounded_quotient(a,b) a /= b
387 : #endif
388 :
389 : #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
390 : #define Big1 0xffffffff
391 :
392 : #ifndef Just_16
393 : /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
394 : * * This makes some inner loops simpler and sometimes saves work
395 : * * during multiplications, but it often seems to make things slightly
396 : * * slower. Hence the default is now to store 32 bits per Long.
397 : * */
398 : #ifndef Pack_32
399 : #define Pack_32
400 : #endif
401 : #endif
402 :
403 : #define Kmax 15
404 :
405 : struct Bigint {
406 : struct Bigint *next;
407 : int k, maxwds, sign, wds;
408 : ULong x[1];
409 : };
410 :
411 : typedef struct Bigint Bigint;
412 :
413 : /* static variables, multithreading fun! */
414 : static Bigint *freelist[Kmax+1];
415 : static Bigint *p5s;
416 :
417 : static void destroy_freelist(void);
418 :
419 : #ifdef ZTS
420 :
421 : static MUTEX_T dtoa_mutex;
422 : static MUTEX_T pow5mult_mutex;
423 :
424 : #define _THREAD_PRIVATE_MUTEX_LOCK(x) tsrm_mutex_lock(x);
425 : #define _THREAD_PRIVATE_MUTEX_UNLOCK(x) tsrm_mutex_unlock(x);
426 :
427 : #else
428 :
429 : #define _THREAD_PRIVATE_MUTEX_LOCK(x)
430 : #define _THREAD_PRIVATE_MUTEX_UNLOCK(x)
431 :
432 : #endif /* ZTS */
433 :
434 : ZEND_API int zend_startup_strtod(void) /* {{{ */
435 220 : {
436 : #ifdef ZTS
437 : dtoa_mutex = tsrm_mutex_alloc();
438 : pow5mult_mutex = tsrm_mutex_alloc();
439 : #endif
440 220 : return 1;
441 : }
442 : /* }}} */
443 : ZEND_API int zend_shutdown_strtod(void) /* {{{ */
444 219 : {
445 219 : destroy_freelist();
446 : #ifdef ZTS
447 : tsrm_mutex_free(dtoa_mutex);
448 : dtoa_mutex = NULL;
449 :
450 : tsrm_mutex_free(pow5mult_mutex);
451 : pow5mult_mutex = NULL;
452 : #endif
453 219 : return 1;
454 : }
455 : /* }}} */
456 :
457 : static Bigint * Balloc(int k)
458 317 : {
459 : int x;
460 : Bigint *rv;
461 :
462 317 : if (k > Kmax) {
463 0 : zend_error(E_ERROR, "Balloc() allocation exceeds list boundary");
464 : }
465 :
466 : _THREAD_PRIVATE_MUTEX_LOCK(dtoa_mutex);
467 317 : if ((rv = freelist[k])) {
468 159 : freelist[k] = rv->next;
469 : } else {
470 158 : x = 1 << k;
471 158 : rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(Long));
472 158 : if (!rv) {
473 : _THREAD_PRIVATE_MUTEX_UNLOCK(dtoa_mutex);
474 0 : zend_error(E_ERROR, "Balloc() failed to allocate memory");
475 : }
476 158 : rv->k = k;
477 158 : rv->maxwds = x;
478 : }
479 : _THREAD_PRIVATE_MUTEX_UNLOCK(dtoa_mutex);
480 317 : rv->sign = rv->wds = 0;
481 317 : return rv;
482 : }
483 :
484 : static void Bfree(Bigint *v)
485 317 : {
486 317 : if (v) {
487 : _THREAD_PRIVATE_MUTEX_LOCK(dtoa_mutex);
488 317 : v->next = freelist[v->k];
489 317 : freelist[v->k] = v;
490 : _THREAD_PRIVATE_MUTEX_UNLOCK(dtoa_mutex);
491 : }
492 317 : }
493 :
494 : #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
495 : y->wds*sizeof(Long) + 2*sizeof(int))
496 :
497 : /* return value is only used as a simple string, so mis-aligned parts
498 : * inside the Bigint are not at risk on strict align architectures
499 : */
500 162 : static char * rv_alloc(int i) {
501 : int j, k, *r;
502 :
503 162 : j = sizeof(ULong);
504 162 : for(k = 0;
505 324 : sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
506 0 : j <<= 1) {
507 0 : k++;
508 : }
509 162 : r = (int*)Balloc(k);
510 162 : *r = k;
511 162 : return (char *)(r+1);
512 : }
513 :
514 :
515 : static char * nrv_alloc(char *s, char **rve, int n)
516 7 : {
517 : char *rv, *t;
518 :
519 7 : t = rv = rv_alloc(n);
520 21 : while((*t = *s++) !=0) {
521 7 : t++;
522 : }
523 7 : if (rve) {
524 0 : *rve = t;
525 : }
526 7 : return rv;
527 : }
528 :
529 : static Bigint * multadd(Bigint *b, int m, int a) /* multiply by m and add a */
530 0 : {
531 : int i, wds;
532 : ULong *x, y;
533 : #ifdef Pack_32
534 : ULong xi, z;
535 : #endif
536 : Bigint *b1;
537 :
538 0 : wds = b->wds;
539 0 : x = b->x;
540 0 : i = 0;
541 : do {
542 : #ifdef Pack_32
543 0 : xi = *x;
544 0 : y = (xi & 0xffff) * m + a;
545 0 : z = (xi >> 16) * m + (y >> 16);
546 0 : a = (int)(z >> 16);
547 0 : *x++ = (z << 16) + (y & 0xffff);
548 : #else
549 : y = *x * m + a;
550 : a = (int)(y >> 16);
551 : *x++ = y & 0xffff;
552 : #endif
553 : }
554 0 : while(++i < wds);
555 0 : if (a) {
556 0 : if (wds >= b->maxwds) {
557 0 : b1 = Balloc(b->k+1);
558 0 : Bcopy(b1, b);
559 0 : Bfree(b);
560 0 : b = b1;
561 : }
562 0 : b->x[wds++] = a;
563 0 : b->wds = wds;
564 : }
565 0 : return b;
566 : }
567 :
568 : static int hi0bits(ULong x)
569 0 : {
570 0 : int k = 0;
571 :
572 0 : if (!(x & 0xffff0000)) {
573 0 : k = 16;
574 0 : x <<= 16;
575 : }
576 0 : if (!(x & 0xff000000)) {
577 0 : k += 8;
578 0 : x <<= 8;
579 : }
580 0 : if (!(x & 0xf0000000)) {
581 0 : k += 4;
582 0 : x <<= 4;
583 : }
584 0 : if (!(x & 0xc0000000)) {
585 0 : k += 2;
586 0 : x <<= 2;
587 : }
588 0 : if (!(x & 0x80000000)) {
589 0 : k++;
590 0 : if (!(x & 0x40000000)) {
591 0 : return 32;
592 : }
593 : }
594 0 : return k;
595 : }
596 :
597 : static int lo0bits(ULong *y)
598 155 : {
599 : int k;
600 155 : ULong x = *y;
601 :
602 155 : if (x & 7) {
603 101 : if (x & 1) {
604 13 : return 0;
605 : }
606 88 : if (x & 2) {
607 88 : *y = x >> 1;
608 88 : return 1;
609 : }
610 0 : *y = x >> 2;
611 0 : return 2;
612 : }
613 54 : k = 0;
614 54 : if (!(x & 0xffff)) {
615 34 : k = 16;
616 34 : x >>= 16;
617 : }
618 54 : if (!(x & 0xff)) {
619 16 : k += 8;
620 16 : x >>= 8;
621 : }
622 54 : if (!(x & 0xf)) {
623 14 : k += 4;
624 14 : x >>= 4;
625 : }
626 54 : if (!(x & 0x3)) {
627 50 : k += 2;
628 50 : x >>= 2;
629 : }
630 54 : if (!(x & 1)) {
631 17 : k++;
632 17 : x >>= 1;
633 17 : if (!(x & 1)) {
634 0 : return 32;
635 : }
636 : }
637 54 : *y = x;
638 54 : return k;
639 : }
640 :
641 : static Bigint * i2b(int i)
642 0 : {
643 : Bigint *b;
644 :
645 0 : b = Balloc(1);
646 0 : b->x[0] = i;
647 0 : b->wds = 1;
648 0 : return b;
649 : }
650 :
651 : static Bigint * mult(Bigint *a, Bigint *b)
652 0 : {
653 : Bigint *c;
654 : int k, wa, wb, wc;
655 : ULong carry, y, z;
656 : ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
657 : #ifdef Pack_32
658 : ULong z2;
659 : #endif
660 :
661 0 : if (a->wds < b->wds) {
662 0 : c = a;
663 0 : a = b;
664 0 : b = c;
665 : }
666 0 : k = a->k;
667 0 : wa = a->wds;
668 0 : wb = b->wds;
669 0 : wc = wa + wb;
670 0 : if (wc > a->maxwds) {
671 0 : k++;
672 : }
673 0 : c = Balloc(k);
674 0 : for(x = c->x, xa = x + wc; x < xa; x++) {
675 0 : *x = 0;
676 : }
677 0 : xa = a->x;
678 0 : xae = xa + wa;
679 0 : xb = b->x;
680 0 : xbe = xb + wb;
681 0 : xc0 = c->x;
682 : #ifdef Pack_32
683 0 : for(; xb < xbe; xb++, xc0++) {
684 0 : if ((y = *xb & 0xffff)) {
685 0 : x = xa;
686 0 : xc = xc0;
687 0 : carry = 0;
688 : do {
689 0 : z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
690 0 : carry = z >> 16;
691 0 : z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
692 0 : carry = z2 >> 16;
693 0 : Storeinc(xc, z2, z);
694 : }
695 0 : while(x < xae);
696 0 : *xc = carry;
697 : }
698 0 : if ((y = *xb >> 16)) {
699 0 : x = xa;
700 0 : xc = xc0;
701 0 : carry = 0;
702 0 : z2 = *xc;
703 : do {
704 0 : z = (*x & 0xffff) * y + (*xc >> 16) + carry;
705 0 : carry = z >> 16;
706 0 : Storeinc(xc, z, z2);
707 0 : z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
708 0 : carry = z2 >> 16;
709 : }
710 0 : while(x < xae);
711 0 : *xc = z2;
712 : }
713 : }
714 : #else
715 : for(; xb < xbe; xc0++) {
716 : if (y = *xb++) {
717 : x = xa;
718 : xc = xc0;
719 : carry = 0;
720 : do {
721 : z = *x++ * y + *xc + carry;
722 : carry = z >> 16;
723 : *xc++ = z & 0xffff;
724 : }
725 : while(x < xae);
726 : *xc = carry;
727 : }
728 : }
729 : #endif
730 0 : for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
731 0 : c->wds = wc;
732 0 : return c;
733 : }
734 :
735 : static Bigint * s2b (CONST char *s, int nd0, int nd, ULong y9)
736 0 : {
737 : Bigint *b;
738 : int i, k;
739 : Long x, y;
740 :
741 0 : x = (nd + 8) / 9;
742 0 : for(k = 0, y = 1; x > y; y <<= 1, k++) ;
743 : #ifdef Pack_32
744 0 : b = Balloc(k);
745 0 : b->x[0] = y9;
746 0 : b->wds = 1;
747 : #else
748 : b = Balloc(k+1);
749 : b->x[0] = y9 & 0xffff;
750 : b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
751 : #endif
752 :
753 0 : i = 9;
754 0 : if (9 < nd0) {
755 0 : s += 9;
756 0 : do b = multadd(b, 10, *s++ - '0');
757 0 : while(++i < nd0);
758 0 : s++;
759 : } else {
760 0 : s += 10;
761 : }
762 0 : for(; i < nd; i++) {
763 0 : b = multadd(b, 10, *s++ - '0');
764 : }
765 0 : return b;
766 : }
767 :
768 : static Bigint * pow5mult(Bigint *b, int k)
769 0 : {
770 : Bigint *b1, *p5, *p51;
771 : int i;
772 : static int p05[3] = { 5, 25, 125 };
773 :
774 : _THREAD_PRIVATE_MUTEX_LOCK(pow5mult_mutex);
775 0 : if ((i = k & 3)) {
776 0 : b = multadd(b, p05[i-1], 0);
777 : }
778 :
779 0 : if (!(k >>= 2)) {
780 : _THREAD_PRIVATE_MUTEX_UNLOCK(pow5mult_mutex);
781 0 : return b;
782 : }
783 0 : if (!(p5 = p5s)) {
784 : /* first time */
785 0 : p5 = p5s = i2b(625);
786 0 : p5->next = 0;
787 : }
788 : for(;;) {
789 0 : if (k & 1) {
790 0 : b1 = mult(b, p5);
791 0 : Bfree(b);
792 0 : b = b1;
793 : }
794 0 : if (!(k >>= 1)) {
795 0 : break;
796 : }
797 0 : if (!(p51 = p5->next)) {
798 0 : if (!(p51 = p5->next)) {
799 0 : p51 = p5->next = mult(p5,p5);
800 0 : p51->next = 0;
801 : }
802 : }
803 0 : p5 = p51;
804 0 : }
805 : _THREAD_PRIVATE_MUTEX_UNLOCK(pow5mult_mutex);
806 0 : return b;
807 : }
808 :
809 :
810 : static Bigint *lshift(Bigint *b, int k)
811 0 : {
812 : int i, k1, n, n1;
813 : Bigint *b1;
814 : ULong *x, *x1, *xe, z;
815 :
816 : #ifdef Pack_32
817 0 : n = k >> 5;
818 : #else
819 : n = k >> 4;
820 : #endif
821 0 : k1 = b->k;
822 0 : n1 = n + b->wds + 1;
823 0 : for(i = b->maxwds; n1 > i; i <<= 1) {
824 0 : k1++;
825 : }
826 0 : b1 = Balloc(k1);
827 0 : x1 = b1->x;
828 0 : for(i = 0; i < n; i++) {
829 0 : *x1++ = 0;
830 : }
831 0 : x = b->x;
832 0 : xe = x + b->wds;
833 : #ifdef Pack_32
834 0 : if (k &= 0x1f) {
835 0 : k1 = 32 - k;
836 0 : z = 0;
837 : do {
838 0 : *x1++ = *x << k | z;
839 0 : z = *x++ >> k1;
840 : }
841 0 : while(x < xe);
842 0 : if ((*x1 = z)) {
843 0 : ++n1;
844 : }
845 : }
846 : #else
847 : if (k &= 0xf) {
848 : k1 = 16 - k;
849 : z = 0;
850 : do {
851 : *x1++ = *x << k & 0xffff | z;
852 : z = *x++ >> k1;
853 : }
854 : while(x < xe);
855 : if (*x1 = z) {
856 : ++n1;
857 : }
858 : }
859 : #endif
860 : else do
861 0 : *x1++ = *x++;
862 0 : while(x < xe);
863 0 : b1->wds = n1 - 1;
864 0 : Bfree(b);
865 0 : return b1;
866 : }
867 :
868 : static int cmp(Bigint *a, Bigint *b)
869 0 : {
870 : ULong *xa, *xa0, *xb, *xb0;
871 : int i, j;
872 :
873 0 : i = a->wds;
874 0 : j = b->wds;
875 : #ifdef DEBUG
876 : if (i > 1 && !a->x[i-1])
877 : Bug("cmp called with a->x[a->wds-1] == 0");
878 : if (j > 1 && !b->x[j-1])
879 : Bug("cmp called with b->x[b->wds-1] == 0");
880 : #endif
881 0 : if (i -= j)
882 0 : return i;
883 0 : xa0 = a->x;
884 0 : xa = xa0 + j;
885 0 : xb0 = b->x;
886 0 : xb = xb0 + j;
887 : for(;;) {
888 0 : if (*--xa != *--xb)
889 0 : return *xa < *xb ? -1 : 1;
890 0 : if (xa <= xa0)
891 0 : break;
892 0 : }
893 0 : return 0;
894 : }
895 :
896 :
897 : static Bigint * diff(Bigint *a, Bigint *b)
898 0 : {
899 : Bigint *c;
900 : int i, wa, wb;
901 : Long borrow, y; /* We need signed shifts here. */
902 : ULong *xa, *xae, *xb, *xbe, *xc;
903 : #ifdef Pack_32
904 : Long z;
905 : #endif
906 :
907 0 : i = cmp(a,b);
908 0 : if (!i) {
909 0 : c = Balloc(0);
910 0 : c->wds = 1;
911 0 : c->x[0] = 0;
912 0 : return c;
913 : }
914 0 : if (i < 0) {
915 0 : c = a;
916 0 : a = b;
917 0 : b = c;
918 0 : i = 1;
919 : } else {
920 0 : i = 0;
921 : }
922 0 : c = Balloc(a->k);
923 0 : c->sign = i;
924 0 : wa = a->wds;
925 0 : xa = a->x;
926 0 : xae = xa + wa;
927 0 : wb = b->wds;
928 0 : xb = b->x;
929 0 : xbe = xb + wb;
930 0 : xc = c->x;
931 0 : borrow = 0;
932 : #ifdef Pack_32
933 : do {
934 0 : y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
935 0 : borrow = y >> 16;
936 : Sign_Extend(borrow, y);
937 0 : z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
938 0 : borrow = z >> 16;
939 : Sign_Extend(borrow, z);
940 0 : Storeinc(xc, z, y);
941 0 : } while(xb < xbe);
942 0 : while(xa < xae) {
943 0 : y = (*xa & 0xffff) + borrow;
944 0 : borrow = y >> 16;
945 : Sign_Extend(borrow, y);
946 0 : z = (*xa++ >> 16) + borrow;
947 0 : borrow = z >> 16;
948 : Sign_Extend(borrow, z);
949 0 : Storeinc(xc, z, y);
950 : }
951 : #else
952 : do {
953 : y = *xa++ - *xb++ + borrow;
954 : borrow = y >> 16;
955 : Sign_Extend(borrow, y);
956 : *xc++ = y & 0xffff;
957 : } while(xb < xbe);
958 : while(xa < xae) {
959 : y = *xa++ + borrow;
960 : borrow = y >> 16;
961 : Sign_Extend(borrow, y);
962 : *xc++ = y & 0xffff;
963 : }
964 : #endif
965 0 : while(!*--xc) {
966 0 : wa--;
967 : }
968 0 : c->wds = wa;
969 0 : return c;
970 : }
971 :
972 : static double ulp (double _x)
973 0 : {
974 : _double x;
975 : register Long L;
976 : _double a;
977 :
978 0 : value(x) = _x;
979 0 : L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
980 : #ifndef Sudden_Underflow
981 0 : if (L > 0) {
982 : #endif
983 : #ifdef IBM
984 : L |= Exp_msk1 >> 4;
985 : #endif
986 0 : word0(a) = L;
987 0 : word1(a) = 0;
988 : #ifndef Sudden_Underflow
989 : }
990 : else {
991 0 : L = -L >> Exp_shift;
992 0 : if (L < Exp_shift) {
993 0 : word0(a) = 0x80000 >> L;
994 0 : word1(a) = 0;
995 : }
996 : else {
997 0 : word0(a) = 0;
998 0 : L -= Exp_shift;
999 0 : word1(a) = L >= 31 ? 1 : 1 << (31 - L);
1000 : }
1001 : }
1002 : #endif
1003 0 : return value(a);
1004 : }
1005 :
1006 : static double
1007 : b2d
1008 : #ifdef KR_headers
1009 : (a, e) Bigint *a; int *e;
1010 : #else
1011 : (Bigint *a, int *e)
1012 : #endif
1013 0 : {
1014 : ULong *xa, *xa0, w, y, z;
1015 : int k;
1016 : _double d;
1017 : #ifdef VAX
1018 : ULong d0, d1;
1019 : #else
1020 : #define d0 word0(d)
1021 : #define d1 word1(d)
1022 : #endif
1023 :
1024 0 : xa0 = a->x;
1025 0 : xa = xa0 + a->wds;
1026 0 : y = *--xa;
1027 : #ifdef DEBUG
1028 : if (!y) Bug("zero y in b2d");
1029 : #endif
1030 0 : k = hi0bits(y);
1031 0 : *e = 32 - k;
1032 : #ifdef Pack_32
1033 0 : if (k < Ebits) {
1034 0 : d0 = Exp_1 | y >> (Ebits - k);
1035 0 : w = xa > xa0 ? *--xa : 0;
1036 0 : d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1037 0 : goto ret_d;
1038 : }
1039 0 : z = xa > xa0 ? *--xa : 0;
1040 0 : if (k -= Ebits) {
1041 0 : d0 = Exp_1 | y << k | z >> (32 - k);
1042 0 : y = xa > xa0 ? *--xa : 0;
1043 0 : d1 = z << k | y >> (32 - k);
1044 : }
1045 : else {
1046 0 : d0 = Exp_1 | y;
1047 0 : d1 = z;
1048 : }
1049 : #else
1050 : if (k < Ebits + 16) {
1051 : z = xa > xa0 ? *--xa : 0;
1052 : d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1053 : w = xa > xa0 ? *--xa : 0;
1054 : y = xa > xa0 ? *--xa : 0;
1055 : d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1056 : goto ret_d;
1057 : }
1058 : z = xa > xa0 ? *--xa : 0;
1059 : w = xa > xa0 ? *--xa : 0;
1060 : k -= Ebits + 16;
1061 : d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1062 : y = xa > xa0 ? *--xa : 0;
1063 : d1 = w << k + 16 | y << k;
1064 : #endif
1065 0 : ret_d:
1066 : #ifdef VAX
1067 : word0(d) = d0 >> 16 | d0 << 16;
1068 : word1(d) = d1 >> 16 | d1 << 16;
1069 : #else
1070 : #undef d0
1071 : #undef d1
1072 : #endif
1073 0 : return value(d);
1074 : }
1075 :
1076 :
1077 : static Bigint * d2b(double _d, int *e, int *bits)
1078 155 : {
1079 : Bigint *b;
1080 : int de, i, k;
1081 : ULong *x, y, z;
1082 : _double d;
1083 : #ifdef VAX
1084 : ULong d0, d1;
1085 : #endif
1086 :
1087 155 : value(d) = _d;
1088 : #ifdef VAX
1089 : d0 = word0(d) >> 16 | word0(d) << 16;
1090 : d1 = word1(d) >> 16 | word1(d) << 16;
1091 : #else
1092 : #define d0 word0(d)
1093 : #define d1 word1(d)
1094 : #endif
1095 :
1096 : #ifdef Pack_32
1097 155 : b = Balloc(1);
1098 : #else
1099 : b = Balloc(2);
1100 : #endif
1101 155 : x = b->x;
1102 :
1103 155 : z = d0 & Frac_mask;
1104 155 : d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1105 : #ifdef Sudden_Underflow
1106 : de = (int)(d0 >> Exp_shift);
1107 : #ifndef IBM
1108 : z |= Exp_msk11;
1109 : #endif
1110 : #else
1111 155 : if ((de = (int)(d0 >> Exp_shift)))
1112 155 : z |= Exp_msk1;
1113 : #endif
1114 : #ifdef Pack_32
1115 155 : if ((y = d1)) {
1116 115 : if ((k = lo0bits(&y))) {
1117 102 : x[0] = y | (z << (32 - k));
1118 102 : z >>= k;
1119 : } else {
1120 13 : x[0] = y;
1121 : }
1122 115 : i = b->wds = (x[1] = z) ? 2 : 1;
1123 : } else {
1124 : #ifdef DEBUG
1125 : if (!z)
1126 : Bug("Zero passed to d2b");
1127 : #endif
1128 40 : k = lo0bits(&z);
1129 40 : x[0] = z;
1130 40 : i = b->wds = 1;
1131 40 : k += 32;
1132 : }
1133 : #else
1134 : if (y = d1) {
1135 : if (k = lo0bits(&y)) {
1136 : if (k >= 16) {
1137 : x[0] = y | z << 32 - k & 0xffff;
1138 : x[1] = z >> k - 16 & 0xffff;
1139 : x[2] = z >> k;
1140 : i = 2;
1141 : } else {
1142 : x[0] = y & 0xffff;
1143 : x[1] = y >> 16 | z << 16 - k & 0xffff;
1144 : x[2] = z >> k & 0xffff;
1145 : x[3] = z >> k+16;
1146 : i = 3;
1147 : }
1148 : } else {
1149 : x[0] = y & 0xffff;
1150 : x[1] = y >> 16;
1151 : x[2] = z & 0xffff;
1152 : x[3] = z >> 16;
1153 : i = 3;
1154 : }
1155 : } else {
1156 : #ifdef DEBUG
1157 : if (!z)
1158 : Bug("Zero passed to d2b");
1159 : #endif
1160 : k = lo0bits(&z);
1161 : if (k >= 16) {
1162 : x[0] = z;
1163 : i = 0;
1164 : } else {
1165 : x[0] = z & 0xffff;
1166 : x[1] = z >> 16;
1167 : i = 1;
1168 : }
1169 : k += 32;
1170 : }
1171 : while(!x[i])
1172 : --i;
1173 : b->wds = i + 1;
1174 : #endif
1175 : #ifndef Sudden_Underflow
1176 155 : if (de) {
1177 : #endif
1178 : #ifdef IBM
1179 : *e = (de - Bias - (P-1) << 2) + k;
1180 : *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1181 : #else
1182 155 : *e = de - Bias - (P-1) + k;
1183 155 : *bits = P - k;
1184 : #endif
1185 : #ifndef Sudden_Underflow
1186 : } else {
1187 0 : *e = de - Bias - (P-1) + 1 + k;
1188 : #ifdef Pack_32
1189 0 : *bits = 32*i - hi0bits(x[i-1]);
1190 : #else
1191 : *bits = (i+2)*16 - hi0bits(x[i]);
1192 : #endif
1193 : }
1194 : #endif
1195 155 : return b;
1196 : }
1197 : #undef d0
1198 : #undef d1
1199 :
1200 :
1201 : static double ratio (Bigint *a, Bigint *b)
1202 0 : {
1203 : _double da, db;
1204 : int k, ka, kb;
1205 :
1206 0 : value(da) = b2d(a, &ka);
1207 0 : value(db) = b2d(b, &kb);
1208 : #ifdef Pack_32
1209 0 : k = ka - kb + 32*(a->wds - b->wds);
1210 : #else
1211 : k = ka - kb + 16*(a->wds - b->wds);
1212 : #endif
1213 : #ifdef IBM
1214 : if (k > 0) {
1215 : word0(da) += (k >> 2)*Exp_msk1;
1216 : if (k &= 3) {
1217 : da *= 1 << k;
1218 : }
1219 : } else {
1220 : k = -k;
1221 : word0(db) += (k >> 2)*Exp_msk1;
1222 : if (k &= 3)
1223 : db *= 1 << k;
1224 : }
1225 : #else
1226 0 : if (k > 0) {
1227 0 : word0(da) += k*Exp_msk1;
1228 : } else {
1229 0 : k = -k;
1230 0 : word0(db) += k*Exp_msk1;
1231 : }
1232 : #endif
1233 0 : return value(da) / value(db);
1234 : }
1235 :
1236 : static CONST double
1237 : tens[] = {
1238 : 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1239 : 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1240 : 1e20, 1e21, 1e22
1241 : #ifdef VAX
1242 : , 1e23, 1e24
1243 : #endif
1244 : };
1245 :
1246 : #ifdef IEEE_Arith
1247 : static CONST double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1248 : static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1249 : #define n_bigtens 5
1250 : #else
1251 : #ifdef IBM
1252 : static CONST double bigtens[] = { 1e16, 1e32, 1e64 };
1253 : static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1254 : #define n_bigtens 3
1255 : #else
1256 : static CONST double bigtens[] = { 1e16, 1e32 };
1257 : static CONST double tinytens[] = { 1e-16, 1e-32 };
1258 : #define n_bigtens 2
1259 : #endif
1260 : #endif
1261 :
1262 :
1263 : static int quorem(Bigint *b, Bigint *S)
1264 0 : {
1265 : int n;
1266 : Long borrow, y;
1267 : ULong carry, q, ys;
1268 : ULong *bx, *bxe, *sx, *sxe;
1269 : #ifdef Pack_32
1270 : Long z;
1271 : ULong si, zs;
1272 : #endif
1273 :
1274 0 : n = S->wds;
1275 : #ifdef DEBUG
1276 : /*debug*/ if (b->wds > n)
1277 : /*debug*/ Bug("oversize b in quorem");
1278 : #endif
1279 0 : if (b->wds < n)
1280 0 : return 0;
1281 0 : sx = S->x;
1282 0 : sxe = sx + --n;
1283 0 : bx = b->x;
1284 0 : bxe = bx + n;
1285 0 : q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1286 : #ifdef DEBUG
1287 : /*debug*/ if (q > 9)
1288 : /*debug*/ Bug("oversized quotient in quorem");
1289 : #endif
1290 0 : if (q) {
1291 0 : borrow = 0;
1292 0 : carry = 0;
1293 : do {
1294 : #ifdef Pack_32
1295 0 : si = *sx++;
1296 0 : ys = (si & 0xffff) * q + carry;
1297 0 : zs = (si >> 16) * q + (ys >> 16);
1298 0 : carry = zs >> 16;
1299 0 : y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1300 0 : borrow = y >> 16;
1301 : Sign_Extend(borrow, y);
1302 0 : z = (*bx >> 16) - (zs & 0xffff) + borrow;
1303 0 : borrow = z >> 16;
1304 : Sign_Extend(borrow, z);
1305 0 : Storeinc(bx, z, y);
1306 : #else
1307 : ys = *sx++ * q + carry;
1308 : carry = ys >> 16;
1309 : y = *bx - (ys & 0xffff) + borrow;
1310 : borrow = y >> 16;
1311 : Sign_Extend(borrow, y);
1312 : *bx++ = y & 0xffff;
1313 : #endif
1314 : }
1315 0 : while(sx <= sxe);
1316 0 : if (!*bxe) {
1317 0 : bx = b->x;
1318 0 : while(--bxe > bx && !*bxe)
1319 0 : --n;
1320 0 : b->wds = n;
1321 : }
1322 : }
1323 0 : if (cmp(b, S) >= 0) {
1324 0 : q++;
1325 0 : borrow = 0;
1326 0 : carry = 0;
1327 0 : bx = b->x;
1328 0 : sx = S->x;
1329 : do {
1330 : #ifdef Pack_32
1331 0 : si = *sx++;
1332 0 : ys = (si & 0xffff) + carry;
1333 0 : zs = (si >> 16) + (ys >> 16);
1334 0 : carry = zs >> 16;
1335 0 : y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1336 0 : borrow = y >> 16;
1337 : Sign_Extend(borrow, y);
1338 0 : z = (*bx >> 16) - (zs & 0xffff) + borrow;
1339 0 : borrow = z >> 16;
1340 : Sign_Extend(borrow, z);
1341 0 : Storeinc(bx, z, y);
1342 : #else
1343 : ys = *sx++ + carry;
1344 : carry = ys >> 16;
1345 : y = *bx - (ys & 0xffff) + borrow;
1346 : borrow = y >> 16;
1347 : Sign_Extend(borrow, y);
1348 : *bx++ = y & 0xffff;
1349 : #endif
1350 : }
1351 0 : while(sx <= sxe);
1352 0 : bx = b->x;
1353 0 : bxe = bx + n;
1354 0 : if (!*bxe) {
1355 0 : while(--bxe > bx && !*bxe)
1356 0 : --n;
1357 0 : b->wds = n;
1358 : }
1359 : }
1360 0 : return q;
1361 : }
1362 :
1363 : static void destroy_freelist(void)
1364 219 : {
1365 : int i;
1366 : Bigint *tmp;
1367 :
1368 : _THREAD_PRIVATE_MUTEX_LOCK(dtoa_mutex);
1369 3723 : for (i = 0; i <= Kmax; i++) {
1370 3504 : Bigint **listp = &freelist[i];
1371 7166 : while ((tmp = *listp) != NULL) {
1372 158 : *listp = tmp->next;
1373 158 : free(tmp);
1374 : }
1375 3504 : freelist[i] = NULL;
1376 : }
1377 : _THREAD_PRIVATE_MUTEX_UNLOCK(dtoa_mutex);
1378 :
1379 219 : }
1380 :
1381 :
1382 : ZEND_API void zend_freedtoa(char *s)
1383 162 : {
1384 162 : Bigint *b = (Bigint *)((int *)s - 1);
1385 162 : b->maxwds = 1 << (b->k = *(int*)b);
1386 162 : Bfree(b);
1387 162 : }
1388 :
1389 : /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1390 : *
1391 : * Inspired by "How to Print Floating-Point Numbers Accurately" by
1392 : * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1393 : *
1394 : * Modifications:
1395 : * 1. Rather than iterating, we use a simple numeric overestimate
1396 : * to determine k = floor(log10(d)). We scale relevant
1397 : * quantities using O(log2(k)) rather than O(k) multiplications.
1398 : * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1399 : * try to generate digits strictly left to right. Instead, we
1400 : * compute with fewer bits and propagate the carry if necessary
1401 : * when rounding the final digit up. This is often faster.
1402 : * 3. Under the assumption that input will be rounded nearest,
1403 : * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1404 : * That is, we allow equality in stopping tests when the
1405 : * round-nearest rule will give the same floating-point value
1406 : * as would satisfaction of the stopping test with strict
1407 : * inequality.
1408 : * 4. We remove common factors of powers of 2 from relevant
1409 : * quantities.
1410 : * 5. When converting floating-point integers less than 1e16,
1411 : * we use floating-point arithmetic rather than resorting
1412 : * to multiple-precision integers.
1413 : * 6. When asked to produce fewer than 15 digits, we first try
1414 : * to get by with floating-point arithmetic; we resort to
1415 : * multiple-precision integer arithmetic only if we cannot
1416 : * guarantee that the floating-point calculation has given
1417 : * the correctly rounded result. For k requested digits and
1418 : * "uniformly" distributed input, the probability is
1419 : * something like 10^(k-15) that we must resort to the Long
1420 : * calculation.
1421 : */
1422 :
1423 : ZEND_API char * zend_dtoa(double _d, int mode, int ndigits, int *decpt, int *sign, char **rve)
1424 162 : {
1425 : /* Arguments ndigits, decpt, sign are similar to those
1426 : of ecvt and fcvt; trailing zeros are suppressed from
1427 : the returned string. If not null, *rve is set to point
1428 : to the end of the return value. If d is +-Infinity or NaN,
1429 : then *decpt is set to 9999.
1430 :
1431 : mode:
1432 : 0 ==> shortest string that yields d when read in
1433 : and rounded to nearest.
1434 : 1 ==> like 0, but with Steele & White stopping rule;
1435 : e.g. with IEEE P754 arithmetic , mode 0 gives
1436 : 1e23 whereas mode 1 gives 9.999999999999999e22.
1437 : 2 ==> max(1,ndigits) significant digits. This gives a
1438 : return value similar to that of ecvt, except
1439 : that trailing zeros are suppressed.
1440 : 3 ==> through ndigits past the decimal point. This
1441 : gives a return value similar to that from fcvt,
1442 : except that trailing zeros are suppressed, and
1443 : ndigits can be negative.
1444 : 4-9 should give the same return values as 2-3, i.e.,
1445 : 4 <= mode <= 9 ==> same return as mode
1446 : 2 + (mode & 1). These modes are mainly for
1447 : debugging; often they run slower but sometimes
1448 : faster than modes 2-3.
1449 : 4,5,8,9 ==> left-to-right digit generation.
1450 : 6-9 ==> don't try fast floating-point estimate
1451 : (if applicable).
1452 :
1453 : Values of mode other than 0-9 are treated as mode 0.
1454 :
1455 : Sufficient space is allocated to the return value
1456 : to hold the suppressed trailing zeros.
1457 : */
1458 :
1459 162 : int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1,
1460 : j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1461 162 : spec_case = 0, try_quick;
1462 : Long L;
1463 : #ifndef Sudden_Underflow
1464 : int denorm;
1465 : ULong x;
1466 : #endif
1467 : Bigint *b, *b1, *delta, *mlo, *mhi, *S, *tmp;
1468 : double ds;
1469 : char *s, *s0;
1470 : _double d, d2, eps;
1471 :
1472 162 : value(d) = _d;
1473 :
1474 162 : if (word0(d) & Sign_bit) {
1475 : /* set sign for everything, including 0's and NaNs */
1476 0 : *sign = 1;
1477 0 : word0(d) &= ~Sign_bit; /* clear sign bit */
1478 : }
1479 : else
1480 162 : *sign = 0;
1481 :
1482 : #if defined(IEEE_Arith) + defined(VAX)
1483 : #ifdef IEEE_Arith
1484 162 : if ((word0(d) & Exp_mask) == Exp_mask)
1485 : #else
1486 : if (word0(d) == 0x8000)
1487 : #endif
1488 : {
1489 : /* Infinity or NaN */
1490 0 : *decpt = 9999;
1491 : #ifdef IEEE_Arith
1492 0 : if (!word1(d) && !(word0(d) & 0xfffff))
1493 0 : return nrv_alloc("Infinity", rve, 8);
1494 : #endif
1495 0 : return nrv_alloc("NaN", rve, 3);
1496 : }
1497 : #endif
1498 : #ifdef IBM
1499 : value(d) += 0; /* normalize */
1500 : #endif
1501 162 : if (!value(d)) {
1502 7 : *decpt = 1;
1503 7 : return nrv_alloc("0", rve, 1);
1504 : }
1505 :
1506 155 : b = d2b(value(d), &be, &bbits);
1507 : #ifdef Sudden_Underflow
1508 : i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1509 : #else
1510 155 : if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
1511 : #endif
1512 155 : value(d2) = value(d);
1513 155 : word0(d2) &= Frac_mask1;
1514 155 : word0(d2) |= Exp_11;
1515 : #ifdef IBM
1516 : if (j = 11 - hi0bits(word0(d2) & Frac_mask))
1517 : value(d2) /= 1 << j;
1518 : #endif
1519 :
1520 : /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1521 : * log10(x) = log(x) / log(10)
1522 : * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1523 : * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1524 : *
1525 : * This suggests computing an approximation k to log10(d) by
1526 : *
1527 : * k = (i - Bias)*0.301029995663981
1528 : * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1529 : *
1530 : * We want k to be too large rather than too small.
1531 : * The error in the first-order Taylor series approximation
1532 : * is in our favor, so we just round up the constant enough
1533 : * to compensate for any error in the multiplication of
1534 : * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1535 : * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1536 : * adding 1e-13 to the constant term more than suffices.
1537 : * Hence we adjust the constant term to 0.1760912590558.
1538 : * (We could get a more accurate k by invoking log10,
1539 : * but this is probably not worthwhile.)
1540 : */
1541 :
1542 155 : i -= Bias;
1543 : #ifdef IBM
1544 : i <<= 2;
1545 : i += j;
1546 : #endif
1547 : #ifndef Sudden_Underflow
1548 155 : denorm = 0;
1549 : }
1550 : else {
1551 : /* d is denormalized */
1552 :
1553 0 : i = bbits + be + (Bias + (P-1) - 1);
1554 0 : x = i > 32 ? (word0(d) << (64 - i)) | (word1(d) >> (i - 32))
1555 : : (word1(d) << (32 - i));
1556 0 : value(d2) = x;
1557 0 : word0(d2) -= 31*Exp_msk1; /* adjust exponent */
1558 0 : i -= (Bias + (P-1) - 1) + 1;
1559 0 : denorm = 1;
1560 : }
1561 : #endif
1562 155 : ds = (value(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1563 155 : k = (int)ds;
1564 155 : if (ds < 0. && ds != k)
1565 14 : k--; /* want k = floor(ds) */
1566 155 : k_check = 1;
1567 155 : if (k >= 0 && k <= Ten_pmax) {
1568 141 : if (value(d) < tens[k])
1569 1 : k--;
1570 141 : k_check = 0;
1571 : }
1572 155 : j = bbits - i - 1;
1573 155 : if (j >= 0) {
1574 148 : b2 = 0;
1575 148 : s2 = j;
1576 : }
1577 : else {
1578 7 : b2 = -j;
1579 7 : s2 = 0;
1580 : }
1581 155 : if (k >= 0) {
1582 141 : b5 = 0;
1583 141 : s5 = k;
1584 141 : s2 += k;
1585 : }
1586 : else {
1587 14 : b2 -= k;
1588 14 : b5 = -k;
1589 14 : s5 = 0;
1590 : }
1591 155 : if (mode < 0 || mode > 9)
1592 0 : mode = 0;
1593 155 : try_quick = 1;
1594 155 : if (mode > 5) {
1595 0 : mode -= 4;
1596 0 : try_quick = 0;
1597 : }
1598 155 : leftright = 1;
1599 155 : switch(mode) {
1600 : case 0:
1601 : case 1:
1602 0 : ilim = ilim1 = -1;
1603 0 : i = 18;
1604 0 : ndigits = 0;
1605 0 : break;
1606 : case 2:
1607 39 : leftright = 0;
1608 : /* no break */
1609 : case 4:
1610 39 : if (ndigits <= 0)
1611 0 : ndigits = 1;
1612 39 : ilim = ilim1 = i = ndigits;
1613 39 : break;
1614 : case 3:
1615 116 : leftright = 0;
1616 : /* no break */
1617 : case 5:
1618 116 : i = ndigits + k + 1;
1619 116 : ilim = i;
1620 116 : ilim1 = i - 1;
1621 116 : if (i <= 0)
1622 0 : i = 1;
1623 : }
1624 155 : s = s0 = rv_alloc(i);
1625 :
1626 155 : if (ilim >= 0 && ilim <= Quick_max && try_quick) {
1627 :
1628 : /* Try to get by with floating-point arithmetic. */
1629 :
1630 155 : i = 0;
1631 155 : value(d2) = value(d);
1632 155 : k0 = k;
1633 155 : ilim0 = ilim;
1634 155 : ieps = 2; /* conservative */
1635 155 : if (k > 0) {
1636 14 : ds = tens[k&0xf];
1637 14 : j = k >> 4;
1638 14 : if (j & Bletch) {
1639 : /* prevent overflows */
1640 0 : j &= Bletch - 1;
1641 0 : value(d) /= bigtens[n_bigtens-1];
1642 0 : ieps++;
1643 : }
1644 14 : for(; j; j >>= 1, i++)
1645 0 : if (j & 1) {
1646 0 : ieps++;
1647 0 : ds *= bigtens[i];
1648 : }
1649 14 : value(d) /= ds;
1650 : }
1651 141 : else if ((j1 = -k)) {
1652 14 : value(d) *= tens[j1 & 0xf];
1653 14 : for(j = j1 >> 4; j; j >>= 1, i++)
1654 0 : if (j & 1) {
1655 0 : ieps++;
1656 0 : value(d) *= bigtens[i];
1657 : }
1658 : }
1659 155 : if (k_check && value(d) < 1. && ilim > 0) {
1660 0 : if (ilim1 <= 0)
1661 0 : goto fast_failed;
1662 0 : ilim = ilim1;
1663 0 : k--;
1664 0 : value(d) *= 10.;
1665 0 : ieps++;
1666 : }
1667 155 : value(eps) = ieps*value(d) + 7.;
1668 155 : word0(eps) -= (P-1)*Exp_msk1;
1669 155 : if (ilim == 0) {
1670 0 : S = mhi = 0;
1671 0 : value(d) -= 5.;
1672 0 : if (value(d) > value(eps))
1673 0 : goto one_digit;
1674 0 : if (value(d) < -value(eps))
1675 0 : goto no_digits;
1676 0 : goto fast_failed;
1677 : }
1678 : #ifndef No_leftright
1679 155 : if (leftright) {
1680 : /* Use Steele & White method of only
1681 : * generating digits needed.
1682 : */
1683 0 : value(eps) = 0.5/tens[ilim-1] - value(eps);
1684 0 : for(i = 0;;) {
1685 0 : L = value(d);
1686 0 : value(d) -= L;
1687 0 : *s++ = '0' + (int)L;
1688 0 : if (value(d) < value(eps))
1689 0 : goto ret1;
1690 0 : if (1. - value(d) < value(eps))
1691 0 : goto bump_up;
1692 0 : if (++i >= ilim)
1693 0 : break;
1694 0 : value(eps) *= 10.;
1695 0 : value(d) *= 10.;
1696 0 : }
1697 : }
1698 : else {
1699 : #endif
1700 : /* Generate ilim digits, then fix them up. */
1701 155 : value(eps) *= tens[ilim-1];
1702 818 : for(i = 1;; i++, value(d) *= 10.) {
1703 818 : L = value(d);
1704 818 : value(d) -= L;
1705 818 : *s++ = '0' + (int)L;
1706 818 : if (i == ilim) {
1707 155 : if (value(d) > 0.5 + value(eps))
1708 45 : goto bump_up;
1709 110 : else if (value(d) < 0.5 - value(eps)) {
1710 428 : while(*--s == '0');
1711 110 : s++;
1712 110 : goto ret1;
1713 : }
1714 0 : break;
1715 : }
1716 663 : }
1717 : #ifndef No_leftright
1718 : }
1719 : #endif
1720 0 : fast_failed:
1721 0 : s = s0;
1722 0 : value(d) = value(d2);
1723 0 : k = k0;
1724 0 : ilim = ilim0;
1725 : }
1726 :
1727 : /* Do we have a "small" integer? */
1728 :
1729 0 : if (be >= 0 && k <= Int_max) {
1730 : /* Yes. */
1731 0 : ds = tens[k];
1732 0 : if (ndigits < 0 && ilim <= 0) {
1733 0 : S = mhi = 0;
1734 0 : if (ilim < 0 || value(d) <= 5*ds)
1735 : goto no_digits;
1736 0 : goto one_digit;
1737 : }
1738 0 : for(i = 1;; i++) {
1739 0 : L = value(d) / ds;
1740 0 : value(d) -= L*ds;
1741 : #ifdef Check_FLT_ROUNDS
1742 : /* If FLT_ROUNDS == 2, L will usually be high by 1 */
1743 : if (value(d) < 0) {
1744 : L--;
1745 : value(d) += ds;
1746 : }
1747 : #endif
1748 0 : *s++ = '0' + (int)L;
1749 0 : if (i == ilim) {
1750 0 : value(d) += value(d);
1751 0 : if (value(d) > ds || (value(d) == ds && (L & 1))) {
1752 45 : bump_up:
1753 159 : while(*--s == '9')
1754 69 : if (s == s0) {
1755 0 : k++;
1756 0 : *s = '0';
1757 0 : break;
1758 : }
1759 45 : ++*s++;
1760 : }
1761 45 : break;
1762 : }
1763 0 : if (!(value(d) *= 10.))
1764 0 : break;
1765 0 : }
1766 45 : goto ret1;
1767 : }
1768 :
1769 0 : m2 = b2;
1770 0 : m5 = b5;
1771 0 : mhi = mlo = 0;
1772 0 : if (leftright) {
1773 0 : if (mode < 2) {
1774 0 : i =
1775 : #ifndef Sudden_Underflow
1776 : denorm ? be + (Bias + (P-1) - 1 + 1) :
1777 : #endif
1778 : #ifdef IBM
1779 : 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
1780 : #else
1781 : 1 + P - bbits;
1782 : #endif
1783 : }
1784 : else {
1785 0 : j = ilim - 1;
1786 0 : if (m5 >= j)
1787 0 : m5 -= j;
1788 : else {
1789 0 : s5 += j -= m5;
1790 0 : b5 += j;
1791 0 : m5 = 0;
1792 : }
1793 0 : if ((i = ilim) < 0) {
1794 0 : m2 -= i;
1795 0 : i = 0;
1796 : }
1797 : }
1798 0 : b2 += i;
1799 0 : s2 += i;
1800 0 : mhi = i2b(1);
1801 : }
1802 0 : if (m2 > 0 && s2 > 0) {
1803 0 : i = m2 < s2 ? m2 : s2;
1804 0 : b2 -= i;
1805 0 : m2 -= i;
1806 0 : s2 -= i;
1807 : }
1808 0 : if (b5 > 0) {
1809 0 : if (leftright) {
1810 0 : if (m5 > 0) {
1811 0 : mhi = pow5mult(mhi, m5);
1812 0 : b1 = mult(mhi, b);
1813 0 : Bfree(b);
1814 0 : b = b1;
1815 : }
1816 0 : if ((j = b5 - m5)) {
1817 0 : b = pow5mult(b, j);
1818 : }
1819 : } else {
1820 0 : b = pow5mult(b, b5);
1821 : }
1822 : }
1823 0 : S = i2b(1);
1824 0 : if (s5 > 0)
1825 0 : S = pow5mult(S, s5);
1826 : /* Check for special case that d is a normalized power of 2. */
1827 :
1828 0 : if (mode < 2) {
1829 0 : if (!word1(d) && !(word0(d) & Bndry_mask)
1830 : #ifndef Sudden_Underflow
1831 : && word0(d) & Exp_mask
1832 : #endif
1833 : ) {
1834 : /* The special case */
1835 0 : b2 += Log2P;
1836 0 : s2 += Log2P;
1837 0 : spec_case = 1;
1838 : } else {
1839 0 : spec_case = 0;
1840 : }
1841 : }
1842 :
1843 : /* Arrange for convenient computation of quotients:
1844 : * shift left if necessary so divisor has 4 leading 0 bits.
1845 : *
1846 : * Perhaps we should just compute leading 28 bits of S once
1847 : * and for all and pass them and a shift to quorem, so it
1848 : * can do shifts and ors to compute the numerator for q.
1849 : */
1850 : #ifdef Pack_32
1851 0 : if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
1852 0 : i = 32 - i;
1853 : #else
1854 : if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf))
1855 : i = 16 - i;
1856 : #endif
1857 0 : if (i > 4) {
1858 0 : i -= 4;
1859 0 : b2 += i;
1860 0 : m2 += i;
1861 0 : s2 += i;
1862 : }
1863 0 : else if (i < 4) {
1864 0 : i += 28;
1865 0 : b2 += i;
1866 0 : m2 += i;
1867 0 : s2 += i;
1868 : }
1869 0 : if (b2 > 0)
1870 0 : b = lshift(b, b2);
1871 0 : if (s2 > 0)
1872 0 : S = lshift(S, s2);
1873 0 : if (k_check) {
1874 0 : if (cmp(b,S) < 0) {
1875 0 : k--;
1876 0 : b = multadd(b, 10, 0); /* we botched the k estimate */
1877 0 : if (leftright)
1878 0 : mhi = multadd(mhi, 10, 0);
1879 0 : ilim = ilim1;
1880 : }
1881 : }
1882 0 : if (ilim <= 0 && mode > 2) {
1883 0 : if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
1884 : /* no digits, fcvt style */
1885 0 : no_digits:
1886 0 : k = -1 - ndigits;
1887 0 : goto ret;
1888 : }
1889 0 : one_digit:
1890 0 : *s++ = '1';
1891 0 : k++;
1892 0 : goto ret;
1893 : }
1894 0 : if (leftright) {
1895 0 : if (m2 > 0)
1896 0 : mhi = lshift(mhi, m2);
1897 :
1898 : /* Compute mlo -- check for special case
1899 : * that d is a normalized power of 2.
1900 : */
1901 :
1902 0 : mlo = mhi;
1903 0 : if (spec_case) {
1904 0 : mhi = Balloc(mhi->k);
1905 0 : Bcopy(mhi, mlo);
1906 0 : mhi = lshift(mhi, Log2P);
1907 : }
1908 :
1909 0 : for(i = 1;;i++) {
1910 0 : dig = quorem(b,S) + '0';
1911 : /* Do we yet have the shortest decimal string
1912 : * that will round to d?
1913 : */
1914 0 : j = cmp(b, mlo);
1915 0 : delta = diff(S, mhi);
1916 0 : j1 = delta->sign ? 1 : cmp(b, delta);
1917 0 : Bfree(delta);
1918 : #ifndef ROUND_BIASED
1919 0 : if (j1 == 0 && !mode && !(word1(d) & 1)) {
1920 0 : if (dig == '9')
1921 0 : goto round_9_up;
1922 0 : if (j > 0)
1923 0 : dig++;
1924 0 : *s++ = dig;
1925 0 : goto ret;
1926 : }
1927 : #endif
1928 0 : if (j < 0 || (j == 0 && !mode
1929 : #ifndef ROUND_BIASED
1930 : && !(word1(d) & 1)
1931 : #endif
1932 : )) {
1933 0 : if (j1 > 0) {
1934 0 : b = lshift(b, 1);
1935 0 : j1 = cmp(b, S);
1936 0 : if ((j1 > 0 || (j1 == 0 && (dig & 1)))
1937 : && dig++ == '9')
1938 0 : goto round_9_up;
1939 : }
1940 0 : *s++ = dig;
1941 0 : goto ret;
1942 : }
1943 0 : if (j1 > 0) {
1944 0 : if (dig == '9') { /* possible if i == 1 */
1945 0 : round_9_up:
1946 0 : *s++ = '9';
1947 0 : goto roundoff;
1948 : }
1949 0 : *s++ = dig + 1;
1950 0 : goto ret;
1951 : }
1952 0 : *s++ = dig;
1953 0 : if (i == ilim)
1954 0 : break;
1955 0 : b = multadd(b, 10, 0);
1956 0 : if (mlo == mhi)
1957 0 : mlo = mhi = multadd(mhi, 10, 0);
1958 : else {
1959 0 : mlo = multadd(mlo, 10, 0);
1960 0 : mhi = multadd(mhi, 10, 0);
1961 : }
1962 0 : }
1963 : }
1964 : else
1965 0 : for(i = 1;; i++) {
1966 0 : *s++ = dig = quorem(b,S) + '0';
1967 0 : if (i >= ilim)
1968 0 : break;
1969 0 : b = multadd(b, 10, 0);
1970 0 : }
1971 :
1972 : /* Round off last digit */
1973 :
1974 0 : b = lshift(b, 1);
1975 0 : j = cmp(b, S);
1976 0 : if (j > 0 || (j == 0 && (dig & 1))) {
1977 0 : roundoff:
1978 0 : while(*--s == '9')
1979 0 : if (s == s0) {
1980 0 : k++;
1981 0 : *s++ = '1';
1982 0 : goto ret;
1983 : }
1984 0 : ++*s++;
1985 : }
1986 : else {
1987 0 : while(*--s == '0');
1988 0 : s++;
1989 : }
1990 0 : ret:
1991 0 : Bfree(S);
1992 0 : if (mhi) {
1993 0 : if (mlo && mlo != mhi)
1994 0 : Bfree(mlo);
1995 0 : Bfree(mhi);
1996 : }
1997 155 : ret1:
1998 :
1999 : _THREAD_PRIVATE_MUTEX_LOCK(pow5mult_mutex);
2000 310 : while (p5s) {
2001 0 : tmp = p5s;
2002 0 : p5s = p5s->next;
2003 0 : free(tmp);
2004 : }
2005 : _THREAD_PRIVATE_MUTEX_UNLOCK(pow5mult_mutex);
2006 :
2007 155 : Bfree(b);
2008 :
2009 155 : if (s == s0) { /* don't return empty string */
2010 0 : *s++ = '0';
2011 0 : k = 0;
2012 : }
2013 155 : *s = 0;
2014 155 : *decpt = k + 1;
2015 155 : if (rve)
2016 116 : *rve = s;
2017 155 : return s0;
2018 : }
2019 :
2020 : ZEND_API double zend_strtod (CONST char *s00, char **se)
2021 267 : {
2022 : int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
2023 : e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
2024 : CONST char *s, *s0, *s1;
2025 : double aadj, aadj1, adj;
2026 : _double rv, rv0;
2027 : Long L;
2028 : ULong y, z;
2029 : Bigint *bb, *bb1, *bd, *bd0, *bs, *delta, *tmp;
2030 : double result;
2031 :
2032 267 : CONST char decimal_point = '.';
2033 :
2034 267 : sign = nz0 = nz = 0;
2035 267 : value(rv) = 0.;
2036 :
2037 :
2038 267 : for(s = s00; isspace((unsigned char) *s); s++)
2039 : ;
2040 :
2041 267 : if (*s == '-') {
2042 0 : sign = 1;
2043 0 : s++;
2044 267 : } else if (*s == '+') {
2045 0 : s++;
2046 : }
2047 :
2048 267 : if (*s == '\0') {
2049 0 : s = s00;
2050 0 : goto ret;
2051 : }
2052 :
2053 267 : if (*s == '0') {
2054 11 : nz0 = 1;
2055 11 : while(*++s == '0') ;
2056 11 : if (!*s)
2057 0 : goto ret;
2058 : }
2059 267 : s0 = s;
2060 267 : y = z = 0;
2061 471 : for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
2062 204 : if (nd < 9)
2063 204 : y = 10*y + c - '0';
2064 0 : else if (nd < 16)
2065 0 : z = 10*z + c - '0';
2066 267 : nd0 = nd;
2067 267 : if (c == decimal_point) {
2068 264 : c = *++s;
2069 264 : if (!nd) {
2070 77 : for(; c == '0'; c = *++s)
2071 1 : nz++;
2072 76 : if (c > '0' && c <= '9') {
2073 66 : s0 = s;
2074 66 : nf += nz;
2075 66 : nz = 0;
2076 66 : goto have_dig;
2077 : }
2078 10 : goto dig_done;
2079 : }
2080 514 : for(; c >= '0' && c <= '9'; c = *++s) {
2081 326 : have_dig:
2082 326 : nz++;
2083 326 : if (c -= '0') {
2084 310 : nf += nz;
2085 312 : for(i = 1; i < nz; i++)
2086 2 : if (nd++ < 9)
2087 2 : y *= 10;
2088 0 : else if (nd <= DBL_DIG + 1)
2089 0 : z *= 10;
2090 310 : if (nd++ < 9)
2091 310 : y = 10*y + c;
2092 0 : else if (nd <= DBL_DIG + 1)
2093 0 : z = 10*z + c;
2094 310 : nz = 0;
2095 : }
2096 : }
2097 : }
2098 267 : dig_done:
2099 267 : e = 0;
2100 267 : if (c == 'e' || c == 'E') {
2101 8 : if (!nd && !nz && !nz0) {
2102 0 : s = s00;
2103 0 : goto ret;
2104 : }
2105 8 : s00 = s;
2106 8 : esign = 0;
2107 8 : switch(c = *++s) {
2108 : case '-':
2109 0 : esign = 1;
2110 : case '+':
2111 0 : c = *++s;
2112 : }
2113 11 : if (c >= '0' && c <= '9') {
2114 6 : while(c == '0')
2115 0 : c = *++s;
2116 6 : if (c > '0' && c <= '9') {
2117 3 : L = c - '0';
2118 3 : s1 = s;
2119 6 : while((c = *++s) >= '0' && c <= '9')
2120 0 : L = 10*L + c - '0';
2121 3 : if (s - s1 > 8 || L > 19999)
2122 : /* Avoid confusion from exponents
2123 : * so large that e might overflow.
2124 : */
2125 0 : e = 19999; /* safe for 16 bit ints */
2126 : else
2127 3 : e = (int)L;
2128 3 : if (esign)
2129 0 : e = -e;
2130 : }
2131 : else
2132 0 : e = 0;
2133 : }
2134 : else
2135 5 : s = s00;
2136 : }
2137 267 : if (!nd) {
2138 10 : if (!nz && !nz0)
2139 0 : s = s00;
2140 10 : goto ret;
2141 : }
2142 257 : e1 = e -= nf;
2143 :
2144 : /* Now we have nd0 digits, starting at s0, followed by a
2145 : * decimal point, followed by nd-nd0 digits. The number we're
2146 : * after is the integer represented by those digits times
2147 : * 10**e */
2148 :
2149 257 : if (!nd0)
2150 66 : nd0 = nd;
2151 257 : k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
2152 257 : value(rv) = y;
2153 257 : if (k > 9)
2154 0 : value(rv) = tens[k - 9] * value(rv) + z;
2155 257 : bd0 = 0;
2156 257 : if (nd <= DBL_DIG
2157 : #ifndef RND_PRODQUOT
2158 : && FLT_ROUNDS == 1
2159 : #endif
2160 : ) {
2161 257 : if (!e)
2162 40 : goto ret;
2163 217 : if (e > 0) {
2164 3 : if (e <= Ten_pmax) {
2165 : #ifdef VAX
2166 : goto vax_ovfl_check;
2167 : #else
2168 3 : /* value(rv) = */ rounded_product(value(rv),
2169 : tens[e]);
2170 3 : goto ret;
2171 : #endif
2172 : }
2173 0 : i = DBL_DIG - nd;
2174 0 : if (e <= Ten_pmax + i) {
2175 : /* A fancier test would sometimes let us do
2176 : * this for larger i values.
2177 : */
2178 0 : e -= i;
2179 0 : value(rv) *= tens[i];
2180 : #ifdef VAX
2181 : /* VAX exponent range is so narrow we must
2182 : * worry about overflow here...
2183 : */
2184 : vax_ovfl_check:
2185 : word0(rv) -= P*Exp_msk1;
2186 : /* value(rv) = */ rounded_product(value(rv),
2187 : tens[e]);
2188 : if ((word0(rv) & Exp_mask)
2189 : > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
2190 : goto ovfl;
2191 : word0(rv) += P*Exp_msk1;
2192 : #else
2193 0 : /* value(rv) = */ rounded_product(value(rv),
2194 : tens[e]);
2195 : #endif
2196 0 : goto ret;
2197 : }
2198 : }
2199 : #ifndef Inaccurate_Divide
2200 214 : else if (e >= -Ten_pmax) {
2201 214 : /* value(rv) = */ rounded_quotient(value(rv),
2202 : tens[-e]);
2203 214 : goto ret;
2204 : }
2205 : #endif
2206 : }
2207 0 : e1 += nd - k;
2208 :
2209 : /* Get starting approximation = rv * 10**e1 */
2210 :
2211 0 : if (e1 > 0) {
2212 0 : if ((i = e1 & 15))
2213 0 : value(rv) *= tens[i];
2214 0 : if (e1 &= ~15) {
2215 0 : if (e1 > DBL_MAX_10_EXP) {
2216 0 : ovfl:
2217 0 : errno = ERANGE;
2218 : #ifndef Bad_float_h
2219 0 : value(rv) = HUGE_VAL;
2220 : #else
2221 : /* Can't trust HUGE_VAL */
2222 : #ifdef IEEE_Arith
2223 : word0(rv) = Exp_mask;
2224 : word1(rv) = 0;
2225 : #else
2226 : word0(rv) = Big0;
2227 : word1(rv) = Big1;
2228 : #endif
2229 : #endif
2230 0 : if (bd0)
2231 0 : goto retfree;
2232 0 : goto ret;
2233 : }
2234 0 : if (e1 >>= 4) {
2235 0 : for(j = 0; e1 > 1; j++, e1 >>= 1)
2236 0 : if (e1 & 1)
2237 0 : value(rv) *= bigtens[j];
2238 : /* The last multiplication could overflow. */
2239 0 : word0(rv) -= P*Exp_msk1;
2240 0 : value(rv) *= bigtens[j];
2241 0 : if ((z = word0(rv) & Exp_mask)
2242 : > Exp_msk1*(DBL_MAX_EXP+Bias-P))
2243 0 : goto ovfl;
2244 0 : if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
2245 : /* set to largest number */
2246 : /* (Can't trust DBL_MAX) */
2247 0 : word0(rv) = Big0;
2248 0 : word1(rv) = Big1;
2249 : }
2250 : else
2251 0 : word0(rv) += P*Exp_msk1;
2252 : }
2253 :
2254 : }
2255 : }
2256 0 : else if (e1 < 0) {
2257 0 : e1 = -e1;
2258 0 : if ((i = e1 & 15))
2259 0 : value(rv) /= tens[i];
2260 0 : if (e1 &= ~15) {
2261 0 : e1 >>= 4;
2262 0 : if (e1 >= 1 << n_bigtens)
2263 0 : goto undfl;
2264 0 : for(j = 0; e1 > 1; j++, e1 >>= 1)
2265 0 : if (e1 & 1)
2266 0 : value(rv) *= tinytens[j];
2267 : /* The last multiplication could underflow. */
2268 0 : value(rv0) = value(rv);
2269 0 : value(rv) *= tinytens[j];
2270 0 : if (!value(rv)) {
2271 0 : value(rv) = 2.*value(rv0);
2272 0 : value(rv) *= tinytens[j];
2273 0 : if (!value(rv)) {
2274 0 : undfl:
2275 0 : value(rv) = 0.;
2276 0 : errno = ERANGE;
2277 0 : if (bd0)
2278 0 : goto retfree;
2279 0 : goto ret;
2280 : }
2281 0 : word0(rv) = Tiny0;
2282 0 : word1(rv) = Tiny1;
2283 : /* The refinement below will clean
2284 : * this approximation up.
2285 : */
2286 : }
2287 : }
2288 : }
2289 :
2290 : /* Now the hard part -- adjusting rv to the correct value.*/
2291 :
2292 : /* Put digits into bd: true value = bd * 10^e */
2293 :
2294 0 : bd0 = s2b(s0, nd0, nd, y);
2295 :
2296 : for(;;) {
2297 0 : bd = Balloc(bd0->k);
2298 0 : Bcopy(bd, bd0);
2299 0 : bb = d2b(value(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
2300 0 : bs = i2b(1);
2301 :
2302 0 : if (e >= 0) {
2303 0 : bb2 = bb5 = 0;
2304 0 : bd2 = bd5 = e;
2305 : }
2306 : else {
2307 0 : bb2 = bb5 = -e;
2308 0 : bd2 = bd5 = 0;
2309 : }
2310 0 : if (bbe >= 0)
2311 0 : bb2 += bbe;
2312 : else
2313 0 : bd2 -= bbe;
2314 0 : bs2 = bb2;
2315 : #ifdef Sudden_Underflow
2316 : #ifdef IBM
2317 : j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
2318 : #else
2319 : j = P + 1 - bbbits;
2320 : #endif
2321 : #else
2322 0 : i = bbe + bbbits - 1; /* logb(rv) */
2323 0 : if (i < Emin) /* denormal */
2324 0 : j = bbe + (P-Emin);
2325 : else
2326 0 : j = P + 1 - bbbits;
2327 : #endif
2328 0 : bb2 += j;
2329 0 : bd2 += j;
2330 0 : i = bb2 < bd2 ? bb2 : bd2;
2331 0 : if (i > bs2)
2332 0 : i = bs2;
2333 0 : if (i > 0) {
2334 0 : bb2 -= i;
2335 0 : bd2 -= i;
2336 0 : bs2 -= i;
2337 : }
2338 0 : if (bb5 > 0) {
2339 0 : bs = pow5mult(bs, bb5);
2340 0 : bb1 = mult(bs, bb);
2341 0 : Bfree(bb);
2342 0 : bb = bb1;
2343 : }
2344 0 : if (bb2 > 0)
2345 0 : bb = lshift(bb, bb2);
2346 0 : if (bd5 > 0)
2347 0 : bd = pow5mult(bd, bd5);
2348 0 : if (bd2 > 0)
2349 0 : bd = lshift(bd, bd2);
2350 0 : if (bs2 > 0)
2351 0 : bs = lshift(bs, bs2);
2352 0 : delta = diff(bb, bd);
2353 0 : dsign = delta->sign;
2354 0 : delta->sign = 0;
2355 0 : i = cmp(delta, bs);
2356 0 : if (i < 0) {
2357 : /* Error is less than half an ulp -- check for
2358 : * special case of mantissa a power of two.
2359 : */
2360 0 : if (dsign || word1(rv) || word0(rv) & Bndry_mask)
2361 : break;
2362 0 : delta = lshift(delta,Log2P);
2363 0 : if (cmp(delta, bs) > 0)
2364 0 : goto drop_down;
2365 0 : break;
2366 : }
2367 0 : if (i == 0) {
2368 : /* exactly half-way between */
2369 0 : if (dsign) {
2370 0 : if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2371 : && word1(rv) == 0xffffffff) {
2372 : /*boundary case -- increment exponent*/
2373 0 : word0(rv) = (word0(rv) & Exp_mask)
2374 : + Exp_msk1
2375 : #ifdef IBM
2376 : | Exp_msk1 >> 4
2377 : #endif
2378 : ;
2379 0 : word1(rv) = 0;
2380 0 : break;
2381 : }
2382 : }
2383 0 : else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2384 0 : drop_down:
2385 : /* boundary case -- decrement exponent */
2386 : #ifdef Sudden_Underflow
2387 : L = word0(rv) & Exp_mask;
2388 : #ifdef IBM
2389 : if (L < Exp_msk1)
2390 : #else
2391 : if (L <= Exp_msk1)
2392 : #endif
2393 : goto undfl;
2394 : L -= Exp_msk1;
2395 : #else
2396 0 : L = (word0(rv) & Exp_mask) - Exp_msk1;
2397 : #endif
2398 0 : word0(rv) = L | Bndry_mask1;
2399 0 : word1(rv) = 0xffffffff;
2400 : #ifdef IBM
2401 : goto cont;
2402 : #else
2403 0 : break;
2404 : #endif
2405 : }
2406 : #ifndef ROUND_BIASED
2407 0 : if (!(word1(rv) & LSB))
2408 0 : break;
2409 : #endif
2410 0 : if (dsign)
2411 0 : value(rv) += ulp(value(rv));
2412 : #ifndef ROUND_BIASED
2413 : else {
2414 0 : value(rv) -= ulp(value(rv));
2415 : #ifndef Sudden_Underflow
2416 0 : if (!value(rv))
2417 0 : goto undfl;
2418 : #endif
2419 : }
2420 : #endif
2421 0 : break;
2422 : }
2423 0 : if ((aadj = ratio(delta, bs)) <= 2.) {
2424 0 : if (dsign)
2425 0 : aadj = aadj1 = 1.;
2426 0 : else if (word1(rv) || word0(rv) & Bndry_mask) {
2427 : #ifndef Sudden_Underflow
2428 0 : if (word1(rv) == Tiny1 && !word0(rv))
2429 0 : goto undfl;
2430 : #endif
2431 0 : aadj = 1.;
2432 0 : aadj1 = -1.;
2433 : }
2434 : else {
2435 : /* special case -- power of FLT_RADIX to be */
2436 : /* rounded down... */
2437 :
2438 0 : if (aadj < 2./FLT_RADIX)
2439 0 : aadj = 1./FLT_RADIX;
2440 : else
2441 0 : aadj *= 0.5;
2442 0 : aadj1 = -aadj;
2443 : }
2444 : }
2445 : else {
2446 0 : aadj *= 0.5;
2447 0 : aadj1 = dsign ? aadj : -aadj;
2448 : #ifdef Check_FLT_ROUNDS
2449 : switch(FLT_ROUNDS) {
2450 : case 2: /* towards +infinity */
2451 : aadj1 -= 0.5;
2452 : break;
2453 : case 0: /* towards 0 */
2454 : case 3: /* towards -infinity */
2455 : aadj1 += 0.5;
2456 : }
2457 : #else
2458 : if (FLT_ROUNDS == 0)
2459 : aadj1 += 0.5;
2460 : #endif
2461 : }
2462 0 : y = word0(rv) & Exp_mask;
2463 :
2464 : /* Check for overflow */
2465 :
2466 0 : if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2467 0 : value(rv0) = value(rv);
2468 0 : word0(rv) -= P*Exp_msk1;
2469 0 : adj = aadj1 * ulp(value(rv));
2470 0 : value(rv) += adj;
2471 0 : if ((word0(rv) & Exp_mask) >=
2472 : Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2473 0 : if (word0(rv0) == Big0 && word1(rv0) == Big1)
2474 0 : goto ovfl;
2475 0 : word0(rv) = Big0;
2476 0 : word1(rv) = Big1;
2477 0 : goto cont;
2478 : }
2479 : else
2480 0 : word0(rv) += P*Exp_msk1;
2481 : }
2482 : else {
2483 : #ifdef Sudden_Underflow
2484 : if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2485 : value(rv0) = value(rv);
2486 : word0(rv) += P*Exp_msk1;
2487 : adj = aadj1 * ulp(value(rv));
2488 : value(rv) += adj;
2489 : #ifdef IBM
2490 : if ((word0(rv) & Exp_mask) < P*Exp_msk1)
2491 : #else
2492 : if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2493 : #endif
2494 : {
2495 : if (word0(rv0) == Tiny0
2496 : && word1(rv0) == Tiny1)
2497 : goto undfl;
2498 : word0(rv) = Tiny0;
2499 : word1(rv) = Tiny1;
2500 : goto cont;
2501 : }
2502 : else
2503 : word0(rv) -= P*Exp_msk1;
2504 : }
2505 : else {
2506 : adj = aadj1 * ulp(value(rv));
2507 : value(rv) += adj;
2508 : }
2509 : #else
2510 : /* Compute adj so that the IEEE rounding rules will
2511 : * correctly round rv + adj in some half-way cases.
2512 : * If rv * ulp(rv) is denormalized (i.e.,
2513 : * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2514 : * trouble from bits lost to denormalization;
2515 : * example: 1.2e-307 .
2516 : */
2517 0 : if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
2518 0 : aadj1 = (double)(int)(aadj + 0.5);
2519 0 : if (!dsign)
2520 0 : aadj1 = -aadj1;
2521 : }
2522 0 : adj = aadj1 * ulp(value(rv));
2523 0 : value(rv) += adj;
2524 : #endif
2525 : }
2526 0 : z = word0(rv) & Exp_mask;
2527 0 : if (y == z) {
2528 : /* Can we stop now? */
2529 0 : L = aadj;
2530 0 : aadj -= L;
2531 : /* The tolerances below are conservative. */
2532 0 : if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2533 0 : if (aadj < .4999999 || aadj > .5000001)
2534 : break;
2535 : }
2536 0 : else if (aadj < .4999999/FLT_RADIX)
2537 0 : break;
2538 : }
2539 0 : cont:
2540 0 : Bfree(bb);
2541 0 : Bfree(bd);
2542 0 : Bfree(bs);
2543 0 : Bfree(delta);
2544 0 : }
2545 0 : retfree:
2546 0 : Bfree(bb);
2547 0 : Bfree(bd);
2548 0 : Bfree(bs);
2549 0 : Bfree(bd0);
2550 0 : Bfree(delta);
2551 267 : ret:
2552 267 : if (se)
2553 177 : *se = (char *)s;
2554 267 : result = sign ? -value(rv) : value(rv);
2555 :
2556 : _THREAD_PRIVATE_MUTEX_LOCK(pow5mult_mutex);
2557 534 : while (p5s) {
2558 0 : tmp = p5s;
2559 0 : p5s = p5s->next;
2560 0 : free(tmp);
2561 : }
2562 : _THREAD_PRIVATE_MUTEX_UNLOCK(pow5mult_mutex);
2563 :
2564 267 : return result;
2565 : }
2566 :
2567 : ZEND_API double zend_hex_strtod(const char *str, char **endptr)
2568 0 : {
2569 0 : const char *s = str;
2570 : char c;
2571 0 : int any = 0;
2572 0 : double value = 0;
2573 :
2574 0 : if (*s == '0' && (s[1] == 'x' || s[1] == 'X')) {
2575 0 : s += 2;
2576 : }
2577 :
2578 0 : while ((c = *s++)) {
2579 0 : if (c >= '0' && c <= '9') {
2580 0 : c -= '0';
2581 0 : } else if (c >= 'A' && c <= 'F') {
2582 0 : c -= 'A' - 10;
2583 0 : } else if (c >= 'a' && c <= 'f') {
2584 0 : c -= 'a' - 10;
2585 : } else {
2586 : break;
2587 : }
2588 :
2589 0 : any = 1;
2590 0 : value = value * 16 + c;
2591 : }
2592 :
2593 0 : if (endptr != NULL) {
2594 0 : *endptr = (char *)(any ? s - 1 : str);
2595 : }
2596 :
2597 0 : return value;
2598 : }
2599 :
2600 : /*
2601 : * Local variables:
2602 : * tab-width: 4
2603 : * c-basic-offset: 4
2604 : * End:
2605 : * vim600: sw=4 ts=4 fdm=marker
2606 : * vim<600: sw=4 ts=4
2607 : */
|
