LTP GCOV extension - code coverage report
Current view: directory - Zend - zend_strtod.c
Test: PHP Code Coverage
Date: 2007-04-10 Instrumented lines: 1041
Code covered: 27.0 % Executed lines: 281
Legend: not executed executed

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

Generated by: LTP GCOV extension version 1.5