From 72e7f011b29998d8a3e15eb5b381ef962af5fe5b Mon Sep 17 00:00:00 2001 From: Karen Arutyunov Date: Fri, 5 Apr 2019 10:30:58 +0300 Subject: Upgrade to 8.0.15 --- mysql/strings/dtoa.c | 2784 -------------------------------------------------- 1 file changed, 2784 deletions(-) delete mode 100644 mysql/strings/dtoa.c (limited to 'mysql/strings/dtoa.c') diff --git a/mysql/strings/dtoa.c b/mysql/strings/dtoa.c deleted file mode 100644 index e05f86a..0000000 --- a/mysql/strings/dtoa.c +++ /dev/null @@ -1,2784 +0,0 @@ -/* Copyright (c) 2007, 2013, Oracle and/or its affiliates. All rights reserved. - - This library is free software; you can redistribute it and/or - modify it under the terms of the GNU Library General Public - License as published by the Free Software Foundation; version 2 - of the License. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program; if not, write to the Free Software - Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA */ - -/**************************************************************** - - This file incorporates work covered by the following copyright and - permission notice: - - The author of this software is David M. Gay. - - Copyright (c) 1991, 2000, 2001 by Lucent Technologies. - - Permission to use, copy, modify, and distribute this software for any - purpose without fee is hereby granted, provided that this entire notice - is included in all copies of any software which is or includes a copy - or modification of this software and in all copies of the supporting - documentation for such software. - - THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED - WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY - REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY - OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. - - ***************************************************************/ - -#include -#include /* for memcpy and NOT_FIXED_DEC */ - -#ifndef EOVERFLOW -#define EOVERFLOW 84 -#endif - -/** - Appears to suffice to not call malloc() in most cases. - @todo - see if it is possible to get rid of malloc(). - this constant is sufficient to avoid malloc() on all inputs I have tried. -*/ -#define DTOA_BUFF_SIZE (460 * sizeof(void *)) - -/* Magic value returned by dtoa() to indicate overflow */ -#define DTOA_OVERFLOW 9999 - -static double my_strtod_int(const char *, char **, int *, char *, size_t); -static char *dtoa(double, int, int, int *, int *, char **, char *, size_t); -static void dtoa_free(char *, char *, size_t); - -/** - @brief - Converts a given floating point number to a zero-terminated string - representation using the 'f' format. - - @details - This function is a wrapper around dtoa() to do the same as - sprintf(to, "%-.*f", precision, x), though the conversion is usually more - precise. The only difference is in handling [-,+]infinity and nan values, - in which case we print '0\0' to the output string and indicate an overflow. - - @param x the input floating point number. - @param precision the number of digits after the decimal point. - All properties of sprintf() apply: - - if the number of significant digits after the decimal - point is less than precision, the resulting string is - right-padded with zeros - - if the precision is 0, no decimal point appears - - if a decimal point appears, at least one digit appears - before it - @param to pointer to the output buffer. The longest string which - my_fcvt() can return is FLOATING_POINT_BUFFER bytes - (including the terminating '\0'). - @param error if not NULL, points to a location where the status of - conversion is stored upon return. - FALSE successful conversion - TRUE the input number is [-,+]infinity or nan. - The output string in this case is always '0'. - @return number of written characters (excluding terminating '\0') -*/ - -size_t my_fcvt(double x, int precision, char *to, my_bool *error) -{ - int decpt, sign, len, i; - char *res, *src, *end, *dst= to; - char buf[DTOA_BUFF_SIZE]; - DBUG_ASSERT(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL); - - res= dtoa(x, 5, precision, &decpt, &sign, &end, buf, sizeof(buf)); - - if (decpt == DTOA_OVERFLOW) - { - dtoa_free(res, buf, sizeof(buf)); - *to++= '0'; - *to= '\0'; - if (error != NULL) - *error= TRUE; - return 1; - } - - src= res; - len= (int)(end - src); - - if (sign) - *dst++= '-'; - - if (decpt <= 0) - { - *dst++= '0'; - *dst++= '.'; - for (i= decpt; i < 0; i++) - *dst++= '0'; - } - - for (i= 1; i <= len; i++) - { - *dst++= *src++; - if (i == decpt && i < len) - *dst++= '.'; - } - while (i++ <= decpt) - *dst++= '0'; - - if (precision > 0) - { - if (len <= decpt) - *dst++= '.'; - - for (i= precision - MY_MAX(0, (len - decpt)); i > 0; i--) - *dst++= '0'; - } - - *dst= '\0'; - if (error != NULL) - *error= FALSE; - - dtoa_free(res, buf, sizeof(buf)); - - return dst - to; -} - -/** - @brief - Converts a given floating point number to a zero-terminated string - representation with a given field width using the 'e' format - (aka scientific notation) or the 'f' one. - - @details - The format is chosen automatically to provide the most number of significant - digits (and thus, precision) with a given field width. In many cases, the - result is similar to that of sprintf(to, "%g", x) with a few notable - differences: - - the conversion is usually more precise than C library functions. - - there is no 'precision' argument. instead, we specify the number of - characters available for conversion (i.e. a field width). - - the result never exceeds the specified field width. If the field is too - short to contain even a rounded decimal representation, my_gcvt() - indicates overflow and truncates the output string to the specified width. - - float-type arguments are handled differently than double ones. For a - float input number (i.e. when the 'type' argument is MY_GCVT_ARG_FLOAT) - we deliberately limit the precision of conversion by FLT_DIG digits to - avoid garbage past the significant digits. - - unlike sprintf(), in cases where the 'e' format is preferred, we don't - zero-pad the exponent to save space for significant digits. The '+' sign - for a positive exponent does not appear for the same reason. - - @param x the input floating point number. - @param type is either MY_GCVT_ARG_FLOAT or MY_GCVT_ARG_DOUBLE. - Specifies the type of the input number (see notes above). - @param width field width in characters. The minimal field width to - hold any number representation (albeit rounded) is 7 - characters ("-Ne-NNN"). - @param to pointer to the output buffer. The result is always - zero-terminated, and the longest returned string is thus - 'width + 1' bytes. - @param error if not NULL, points to a location where the status of - conversion is stored upon return. - FALSE successful conversion - TRUE the input number is [-,+]infinity or nan. - The output string in this case is always '0'. - @return number of written characters (excluding terminating '\0') - - @todo - Check if it is possible and makes sense to do our own rounding on top of - dtoa() instead of calling dtoa() twice in (rare) cases when the resulting - string representation does not fit in the specified field width and we want - to re-round the input number with fewer significant digits. Examples: - - my_gcvt(-9e-3, ..., 4, ...); - my_gcvt(-9e-3, ..., 2, ...); - my_gcvt(1.87e-3, ..., 4, ...); - my_gcvt(55, ..., 1, ...); - - We do our best to minimize such cases by: - - - passing to dtoa() the field width as the number of significant digits - - - removing the sign of the number early (and decreasing the width before - passing it to dtoa()) - - - choosing the proper format to preserve the most number of significant - digits. -*/ - -size_t my_gcvt(double x, my_gcvt_arg_type type, int width, char *to, - my_bool *error) -{ - int decpt, sign, len, exp_len; - char *res, *src, *end, *dst= to, *dend= dst + width; - char buf[DTOA_BUFF_SIZE]; - my_bool have_space, force_e_format; - DBUG_ASSERT(width > 0 && to != NULL); - - /* We want to remove '-' from equations early */ - if (x < 0.) - width--; - - res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? width : MY_MIN(width, FLT_DIG), - &decpt, &sign, &end, buf, sizeof(buf)); - if (decpt == DTOA_OVERFLOW) - { - dtoa_free(res, buf, sizeof(buf)); - *to++= '0'; - *to= '\0'; - if (error != NULL) - *error= TRUE; - return 1; - } - - if (error != NULL) - *error= FALSE; - - src= res; - len= (int)(end - res); - - /* - Number of digits in the exponent from the 'e' conversion. - The sign of the exponent is taken into account separetely, we don't need - to count it here. - */ - exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9); - - /* - Do we have enough space for all digits in the 'f' format? - Let 'len' be the number of significant digits returned by dtoa, - and F be the length of the resulting decimal representation. - Consider the following cases: - 1. decpt <= 0, i.e. we have "0.NNN" => F = len - decpt + 2 - 2. 0 < decpt < len, i.e. we have "NNN.NNN" => F = len + 1 - 3. len <= decpt, i.e. we have "NNN00" => F = decpt - */ - have_space= (decpt <= 0 ? len - decpt + 2 : - decpt > 0 && decpt < len ? len + 1 : - decpt) <= width; - /* - The following is true when no significant digits can be placed with the - specified field width using the 'f' format, and the 'e' format - will not be truncated. - */ - force_e_format= (decpt <= 0 && width <= 2 - decpt && width >= 3 + exp_len); - /* - Assume that we don't have enough space to place all significant digits in - the 'f' format. We have to choose between the 'e' format and the 'f' one - to keep as many significant digits as possible. - Let E and F be the lengths of decimal representaion in the 'e' and 'f' - formats, respectively. We want to use the 'f' format if, and only if F <= E. - Consider the following cases: - 1. decpt <= 0. - F = len - decpt + 2 (see above) - E = len + (len > 1) + 1 + 1 (decpt <= -99) + (decpt <= -9) + 1 - ("N.NNe-MMM") - (F <= E) <=> (len == 1 && decpt >= -1) || (len > 1 && decpt >= -2) - We also need to ensure that if the 'f' format is chosen, - the field width allows us to place at least one significant digit - (i.e. width > 2 - decpt). If not, we prefer the 'e' format. - 2. 0 < decpt < len - F = len + 1 (see above) - E = len + 1 + 1 + ... ("N.NNeMMM") - F is always less than E. - 3. len <= decpt <= width - In this case we have enough space to represent the number in the 'f' - format, so we prefer it with some exceptions. - 4. width < decpt - The number cannot be represented in the 'f' format at all, always use - the 'e' 'one. - */ - if ((have_space || - /* - Not enough space, let's see if the 'f' format provides the most number - of significant digits. - */ - ((decpt <= width && (decpt >= -1 || (decpt == -2 && - (len > 1 || !force_e_format)))) && - !force_e_format)) && - - /* - Use the 'e' format in some cases even if we have enough space for the - 'f' one. See comment for MAX_DECPT_FOR_F_FORMAT. - */ - (!have_space || (decpt >= -MAX_DECPT_FOR_F_FORMAT + 1 && - (decpt <= MAX_DECPT_FOR_F_FORMAT || len > decpt)))) - { - /* 'f' format */ - int i; - - width-= (decpt < len) + (decpt <= 0 ? 1 - decpt : 0); - - /* Do we have to truncate any digits? */ - if (width < len) - { - if (width < decpt) - { - if (error != NULL) - *error= TRUE; - width= decpt; - } - - /* - We want to truncate (len - width) least significant digits after the - decimal point. For this we are calling dtoa with mode=5, passing the - number of significant digits = (len-decpt) - (len-width) = width-decpt - */ - dtoa_free(res, buf, sizeof(buf)); - res= dtoa(x, 5, width - decpt, &decpt, &sign, &end, buf, sizeof(buf)); - src= res; - len= (int)(end - res); - } - - if (len == 0) - { - /* Underflow. Just print '0' and exit */ - *dst++= '0'; - goto end; - } - - /* - At this point we are sure we have enough space to put all digits - returned by dtoa - */ - if (sign && dst < dend) - *dst++= '-'; - if (decpt <= 0) - { - if (dst < dend) - *dst++= '0'; - if (len > 0 && dst < dend) - *dst++= '.'; - for (; decpt < 0 && dst < dend; decpt++) - *dst++= '0'; - } - - for (i= 1; i <= len && dst < dend; i++) - { - *dst++= *src++; - if (i == decpt && i < len && dst < dend) - *dst++= '.'; - } - while (i++ <= decpt && dst < dend) - *dst++= '0'; - } - else - { - /* 'e' format */ - int decpt_sign= 0; - - if (--decpt < 0) - { - decpt= -decpt; - width--; - decpt_sign= 1; - } - width-= 1 + exp_len; /* eNNN */ - - if (len > 1) - width--; - - if (width <= 0) - { - /* Overflow */ - if (error != NULL) - *error= TRUE; - width= 0; - } - - /* Do we have to truncate any digits? */ - if (width < len) - { - /* Yes, re-convert with a smaller width */ - dtoa_free(res, buf, sizeof(buf)); - res= dtoa(x, 4, width, &decpt, &sign, &end, buf, sizeof(buf)); - src= res; - len= (int)(end - res); - if (--decpt < 0) - decpt= -decpt; - } - /* - At this point we are sure we have enough space to put all digits - returned by dtoa - */ - if (sign && dst < dend) - *dst++= '-'; - if (dst < dend) - *dst++= *src++; - if (len > 1 && dst < dend) - { - *dst++= '.'; - while (src < end && dst < dend) - *dst++= *src++; - } - if (dst < dend) - *dst++= 'e'; - if (decpt_sign && dst < dend) - *dst++= '-'; - - if (decpt >= 100 && dst < dend) - { - *dst++= decpt / 100 + '0'; - decpt%= 100; - if (dst < dend) - *dst++= decpt / 10 + '0'; - } - else if (decpt >= 10 && dst < dend) - *dst++= decpt / 10 + '0'; - if (dst < dend) - *dst++= decpt % 10 + '0'; - - } - -end: - dtoa_free(res, buf, sizeof(buf)); - *dst= '\0'; - - return dst - to; -} - -/** - @brief - Converts string to double (string does not have to be zero-terminated) - - @details - This is a wrapper around dtoa's version of strtod(). - - @param str input string - @param end address of a pointer to the first character after the input - string. Upon return the pointer is set to point to the first - rejected character. - @param error Upon return is set to EOVERFLOW in case of underflow or - overflow. - - @return The resulting double value. In case of underflow, 0.0 is - returned. In case overflow, signed DBL_MAX is returned. -*/ - -double my_strtod(const char *str, char **end, int *error) -{ - char buf[DTOA_BUFF_SIZE]; - double res; - DBUG_ASSERT(end != NULL && ((str != NULL && *end != NULL) || - (str == NULL && *end == NULL)) && - error != NULL); - - res= my_strtod_int(str, end, error, buf, sizeof(buf)); - return (*error == 0) ? res : (res < 0 ? -DBL_MAX : DBL_MAX); -} - - -double my_atof(const char *nptr) -{ - int error; - const char *end= nptr+65535; /* Should be enough */ - return (my_strtod(nptr, (char**) &end, &error)); -} - - -/**************************************************************** - * - * The author of this software is David M. Gay. - * - * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. - * - * Permission to use, copy, modify, and distribute this software for any - * purpose without fee is hereby granted, provided that this entire notice - * is included in all copies of any software which is or includes a copy - * or modification of this software and in all copies of the supporting - * documentation for such software. - * - * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED - * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY - * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY - * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. - * - ***************************************************************/ -/* Please send bug reports to David M. Gay (dmg at acm dot org, - * with " at " changed at "@" and " dot " changed to "."). */ - -/* - Original copy of the software is located at http://www.netlib.org/fp/dtoa.c - It was adjusted to serve MySQL server needs: - * strtod() was modified to not expect a zero-terminated string. - It now honors 'se' (end of string) argument as the input parameter, - not just as the output one. - * in dtoa(), in case of overflow/underflow/NaN result string now contains "0"; - decpt is set to DTOA_OVERFLOW to indicate overflow. - * support for VAX, IBM mainframe and 16-bit hardware removed - * we always assume that 64-bit integer type is available - * support for Kernigan-Ritchie style headers (pre-ANSI compilers) - removed - * all gcc warnings ironed out - * we always assume multithreaded environment, so we had to change - memory allocation procedures to use stack in most cases; - malloc is used as the last resort. - * pow5mult rewritten to use pre-calculated pow5 list instead of - the one generated on the fly. -*/ - - -/* - On a machine with IEEE extended-precision registers, it is - necessary to specify double-precision (53-bit) rounding precision - before invoking strtod or dtoa. If the machine uses (the equivalent - of) Intel 80x87 arithmetic, the call - _control87(PC_53, MCW_PC); - does this with many compilers. Whether this or another call is - appropriate depends on the compiler; for this to work, it may be - necessary to #include "float.h" or another system-dependent header - file. -*/ - -/* - #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 - and dtoa should round accordingly. - #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3 - and Honor_FLT_ROUNDS is not #defined. - - TODO: check if we can get rid of the above two -*/ - -typedef int32 Long; -typedef uint32 ULong; -typedef int64 LLong; -typedef uint64 ULLong; - -typedef union { double d; ULong L[2]; } U; - -#if defined(WORDS_BIGENDIAN) || (defined(__FLOAT_WORD_ORDER) && \ - (__FLOAT_WORD_ORDER == __BIG_ENDIAN)) -#define word0(x) (x)->L[0] -#define word1(x) (x)->L[1] -#else -#define word0(x) (x)->L[1] -#define word1(x) (x)->L[0] -#endif - -#define dval(x) (x)->d - -/* #define P DBL_MANT_DIG */ -/* Ten_pmax= floor(P*log(2)/log(5)) */ -/* Bletch= (highest power of 2 < DBL_MAX_10_EXP) / 16 */ -/* Quick_max= floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ -/* Int_max= floor(P*log(FLT_RADIX)/log(10) - 1) */ - -#define Exp_shift 20 -#define Exp_shift1 20 -#define Exp_msk1 0x100000 -#define Exp_mask 0x7ff00000 -#define P 53 -#define Bias 1023 -#define Emin (-1022) -#define Exp_1 0x3ff00000 -#define Exp_11 0x3ff00000 -#define Ebits 11 -#define Frac_mask 0xfffff -#define Frac_mask1 0xfffff -#define Ten_pmax 22 -#define Bletch 0x10 -#define Bndry_mask 0xfffff -#define Bndry_mask1 0xfffff -#define LSB 1 -#define Sign_bit 0x80000000 -#define Log2P 1 -#define Tiny1 1 -#define Quick_max 14 -#define Int_max 14 - -#ifndef Flt_Rounds -#ifdef FLT_ROUNDS -#define Flt_Rounds FLT_ROUNDS -#else -#define Flt_Rounds 1 -#endif -#endif /*Flt_Rounds*/ - -#ifdef Honor_FLT_ROUNDS -#define Rounding rounding -#undef Check_FLT_ROUNDS -#define Check_FLT_ROUNDS -#else -#define Rounding Flt_Rounds -#endif - -#define rounded_product(a,b) a*= b -#define rounded_quotient(a,b) a/= b - -#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) -#define Big1 0xffffffff -#define FFFFFFFF 0xffffffffUL - -/* This is tested to be enough for dtoa */ - -#define Kmax 15 - -#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ - 2*sizeof(int) + y->wds*sizeof(ULong)) - -/* Arbitrary-length integer */ - -typedef struct Bigint -{ - union { - ULong *x; /* points right after this Bigint object */ - struct Bigint *next; /* to maintain free lists */ - } p; - int k; /* 2^k = maxwds */ - int maxwds; /* maximum length in 32-bit words */ - int sign; /* not zero if number is negative */ - int wds; /* current length in 32-bit words */ -} Bigint; - - -/* A simple stack-memory based allocator for Bigints */ - -typedef struct Stack_alloc -{ - char *begin; - char *free; - char *end; - /* - Having list of free blocks lets us reduce maximum required amount - of memory from ~4000 bytes to < 1680 (tested on x86). - */ - Bigint *freelist[Kmax+1]; -} Stack_alloc; - - -/* - Try to allocate object on stack, and resort to malloc if all - stack memory is used. Ensure allocated objects to be aligned by the pointer - size in order to not break the alignment rules when storing a pointer to a - Bigint. -*/ - -static Bigint *Balloc(int k, Stack_alloc *alloc) -{ - Bigint *rv; - DBUG_ASSERT(k <= Kmax); - if (k <= Kmax && alloc->freelist[k]) - { - rv= alloc->freelist[k]; - alloc->freelist[k]= rv->p.next; - } - else - { - int x, len; - - x= 1 << k; - len= MY_ALIGN(sizeof(Bigint) + x * sizeof(ULong), SIZEOF_CHARP); - - if (alloc->free + len <= alloc->end) - { - rv= (Bigint*) alloc->free; - alloc->free+= len; - } - else - rv= (Bigint*) malloc(len); - - rv->k= k; - rv->maxwds= x; - } - rv->sign= rv->wds= 0; - rv->p.x= (ULong*) (rv + 1); - return rv; -} - - -/* - If object was allocated on stack, try putting it to the free - list. Otherwise call free(). -*/ - -static void Bfree(Bigint *v, Stack_alloc *alloc) -{ - char *gptr= (char*) v; /* generic pointer */ - if (gptr < alloc->begin || gptr >= alloc->end) - free(gptr); - else if (v->k <= Kmax) - { - /* - Maintain free lists only for stack objects: this way we don't - have to bother with freeing lists in the end of dtoa; - heap should not be used normally anyway. - */ - v->p.next= alloc->freelist[v->k]; - alloc->freelist[v->k]= v; - } -} - - -/* - This is to place return value of dtoa in: tries to use stack - as well, but passes by free lists management and just aligns len by - the pointer size in order to not break the alignment rules when storing a - pointer to a Bigint. -*/ - -static char *dtoa_alloc(int i, Stack_alloc *alloc) -{ - char *rv; - int aligned_size= MY_ALIGN(i, SIZEOF_CHARP); - if (alloc->free + aligned_size <= alloc->end) - { - rv= alloc->free; - alloc->free+= aligned_size; - } - else - rv= malloc(i); - return rv; -} - - -/* - dtoa_free() must be used to free values s returned by dtoa() - This is the counterpart of dtoa_alloc() -*/ - -static void dtoa_free(char *gptr, char *buf, size_t buf_size) -{ - if (gptr < buf || gptr >= buf + buf_size) - free(gptr); -} - - -/* Bigint arithmetic functions */ - -/* Multiply by m and add a */ - -static Bigint *multadd(Bigint *b, int m, int a, Stack_alloc *alloc) -{ - int i, wds; - ULong *x; - ULLong carry, y; - Bigint *b1; - - wds= b->wds; - x= b->p.x; - i= 0; - carry= a; - do - { - y= *x * (ULLong)m + carry; - carry= y >> 32; - *x++= (ULong)(y & FFFFFFFF); - } - while (++i < wds); - if (carry) - { - if (wds >= b->maxwds) - { - b1= Balloc(b->k+1, alloc); - Bcopy(b1, b); - Bfree(b, alloc); - b= b1; - } - b->p.x[wds++]= (ULong) carry; - b->wds= wds; - } - return b; -} - -/** - Converts a string to Bigint. - - Now we have nd0 digits, starting at s, followed by a - decimal point, followed by nd-nd0 digits. - Unless nd0 == nd, in which case we have a number of the form: - ".xxxxxx" or "xxxxxx." - - @param s Input string, already partially parsed by my_strtod_int(). - @param nd0 Number of digits before decimal point. - @param nd Total number of digits. - @param y9 Pre-computed value of the first nine digits. - @param alloc Stack allocator for Bigints. - */ -static Bigint *s2b(const char *s, int nd0, int nd, ULong y9, Stack_alloc *alloc) -{ - Bigint *b; - int i, k; - Long x, y; - - x= (nd + 8) / 9; - for (k= 0, y= 1; x > y; y <<= 1, k++) ; - b= Balloc(k, alloc); - b->p.x[0]= y9; - b->wds= 1; - - i= 9; - if (9 < nd0) - { - s+= 9; - do - b= multadd(b, 10, *s++ - '0', alloc); - while (++i < nd0); - s++; /* skip '.' */ - } - else - s+= 10; - /* now do the fractional part */ - for(; i < nd; i++) - b= multadd(b, 10, *s++ - '0', alloc); - return b; -} - - -static int hi0bits(ULong x) -{ - int k= 0; - - if (!(x & 0xffff0000)) - { - k= 16; - x<<= 16; - } - if (!(x & 0xff000000)) - { - k+= 8; - x<<= 8; - } - if (!(x & 0xf0000000)) - { - k+= 4; - x<<= 4; - } - if (!(x & 0xc0000000)) - { - k+= 2; - x<<= 2; - } - if (!(x & 0x80000000)) - { - k++; - if (!(x & 0x40000000)) - return 32; - } - return k; -} - - -static int lo0bits(ULong *y) -{ - int k; - ULong x= *y; - - if (x & 7) - { - if (x & 1) - return 0; - if (x & 2) - { - *y= x >> 1; - return 1; - } - *y= x >> 2; - return 2; - } - k= 0; - if (!(x & 0xffff)) - { - k= 16; - x>>= 16; - } - if (!(x & 0xff)) - { - k+= 8; - x>>= 8; - } - if (!(x & 0xf)) - { - k+= 4; - x>>= 4; - } - if (!(x & 0x3)) - { - k+= 2; - x>>= 2; - } - if (!(x & 1)) - { - k++; - x>>= 1; - if (!x) - return 32; - } - *y= x; - return k; -} - - -/* Convert integer to Bigint number */ - -static Bigint *i2b(int i, Stack_alloc *alloc) -{ - Bigint *b; - - b= Balloc(1, alloc); - b->p.x[0]= i; - b->wds= 1; - return b; -} - - -/* Multiply two Bigint numbers */ - -static Bigint *mult(Bigint *a, Bigint *b, Stack_alloc *alloc) -{ - Bigint *c; - int k, wa, wb, wc; - ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; - ULong y; - ULLong carry, z; - - if (a->wds < b->wds) - { - c= a; - a= b; - b= c; - } - k= a->k; - wa= a->wds; - wb= b->wds; - wc= wa + wb; - if (wc > a->maxwds) - k++; - c= Balloc(k, alloc); - for (x= c->p.x, xa= x + wc; x < xa; x++) - *x= 0; - xa= a->p.x; - xae= xa + wa; - xb= b->p.x; - xbe= xb + wb; - xc0= c->p.x; - for (; xb < xbe; xc0++) - { - if ((y= *xb++)) - { - x= xa; - xc= xc0; - carry= 0; - do - { - z= *x++ * (ULLong)y + *xc + carry; - carry= z >> 32; - *xc++= (ULong) (z & FFFFFFFF); - } - while (x < xae); - *xc= (ULong) carry; - } - } - for (xc0= c->p.x, xc= xc0 + wc; wc > 0 && !*--xc; --wc) ; - c->wds= wc; - return c; -} - - -/* - Precalculated array of powers of 5: tested to be enough for - vasting majority of dtoa_r cases. -*/ - -static ULong powers5[]= -{ - 625UL, - - 390625UL, - - 2264035265UL, 35UL, - - 2242703233UL, 762134875UL, 1262UL, - - 3211403009UL, 1849224548UL, 3668416493UL, 3913284084UL, 1593091UL, - - 781532673UL, 64985353UL, 253049085UL, 594863151UL, 3553621484UL, - 3288652808UL, 3167596762UL, 2788392729UL, 3911132675UL, 590UL, - - 2553183233UL, 3201533787UL, 3638140786UL, 303378311UL, 1809731782UL, - 3477761648UL, 3583367183UL, 649228654UL, 2915460784UL, 487929380UL, - 1011012442UL, 1677677582UL, 3428152256UL, 1710878487UL, 1438394610UL, - 2161952759UL, 4100910556UL, 1608314830UL, 349175UL -}; - - -static Bigint p5_a[]= -{ - /* { x } - k - maxwds - sign - wds */ - { { powers5 }, 1, 1, 0, 1 }, - { { powers5 + 1 }, 1, 1, 0, 1 }, - { { powers5 + 2 }, 1, 2, 0, 2 }, - { { powers5 + 4 }, 2, 3, 0, 3 }, - { { powers5 + 7 }, 3, 5, 0, 5 }, - { { powers5 + 12 }, 4, 10, 0, 10 }, - { { powers5 + 22 }, 5, 19, 0, 19 } -}; - -#define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1) - -static Bigint *pow5mult(Bigint *b, int k, Stack_alloc *alloc) -{ - Bigint *b1, *p5, *p51=NULL; - int i; - static int p05[3]= { 5, 25, 125 }; - my_bool overflow= FALSE; - - if ((i= k & 3)) - b= multadd(b, p05[i-1], 0, alloc); - - if (!(k>>= 2)) - return b; - p5= p5_a; - for (;;) - { - if (k & 1) - { - b1= mult(b, p5, alloc); - Bfree(b, alloc); - b= b1; - } - if (!(k>>= 1)) - break; - /* Calculate next power of 5 */ - if (overflow) - { - p51= mult(p5, p5, alloc); - Bfree(p5, alloc); - p5= p51; - } - else if (p5 < p5_a + P5A_MAX) - ++p5; - else if (p5 == p5_a + P5A_MAX) - { - p5= mult(p5, p5, alloc); - overflow= TRUE; - } - } - if (p51) - Bfree(p51, alloc); - return b; -} - - -static Bigint *lshift(Bigint *b, int k, Stack_alloc *alloc) -{ - int i, k1, n, n1; - Bigint *b1; - ULong *x, *x1, *xe, z; - - n= k >> 5; - k1= b->k; - n1= n + b->wds + 1; - for (i= b->maxwds; n1 > i; i<<= 1) - k1++; - b1= Balloc(k1, alloc); - x1= b1->p.x; - for (i= 0; i < n; i++) - *x1++= 0; - x= b->p.x; - xe= x + b->wds; - if (k&= 0x1f) - { - k1= 32 - k; - z= 0; - do - { - *x1++= *x << k | z; - z= *x++ >> k1; - } - while (x < xe); - if ((*x1= z)) - ++n1; - } - else - do - *x1++= *x++; - while (x < xe); - b1->wds= n1 - 1; - Bfree(b, alloc); - return b1; -} - - -static int cmp(Bigint *a, Bigint *b) -{ - ULong *xa, *xa0, *xb, *xb0; - int i, j; - - i= a->wds; - j= b->wds; - if (i-= j) - return i; - xa0= a->p.x; - xa= xa0 + j; - xb0= b->p.x; - xb= xb0 + j; - for (;;) - { - if (*--xa != *--xb) - return *xa < *xb ? -1 : 1; - if (xa <= xa0) - break; - } - return 0; -} - - -static Bigint *diff(Bigint *a, Bigint *b, Stack_alloc *alloc) -{ - Bigint *c; - int i, wa, wb; - ULong *xa, *xae, *xb, *xbe, *xc; - ULLong borrow, y; - - i= cmp(a,b); - if (!i) - { - c= Balloc(0, alloc); - c->wds= 1; - c->p.x[0]= 0; - return c; - } - if (i < 0) - { - c= a; - a= b; - b= c; - i= 1; - } - else - i= 0; - c= Balloc(a->k, alloc); - c->sign= i; - wa= a->wds; - xa= a->p.x; - xae= xa + wa; - wb= b->wds; - xb= b->p.x; - xbe= xb + wb; - xc= c->p.x; - borrow= 0; - do - { - y= (ULLong)*xa++ - *xb++ - borrow; - borrow= y >> 32 & (ULong)1; - *xc++= (ULong) (y & FFFFFFFF); - } - while (xb < xbe); - while (xa < xae) - { - y= *xa++ - borrow; - borrow= y >> 32 & (ULong)1; - *xc++= (ULong) (y & FFFFFFFF); - } - while (!*--xc) - wa--; - c->wds= wa; - return c; -} - - -static double ulp(U *x) -{ - Long L; - U u; - - L= (word0(x) & Exp_mask) - (P - 1)*Exp_msk1; - word0(&u) = L; - word1(&u) = 0; - return dval(&u); -} - - -static double b2d(Bigint *a, int *e) -{ - ULong *xa, *xa0, w, y, z; - int k; - U d; -#define d0 word0(&d) -#define d1 word1(&d) - - xa0= a->p.x; - xa= xa0 + a->wds; - y= *--xa; - k= hi0bits(y); - *e= 32 - k; - if (k < Ebits) - { - d0= Exp_1 | y >> (Ebits - k); - w= xa > xa0 ? *--xa : 0; - d1= y << ((32-Ebits) + k) | w >> (Ebits - k); - goto ret_d; - } - z= xa > xa0 ? *--xa : 0; - if (k-= Ebits) - { - d0= Exp_1 | y << k | z >> (32 - k); - y= xa > xa0 ? *--xa : 0; - d1= z << k | y >> (32 - k); - } - else - { - d0= Exp_1 | y; - d1= z; - } - ret_d: -#undef d0 -#undef d1 - return dval(&d); -} - - -static Bigint *d2b(U *d, int *e, int *bits, Stack_alloc *alloc) -{ - Bigint *b; - int de, k; - ULong *x, y, z; - int i; -#define d0 word0(d) -#define d1 word1(d) - - b= Balloc(1, alloc); - x= b->p.x; - - z= d0 & Frac_mask; - d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ - if ((de= (int)(d0 >> Exp_shift))) - z|= Exp_msk1; - if ((y= d1)) - { - if ((k= lo0bits(&y))) - { - x[0]= y | z << (32 - k); - z>>= k; - } - else - x[0]= y; - i= b->wds= (x[1]= z) ? 2 : 1; - } - else - { - k= lo0bits(&z); - x[0]= z; - i= b->wds= 1; - k+= 32; - } - if (de) - { - *e= de - Bias - (P-1) + k; - *bits= P - k; - } - else - { - *e= de - Bias - (P-1) + 1 + k; - *bits= 32*i - hi0bits(x[i-1]); - } - return b; -#undef d0 -#undef d1 -} - - -static double ratio(Bigint *a, Bigint *b) -{ - U da, db; - int k, ka, kb; - - dval(&da)= b2d(a, &ka); - dval(&db)= b2d(b, &kb); - k= ka - kb + 32*(a->wds - b->wds); - if (k > 0) - word0(&da)+= k*Exp_msk1; - else - { - k= -k; - word0(&db)+= k*Exp_msk1; - } - return dval(&da) / dval(&db); -} - -static const double tens[] = -{ - 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, - 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, - 1e20, 1e21, 1e22 -}; - -static const double bigtens[]= { 1e16, 1e32, 1e64, 1e128, 1e256 }; -static const double tinytens[]= -{ 1e-16, 1e-32, 1e-64, 1e-128, - 9007199254740992.*9007199254740992.e-256 /* = 2^106 * 1e-53 */ -}; -/* - The factor of 2^53 in tinytens[4] helps us avoid setting the underflow - flag unnecessarily. It leads to a song and dance at the end of strtod. -*/ -#define Scale_Bit 0x10 -#define n_bigtens 5 - -/* - strtod for IEEE--arithmetic machines. - - This strtod returns a nearest machine number to the input decimal - string (or sets errno to EOVERFLOW). Ties are broken by the IEEE round-even - rule. - - Inspired loosely by William D. Clinger's paper "How to Read Floating - Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. - - Modifications: - - 1. We only require IEEE (not IEEE double-extended). - 2. We get by with floating-point arithmetic in a case that - Clinger missed -- when we're computing d * 10^n - for a small integer d and the integer n is not too - much larger than 22 (the maximum integer k for which - we can represent 10^k exactly), we may be able to - compute (d*10^k) * 10^(e-k) with just one roundoff. - 3. Rather than a bit-at-a-time adjustment of the binary - result in the hard case, we use floating-point - arithmetic to determine the adjustment to within - one bit; only in really hard cases do we need to - compute a second residual. - 4. Because of 3., we don't need a large table of powers of 10 - for ten-to-e (just some small tables, e.g. of 10^k - for 0 <= k <= 22). -*/ - -static double my_strtod_int(const char *s00, char **se, int *error, char *buf, size_t buf_size) -{ - int scale; - int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c= 0, dsign, - e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; - const char *s, *s0, *s1, *end = *se; - double aadj, aadj1; - U aadj2, adj, rv, rv0; - Long L; - ULong y, z; - Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; -#ifdef Honor_FLT_ROUNDS - int rounding; -#endif - Stack_alloc alloc; - - *error= 0; - - alloc.begin= alloc.free= buf; - alloc.end= buf + buf_size; - memset(alloc.freelist, 0, sizeof(alloc.freelist)); - - sign= nz0= nz= 0; - dval(&rv)= 0.; - for (s= s00; s < end; s++) - switch (*s) { - case '-': - sign= 1; - /* no break */ - case '+': - s++; - goto break2; - case '\t': - case '\n': - case '\v': - case '\f': - case '\r': - case ' ': - continue; - default: - goto break2; - } - break2: - if (s >= end) - goto ret0; - - if (*s == '0') - { - nz0= 1; - while (++s < end && *s == '0') ; - if (s >= end) - goto ret; - } - s0= s; - y= z= 0; - for (nd= nf= 0; s < end && (c= *s) >= '0' && c <= '9'; nd++, s++) - if (nd < 9) - y= 10*y + c - '0'; - else if (nd < 16) - z= 10*z + c - '0'; - nd0= nd; - if (s < end && c == '.') - { - c= *++s; - if (!nd) - { - for (; s < end; ++s) - { - c= *s; - if (c != '0') - break; - nz++; - } - if (s < end && c > '0' && c <= '9') - { - s0= s; - nf+= nz; - nz= 0; - } - else - goto dig_done; - } - for (; s < end; ++s) - { - c= *s; - if (c < '0' || c > '9') - break; - /* - Here we are parsing the fractional part. - We can stop counting digits after a while: the extra digits - will not contribute to the actual result produced by s2b(). - We have to continue scanning, in case there is an exponent part. - */ - if (nd < 2 * DBL_DIG) - { - nz++; - if (c-= '0') - { - nf+= nz; - for (i= 1; i < nz; i++) - if (nd++ < 9) - y*= 10; - else if (nd <= DBL_DIG + 1) - z*= 10; - if (nd++ < 9) - y= 10*y + c; - else if (nd <= DBL_DIG + 1) - z= 10*z + c; - nz= 0; - } - } - } - } - dig_done: - e= 0; - if (s < end && (c == 'e' || c == 'E')) - { - if (!nd && !nz && !nz0) - goto ret0; - s00= s; - esign= 0; - if (++s < end) - switch (c= *s) { - case '-': - esign= 1; - case '+': - if (++s < end) - c= *s; - } - if (s < end && c >= '0' && c <= '9') - { - while (s < end && c == '0') - c= *++s; - if (s < end && c > '0' && c <= '9') { - L= c - '0'; - s1= s; - while (++s < end && (c= *s) >= '0' && c <= '9') - L= 10*L + c - '0'; - if (s - s1 > 8 || L > 19999) - /* Avoid confusion from exponents - * so large that e might overflow. - */ - e= 19999; /* safe for 16 bit ints */ - else - e= (int)L; - if (esign) - e= -e; - } - else - e= 0; - } - else - s= s00; - } - if (!nd) - { - if (!nz && !nz0) - { - ret0: - s= s00; - sign= 0; - } - goto ret; - } - e1= e -= nf; - - /* - Now we have nd0 digits, starting at s0, followed by a - decimal point, followed by nd-nd0 digits. The number we're - after is the integer represented by those digits times - 10**e - */ - - if (!nd0) - nd0= nd; - k= nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; - dval(&rv)= y; - if (k > 9) - { - dval(&rv)= tens[k - 9] * dval(&rv) + z; - } - bd0= 0; - if (nd <= DBL_DIG -#ifndef Honor_FLT_ROUNDS - && Flt_Rounds == 1 -#endif - ) - { - if (!e) - goto ret; - if (e > 0) - { - if (e <= Ten_pmax) - { -#ifdef Honor_FLT_ROUNDS - /* round correctly FLT_ROUNDS = 2 or 3 */ - if (sign) - { - rv.d= -rv.d; - sign= 0; - } -#endif - /* rv = */ rounded_product(dval(&rv), tens[e]); - goto ret; - } - i= DBL_DIG - nd; - if (e <= Ten_pmax + i) - { - /* - A fancier test would sometimes let us do - this for larger i values. - */ -#ifdef Honor_FLT_ROUNDS - /* round correctly FLT_ROUNDS = 2 or 3 */ - if (sign) - { - rv.d= -rv.d; - sign= 0; - } -#endif - e-= i; - dval(&rv)*= tens[i]; - /* rv = */ rounded_product(dval(&rv), tens[e]); - goto ret; - } - } -#ifndef Inaccurate_Divide - else if (e >= -Ten_pmax) - { -#ifdef Honor_FLT_ROUNDS - /* round correctly FLT_ROUNDS = 2 or 3 */ - if (sign) - { - rv.d= -rv.d; - sign= 0; - } -#endif - /* rv = */ rounded_quotient(dval(&rv), tens[-e]); - goto ret; - } -#endif - } - e1+= nd - k; - - scale= 0; -#ifdef Honor_FLT_ROUNDS - if ((rounding= Flt_Rounds) >= 2) - { - if (sign) - rounding= rounding == 2 ? 0 : 2; - else - if (rounding != 2) - rounding= 0; - } -#endif - - /* Get starting approximation = rv * 10**e1 */ - - if (e1 > 0) - { - if ((i= e1 & 15)) - dval(&rv)*= tens[i]; - if (e1&= ~15) - { - if (e1 > DBL_MAX_10_EXP) - { - ovfl: - *error= EOVERFLOW; - /* Can't trust HUGE_VAL */ -#ifdef Honor_FLT_ROUNDS - switch (rounding) - { - case 0: /* toward 0 */ - case 3: /* toward -infinity */ - word0(&rv)= Big0; - word1(&rv)= Big1; - break; - default: - word0(&rv)= Exp_mask; - word1(&rv)= 0; - } -#else /*Honor_FLT_ROUNDS*/ - word0(&rv)= Exp_mask; - word1(&rv)= 0; -#endif /*Honor_FLT_ROUNDS*/ - if (bd0) - goto retfree; - goto ret; - } - e1>>= 4; - for(j= 0; e1 > 1; j++, e1>>= 1) - if (e1 & 1) - dval(&rv)*= bigtens[j]; - /* The last multiplication could overflow. */ - word0(&rv)-= P*Exp_msk1; - dval(&rv)*= bigtens[j]; - if ((z= word0(&rv) & Exp_mask) > Exp_msk1 * (DBL_MAX_EXP + Bias - P)) - goto ovfl; - if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) - { - /* set to largest number (Can't trust DBL_MAX) */ - word0(&rv)= Big0; - word1(&rv)= Big1; - } - else - word0(&rv)+= P*Exp_msk1; - } - } - else if (e1 < 0) - { - e1= -e1; - if ((i= e1 & 15)) - dval(&rv)/= tens[i]; - if ((e1>>= 4)) - { - if (e1 >= 1 << n_bigtens) - goto undfl; - if (e1 & Scale_Bit) - scale= 2 * P; - for(j= 0; e1 > 0; j++, e1>>= 1) - if (e1 & 1) - dval(&rv)*= tinytens[j]; - if (scale && (j = 2 * P + 1 - ((word0(&rv) & Exp_mask) >> Exp_shift)) > 0) - { - /* scaled rv is denormal; zap j low bits */ - if (j >= 32) - { - word1(&rv)= 0; - if (j >= 53) - word0(&rv)= (P + 2) * Exp_msk1; - else - word0(&rv)&= 0xffffffff << (j - 32); - } - else - word1(&rv)&= 0xffffffff << j; - } - if (!dval(&rv)) - { - undfl: - dval(&rv)= 0.; - if (bd0) - goto retfree; - goto ret; - } - } - } - - /* Now the hard part -- adjusting rv to the correct value.*/ - - /* Put digits into bd: true value = bd * 10^e */ - - bd0= s2b(s0, nd0, nd, y, &alloc); - - for(;;) - { - bd= Balloc(bd0->k, &alloc); - Bcopy(bd, bd0); - bb= d2b(&rv, &bbe, &bbbits, &alloc); /* rv = bb * 2^bbe */ - bs= i2b(1, &alloc); - - if (e >= 0) - { - bb2= bb5= 0; - bd2= bd5= e; - } - else - { - bb2= bb5= -e; - bd2= bd5= 0; - } - if (bbe >= 0) - bb2+= bbe; - else - bd2-= bbe; - bs2= bb2; -#ifdef Honor_FLT_ROUNDS - if (rounding != 1) - bs2++; -#endif - j= bbe - scale; - i= j + bbbits - 1; /* logb(rv) */ - if (i < Emin) /* denormal */ - j+= P - Emin; - else - j= P + 1 - bbbits; - bb2+= j; - bd2+= j; - bd2+= scale; - i= bb2 < bd2 ? bb2 : bd2; - if (i > bs2) - i= bs2; - if (i > 0) - { - bb2-= i; - bd2-= i; - bs2-= i; - } - if (bb5 > 0) - { - bs= pow5mult(bs, bb5, &alloc); - bb1= mult(bs, bb, &alloc); - Bfree(bb, &alloc); - bb= bb1; - } - if (bb2 > 0) - bb= lshift(bb, bb2, &alloc); - if (bd5 > 0) - bd= pow5mult(bd, bd5, &alloc); - if (bd2 > 0) - bd= lshift(bd, bd2, &alloc); - if (bs2 > 0) - bs= lshift(bs, bs2, &alloc); - delta= diff(bb, bd, &alloc); - dsign= delta->sign; - delta->sign= 0; - i= cmp(delta, bs); -#ifdef Honor_FLT_ROUNDS - if (rounding != 1) - { - if (i < 0) - { - /* Error is less than an ulp */ - if (!delta->p.x[0] && delta->wds <= 1) - { - /* exact */ - break; - } - if (rounding) - { - if (dsign) - { - adj.d= 1.; - goto apply_adj; - } - } - else if (!dsign) - { - adj.d= -1.; - if (!word1(&rv) && !(word0(&rv) & Frac_mask)) - { - y= word0(&rv) & Exp_mask; - if (!scale || y > 2*P*Exp_msk1) - { - delta= lshift(delta, Log2P, &alloc); - if (cmp(delta, bs) <= 0) - adj.d= -0.5; - } - } - apply_adj: - if (scale && (y= word0(&rv) & Exp_mask) <= 2 * P * Exp_msk1) - word0(&adj)+= (2 * P + 1) * Exp_msk1 - y; - dval(&rv)+= adj.d * ulp(&rv); - } - break; - } - adj.d= ratio(delta, bs); - if (adj.d < 1.) - adj.d= 1.; - if (adj.d <= 0x7ffffffe) - { - /* adj = rounding ? ceil(adj) : floor(adj); */ - y= adj.d; - if (y != adj.d) - { - if (!((rounding >> 1) ^ dsign)) - y++; - adj.d= y; - } - } - if (scale && (y= word0(&rv) & Exp_mask) <= 2 * P * Exp_msk1) - word0(&adj)+= (2 * P + 1) * Exp_msk1 - y; - adj.d*= ulp(&rv); - if (dsign) - dval(&rv)+= adj.d; - else - dval(&rv)-= adj.d; - goto cont; - } -#endif /*Honor_FLT_ROUNDS*/ - - if (i < 0) - { - /* - Error is less than half an ulp -- check for special case of mantissa - a power of two. - */ - if (dsign || word1(&rv) || word0(&rv) & Bndry_mask || - (word0(&rv) & Exp_mask) <= (2 * P + 1) * Exp_msk1) - { - break; - } - if (!delta->p.x[0] && delta->wds <= 1) - { - /* exact result */ - break; - } - delta= lshift(delta, Log2P, &alloc); - if (cmp(delta, bs) > 0) - goto drop_down; - break; - } - if (i == 0) - { - /* exactly half-way between */ - if (dsign) - { - if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 && - word1(&rv) == - ((scale && (y = word0(&rv) & Exp_mask) <= 2 * P * Exp_msk1) ? - (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : - 0xffffffff)) - { - /*boundary case -- increment exponent*/ - word0(&rv)= (word0(&rv) & Exp_mask) + Exp_msk1; - word1(&rv) = 0; - dsign = 0; - break; - } - } - else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) - { - drop_down: - /* boundary case -- decrement exponent */ - if (scale) - { - L= word0(&rv) & Exp_mask; - if (L <= (2 *P + 1) * Exp_msk1) - { - if (L > (P + 2) * Exp_msk1) - /* round even ==> accept rv */ - break; - /* rv = smallest denormal */ - goto undfl; - } - } - L= (word0(&rv) & Exp_mask) - Exp_msk1; - word0(&rv)= L | Bndry_mask1; - word1(&rv)= 0xffffffff; - break; - } - if (!(word1(&rv) & LSB)) - break; - if (dsign) - dval(&rv)+= ulp(&rv); - else - { - dval(&rv)-= ulp(&rv); - if (!dval(&rv)) - goto undfl; - } - dsign= 1 - dsign; - break; - } - if ((aadj= ratio(delta, bs)) <= 2.) - { - if (dsign) - aadj= aadj1= 1.; - else if (word1(&rv) || word0(&rv) & Bndry_mask) - { - if (word1(&rv) == Tiny1 && !word0(&rv)) - goto undfl; - aadj= 1.; - aadj1= -1.; - } - else - { - /* special case -- power of FLT_RADIX to be rounded down... */ - if (aadj < 2. / FLT_RADIX) - aadj= 1. / FLT_RADIX; - else - aadj*= 0.5; - aadj1= -aadj; - } - } - else - { - aadj*= 0.5; - aadj1= dsign ? aadj : -aadj; -#ifdef Check_FLT_ROUNDS - switch (Rounding) - { - case 2: /* towards +infinity */ - aadj1-= 0.5; - break; - case 0: /* towards 0 */ - case 3: /* towards -infinity */ - aadj1+= 0.5; - } -#else - if (Flt_Rounds == 0) - aadj1+= 0.5; -#endif /*Check_FLT_ROUNDS*/ - } - y= word0(&rv) & Exp_mask; - - /* Check for overflow */ - - if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1)) - { - dval(&rv0)= dval(&rv); - word0(&rv)-= P * Exp_msk1; - adj.d= aadj1 * ulp(&rv); - dval(&rv)+= adj.d; - if ((word0(&rv) & Exp_mask) >= Exp_msk1 * (DBL_MAX_EXP + Bias - P)) - { - if (word0(&rv0) == Big0 && word1(&rv0) == Big1) - goto ovfl; - word0(&rv)= Big0; - word1(&rv)= Big1; - goto cont; - } - else - word0(&rv)+= P * Exp_msk1; - } - else - { - if (scale && y <= 2 * P * Exp_msk1) - { - if (aadj <= 0x7fffffff) - { - if ((z= (ULong) aadj) <= 0) - z= 1; - aadj= z; - aadj1= dsign ? aadj : -aadj; - } - dval(&aadj2) = aadj1; - word0(&aadj2)+= (2 * P + 1) * Exp_msk1 - y; - aadj1= dval(&aadj2); - adj.d= aadj1 * ulp(&rv); - dval(&rv)+= adj.d; - if (rv.d == 0.) - goto undfl; - } - else - { - adj.d= aadj1 * ulp(&rv); - dval(&rv)+= adj.d; - } - } - z= word0(&rv) & Exp_mask; - if (!scale) - if (y == z) - { - /* Can we stop now? */ - L= (Long)aadj; - aadj-= L; - /* The tolerances below are conservative. */ - if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) - { - if (aadj < .4999999 || aadj > .5000001) - break; - } - else if (aadj < .4999999 / FLT_RADIX) - break; - } - cont: - Bfree(bb, &alloc); - Bfree(bd, &alloc); - Bfree(bs, &alloc); - Bfree(delta, &alloc); - } - if (scale) - { - word0(&rv0)= Exp_1 - 2 * P * Exp_msk1; - word1(&rv0)= 0; - dval(&rv)*= dval(&rv0); - } - retfree: - Bfree(bb, &alloc); - Bfree(bd, &alloc); - Bfree(bs, &alloc); - Bfree(bd0, &alloc); - Bfree(delta, &alloc); - ret: - *se= (char *)s; - return sign ? -dval(&rv) : dval(&rv); -} - - -static int quorem(Bigint *b, Bigint *S) -{ - int n; - ULong *bx, *bxe, q, *sx, *sxe; - ULLong borrow, carry, y, ys; - - n= S->wds; - if (b->wds < n) - return 0; - sx= S->p.x; - sxe= sx + --n; - bx= b->p.x; - bxe= bx + n; - q= *bxe / (*sxe + 1); /* ensure q <= true quotient */ - if (q) - { - borrow= 0; - carry= 0; - do - { - ys= *sx++ * (ULLong)q + carry; - carry= ys >> 32; - y= *bx - (ys & FFFFFFFF) - borrow; - borrow= y >> 32 & (ULong)1; - *bx++= (ULong) (y & FFFFFFFF); - } - while (sx <= sxe); - if (!*bxe) - { - bx= b->p.x; - while (--bxe > bx && !*bxe) - --n; - b->wds= n; - } - } - if (cmp(b, S) >= 0) - { - q++; - borrow= 0; - carry= 0; - bx= b->p.x; - sx= S->p.x; - do - { - ys= *sx++ + carry; - carry= ys >> 32; - y= *bx - (ys & FFFFFFFF) - borrow; - borrow= y >> 32 & (ULong)1; - *bx++= (ULong) (y & FFFFFFFF); - } - while (sx <= sxe); - bx= b->p.x; - bxe= bx + n; - if (!*bxe) - { - while (--bxe > bx && !*bxe) - --n; - b->wds= n; - } - } - return q; -} - - -/* - dtoa for IEEE arithmetic (dmg): convert double to ASCII string. - - Inspired by "How to Print Floating-Point Numbers Accurately" by - Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. - - Modifications: - 1. Rather than iterating, we use a simple numeric overestimate - to determine k= floor(log10(d)). We scale relevant - quantities using O(log2(k)) rather than O(k) multiplications. - 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't - try to generate digits strictly left to right. Instead, we - compute with fewer bits and propagate the carry if necessary - when rounding the final digit up. This is often faster. - 3. Under the assumption that input will be rounded nearest, - mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. - That is, we allow equality in stopping tests when the - round-nearest rule will give the same floating-point value - as would satisfaction of the stopping test with strict - inequality. - 4. We remove common factors of powers of 2 from relevant - quantities. - 5. When converting floating-point integers less than 1e16, - we use floating-point arithmetic rather than resorting - to multiple-precision integers. - 6. When asked to produce fewer than 15 digits, we first try - to get by with floating-point arithmetic; we resort to - multiple-precision integer arithmetic only if we cannot - guarantee that the floating-point calculation has given - the correctly rounded result. For k requested digits and - "uniformly" distributed input, the probability is - something like 10^(k-15) that we must resort to the Long - calculation. - */ - -static char *dtoa(double dd, int mode, int ndigits, int *decpt, int *sign, - char **rve, char *buf, size_t buf_size) -{ - /* - Arguments ndigits, decpt, sign are similar to those - of ecvt and fcvt; trailing zeros are suppressed from - the returned string. If not null, *rve is set to point - to the end of the return value. If d is +-Infinity or NaN, - then *decpt is set to DTOA_OVERFLOW. - - mode: - 0 ==> shortest string that yields d when read in - and rounded to nearest. - 1 ==> like 0, but with Steele & White stopping rule; - e.g. with IEEE P754 arithmetic , mode 0 gives - 1e23 whereas mode 1 gives 9.999999999999999e22. - 2 ==> max(1,ndigits) significant digits. This gives a - return value similar to that of ecvt, except - that trailing zeros are suppressed. - 3 ==> through ndigits past the decimal point. This - gives a return value similar to that from fcvt, - except that trailing zeros are suppressed, and - ndigits can be negative. - 4,5 ==> similar to 2 and 3, respectively, but (in - round-nearest mode) with the tests of mode 0 to - possibly return a shorter string that rounds to d. - With IEEE arithmetic and compilation with - -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same - as modes 2 and 3 when FLT_ROUNDS != 1. - 6-9 ==> Debugging modes similar to mode - 4: don't try - fast floating-point estimate (if applicable). - - Values of mode other than 0-9 are treated as mode 0. - - Sufficient space is allocated to the return value - to hold the suppressed trailing zeros. - */ - - int bbits, b2, b5, be, dig, i, ieps, ilim= 0, ilim0, - ilim1= 0, j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, - spec_case, try_quick; - Long L; - int denorm; - ULong x; - Bigint *b, *b1, *delta, *mlo, *mhi, *S; - U d2, eps, u; - double ds; - char *s, *s0; -#ifdef Honor_FLT_ROUNDS - int rounding; -#endif - Stack_alloc alloc; - - alloc.begin= alloc.free= buf; - alloc.end= buf + buf_size; - memset(alloc.freelist, 0, sizeof(alloc.freelist)); - - u.d= dd; - if (word0(&u) & Sign_bit) - { - /* set sign for everything, including 0's and NaNs */ - *sign= 1; - word0(&u) &= ~Sign_bit; /* clear sign bit */ - } - else - *sign= 0; - - /* If infinity, set decpt to DTOA_OVERFLOW, if 0 set it to 1 */ - if (((word0(&u) & Exp_mask) == Exp_mask && (*decpt= DTOA_OVERFLOW)) || - (!dval(&u) && (*decpt= 1))) - { - /* Infinity, NaN, 0 */ - char *res= (char*) dtoa_alloc(2, &alloc); - res[0]= '0'; - res[1]= '\0'; - if (rve) - *rve= res + 1; - return res; - } - -#ifdef Honor_FLT_ROUNDS - if ((rounding= Flt_Rounds) >= 2) - { - if (*sign) - rounding= rounding == 2 ? 0 : 2; - else - if (rounding != 2) - rounding= 0; - } -#endif - - b= d2b(&u, &be, &bbits, &alloc); - if ((i= (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) - { - dval(&d2)= dval(&u); - word0(&d2) &= Frac_mask1; - word0(&d2) |= Exp_11; - - /* - log(x) ~=~ log(1.5) + (x-1.5)/1.5 - log10(x) = log(x) / log(10) - ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) - log10(d)= (i-Bias)*log(2)/log(10) + log10(d2) - - This suggests computing an approximation k to log10(d) by - - k= (i - Bias)*0.301029995663981 - + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); - - We want k to be too large rather than too small. - The error in the first-order Taylor series approximation - is in our favor, so we just round up the constant enough - to compensate for any error in the multiplication of - (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, - and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, - adding 1e-13 to the constant term more than suffices. - Hence we adjust the constant term to 0.1760912590558. - (We could get a more accurate k by invoking log10, - but this is probably not worthwhile.) - */ - - i-= Bias; - denorm= 0; - } - else - { - /* d is denormalized */ - - i= bbits + be + (Bias + (P-1) - 1); - x= i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32) - : word1(&u) << (32 - i); - dval(&d2)= x; - word0(&d2)-= 31*Exp_msk1; /* adjust exponent */ - i-= (Bias + (P-1) - 1) + 1; - denorm= 1; - } - ds= (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; - k= (int)ds; - if (ds < 0. && ds != k) - k--; /* want k= floor(ds) */ - k_check= 1; - if (k >= 0 && k <= Ten_pmax) - { - if (dval(&u) < tens[k]) - k--; - k_check= 0; - } - j= bbits - i - 1; - if (j >= 0) - { - b2= 0; - s2= j; - } - else - { - b2= -j; - s2= 0; - } - if (k >= 0) - { - b5= 0; - s5= k; - s2+= k; - } - else - { - b2-= k; - b5= -k; - s5= 0; - } - if (mode < 0 || mode > 9) - mode= 0; - -#ifdef Check_FLT_ROUNDS - try_quick= Rounding == 1; -#else - try_quick= 1; -#endif - - if (mode > 5) - { - mode-= 4; - try_quick= 0; - } - leftright= 1; - switch (mode) { - case 0: - case 1: - ilim= ilim1= -1; - i= 18; - ndigits= 0; - break; - case 2: - leftright= 0; - /* no break */ - case 4: - if (ndigits <= 0) - ndigits= 1; - ilim= ilim1= i= ndigits; - break; - case 3: - leftright= 0; - /* no break */ - case 5: - i= ndigits + k + 1; - ilim= i; - ilim1= i - 1; - if (i <= 0) - i= 1; - } - s= s0= dtoa_alloc(i, &alloc); - -#ifdef Honor_FLT_ROUNDS - if (mode > 1 && rounding != 1) - leftright= 0; -#endif - - if (ilim >= 0 && ilim <= Quick_max && try_quick) - { - /* Try to get by with floating-point arithmetic. */ - i= 0; - dval(&d2)= dval(&u); - k0= k; - ilim0= ilim; - ieps= 2; /* conservative */ - if (k > 0) - { - ds= tens[k&0xf]; - j= k >> 4; - if (j & Bletch) - { - /* prevent overflows */ - j&= Bletch - 1; - dval(&u)/= bigtens[n_bigtens-1]; - ieps++; - } - for (; j; j>>= 1, i++) - { - if (j & 1) - { - ieps++; - ds*= bigtens[i]; - } - } - dval(&u)/= ds; - } - else if ((j1= -k)) - { - dval(&u)*= tens[j1 & 0xf]; - for (j= j1 >> 4; j; j>>= 1, i++) - { - if (j & 1) - { - ieps++; - dval(&u)*= bigtens[i]; - } - } - } - if (k_check && dval(&u) < 1. && ilim > 0) - { - if (ilim1 <= 0) - goto fast_failed; - ilim= ilim1; - k--; - dval(&u)*= 10.; - ieps++; - } - dval(&eps)= ieps*dval(&u) + 7.; - word0(&eps)-= (P-1)*Exp_msk1; - if (ilim == 0) - { - S= mhi= 0; - dval(&u)-= 5.; - if (dval(&u) > dval(&eps)) - goto one_digit; - if (dval(&u) < -dval(&eps)) - goto no_digits; - goto fast_failed; - } - if (leftright) - { - /* Use Steele & White method of only generating digits needed. */ - dval(&eps)= 0.5/tens[ilim-1] - dval(&eps); - for (i= 0;;) - { - L= (Long) dval(&u); - dval(&u)-= L; - *s++= '0' + (int)L; - if (dval(&u) < dval(&eps)) - goto ret1; - if (1. - dval(&u) < dval(&eps)) - goto bump_up; - if (++i >= ilim) - break; - dval(&eps)*= 10.; - dval(&u)*= 10.; - } - } - else - { - /* Generate ilim digits, then fix them up. */ - dval(&eps)*= tens[ilim-1]; - for (i= 1;; i++, dval(&u)*= 10.) - { - L= (Long)(dval(&u)); - if (!(dval(&u)-= L)) - ilim= i; - *s++= '0' + (int)L; - if (i == ilim) - { - if (dval(&u) > 0.5 + dval(&eps)) - goto bump_up; - else if (dval(&u) < 0.5 - dval(&eps)) - { - while (*--s == '0'); - s++; - goto ret1; - } - break; - } - } - } - fast_failed: - s= s0; - dval(&u)= dval(&d2); - k= k0; - ilim= ilim0; - } - - /* Do we have a "small" integer? */ - - if (be >= 0 && k <= Int_max) - { - /* Yes. */ - ds= tens[k]; - if (ndigits < 0 && ilim <= 0) - { - S= mhi= 0; - if (ilim < 0 || dval(&u) <= 5*ds) - goto no_digits; - goto one_digit; - } - for (i= 1;; i++, dval(&u)*= 10.) - { - L= (Long)(dval(&u) / ds); - dval(&u)-= L*ds; -#ifdef Check_FLT_ROUNDS - /* If FLT_ROUNDS == 2, L will usually be high by 1 */ - if (dval(&u) < 0) - { - L--; - dval(&u)+= ds; - } -#endif - *s++= '0' + (int)L; - if (!dval(&u)) - { - break; - } - if (i == ilim) - { -#ifdef Honor_FLT_ROUNDS - if (mode > 1) - { - switch (rounding) { - case 0: goto ret1; - case 2: goto bump_up; - } - } -#endif - dval(&u)+= dval(&u); - if (dval(&u) > ds || (dval(&u) == ds && L & 1)) - { -bump_up: - while (*--s == '9') - if (s == s0) - { - k++; - *s= '0'; - break; - } - ++*s++; - } - break; - } - } - goto ret1; - } - - m2= b2; - m5= b5; - mhi= mlo= 0; - if (leftright) - { - i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits; - b2+= i; - s2+= i; - mhi= i2b(1, &alloc); - } - if (m2 > 0 && s2 > 0) - { - i= m2 < s2 ? m2 : s2; - b2-= i; - m2-= i; - s2-= i; - } - if (b5 > 0) - { - if (leftright) - { - if (m5 > 0) - { - mhi= pow5mult(mhi, m5, &alloc); - b1= mult(mhi, b, &alloc); - Bfree(b, &alloc); - b= b1; - } - if ((j= b5 - m5)) - b= pow5mult(b, j, &alloc); - } - else - b= pow5mult(b, b5, &alloc); - } - S= i2b(1, &alloc); - if (s5 > 0) - S= pow5mult(S, s5, &alloc); - - /* Check for special case that d is a normalized power of 2. */ - - spec_case= 0; - if ((mode < 2 || leftright) -#ifdef Honor_FLT_ROUNDS - && rounding == 1 -#endif - ) - { - if (!word1(&u) && !(word0(&u) & Bndry_mask) && - word0(&u) & (Exp_mask & ~Exp_msk1) - ) - { - /* The special case */ - b2+= Log2P; - s2+= Log2P; - spec_case= 1; - } - } - - /* - Arrange for convenient computation of quotients: - shift left if necessary so divisor has 4 leading 0 bits. - - Perhaps we should just compute leading 28 bits of S once - a nd for all and pass them and a shift to quorem, so it - can do shifts and ors to compute the numerator for q. - */ - if ((i= ((s5 ? 32 - hi0bits(S->p.x[S->wds-1]) : 1) + s2) & 0x1f)) - i= 32 - i; - if (i > 4) - { - i-= 4; - b2+= i; - m2+= i; - s2+= i; - } - else if (i < 4) - { - i+= 28; - b2+= i; - m2+= i; - s2+= i; - } - if (b2 > 0) - b= lshift(b, b2, &alloc); - if (s2 > 0) - S= lshift(S, s2, &alloc); - if (k_check) - { - if (cmp(b,S) < 0) - { - k--; - /* we botched the k estimate */ - b= multadd(b, 10, 0, &alloc); - if (leftright) - mhi= multadd(mhi, 10, 0, &alloc); - ilim= ilim1; - } - } - if (ilim <= 0 && (mode == 3 || mode == 5)) - { - if (ilim < 0 || cmp(b,S= multadd(S,5,0, &alloc)) <= 0) - { - /* no digits, fcvt style */ -no_digits: - k= -1 - ndigits; - goto ret; - } -one_digit: - *s++= '1'; - k++; - goto ret; - } - if (leftright) - { - if (m2 > 0) - mhi= lshift(mhi, m2, &alloc); - - /* - Compute mlo -- check for special case that d is a normalized power of 2. - */ - - mlo= mhi; - if (spec_case) - { - mhi= Balloc(mhi->k, &alloc); - Bcopy(mhi, mlo); - mhi= lshift(mhi, Log2P, &alloc); - } - - for (i= 1;;i++) - { - dig= quorem(b,S) + '0'; - /* Do we yet have the shortest decimal string that will round to d? */ - j= cmp(b, mlo); - delta= diff(S, mhi, &alloc); - j1= delta->sign ? 1 : cmp(b, delta); - Bfree(delta, &alloc); - if (j1 == 0 && mode != 1 && !(word1(&u) & 1) -#ifdef Honor_FLT_ROUNDS - && rounding >= 1 -#endif - ) - { - if (dig == '9') - goto round_9_up; - if (j > 0) - dig++; - *s++= dig; - goto ret; - } - if (j < 0 || (j == 0 && mode != 1 && !(word1(&u) & 1))) - { - if (!b->p.x[0] && b->wds <= 1) - { - goto accept_dig; - } -#ifdef Honor_FLT_ROUNDS - if (mode > 1) - switch (rounding) { - case 0: goto accept_dig; - case 2: goto keep_dig; - } -#endif /*Honor_FLT_ROUNDS*/ - if (j1 > 0) - { - b= lshift(b, 1, &alloc); - j1= cmp(b, S); - if ((j1 > 0 || (j1 == 0 && dig & 1)) - && dig++ == '9') - goto round_9_up; - } -accept_dig: - *s++= dig; - goto ret; - } - if (j1 > 0) - { -#ifdef Honor_FLT_ROUNDS - if (!rounding) - goto accept_dig; -#endif - if (dig == '9') - { /* possible if i == 1 */ -round_9_up: - *s++= '9'; - goto roundoff; - } - *s++= dig + 1; - goto ret; - } -#ifdef Honor_FLT_ROUNDS -keep_dig: -#endif - *s++= dig; - if (i == ilim) - break; - b= multadd(b, 10, 0, &alloc); - if (mlo == mhi) - mlo= mhi= multadd(mhi, 10, 0, &alloc); - else - { - mlo= multadd(mlo, 10, 0, &alloc); - mhi= multadd(mhi, 10, 0, &alloc); - } - } - } - else - for (i= 1;; i++) - { - *s++= dig= quorem(b,S) + '0'; - if (!b->p.x[0] && b->wds <= 1) - { - goto ret; - } - if (i >= ilim) - break; - b= multadd(b, 10, 0, &alloc); - } - - /* Round off last digit */ - -#ifdef Honor_FLT_ROUNDS - switch (rounding) { - case 0: goto trimzeros; - case 2: goto roundoff; - } -#endif - b= lshift(b, 1, &alloc); - j= cmp(b, S); - if (j > 0 || (j == 0 && dig & 1)) - { -roundoff: - while (*--s == '9') - if (s == s0) - { - k++; - *s++= '1'; - goto ret; - } - ++*s++; - } - else - { -#ifdef Honor_FLT_ROUNDS -trimzeros: -#endif - while (*--s == '0'); - s++; - } -ret: - Bfree(S, &alloc); - if (mhi) - { - if (mlo && mlo != mhi) - Bfree(mlo, &alloc); - Bfree(mhi, &alloc); - } -ret1: - Bfree(b, &alloc); - *s= 0; - *decpt= k + 1; - if (rve) - *rve= s; - return s0; -} -- cgit v1.1