summaryrefslogtreecommitdiff
path: root/mysql/strings/dtoa.c
diff options
context:
space:
mode:
Diffstat (limited to 'mysql/strings/dtoa.c')
-rw-r--r--mysql/strings/dtoa.c2784
1 files changed, 0 insertions, 2784 deletions
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 <my_global.h>
-#include <m_string.h> /* 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;
-}
generated by cgit v1.1 at 2024-04-26 17:09:16 +0000