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