root/Zend/zend_strtod.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
  1. Bug
  2. zend_startup_strtod
  3. zend_shutdown_strtod
  4. Balloc
  5. Bfree
  6. rv_alloc
  7. nrv_alloc
  8. multadd
  9. hi0bits
  10. lo0bits
  11. i2b
  12. mult
  13. s2b
  14. pow5mult
  15. lshift
  16. cmp
  17. diff
  18. ulp
  19. d2b
  20. ratio
  21. quorem
  22. destroy_freelist
  23. zend_freedtoa
  24. zend_dtoa
  25. zend_strtod
  26. zend_hex_strtod
  27. zend_oct_strtod
  28. zend_bin_strtod

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

/* [<][>][^][v][top][bottom][index][help] */