Drizzled Public API Documentation

dtoa.cc
00001 /* Copyright (C) 2007 MySQL AB
00002    This program is free software; you can redistribute it and/or modify
00003    it under the terms of the GNU General Public License as published by
00004    the Free Software Foundation; either version 2 of the License, or
00005    (at your option) any later version.
00006 
00007    This program is distributed in the hope that it will be useful,
00008    but WITHOUT ANY WARRANTY; without even the implied warranty of
00009    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00010    GNU General Public License for more details.
00011 
00012    You should have received a copy of the GNU General Public License
00013    along with this program; if not, write to the Free Software
00014    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA */
00015 
00016 /****************************************************************
00017 
00018   This file incorporates work covered by the following copyright and
00019   permission notice:
00020 
00021   The author of this software is David M. Gay.
00022 
00023   Copyright (C) 1991, 2000, 2001 by Lucent Technologies.
00024 
00025   Permission to use, copy, modify, and distribute this software for any
00026   purpose without fee is hereby granted, provided that this entire notice
00027   is included in all copies of any software which is or includes a copy
00028   or modification of this software and in all copies of the supporting
00029   documentation for such software.
00030 
00031   THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
00032   WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
00033   REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
00034   OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
00035 
00036  ***************************************************************/
00037 
00038 #include <config.h>
00039 
00040 #include <drizzled/internal/m_string.h>  /* for memcpy and NOT_FIXED_DEC */
00041 
00042 #include <float.h>
00043 
00044 #include <cstdlib>
00045 #include <cerrno>
00046 #include <algorithm>
00047 
00048 using namespace std;
00049 
00050 namespace drizzled
00051 {
00052 namespace internal
00053 {
00054 
00055 /* Magic value returned by dtoa() to indicate overflow */
00056 #define DTOA_OVERFLOW 9999
00057 
00058 static double my_strtod_int(const char *, char **, int *);
00059 static char *dtoa(double, int, int, int *, int *, char **);
00060 
00092 size_t my_fcvt(double x, int precision, char *to, bool *error)
00093 {
00094   int decpt, sign, len, i;
00095   char *res, *src, *end, *dst= to;
00096   assert(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
00097 
00098   res= dtoa(x, 5, precision, &decpt, &sign, &end);
00099 
00100   if (decpt == DTOA_OVERFLOW)
00101   {
00102     free(res);
00103     *to++= '0';
00104     *to= '\0';
00105     if (error != NULL)
00106       *error= true;
00107     return 1;
00108   }
00109 
00110   src= res;
00111   len= end - src;
00112 
00113   if (sign)
00114     *dst++= '-';
00115 
00116   if (decpt <= 0)
00117   {
00118     *dst++= '0';
00119     *dst++= '.';
00120     for (i= decpt; i < 0; i++)
00121       *dst++= '0';
00122   }
00123 
00124   for (i= 1; i <= len; i++)
00125   {
00126     *dst++= *src++;
00127     if (i == decpt && i < len)
00128       *dst++= '.';
00129   }
00130   while (i++ <= decpt)
00131     *dst++= '0';
00132 
00133   if (precision > 0)
00134   {
00135     if (len <= decpt)
00136       *dst++= '.';
00137 
00138     for (i= precision - max(0, (len - decpt)); i > 0; i--)
00139       *dst++= '0';
00140   }
00141 
00142   *dst= '\0';
00143   if (error != NULL)
00144     *error= false;
00145 
00146   free(res);
00147 
00148   return dst - to;
00149 }
00150 
00214 size_t my_gcvt(double x, my_gcvt_arg_type type, int width, char *to,
00215                bool *error)
00216 {
00217   int decpt, sign, len, exp_len;
00218   char *res, *src, *end, *dst= to, *dend= dst + width;
00219   bool have_space, force_e_format;
00220   assert(width > 0 && to != NULL);
00221 
00222   /* We want to remove '-' from equations early */
00223   if (x < 0.)
00224     width--;
00225 
00226   res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? min(width, DBL_DIG) : 
00227             min(width, FLT_DIG), &decpt, &sign, &end);
00228 
00229   if (decpt == DTOA_OVERFLOW)
00230   {
00231     free(res);
00232     *to++= '0';
00233     *to= '\0';
00234     if (error != NULL)
00235       *error= true;
00236     return 1;
00237   }
00238 
00239   if (error != NULL)
00240     *error= false;
00241 
00242   src= res;
00243   len= end - res;
00244 
00245   /*
00246     Number of digits in the exponent from the 'e' conversion.
00247      The sign of the exponent is taken into account separetely, we don't need
00248      to count it here.
00249    */
00250   exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
00251 
00252   /*
00253      Do we have enough space for all digits in the 'f' format?
00254      Let 'len' be the number of significant digits returned by dtoa,
00255      and F be the length of the resulting decimal representation.
00256      Consider the following cases:
00257      1. decpt <= 0, i.e. we have "0.NNN" => F = len - decpt + 2
00258      2. 0 < decpt < len, i.e. we have "NNN.NNN" => F = len + 1
00259      3. len <= decpt, i.e. we have "NNN00" => F = decpt
00260   */
00261   have_space= (decpt <= 0 ? len - decpt + 2 :
00262                decpt > 0 && decpt < len ? len + 1 :
00263                decpt) <= width;
00264   /*
00265     The following is true when no significant digits can be placed with the
00266     specified field width using the 'f' format, and the 'e' format
00267     will not be truncated.
00268   */
00269   force_e_format= (decpt <= 0 && width <= 2 - decpt && width >= 3 + exp_len);
00270   /*
00271     Assume that we don't have enough space to place all significant digits in
00272     the 'f' format. We have to choose between the 'e' format and the 'f' one
00273     to keep as many significant digits as possible.
00274     Let E and F be the lengths of decimal representaion in the 'e' and 'f'
00275     formats, respectively. We want to use the 'f' format if, and only if F <= E.
00276     Consider the following cases:
00277     1. decpt <= 0.
00278        F = len - decpt + 2 (see above)
00279        E = len + (len > 1) + 1 + 1 (decpt <= -99) + (decpt <= -9) + 1
00280        ("N.NNe-MMM")
00281        (F <= E) <=> (len == 1 && decpt >= -1) || (len > 1 && decpt >= -2)
00282        We also need to ensure that if the 'f' format is chosen,
00283        the field width allows us to place at least one significant digit
00284        (i.e. width > 2 - decpt). If not, we prefer the 'e' format.
00285     2. 0 < decpt < len
00286        F = len + 1 (see above)
00287        E = len + 1 + 1 + ... ("N.NNeMMM")
00288        F is always less than E.
00289     3. len <= decpt <= width
00290        In this case we have enough space to represent the number in the 'f'
00291        format, so we prefer it with some exceptions.
00292     4. width < decpt
00293        The number cannot be represented in the 'f' format at all, always use
00294        the 'e' 'one.
00295   */
00296   if ((have_space ||
00297       /*
00298         Not enough space, let's see if the 'f' format provides the most number
00299         of significant digits.
00300       */
00301        ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
00302                                             (len > 1 || !force_e_format)))) &&
00303          !force_e_format)) &&
00304 
00305        /*
00306          Use the 'e' format in some cases even if we have enough space for the
00307          'f' one. See comment for MAX_DECPT_FOR_F_FORMAT.
00308        */
00309       (!have_space || (decpt >= -MAX_DECPT_FOR_F_FORMAT + 1 &&
00310                        (decpt <= MAX_DECPT_FOR_F_FORMAT || len > decpt))))
00311   {
00312     /* 'f' format */
00313     int i;
00314 
00315     width-= (decpt < len) + (decpt <= 0 ? 1 - decpt : 0);
00316 
00317     /* Do we have to truncate any digits? */
00318     if (width < len)
00319     {
00320       if (width < decpt)
00321       {
00322         if (error != NULL)
00323           *error= true;
00324         width= decpt;
00325       }
00326 
00327       /*
00328         We want to truncate (len - width) least significant digits after the
00329         decimal point. For this we are calling dtoa with mode=5, passing the
00330         number of significant digits = (len-decpt) - (len-width) = width-decpt
00331       */
00332       free(res);
00333       res= dtoa(x, 5, width - decpt, &decpt, &sign, &end);
00334       src= res;
00335       len= end - res;
00336     }
00337 
00338     if (len == 0)
00339     {
00340       /* Underflow. Just print '0' and exit */
00341       *dst++= '0';
00342       goto end;
00343     }
00344 
00345     /*
00346       At this point we are sure we have enough space to put all digits
00347       returned by dtoa
00348     */
00349     if (sign && dst < dend)
00350       *dst++= '-';
00351     if (decpt <= 0)
00352     {
00353       if (dst < dend)
00354         *dst++= '0';
00355       if (len > 0 && dst < dend)
00356         *dst++= '.';
00357       for (; decpt < 0 && dst < dend; decpt++)
00358         *dst++= '0';
00359     }
00360 
00361     for (i= 1; i <= len && dst < dend; i++)
00362     {
00363       *dst++= *src++;
00364       if (i == decpt && i < len && dst < dend)
00365         *dst++= '.';
00366     }
00367     while (i++ <= decpt && dst < dend)
00368       *dst++= '0';
00369   }
00370   else
00371   {
00372     /* 'e' format */
00373     int decpt_sign= 0;
00374 
00375     if (--decpt < 0)
00376     {
00377       decpt= -decpt;
00378       width--;
00379       decpt_sign= 1;
00380     }
00381     width-= 1 + exp_len; /* eNNN */
00382 
00383     if (len > 1)
00384       width--;
00385 
00386     if (width <= 0)
00387     {
00388       /* Overflow */
00389       if (error != NULL)
00390         *error= true;
00391       width= 0;
00392     }
00393 
00394     /* Do we have to truncate any digits? */
00395     if (width < len)
00396     {
00397       /* Yes, re-convert with a smaller width */
00398       free(res);
00399       res= dtoa(x, 4, width, &decpt, &sign, &end);
00400       src= res;
00401       len= end - res;
00402       if (--decpt < 0)
00403         decpt= -decpt;
00404     }
00405     /*
00406       At this point we are sure we have enough space to put all digits
00407       returned by dtoa
00408     */
00409     if (sign && dst < dend)
00410       *dst++= '-';
00411     if (dst < dend)
00412       *dst++= *src++;
00413     if (len > 1 && dst < dend)
00414     {
00415       *dst++= '.';
00416       while (src < end && dst < dend)
00417         *dst++= *src++;
00418     }
00419     if (dst < dend)
00420       *dst++= 'e';
00421     if (decpt_sign && dst < dend)
00422       *dst++= '-';
00423 
00424     if (decpt >= 100 && dst < dend)
00425     {
00426       *dst++= decpt / 100 + '0';
00427       decpt%= 100;
00428       if (dst < dend)
00429         *dst++= decpt / 10 + '0';
00430     }
00431     else if (decpt >= 10 && dst < dend)
00432       *dst++= decpt / 10 + '0';
00433     if (dst < dend)
00434       *dst++= decpt % 10 + '0';
00435 
00436   }
00437 
00438 end:
00439   free(res);
00440   *dst= '\0';
00441 
00442   return dst - to;
00443 }
00444 
00463 double my_strtod(const char *str, char **end, int *error)
00464 {
00465   double res;
00466   assert(str != NULL && end != NULL && *end != NULL && error != NULL);
00467 
00468   res= my_strtod_int(str, end, error);
00469   return (*error == 0) ? res : (res < 0 ? -DBL_MAX : DBL_MAX);
00470 }
00471 
00472 
00473 double my_atof(const char *nptr)
00474 {
00475   int error;
00476   const char *end= nptr+65535;                  /* Should be enough */
00477   return (my_strtod(nptr, (char**) &end, &error));
00478 }
00479 
00480 
00481 /****************************************************************
00482  *
00483  * The author of this software is David M. Gay.
00484  *
00485  * Copyright (C) 1991, 2000, 2001 by Lucent Technologies.
00486  *
00487  * Permission to use, copy, modify, and distribute this software for any
00488  * purpose without fee is hereby granted, provided that this entire notice
00489  * is included in all copies of any software which is or includes a copy
00490  * or modification of this software and in all copies of the supporting
00491  * documentation for such software.
00492  *
00493  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
00494  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
00495  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
00496  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
00497  *
00498  ***************************************************************/
00499 /* Please send bug reports to David M. Gay (dmg at acm dot org,
00500  * with " at " changed at "@" and " dot " changed to ".").      */
00501 
00502 /*
00503   Original copy of the software is located at http://www.netlib.org/fp/dtoa.c
00504   It was adjusted to serve MySQL server needs:
00505   * strtod() was modified to not expect a zero-terminated string.
00506     It now honors 'se' (end of string) argument as the input parameter,
00507     not just as the output one.
00508   * in dtoa(), in case of overflow/underflow/NaN result string now contains "0";
00509     decpt is set to DTOA_OVERFLOW to indicate overflow.
00510   * support for VAX, IBM mainframe and 16-bit hardware removed
00511   * we always assume that 64-bit integer type is available
00512   * support for Kernigan-Ritchie style headers (pre-ANSI compilers)
00513     removed
00514   * all gcc warnings ironed out
00515   * we always assume multithreaded environment, so we had to change
00516     memory allocation procedures to use stack in most cases;
00517     malloc is used as the last resort.
00518   * pow5mult rewritten to use pre-calculated pow5 list instead of
00519     the one generated on the fly.
00520 */
00521 
00522 
00523 /*
00524   On a machine with IEEE extended-precision registers, it is
00525   necessary to specify double-precision (53-bit) rounding precision
00526   before invoking strtod or dtoa.  If the machine uses (the equivalent
00527   of) Intel 80x87 arithmetic, the call
00528        _control87(PC_53, MCW_PC);
00529   does this with many compilers.  Whether this or another call is
00530   appropriate depends on the compiler; for this to work, it may be
00531   necessary to #include "float.h" or another system-dependent header
00532   file.
00533 */
00534 
00535 typedef int32_t Long;
00536 typedef uint32_t ULong;
00537 typedef int64_t LLong;
00538 typedef uint64_t ULLong;
00539 
00540 typedef union { double d; ULong L[2]; } U;
00541 
00542 #if defined(WORDS_BIGENDIAN) || (defined(__FLOAT_WORD_ORDER) &&        \
00543                                  (__FLOAT_WORD_ORDER == __BIG_ENDIAN))
00544 #define word0(x) ((U*)&x)->L[0]
00545 #define word1(x) ((U*)&x)->L[1]
00546 #else
00547 #define word0(x) ((U*)&x)->L[1]
00548 #define word1(x) ((U*)&x)->L[0]
00549 #endif
00550 
00551 #define dval(x) ((U*)&x)->d
00552 
00553 /* #define P DBL_MANT_DIG */
00554 /* Ten_pmax= floor(P*log(2)/log(5)) */
00555 /* Bletch= (highest power of 2 < DBL_MAX_10_EXP) / 16 */
00556 /* Quick_max= floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
00557 /* Int_max= floor(P*log(FLT_RADIX)/log(10) - 1) */
00558 
00559 #define Exp_shift  20
00560 #define Exp_shift1 20
00561 #define Exp_msk1    0x100000
00562 #define Exp_mask  0x7ff00000
00563 #define P 53
00564 #define Bias 1023
00565 #define Emin (-1022)
00566 #define Exp_1  0x3ff00000
00567 #define Exp_11 0x3ff00000
00568 #define Ebits 11
00569 #define Frac_mask  0xfffff
00570 #define Frac_mask1 0xfffff
00571 #define Ten_pmax 22
00572 #define Bletch 0x10
00573 #define Bndry_mask  0xfffff
00574 #define Bndry_mask1 0xfffff
00575 #define LSB 1
00576 #define Sign_bit 0x80000000
00577 #define Log2P 1
00578 #define Tiny1 1
00579 #define Quick_max 14
00580 #define Int_max 14
00581 
00582 #ifndef Flt_Rounds
00583 #ifdef FLT_ROUNDS
00584 #define Flt_Rounds FLT_ROUNDS
00585 #else
00586 #define Flt_Rounds 1
00587 #endif
00588 #endif /*Flt_Rounds*/
00589 
00590 #define rounded_product(a,b) a*= b
00591 #define rounded_quotient(a,b) a/= b
00592 
00593 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
00594 #define Big1 0xffffffff
00595 #define FFFFFFFF 0xffffffffUL
00596 
00597 /* Arbitrary-length integer */
00598 
00599 typedef struct Bigint
00600 {
00601   union {
00602     ULong *x;              /* points right after this Bigint object */
00603     struct Bigint *next;   /* to maintain free lists */
00604   } p;
00605   int k;                   /* 2^k = maxwds */
00606   int maxwds;              /* maximum length in 32-bit words */
00607   int sign;                /* not zero if number is negative */
00608   int wds;                 /* current length in 32-bit words */
00609 } Bigint;
00610 
00611 static Bigint *Bcopy(Bigint* dst, Bigint* src)
00612 {
00613   dst->sign= src->sign;
00614   dst->wds= src->wds;
00615 
00616   assert(dst->maxwds >= src->wds);
00617 
00618   memcpy(dst->p.x, src->p.x, src->wds*sizeof(ULong));
00619 
00620   return dst;
00621 }
00622 
00623 static Bigint *Balloc(int k)
00624 {
00625   Bigint *rv;
00626 
00627   /* TODO: some malloc failure checking */
00628 
00629   int x= 1 << k;
00630   rv= (Bigint*) malloc(sizeof(Bigint));
00631 
00632   rv->p.x= (ULong*)malloc(x * sizeof(ULong));
00633 
00634   rv->k= k;
00635   rv->maxwds= x;
00636 
00637   rv->sign= rv->wds= 0;
00638 
00639   return rv;
00640 }
00641 
00642 
00643 /*
00644   If object was allocated on stack, try putting it to the free
00645   list. Otherwise call free().
00646 */
00647 
00648 static void Bfree(Bigint *v)
00649 {
00650   if(!v)
00651     return;
00652   free(v->p.x);
00653   free(v);
00654 }
00655 
00656 /* Bigint arithmetic functions */
00657 
00658 /* Multiply by m and add a */
00659 
00660 static Bigint *multadd(Bigint *b, int m, int a)
00661 {
00662   int i, wds;
00663   ULong *x;
00664   ULLong carry, y;
00665   Bigint *b1;
00666 
00667   wds= b->wds;
00668   x= b->p.x;
00669   i= 0;
00670   carry= a;
00671   do
00672   {
00673     y= *x * (ULLong)m + carry;
00674     carry= y >> 32;
00675     *x++= (ULong)(y & FFFFFFFF);
00676   }
00677   while (++i < wds);
00678   if (carry)
00679   {
00680     if (wds >= b->maxwds)
00681     {
00682       b1= Balloc(b->k+1);
00683       Bcopy(b1, b);
00684       Bfree(b);
00685       b= b1;
00686     }
00687     b->p.x[wds++]= (ULong) carry;
00688     b->wds= wds;
00689   }
00690   return b;
00691 }
00692 
00693 
00694 static Bigint *s2b(const char *s, int nd0, int nd, ULong y9)
00695 {
00696   Bigint *b;
00697   int i, k;
00698   Long x, y;
00699 
00700   x= (nd + 8) / 9;
00701   for (k= 0, y= 1; x > y; y <<= 1, k++) ;
00702   b= Balloc(k);
00703   b->p.x[0]= y9;
00704   b->wds= 1;
00705 
00706   i= 9;
00707   if (9 < nd0)
00708   {
00709     s+= 9;
00710     do
00711       b= multadd(b, 10, *s++ - '0');
00712     while (++i < nd0);
00713     s++;
00714   }
00715   else
00716     s+= 10;
00717   for(; i < nd; i++)
00718     b= multadd(b, 10, *s++ - '0');
00719   return b;
00720 }
00721 
00722 
00723 static int hi0bits(ULong x)
00724 {
00725   int k= 0;
00726 
00727   if (!(x & 0xffff0000))
00728   {
00729     k= 16;
00730     x<<= 16;
00731   }
00732   if (!(x & 0xff000000))
00733   {
00734     k+= 8;
00735     x<<= 8;
00736   }
00737   if (!(x & 0xf0000000))
00738   {
00739     k+= 4;
00740     x<<= 4;
00741   }
00742   if (!(x & 0xc0000000))
00743   {
00744     k+= 2;
00745     x<<= 2;
00746   }
00747   if (!(x & 0x80000000))
00748   {
00749     k++;
00750     if (!(x & 0x40000000))
00751       return 32;
00752   }
00753   return k;
00754 }
00755 
00756 
00757 static int lo0bits(ULong *y)
00758 {
00759   int k;
00760   ULong x= *y;
00761 
00762   if (x & 7)
00763   {
00764     if (x & 1)
00765       return 0;
00766     if (x & 2)
00767     {
00768       *y= x >> 1;
00769       return 1;
00770     }
00771     *y= x >> 2;
00772     return 2;
00773   }
00774   k= 0;
00775   if (!(x & 0xffff))
00776   {
00777     k= 16;
00778     x>>= 16;
00779   }
00780   if (!(x & 0xff))
00781   {
00782     k+= 8;
00783     x>>= 8;
00784   }
00785   if (!(x & 0xf))
00786   {
00787     k+= 4;
00788     x>>= 4;
00789   }
00790   if (!(x & 0x3))
00791   {
00792     k+= 2;
00793     x>>= 2;
00794   }
00795   if (!(x & 1))
00796   {
00797     k++;
00798     x>>= 1;
00799     if (!x)
00800       return 32;
00801   }
00802   *y= x;
00803   return k;
00804 }
00805 
00806 
00807 /* Convert integer to Bigint number */
00808 
00809 static Bigint *i2b(int i)
00810 {
00811   Bigint *b;
00812 
00813   b= Balloc(1);
00814   b->p.x[0]= i;
00815   b->wds= 1;
00816   return b;
00817 }
00818 
00819 
00820 /* Multiply two Bigint numbers */
00821 
00822 static Bigint *mult(Bigint *a, Bigint *b)
00823 {
00824   Bigint *c;
00825   int k, wa, wb, wc;
00826   ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
00827   ULong y;
00828   ULLong carry, z;
00829 
00830   if (a->wds < b->wds)
00831   {
00832     c= a;
00833     a= b;
00834     b= c;
00835   }
00836   k= a->k;
00837   wa= a->wds;
00838   wb= b->wds;
00839   wc= wa + wb;
00840   if (wc > a->maxwds)
00841     k++;
00842   c= Balloc(k);
00843   for (x= c->p.x, xa= x + wc; x < xa; x++)
00844     *x= 0;
00845   xa= a->p.x;
00846   xae= xa + wa;
00847   xb= b->p.x;
00848   xbe= xb + wb;
00849   xc0= c->p.x;
00850   for (; xb < xbe; xc0++)
00851   {
00852     if ((y= *xb++))
00853     {
00854       x= xa;
00855       xc= xc0;
00856       carry= 0;
00857       do
00858       {
00859         z= *x++ * (ULLong)y + *xc + carry;
00860         carry= z >> 32;
00861         *xc++= (ULong) (z & FFFFFFFF);
00862       }
00863       while (x < xae);
00864       *xc= (ULong) carry;
00865     }
00866   }
00867   for (xc0= c->p.x, xc= xc0 + wc; wc > 0 && !*--xc; --wc) ;
00868   c->wds= wc;
00869   return c;
00870 }
00871 
00872 
00873 /*
00874   Precalculated array of powers of 5: tested to be enough for
00875   vasting majority of dtoa_r cases.
00876 */
00877 
00878 static ULong powers5[]=
00879 {
00880   625UL,
00881 
00882   390625UL,
00883 
00884   2264035265UL, 35UL,
00885 
00886   2242703233UL, 762134875UL,  1262UL,
00887 
00888   3211403009UL, 1849224548UL, 3668416493UL, 3913284084UL, 1593091UL,
00889 
00890   781532673UL,  64985353UL,   253049085UL,  594863151UL,  3553621484UL,
00891   3288652808UL, 3167596762UL, 2788392729UL, 3911132675UL, 590UL,
00892 
00893   2553183233UL, 3201533787UL, 3638140786UL, 303378311UL, 1809731782UL,
00894   3477761648UL, 3583367183UL, 649228654UL, 2915460784UL, 487929380UL,
00895   1011012442UL, 1677677582UL, 3428152256UL, 1710878487UL, 1438394610UL,
00896   2161952759UL, 4100910556UL, 1608314830UL, 349175UL
00897 };
00898 
00899 
00900 static Bigint p5_a[]=
00901 {
00902   /*  { x } - k - maxwds - sign - wds */
00903   { { powers5 }, 1, 1, 0, 1 },
00904   { { powers5 + 1 }, 1, 1, 0, 1 },
00905   { { powers5 + 2 }, 1, 2, 0, 2 },
00906   { { powers5 + 4 }, 2, 3, 0, 3 },
00907   { { powers5 + 7 }, 3, 5, 0, 5 },
00908   { { powers5 + 12 }, 4, 10, 0, 10 },
00909   { { powers5 + 22 }, 5, 19, 0, 19 }
00910 };
00911 
00912 #define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1)
00913 
00914 static Bigint *pow5mult(Bigint *b, int k)
00915 {
00916   Bigint *b1, *p5, *p51;
00917   int i;
00918   static int p05[3]= { 5, 25, 125 };
00919 
00920   if ((i= k & 3))
00921     b= multadd(b, p05[i-1], 0);
00922 
00923   if (!(k>>= 2))
00924     return b;
00925   p5= p5_a;
00926   for (;;)
00927   {
00928     if (k & 1)
00929     {
00930       b1= mult(b, p5);
00931       Bfree(b);
00932       b= b1;
00933     }
00934     if (!(k>>= 1))
00935       break;
00936     /* Calculate next power of 5 */
00937     if (p5 < p5_a + P5A_MAX)
00938       ++p5;
00939     else if (p5 == p5_a + P5A_MAX)
00940       p5= mult(p5, p5);
00941     else
00942     {
00943       p51= mult(p5, p5);
00944       Bfree(p5);
00945       p5= p51;
00946     }
00947   }
00948   return b;
00949 }
00950 
00951 
00952 static Bigint *lshift(Bigint *b, int k)
00953 {
00954   int i, k1, n, n1;
00955   Bigint *b1;
00956   ULong *x, *x1, *xe, z;
00957 
00958   n= k >> 5;
00959   k1= b->k;
00960   n1= n + b->wds + 1;
00961   for (i= b->maxwds; n1 > i; i<<= 1)
00962     k1++;
00963   b1= Balloc(k1);
00964   x1= b1->p.x;
00965   for (i= 0; i < n; i++)
00966     *x1++= 0;
00967   x= b->p.x;
00968   xe= x + b->wds;
00969   if (k&= 0x1f)
00970   {
00971     k1= 32 - k;
00972     z= 0;
00973     do
00974     {
00975       *x1++= *x << k | z;
00976       z= *x++ >> k1;
00977     }
00978     while (x < xe);
00979     if ((*x1= z))
00980       ++n1;
00981   }
00982   else
00983     do
00984       *x1++= *x++;
00985     while (x < xe);
00986   b1->wds= n1 - 1;
00987   Bfree(b);
00988   return b1;
00989 }
00990 
00991 
00992 static int cmp(Bigint *a, Bigint *b)
00993 {
00994   ULong *xa, *xa0, *xb, *xb0;
00995   int i, j;
00996 
00997   i= a->wds;
00998   j= b->wds;
00999   if (i-= j)
01000     return i;
01001   xa0= a->p.x;
01002   xa= xa0 + j;
01003   xb0= b->p.x;
01004   xb= xb0 + j;
01005   for (;;)
01006   {
01007     if (*--xa != *--xb)
01008       return *xa < *xb ? -1 : 1;
01009     if (xa <= xa0)
01010       break;
01011   }
01012   return 0;
01013 }
01014 
01015 
01016 static Bigint *diff(Bigint *a, Bigint *b)
01017 {
01018   Bigint *c;
01019   int i, wa, wb;
01020   ULong *xa, *xae, *xb, *xbe, *xc;
01021   ULLong borrow, y;
01022 
01023   i= cmp(a,b);
01024   if (!i)
01025   {
01026     c= Balloc(0);
01027     c->wds= 1;
01028     c->p.x[0]= 0;
01029     return c;
01030   }
01031   if (i < 0)
01032   {
01033     c= a;
01034     a= b;
01035     b= c;
01036     i= 1;
01037   }
01038   else
01039     i= 0;
01040   c= Balloc(a->k);
01041   c->sign= i;
01042   wa= a->wds;
01043   xa= a->p.x;
01044   xae= xa + wa;
01045   wb= b->wds;
01046   xb= b->p.x;
01047   xbe= xb + wb;
01048   xc= c->p.x;
01049   borrow= 0;
01050   do
01051   {
01052     y= (ULLong)*xa++ - *xb++ - borrow;
01053     borrow= y >> 32 & (ULong)1;
01054     *xc++= (ULong) (y & FFFFFFFF);
01055   }
01056   while (xb < xbe);
01057   while (xa < xae)
01058   {
01059     y= *xa++ - borrow;
01060     borrow= y >> 32 & (ULong)1;
01061     *xc++= (ULong) (y & FFFFFFFF);
01062   }
01063   while (!*--xc)
01064     wa--;
01065   c->wds= wa;
01066   return c;
01067 }
01068 
01069 
01070 static double ulp(double x)
01071 {
01072   Long L;
01073   double a;
01074 
01075   L= (word0(x) & Exp_mask) - (P - 1)*Exp_msk1;
01076   word0(a) = L;
01077   word1(a) = 0;
01078   return dval(a);
01079 }
01080 
01081 
01082 static double b2d(Bigint *a, int *e)
01083 {
01084   ULong *xa, *xa0, w, y, z;
01085   int k;
01086   double d;
01087 #define d0 word0(d)
01088 #define d1 word1(d)
01089 
01090   xa0= a->p.x;
01091   xa= xa0 + a->wds;
01092   y= *--xa;
01093   k= hi0bits(y);
01094   *e= 32 - k;
01095   if (k < Ebits)
01096   {
01097     d0= Exp_1 | y >> (Ebits - k);
01098     w= xa > xa0 ? *--xa : 0;
01099     d1= y << ((32-Ebits) + k) | w >> (Ebits - k);
01100     goto ret_d;
01101   }
01102   z= xa > xa0 ? *--xa : 0;
01103   if (k-= Ebits)
01104   {
01105     d0= Exp_1 | y << k | z >> (32 - k);
01106     y= xa > xa0 ? *--xa : 0;
01107     d1= z << k | y >> (32 - k);
01108   }
01109   else
01110   {
01111     d0= Exp_1 | y;
01112     d1= z;
01113   }
01114  ret_d:
01115 #undef d0
01116 #undef d1
01117   return dval(d);
01118 }
01119 
01120 
01121 static Bigint *d2b(double d, int *e, int *bits)
01122 {
01123   Bigint *b;
01124   int de, k;
01125   ULong *x, y, z;
01126   int i;
01127 #define d0 word0(d)
01128 #define d1 word1(d)
01129 
01130   b= Balloc(1);
01131   x= b->p.x;
01132 
01133   z= d0 & Frac_mask;
01134   d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
01135   if ((de= (int)(d0 >> Exp_shift)))
01136     z|= Exp_msk1;
01137   if ((y= d1))
01138   {
01139     if ((k= lo0bits(&y)))
01140     {
01141       x[0]= y | z << (32 - k);
01142       z>>= k;
01143     }
01144     else
01145       x[0]= y;
01146     i= b->wds= (x[1]= z) ? 2 : 1;
01147   }
01148   else
01149   {
01150     k= lo0bits(&z);
01151     x[0]= z;
01152     i= b->wds= 1;
01153     k+= 32;
01154   }
01155   if (de)
01156   {
01157     *e= de - Bias - (P-1) + k;
01158     *bits= P - k;
01159   }
01160   else
01161   {
01162     *e= de - Bias - (P-1) + 1 + k;
01163     *bits= 32*i - hi0bits(x[i-1]);
01164   }
01165   return b;
01166 #undef d0
01167 #undef d1
01168 }
01169 
01170 
01171 static double ratio(Bigint *a, Bigint *b)
01172 {
01173   double da, db;
01174   int k, ka, kb;
01175 
01176   dval(da)= b2d(a, &ka);
01177   dval(db)= b2d(b, &kb);
01178   k= ka - kb + 32*(a->wds - b->wds);
01179   if (k > 0)
01180     word0(da)+= k*Exp_msk1;
01181   else
01182   {
01183     k= -k;
01184     word0(db)+= k*Exp_msk1;
01185   }
01186   return dval(da) / dval(db);
01187 }
01188 
01189 static const double tens[] =
01190 {
01191   1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
01192   1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
01193   1e20, 1e21, 1e22
01194 };
01195 
01196 static const double bigtens[]= { 1e16, 1e32, 1e64, 1e128, 1e256 };
01197 static const double tinytens[]=
01198 { 1e-16, 1e-32, 1e-64, 1e-128,
01199   9007199254740992.*9007199254740992.e-256 /* = 2^106 * 1e-53 */
01200 };
01201 /*
01202   The factor of 2^53 in tinytens[4] helps us avoid setting the underflow
01203   flag unnecessarily.  It leads to a song and dance at the end of strtod.
01204 */
01205 #define Scale_Bit 0x10
01206 #define n_bigtens 5
01207 
01208 /*
01209   strtod for IEEE--arithmetic machines.
01210 
01211   This strtod returns a nearest machine number to the input decimal
01212   string (or sets errno to EOVERFLOW). Ties are broken by the IEEE round-even
01213   rule.
01214 
01215   Inspired loosely by William D. Clinger's paper "How to Read Floating
01216   Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
01217 
01218   Modifications:
01219 
01220    1. We only require IEEE (not IEEE double-extended).
01221    2. We get by with floating-point arithmetic in a case that
01222      Clinger missed -- when we're computing d * 10^n
01223      for a small integer d and the integer n is not too
01224      much larger than 22 (the maximum integer k for which
01225      we can represent 10^k exactly), we may be able to
01226      compute (d*10^k) * 10^(e-k) with just one roundoff.
01227    3. Rather than a bit-at-a-time adjustment of the binary
01228      result in the hard case, we use floating-point
01229      arithmetic to determine the adjustment to within
01230      one bit; only in really hard cases do we need to
01231      compute a second residual.
01232    4. Because of 3., we don't need a large table of powers of 10
01233      for ten-to-e (just some small tables, e.g. of 10^k
01234      for 0 <= k <= 22).
01235 */
01236 
01237 static double my_strtod_int(const char *s00, char **se, int *error)
01238 {
01239   int scale;
01240   int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
01241      e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
01242   const char *s, *s0, *s1, *end = *se;
01243   double aadj, aadj1, adj, rv, rv0;
01244   Long L;
01245   ULong y, z;
01246   Bigint *bb= NULL, *bb1= NULL, *bd= NULL, *bd0= NULL, *bs= NULL, *delta= NULL;
01247 #ifdef SET_INEXACT
01248   int inexact, oldinexact;
01249 #endif
01250 
01251   c= 0;
01252   *error= 0;
01253 
01254   sign= nz0= nz= 0;
01255   dval(rv)= 0.;
01256   for (s= s00; s < end; s++)
01257     switch (*s) {
01258     case '-':
01259       sign= 1;
01260       /* no break */
01261     case '+':
01262       s++;
01263       goto break2;
01264     case '\t':
01265     case '\n':
01266     case '\v':
01267     case '\f':
01268     case '\r':
01269     case ' ':
01270       continue;
01271     default:
01272       goto break2;
01273     }
01274  break2:
01275   if (s >= end)
01276     goto ret0;
01277 
01278   if (*s == '0')
01279   {
01280     nz0= 1;
01281     while (++s < end && *s == '0') ;
01282     if (s >= end)
01283       goto ret;
01284   }
01285   s0= s;
01286   y= z= 0;
01287   for (nd= nf= 0; s < end && (c= *s) >= '0' && c <= '9'; nd++, s++)
01288     if (nd < 9)
01289       y= 10*y + c - '0';
01290     else if (nd < 16)
01291       z= 10*z + c - '0';
01292   nd0= nd;
01293   if (s < end - 1 && c == '.')
01294   {
01295     c= *++s;
01296     if (!nd)
01297     {
01298       for (; s < end && c == '0'; c= *++s)
01299         nz++;
01300       if (s < end && c > '0' && c <= '9')
01301       {
01302         s0= s;
01303         nf+= nz;
01304         nz= 0;
01305         goto have_dig;
01306       }
01307       goto dig_done;
01308     }
01309     for (; s < end && c >= '0' && c <= '9'; c = *++s)
01310     {
01311  have_dig:
01312       nz++;
01313       if (c-= '0')
01314       {
01315         nf+= nz;
01316         for (i= 1; i < nz; i++)
01317           if (nd++ < 9)
01318             y*= 10;
01319           else if (nd <= DBL_DIG + 1)
01320             z*= 10;
01321         if (nd++ < 9)
01322           y= 10*y + c;
01323         else if (nd <= DBL_DIG + 1)
01324           z= 10*z + c;
01325         nz= 0;
01326       }
01327     }
01328   }
01329  dig_done:
01330   e= 0;
01331   if (s < end && (c == 'e' || c == 'E'))
01332   {
01333     if (!nd && !nz && !nz0)
01334       goto ret0;
01335     s00= s;
01336     esign= 0;
01337     if (++s < end)
01338       switch (c= *s) {
01339       case '-':
01340         esign= 1;
01341       case '+':
01342         c= *++s;
01343       }
01344     if (s < end && c >= '0' && c <= '9')
01345     {
01346       while (s < end && c == '0')
01347         c= *++s;
01348       if (s < end && c > '0' && c <= '9') {
01349         L= c - '0';
01350         s1= s;
01351         while (++s < end && (c= *s) >= '0' && c <= '9')
01352           L= 10*L + c - '0';
01353         if (s - s1 > 8 || L > 19999)
01354           /* Avoid confusion from exponents
01355            * so large that e might overflow.
01356            */
01357           e= 19999; /* safe for 16 bit ints */
01358         else
01359           e= (int)L;
01360         if (esign)
01361           e= -e;
01362       }
01363       else
01364         e= 0;
01365     }
01366     else
01367       s= s00;
01368   }
01369   if (!nd)
01370   {
01371     if (!nz && !nz0)
01372     {
01373  ret0:
01374       s= s00;
01375       sign= 0;
01376     }
01377     goto ret;
01378   }
01379   e1= e -= nf;
01380 
01381   /*
01382     Now we have nd0 digits, starting at s0, followed by a
01383     decimal point, followed by nd-nd0 digits.  The number we're
01384     after is the integer represented by those digits times
01385     10**e
01386    */
01387 
01388   if (!nd0)
01389     nd0= nd;
01390   k= nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
01391   dval(rv)= y;
01392   if (k > 9)
01393   {
01394 #ifdef SET_INEXACT
01395     if (k > DBL_DIG)
01396       oldinexact = get_inexact();
01397 #endif
01398     dval(rv)= tens[k - 9] * dval(rv) + z;
01399   }
01400   bd0= 0;
01401   if (nd <= DBL_DIG
01402       )
01403   {
01404     if (!e)
01405       goto ret;
01406     if (e > 0)
01407     {
01408       if (e <= Ten_pmax)
01409       {
01410         /* rv = */ rounded_product(dval(rv), tens[e]);
01411         goto ret;
01412       }
01413       i= DBL_DIG - nd;
01414       if (e <= Ten_pmax + i)
01415       {
01416         /*
01417           A fancier test would sometimes let us do
01418           this for larger i values.
01419         */
01420         e-= i;
01421         dval(rv)*= tens[i];
01422         /* rv = */ rounded_product(dval(rv), tens[e]);
01423         goto ret;
01424       }
01425     }
01426 #ifndef Inaccurate_Divide
01427     else if (e >= -Ten_pmax)
01428     {
01429       /* rv = */ rounded_quotient(dval(rv), tens[-e]);
01430       goto ret;
01431     }
01432 #endif
01433   }
01434   e1+= nd - k;
01435 
01436 #ifdef SET_INEXACT
01437   inexact= 1;
01438   if (k <= DBL_DIG)
01439     oldinexact= get_inexact();
01440 #endif
01441   scale= 0;
01442 
01443   /* Get starting approximation = rv * 10**e1 */
01444 
01445   if (e1 > 0)
01446   {
01447     if ((i= e1 & 15))
01448       dval(rv)*= tens[i];
01449     if (e1&= ~15)
01450     {
01451       if (e1 > DBL_MAX_10_EXP)
01452       {
01453  ovfl:
01454         *error= EOVERFLOW;
01455         /* Can't trust HUGE_VAL */
01456         word0(rv)= Exp_mask;
01457         word1(rv)= 0;
01458 #ifdef SET_INEXACT
01459         /* set overflow bit */
01460         dval(rv0)= 1e300;
01461         dval(rv0)*= dval(rv0);
01462 #endif
01463         if (bd0)
01464           goto retfree;
01465         goto ret;
01466       }
01467       e1>>= 4;
01468       for(j= 0; e1 > 1; j++, e1>>= 1)
01469         if (e1 & 1)
01470           dval(rv)*= bigtens[j];
01471     /* The last multiplication could overflow. */
01472       word0(rv)-= P*Exp_msk1;
01473       dval(rv)*= bigtens[j];
01474       if ((z= word0(rv) & Exp_mask) > Exp_msk1 * (DBL_MAX_EXP + Bias - P))
01475         goto ovfl;
01476       if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P))
01477       {
01478         /* set to largest number (Can't trust DBL_MAX) */
01479         word0(rv)= Big0;
01480         word1(rv)= Big1;
01481       }
01482       else
01483         word0(rv)+= P*Exp_msk1;
01484     }
01485   }
01486   else if (e1 < 0)
01487   {
01488     e1= -e1;
01489     if ((i= e1 & 15))
01490       dval(rv)/= tens[i];
01491     if ((e1>>= 4))
01492     {
01493       if (e1 >= 1 << n_bigtens)
01494         goto undfl;
01495       if (e1 & Scale_Bit)
01496         scale= 2 * P;
01497       for(j= 0; e1 > 0; j++, e1>>= 1)
01498         if (e1 & 1)
01499           dval(rv)*= tinytens[j];
01500       if (scale && (j = 2 * P + 1 - ((word0(rv) & Exp_mask) >> Exp_shift)) > 0)
01501       {
01502         /* scaled rv is denormal; zap j low bits */
01503         if (j >= 32)
01504         {
01505           word1(rv)= 0;
01506           if (j >= 53)
01507             word0(rv)= (P + 2) * Exp_msk1;
01508           else
01509             word0(rv)&= 0xffffffff << (j - 32);
01510         }
01511         else
01512           word1(rv)&= 0xffffffff << j;
01513       }
01514       if (!dval(rv))
01515       {
01516  undfl:
01517           dval(rv)= 0.;
01518           if (bd0)
01519             goto retfree;
01520           goto ret;
01521       }
01522     }
01523   }
01524 
01525   /* Now the hard part -- adjusting rv to the correct value.*/
01526 
01527   /* Put digits into bd: true value = bd * 10^e */
01528 
01529   bd0= s2b(s0, nd0, nd, y);
01530 
01531   for(;;)
01532   {
01533     bd= Balloc(bd0->k);
01534     Bcopy(bd, bd0);
01535     bb= d2b(dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
01536     bs= i2b(1);
01537 
01538     if (e >= 0)
01539     {
01540       bb2= bb5= 0;
01541       bd2= bd5= e;
01542     }
01543     else
01544     {
01545       bb2= bb5= -e;
01546       bd2= bd5= 0;
01547     }
01548     if (bbe >= 0)
01549       bb2+= bbe;
01550     else
01551       bd2-= bbe;
01552     bs2= bb2;
01553     j= bbe - scale;
01554     i= j + bbbits - 1;  /* logb(rv) */
01555     if (i < Emin)  /* denormal */
01556       j+= P - Emin;
01557     else
01558       j= P + 1 - bbbits;
01559     bb2+= j;
01560     bd2+= j;
01561     bd2+= scale;
01562     i= bb2 < bd2 ? bb2 : bd2;
01563     if (i > bs2)
01564       i= bs2;
01565     if (i > 0)
01566     {
01567       bb2-= i;
01568       bd2-= i;
01569       bs2-= i;
01570     }
01571     if (bb5 > 0)
01572     {
01573       bs= pow5mult(bs, bb5);
01574       bb1= mult(bs, bb);
01575       Bfree(bb);
01576       bb= bb1;
01577     }
01578     if (bb2 > 0)
01579       bb= lshift(bb, bb2);
01580     if (bd5 > 0)
01581       bd= pow5mult(bd, bd5);
01582     if (bd2 > 0)
01583       bd= lshift(bd, bd2);
01584     if (bs2 > 0)
01585       bs= lshift(bs, bs2);
01586     delta= diff(bb, bd);
01587     dsign= delta->sign;
01588     delta->sign= 0;
01589     i= cmp(delta, bs);
01590 
01591     if (i < 0)
01592     {
01593       /*
01594         Error is less than half an ulp -- check for special case of mantissa
01595         a power of two.
01596       */
01597       if (dsign || word1(rv) || word0(rv) & Bndry_mask ||
01598           (word0(rv) & Exp_mask) <= (2 * P + 1) * Exp_msk1)
01599       {
01600 #ifdef SET_INEXACT
01601         if (!delta->x[0] && delta->wds <= 1)
01602           inexact= 0;
01603 #endif
01604         break;
01605       }
01606       if (!delta->p.x[0] && delta->wds <= 1)
01607       {
01608         /* exact result */
01609 #ifdef SET_INEXACT
01610         inexact= 0;
01611 #endif
01612         break;
01613       }
01614       delta= lshift(delta, Log2P);
01615       if (cmp(delta, bs) > 0)
01616         goto drop_down;
01617       break;
01618     }
01619     if (i == 0)
01620     {
01621       /* exactly half-way between */
01622       if (dsign)
01623       {
01624         if ((word0(rv) & Bndry_mask1) == Bndry_mask1 &&
01625             word1(rv) ==
01626             ((scale && (y = word0(rv) & Exp_mask) <= 2 * P * Exp_msk1) ?
01627              (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
01628              0xffffffff))
01629         {
01630           /*boundary case -- increment exponent*/
01631           word0(rv)= (word0(rv) & Exp_mask) + Exp_msk1;
01632           word1(rv) = 0;
01633           dsign = 0;
01634           break;
01635         }
01636       }
01637       else if (!(word0(rv) & Bndry_mask) && !word1(rv))
01638       {
01639  drop_down:
01640         /* boundary case -- decrement exponent */
01641         if (scale)
01642         {
01643           L= word0(rv) & Exp_mask;
01644           if (L <= (2 *P + 1) * Exp_msk1)
01645           {
01646             if (L > (P + 2) * Exp_msk1)
01647               /* round even ==> accept rv */
01648               break;
01649             /* rv = smallest denormal */
01650             goto undfl;
01651           }
01652         }
01653         L= (word0(rv) & Exp_mask) - Exp_msk1;
01654         word0(rv)= L | Bndry_mask1;
01655         word1(rv)= 0xffffffff;
01656         break;
01657       }
01658       if (!(word1(rv) & LSB))
01659         break;
01660       if (dsign)
01661         dval(rv)+= ulp(dval(rv));
01662       else
01663       {
01664         dval(rv)-= ulp(dval(rv));
01665         if (!dval(rv))
01666           goto undfl;
01667       }
01668       dsign= 1 - dsign;
01669       break;
01670     }
01671     if ((aadj= ratio(delta, bs)) <= 2.)
01672     {
01673       if (dsign)
01674         aadj= aadj1= 1.;
01675       else if (word1(rv) || word0(rv) & Bndry_mask)
01676       {
01677         if (word1(rv) == Tiny1 && !word0(rv))
01678           goto undfl;
01679         aadj= 1.;
01680         aadj1= -1.;
01681       }
01682       else
01683       {
01684         /* special case -- power of FLT_RADIX to be rounded down... */
01685         if (aadj < 2. / FLT_RADIX)
01686           aadj= 1. / FLT_RADIX;
01687         else
01688           aadj*= 0.5;
01689         aadj1= -aadj;
01690       }
01691     }
01692     else
01693     {
01694       aadj*= 0.5;
01695       aadj1= dsign ? aadj : -aadj;
01696       if (Flt_Rounds == 0)
01697         aadj1+= 0.5;
01698     }
01699     y= word0(rv) & Exp_mask;
01700 
01701     /* Check for overflow */
01702 
01703     if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1))
01704     {
01705       dval(rv0)= dval(rv);
01706       word0(rv)-= P * Exp_msk1;
01707       adj= aadj1 * ulp(dval(rv));
01708       dval(rv)+= adj;
01709       if ((word0(rv) & Exp_mask) >= Exp_msk1 * (DBL_MAX_EXP + Bias - P))
01710       {
01711         if (word0(rv0) == Big0 && word1(rv0) == Big1)
01712           goto ovfl;
01713         word0(rv)= Big0;
01714         word1(rv)= Big1;
01715         goto cont;
01716       }
01717       else
01718         word0(rv)+= P * Exp_msk1;
01719     }
01720     else
01721     {
01722       if (scale && y <= 2 * P * Exp_msk1)
01723       {
01724         if (aadj <= 0x7fffffff)
01725         {
01726           if ((z= (ULong) aadj) <= 0)
01727             z= 1;
01728           aadj= z;
01729           aadj1= dsign ? aadj : -aadj;
01730         }
01731         word0(aadj1)+= (2 * P + 1) * Exp_msk1 - y;
01732       }
01733       adj = aadj1 * ulp(dval(rv));
01734       dval(rv) += adj;
01735     }
01736     z= word0(rv) & Exp_mask;
01737 #ifndef SET_INEXACT
01738     if (!scale)
01739       if (y == z)
01740       {
01741         /* Can we stop now? */
01742         L= (Long)aadj;
01743         aadj-= L;
01744         /* The tolerances below are conservative. */
01745         if (dsign || word1(rv) || word0(rv) & Bndry_mask)
01746         {
01747           if (aadj < .4999999 || aadj > .5000001)
01748             break;
01749         }
01750         else if (aadj < .4999999 / FLT_RADIX)
01751           break;
01752       }
01753 #endif
01754  cont:
01755     Bfree(bb);
01756     Bfree(bd);
01757     Bfree(bs);
01758     Bfree(delta);
01759   }
01760 #ifdef SET_INEXACT
01761   if (inexact)
01762   {
01763     if (!oldinexact)
01764     {
01765       word0(rv0)= Exp_1 + (70 << Exp_shift);
01766       word1(rv0)= 0;
01767       dval(rv0)+= 1.;
01768     }
01769   }
01770   else if (!oldinexact)
01771     clear_inexact();
01772 #endif
01773   if (scale)
01774   {
01775     word0(rv0)= Exp_1 - 2 * P * Exp_msk1;
01776     word1(rv0)= 0;
01777     dval(rv)*= dval(rv0);
01778   }
01779 #ifdef SET_INEXACT
01780   if (inexact && !(word0(rv) & Exp_mask))
01781   {
01782     /* set underflow bit */
01783     dval(rv0)= 1e-300;
01784     dval(rv0)*= dval(rv0);
01785   }
01786 #endif
01787  retfree:
01788   Bfree(bb);
01789   Bfree(bd);
01790   Bfree(bs);
01791   Bfree(bd0);
01792   Bfree(delta);
01793  ret:
01794   *se= (char *)s;
01795   return sign ? -dval(rv) : dval(rv);
01796 }
01797 
01798 
01799 static int quorem(Bigint *b, Bigint *S)
01800 {
01801   int n;
01802   ULong *bx, *bxe, q, *sx, *sxe;
01803   ULLong borrow, carry, y, ys;
01804 
01805   n= S->wds;
01806   if (b->wds < n)
01807     return 0;
01808   sx= S->p.x;
01809   sxe= sx + --n;
01810   bx= b->p.x;
01811   bxe= bx + n;
01812   q= *bxe / (*sxe + 1);  /* ensure q <= true quotient */
01813   if (q)
01814   {
01815     borrow= 0;
01816     carry= 0;
01817     do
01818     {
01819       ys= *sx++ * (ULLong)q + carry;
01820       carry= ys >> 32;
01821       y= *bx - (ys & FFFFFFFF) - borrow;
01822       borrow= y >> 32 & (ULong)1;
01823       *bx++= (ULong) (y & FFFFFFFF);
01824     }
01825     while (sx <= sxe);
01826     if (!*bxe)
01827     {
01828       bx= b->p.x;
01829       while (--bxe > bx && !*bxe)
01830         --n;
01831       b->wds= n;
01832     }
01833   }
01834   if (cmp(b, S) >= 0)
01835   {
01836     q++;
01837     borrow= 0;
01838     carry= 0;
01839     bx= b->p.x;
01840     sx= S->p.x;
01841     do
01842     {
01843       ys= *sx++ + carry;
01844       carry= ys >> 32;
01845       y= *bx - (ys & FFFFFFFF) - borrow;
01846       borrow= y >> 32 & (ULong)1;
01847       *bx++= (ULong) (y & FFFFFFFF);
01848     }
01849     while (sx <= sxe);
01850     bx= b->p.x;
01851     bxe= bx + n;
01852     if (!*bxe)
01853     {
01854       while (--bxe > bx && !*bxe)
01855         --n;
01856       b->wds= n;
01857     }
01858   }
01859   return q;
01860 }
01861 
01862 
01863 /*
01864    dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
01865 
01866    Inspired by "How to Print Floating-Point Numbers Accurately" by
01867    Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
01868 
01869    Modifications:
01870         1. Rather than iterating, we use a simple numeric overestimate
01871            to determine k= floor(log10(d)).  We scale relevant
01872            quantities using O(log2(k)) rather than O(k) multiplications.
01873         2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
01874            try to generate digits strictly left to right.  Instead, we
01875            compute with fewer bits and propagate the carry if necessary
01876            when rounding the final digit up.  This is often faster.
01877         3. Under the assumption that input will be rounded nearest,
01878            mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
01879            That is, we allow equality in stopping tests when the
01880            round-nearest rule will give the same floating-point value
01881            as would satisfaction of the stopping test with strict
01882            inequality.
01883         4. We remove common factors of powers of 2 from relevant
01884            quantities.
01885         5. When converting floating-point integers less than 1e16,
01886            we use floating-point arithmetic rather than resorting
01887            to multiple-precision integers.
01888         6. When asked to produce fewer than 15 digits, we first try
01889            to get by with floating-point arithmetic; we resort to
01890            multiple-precision integer arithmetic only if we cannot
01891            guarantee that the floating-point calculation has given
01892            the correctly rounded result.  For k requested digits and
01893            "uniformly" distributed input, the probability is
01894            something like 10^(k-15) that we must resort to the Long
01895            calculation.
01896  */
01897 
01898 static char *dtoa(double d, int mode, int ndigits, int *decpt, int *sign,
01899                   char **rve)
01900 {
01901   /*
01902     Arguments ndigits, decpt, sign are similar to those
01903     of ecvt and fcvt; trailing zeros are suppressed from
01904     the returned string.  If not null, *rve is set to point
01905     to the end of the return value.  If d is +-Infinity or NaN,
01906     then *decpt is set to DTOA_OVERFLOW.
01907 
01908     mode:
01909           0 ==> shortest string that yields d when read in
01910                 and rounded to nearest.
01911 
01912           1 ==> like 0, but with Steele & White stopping rule;
01913                 e.g. with IEEE P754 arithmetic , mode 0 gives
01914                 1e23 whereas mode 1 gives 9.999999999999999e22.
01915           2 ==> cmax(1,ndigits) significant digits.  This gives a
01916                 return value similar to that of ecvt, except
01917                 that trailing zeros are suppressed.
01918           3 ==> through ndigits past the decimal point.  This
01919                 gives a return value similar to that from fcvt,
01920                 except that trailing zeros are suppressed, and
01921                 ndigits can be negative.
01922           4,5 ==> similar to 2 and 3, respectively, but (in
01923                 round-nearest mode) with the tests of mode 0 to
01924                 possibly return a shorter string that rounds to d.
01925           6-9 ==> Debugging modes similar to mode - 4:  don't try
01926                 fast floating-point estimate (if applicable).
01927 
01928       Values of mode other than 0-9 are treated as mode 0.
01929 
01930     Sufficient space is allocated to the return value
01931     to hold the suppressed trailing zeros.
01932   */
01933 
01934   int bbits, b2, b5, be, dig, i, ieps, ilim=0, ilim0, ilim1= 0,
01935     j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
01936     spec_case, try_quick;
01937   Long L;
01938   int denorm;
01939   ULong x;
01940   Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
01941   double d2, ds, eps;
01942   char *s, *s0;
01943 
01944   if (word0(d) & Sign_bit)
01945   {
01946     /* set sign for everything, including 0's and NaNs */
01947     *sign= 1;
01948     word0(d) &= ~Sign_bit;  /* clear sign bit */
01949   }
01950   else
01951     *sign= 0;
01952 
01953   /* If infinity, set decpt to DTOA_OVERFLOW, if 0 set it to 1 */
01954   if ((((word0(d) & Exp_mask) == Exp_mask) && ((*decpt= DTOA_OVERFLOW) != 0)) ||
01955       (!dval(d) && ((*decpt= 1) != 0)))
01956   {
01957     /* Infinity, NaN, 0 */
01958     char *res= (char*) malloc(2);
01959     res[0]= '0';
01960     res[1]= '\0';
01961     if (rve)
01962       *rve= res + 1;
01963     return res;
01964   }
01965 
01966 
01967   b= d2b(dval(d), &be, &bbits);
01968   if ((i= (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
01969   {
01970     dval(d2)= dval(d);
01971     word0(d2) &= Frac_mask1;
01972     word0(d2) |= Exp_11;
01973 
01974     /*
01975       log(x)       ~=~ log(1.5) + (x-1.5)/1.5
01976       log10(x)      =  log(x) / log(10)
01977                    ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
01978       log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
01979 
01980       This suggests computing an approximation k to log10(d) by
01981 
01982       k= (i - Bias)*0.301029995663981
01983            + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
01984 
01985       We want k to be too large rather than too small.
01986       The error in the first-order Taylor series approximation
01987       is in our favor, so we just round up the constant enough
01988       to compensate for any error in the multiplication of
01989       (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
01990       and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
01991       adding 1e-13 to the constant term more than suffices.
01992       Hence we adjust the constant term to 0.1760912590558.
01993       (We could get a more accurate k by invoking log10,
01994        but this is probably not worthwhile.)
01995     */
01996 
01997     i-= Bias;
01998     denorm= 0;
01999   }
02000   else
02001   {
02002     /* d is denormalized */
02003 
02004     i= bbits + be + (Bias + (P-1) - 1);
02005     x= i > 32  ? word0(d) << (64 - i) | word1(d) >> (i - 32)
02006       : word1(d) << (32 - i);
02007     dval(d2)= x;
02008     word0(d2)-= 31*Exp_msk1; /* adjust exponent */
02009     i-= (Bias + (P-1) - 1) + 1;
02010     denorm= 1;
02011   }
02012   ds= (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
02013   k= (int)ds;
02014   if (ds < 0. && ds != k)
02015     k--;    /* want k= floor(ds) */
02016   k_check= 1;
02017   if (k >= 0 && k <= Ten_pmax)
02018   {
02019     if (dval(d) < tens[k])
02020       k--;
02021     k_check= 0;
02022   }
02023   j= bbits - i - 1;
02024   if (j >= 0)
02025   {
02026     b2= 0;
02027     s2= j;
02028   }
02029   else
02030   {
02031     b2= -j;
02032     s2= 0;
02033   }
02034   if (k >= 0)
02035   {
02036     b5= 0;
02037     s5= k;
02038     s2+= k;
02039   }
02040   else
02041   {
02042     b2-= k;
02043     b5= -k;
02044     s5= 0;
02045   }
02046   if (mode < 0 || mode > 9)
02047     mode= 0;
02048 
02049   try_quick= 1;
02050 
02051   if (mode > 5)
02052   {
02053     mode-= 4;
02054     try_quick= 0;
02055   }
02056   leftright= 1;
02057   switch (mode) {
02058   case 0:
02059   case 1:
02060     ilim= ilim1= -1;
02061     i= 18;
02062     ndigits= 0;
02063     break;
02064   case 2:
02065     leftright= 0;
02066     /* no break */
02067   case 4:
02068     if (ndigits <= 0)
02069       ndigits= 1;
02070     ilim= ilim1= i= ndigits;
02071     break;
02072   case 3:
02073     leftright= 0;
02074     /* no break */
02075   case 5:
02076     i= ndigits + k + 1;
02077     ilim= i;
02078     ilim1= i - 1;
02079     if (i <= 0)
02080       i= 1;
02081   }
02082   s= s0= (char*)malloc(i+1); /* +1 for trailing '\0' appended
02083         at end of function */
02084 
02085   if (ilim >= 0 && ilim <= Quick_max && try_quick)
02086   {
02087     /* Try to get by with floating-point arithmetic. */
02088     i= 0;
02089     dval(d2)= dval(d);
02090     k0= k;
02091     ilim0= ilim;
02092     ieps= 2; /* conservative */
02093     if (k > 0)
02094     {
02095       ds= tens[k&0xf];
02096       j= k >> 4;
02097       if (j & Bletch)
02098       {
02099         /* prevent overflows */
02100         j&= Bletch - 1;
02101         dval(d)/= bigtens[n_bigtens-1];
02102         ieps++;
02103       }
02104       for (; j; j>>= 1, i++)
02105       {
02106         if (j & 1)
02107         {
02108           ieps++;
02109           ds*= bigtens[i];
02110         }
02111       }
02112       dval(d)/= ds;
02113     }
02114     else if ((j1= -k))
02115     {
02116       dval(d)*= tens[j1 & 0xf];
02117       for (j= j1 >> 4; j; j>>= 1, i++)
02118       {
02119         if (j & 1)
02120         {
02121           ieps++;
02122           dval(d)*= bigtens[i];
02123         }
02124       }
02125     }
02126     if (k_check && dval(d) < 1. && ilim > 0)
02127     {
02128       if (ilim1 <= 0)
02129         goto fast_failed;
02130       ilim= ilim1;
02131       k--;
02132       dval(d)*= 10.;
02133       ieps++;
02134     }
02135     dval(eps)= ieps*dval(d) + 7.;
02136     word0(eps)-= (P-1)*Exp_msk1;
02137     if (ilim == 0)
02138     {
02139       S= mhi= 0;
02140       dval(d)-= 5.;
02141       if (dval(d) > dval(eps))
02142         goto one_digit;
02143       if (dval(d) < -dval(eps))
02144         goto no_digits;
02145       goto fast_failed;
02146     }
02147     if (leftright)
02148     {
02149       /* Use Steele & White method of only generating digits needed. */
02150       dval(eps)= 0.5/tens[ilim-1] - dval(eps);
02151       for (i= 0;;)
02152       {
02153         L= (Long) dval(d);
02154         dval(d)-= L;
02155         *s++= '0' + (int)L;
02156         if (dval(d) < dval(eps))
02157           goto ret1;
02158         if (1. - dval(d) < dval(eps))
02159           goto bump_up;
02160         if (++i >= ilim)
02161           break;
02162         dval(eps)*= 10.;
02163         dval(d)*= 10.;
02164       }
02165     }
02166     else
02167     {
02168       /* Generate ilim digits, then fix them up. */
02169       dval(eps)*= tens[ilim-1];
02170       for (i= 1;; i++, dval(d)*= 10.)
02171       {
02172         L= (Long)(dval(d));
02173         if (!(dval(d)-= L))
02174           ilim= i;
02175         *s++= '0' + (int)L;
02176         if (i == ilim)
02177         {
02178           if (dval(d) > 0.5 + dval(eps))
02179             goto bump_up;
02180           else if (dval(d) < 0.5 - dval(eps))
02181           {
02182             while (*--s == '0') {}
02183             s++;
02184             goto ret1;
02185           }
02186           break;
02187         }
02188       }
02189     }
02190   fast_failed:
02191     s= s0;
02192     dval(d)= dval(d2);
02193     k= k0;
02194     ilim= ilim0;
02195   }
02196 
02197   /* Do we have a "small" integer? */
02198 
02199   if (be >= 0 && k <= Int_max)
02200   {
02201     /* Yes. */
02202     ds= tens[k];
02203     if (ndigits < 0 && ilim <= 0)
02204     {
02205       S= mhi= 0;
02206       if (ilim < 0 || dval(d) <= 5*ds)
02207         goto no_digits;
02208       goto one_digit;
02209     }
02210     for (i= 1;; i++, dval(d)*= 10.)
02211     {
02212       L= (Long)(dval(d) / ds);
02213       dval(d)-= L*ds;
02214       *s++= '0' + (int)L;
02215       if (!dval(d))
02216       {
02217         break;
02218       }
02219       if (i == ilim)
02220       {
02221         dval(d)+= dval(d);
02222         if (dval(d) > ds || (dval(d) == ds && L & 1))
02223         {
02224 bump_up:
02225           while (*--s == '9')
02226             if (s == s0)
02227             {
02228               k++;
02229               *s= '0';
02230               break;
02231             }
02232           ++*s++;
02233         }
02234         break;
02235       }
02236     }
02237     goto ret1;
02238   }
02239 
02240   m2= b2;
02241   m5= b5;
02242   mhi= mlo= 0;
02243   if (leftright)
02244   {
02245     i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
02246     b2+= i;
02247     s2+= i;
02248     mhi= i2b(1);
02249   }
02250   if (m2 > 0 && s2 > 0)
02251   {
02252     i= m2 < s2 ? m2 : s2;
02253     b2-= i;
02254     m2-= i;
02255     s2-= i;
02256   }
02257   if (b5 > 0)
02258   {
02259     if (leftright)
02260     {
02261       if (m5 > 0)
02262       {
02263         mhi= pow5mult(mhi, m5);
02264         b1= mult(mhi, b);
02265         Bfree(b);
02266         b= b1;
02267       }
02268       if ((j= b5 - m5))
02269         b= pow5mult(b, j);
02270     }
02271     else
02272       b= pow5mult(b, b5);
02273   }
02274   S= i2b(1);
02275   if (s5 > 0)
02276     S= pow5mult(S, s5);
02277 
02278   /* Check for special case that d is a normalized power of 2. */
02279 
02280   spec_case= 0;
02281   if ((mode < 2 || leftright)
02282      )
02283   {
02284     if (!word1(d) && !(word0(d) & Bndry_mask) &&
02285         word0(d) & (Exp_mask & ~Exp_msk1)
02286        )
02287     {
02288       /* The special case */
02289       b2+= Log2P;
02290       s2+= Log2P;
02291       spec_case= 1;
02292     }
02293   }
02294 
02295   /*
02296     Arrange for convenient computation of quotients:
02297     shift left if necessary so divisor has 4 leading 0 bits.
02298 
02299     Perhaps we should just compute leading 28 bits of S once
02300     a nd for all and pass them and a shift to quorem, so it
02301     can do shifts and ors to compute the numerator for q.
02302   */
02303   if ((i= ((s5 ? 32 - hi0bits(S->p.x[S->wds-1]) : 1) + s2) & 0x1f))
02304     i= 32 - i;
02305   if (i > 4)
02306   {
02307     i-= 4;
02308     b2+= i;
02309     m2+= i;
02310     s2+= i;
02311   }
02312   else if (i < 4)
02313   {
02314     i+= 28;
02315     b2+= i;
02316     m2+= i;
02317     s2+= i;
02318   }
02319   if (b2 > 0)
02320     b= lshift(b, b2);
02321   if (s2 > 0)
02322     S= lshift(S, s2);
02323   if (k_check)
02324   {
02325     if (cmp(b,S) < 0)
02326     {
02327       k--;
02328       /* we botched the k estimate */
02329       b= multadd(b, 10, 0);
02330       if (leftright)
02331         mhi= multadd(mhi, 10, 0);
02332       ilim= ilim1;
02333     }
02334   }
02335   if (ilim <= 0 && (mode == 3 || mode == 5))
02336   {
02337     if (ilim < 0 || cmp(b,S= multadd(S,5,0)) <= 0)
02338     {
02339       /* no digits, fcvt style */
02340 no_digits:
02341       k= -1 - ndigits;
02342       goto ret;
02343     }
02344 one_digit:
02345     *s++= '1';
02346     k++;
02347     goto ret;
02348   }
02349   if (leftright)
02350   {
02351     if (m2 > 0)
02352       mhi= lshift(mhi, m2);
02353 
02354     /*
02355       Compute mlo -- check for special case that d is a normalized power of 2.
02356     */
02357 
02358     mlo= mhi;
02359     if (spec_case)
02360     {
02361       mhi= Balloc(mhi->k);
02362       Bcopy(mhi, mlo);
02363       mhi= lshift(mhi, Log2P);
02364     }
02365 
02366     for (i= 1;;i++)
02367     {
02368       dig= quorem(b,S) + '0';
02369       /* Do we yet have the shortest decimal string that will round to d? */
02370       j= cmp(b, mlo);
02371       delta= diff(S, mhi);
02372       j1= delta->sign ? 1 : cmp(b, delta);
02373       Bfree(delta);
02374       if (j1 == 0 && mode != 1 && !(word1(d) & 1)
02375          )
02376       {
02377         if (dig == '9')
02378           goto round_9_up;
02379         if (j > 0)
02380           dig++;
02381         *s++= dig;
02382         goto ret;
02383       }
02384       if (j < 0 || (j == 0 && mode != 1 && !(word1(d) & 1)))
02385       {
02386         if (!b->p.x[0] && b->wds <= 1)
02387         {
02388           goto accept_dig;
02389         }
02390         if (j1 > 0)
02391         {
02392           b= lshift(b, 1);
02393           j1= cmp(b, S);
02394           if ((j1 > 0 || (j1 == 0 && dig & 1))
02395               && dig++ == '9')
02396             goto round_9_up;
02397         }
02398 accept_dig:
02399         *s++= dig;
02400         goto ret;
02401       }
02402       if (j1 > 0)
02403       {
02404         if (dig == '9')
02405         { /* possible if i == 1 */
02406 round_9_up:
02407           *s++= '9';
02408           goto roundoff;
02409         }
02410         *s++= dig + 1;
02411         goto ret;
02412       }
02413       *s++= dig;
02414       if (i == ilim)
02415         break;
02416       b= multadd(b, 10, 0);
02417       if (mlo == mhi)
02418         mlo= mhi= multadd(mhi, 10, 0);
02419       else
02420       {
02421         mlo= multadd(mlo, 10, 0);
02422         mhi= multadd(mhi, 10, 0);
02423       }
02424     }
02425   }
02426   else
02427     for (i= 1;; i++)
02428     {
02429       *s++= dig= quorem(b,S) + '0';
02430       if (!b->p.x[0] && b->wds <= 1)
02431       {
02432         goto ret;
02433       }
02434       if (i >= ilim)
02435         break;
02436       b= multadd(b, 10, 0);
02437     }
02438 
02439   /* Round off last digit */
02440 
02441   b= lshift(b, 1);
02442   j= cmp(b, S);
02443   if (j > 0 || (j == 0 && dig & 1))
02444   {
02445 roundoff:
02446     while (*--s == '9')
02447       if (s == s0)
02448       {
02449         k++;
02450         *s++= '1';
02451         goto ret;
02452       }
02453     ++*s++;
02454   }
02455   else
02456   {
02457     while (*--s == '0') {}
02458     s++;
02459   }
02460 ret:
02461   Bfree(S);
02462   if (mhi)
02463   {
02464     if (mlo && mlo != mhi)
02465       Bfree(mlo);
02466     Bfree(mhi);
02467   }
02468 ret1:
02469   Bfree(b);
02470   *s= 0;
02471   *decpt= k + 1;
02472   if (rve)
02473     *rve= s;
02474   return s0;
02475 }
02476 
02477 } /* namespace internal */
02478 } /* namespace drizzled */