File:  [ELWIX - Embedded LightWeight unIX -] / embedaddon / php / Zend / zend_strtod.c
Revision 1.1.1.3 (vendor branch): download - view: text, annotated - select for diffs - revision graph
Sun Jun 15 20:04:04 2014 UTC (10 years, 4 months ago) by misho
Branches: php, MAIN
CVS tags: v5_4_29, HEAD
php 5.4.29

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

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>