Drizzled Public API Documentation

dtoa.cc
1 /* Copyright (C) 2007 MySQL AB
2  This program is free software; you can redistribute it and/or modify
3  it under the terms of the GNU General Public License as published by
4  the Free Software Foundation; either version 2 of the License, or
5  (at your option) any later version.
6 
7  This program is distributed in the hope that it will be useful,
8  but WITHOUT ANY WARRANTY; without even the implied warranty of
9  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  GNU General Public License for more details.
11 
12  You should have received a copy of the GNU General Public License
13  along with this program; if not, write to the Free Software
14  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA */
15 
16 /****************************************************************
17 
18  This file incorporates work covered by the following copyright and
19  permission notice:
20 
21  The author of this software is David M. Gay.
22 
23  Copyright (C) 1991, 2000, 2001 by Lucent Technologies.
24 
25  Permission to use, copy, modify, and distribute this software for any
26  purpose without fee is hereby granted, provided that this entire notice
27  is included in all copies of any software which is or includes a copy
28  or modification of this software and in all copies of the supporting
29  documentation for such software.
30 
31  THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
32  WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
33  REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
34  OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
35 
36  ***************************************************************/
37 
38 #include <config.h>
39 
40 #include <drizzled/internal/m_string.h> /* for memcpy and NOT_FIXED_DEC */
41 
42 #include <float.h>
43 
44 #include <cstdlib>
45 #include <cerrno>
46 #include <algorithm>
47 
48 using namespace std;
49 
50 namespace drizzled
51 {
52 namespace internal
53 {
54 
55 /* Magic value returned by dtoa() to indicate overflow */
56 #define DTOA_OVERFLOW 9999
57 
58 static double my_strtod_int(const char *, char **, int *);
59 static char *dtoa(double, int, int, int *, int *, char **);
60 
92 size_t my_fcvt(double x, int precision, char *to, bool *error)
93 {
94  int decpt, sign, len, i;
95  char *res, *src, *end, *dst= to;
96  assert(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
97 
98  res= dtoa(x, 5, precision, &decpt, &sign, &end);
99 
100  if (decpt == DTOA_OVERFLOW)
101  {
102  free(res);
103  *to++= '0';
104  *to= '\0';
105  if (error != NULL)
106  *error= true;
107  return 1;
108  }
109 
110  src= res;
111  len= end - src;
112 
113  if (sign)
114  *dst++= '-';
115 
116  if (decpt <= 0)
117  {
118  *dst++= '0';
119  *dst++= '.';
120  for (i= decpt; i < 0; i++)
121  *dst++= '0';
122  }
123 
124  for (i= 1; i <= len; i++)
125  {
126  *dst++= *src++;
127  if (i == decpt && i < len)
128  *dst++= '.';
129  }
130  while (i++ <= decpt)
131  *dst++= '0';
132 
133  if (precision > 0)
134  {
135  if (len <= decpt)
136  *dst++= '.';
137 
138  for (i= precision - max(0, (len - decpt)); i > 0; i--)
139  *dst++= '0';
140  }
141 
142  *dst= '\0';
143  if (error != NULL)
144  *error= false;
145 
146  free(res);
147 
148  return dst - to;
149 }
150 
214 size_t my_gcvt(double x, my_gcvt_arg_type type, int width, char *to,
215  bool *error)
216 {
217  int decpt, sign, len, exp_len;
218  char *res, *src, *end, *dst= to, *dend= dst + width;
219  bool have_space, force_e_format;
220  assert(width > 0 && to != NULL);
221 
222  /* We want to remove '-' from equations early */
223  if (x < 0.)
224  width--;
225 
226  res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? min(width, DBL_DIG) :
227  min(width, FLT_DIG), &decpt, &sign, &end);
228 
229  if (decpt == DTOA_OVERFLOW)
230  {
231  free(res);
232  *to++= '0';
233  *to= '\0';
234  if (error != NULL)
235  *error= true;
236  return 1;
237  }
238 
239  if (error != NULL)
240  *error= false;
241 
242  src= res;
243  len= end - res;
244 
245  /*
246  Number of digits in the exponent from the 'e' conversion.
247  The sign of the exponent is taken into account separetely, we don't need
248  to count it here.
249  */
250  exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
251 
252  /*
253  Do we have enough space for all digits in the 'f' format?
254  Let 'len' be the number of significant digits returned by dtoa,
255  and F be the length of the resulting decimal representation.
256  Consider the following cases:
257  1. decpt <= 0, i.e. we have "0.NNN" => F = len - decpt + 2
258  2. 0 < decpt < len, i.e. we have "NNN.NNN" => F = len + 1
259  3. len <= decpt, i.e. we have "NNN00" => F = decpt
260  */
261  have_space= (decpt <= 0 ? len - decpt + 2 :
262  decpt > 0 && decpt < len ? len + 1 :
263  decpt) <= width;
264  /*
265  The following is true when no significant digits can be placed with the
266  specified field width using the 'f' format, and the 'e' format
267  will not be truncated.
268  */
269  force_e_format= (decpt <= 0 && width <= 2 - decpt && width >= 3 + exp_len);
270  /*
271  Assume that we don't have enough space to place all significant digits in
272  the 'f' format. We have to choose between the 'e' format and the 'f' one
273  to keep as many significant digits as possible.
274  Let E and F be the lengths of decimal representaion in the 'e' and 'f'
275  formats, respectively. We want to use the 'f' format if, and only if F <= E.
276  Consider the following cases:
277  1. decpt <= 0.
278  F = len - decpt + 2 (see above)
279  E = len + (len > 1) + 1 + 1 (decpt <= -99) + (decpt <= -9) + 1
280  ("N.NNe-MMM")
281  (F <= E) <=> (len == 1 && decpt >= -1) || (len > 1 && decpt >= -2)
282  We also need to ensure that if the 'f' format is chosen,
283  the field width allows us to place at least one significant digit
284  (i.e. width > 2 - decpt). If not, we prefer the 'e' format.
285  2. 0 < decpt < len
286  F = len + 1 (see above)
287  E = len + 1 + 1 + ... ("N.NNeMMM")
288  F is always less than E.
289  3. len <= decpt <= width
290  In this case we have enough space to represent the number in the 'f'
291  format, so we prefer it with some exceptions.
292  4. width < decpt
293  The number cannot be represented in the 'f' format at all, always use
294  the 'e' 'one.
295  */
296  if ((have_space ||
297  /*
298  Not enough space, let's see if the 'f' format provides the most number
299  of significant digits.
300  */
301  ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
302  (len > 1 || !force_e_format)))) &&
303  !force_e_format)) &&
304 
305  /*
306  Use the 'e' format in some cases even if we have enough space for the
307  'f' one. See comment for MAX_DECPT_FOR_F_FORMAT.
308  */
309  (!have_space || (decpt >= -MAX_DECPT_FOR_F_FORMAT + 1 &&
310  (decpt <= MAX_DECPT_FOR_F_FORMAT || len > decpt))))
311  {
312  /* 'f' format */
313  int i;
314 
315  width-= (decpt < len) + (decpt <= 0 ? 1 - decpt : 0);
316 
317  /* Do we have to truncate any digits? */
318  if (width < len)
319  {
320  if (width < decpt)
321  {
322  if (error != NULL)
323  *error= true;
324  width= decpt;
325  }
326 
327  /*
328  We want to truncate (len - width) least significant digits after the
329  decimal point. For this we are calling dtoa with mode=5, passing the
330  number of significant digits = (len-decpt) - (len-width) = width-decpt
331  */
332  free(res);
333  res= dtoa(x, 5, width - decpt, &decpt, &sign, &end);
334  src= res;
335  len= end - res;
336  }
337 
338  if (len == 0)
339  {
340  /* Underflow. Just print '0' and exit */
341  *dst++= '0';
342  goto end;
343  }
344 
345  /*
346  At this point we are sure we have enough space to put all digits
347  returned by dtoa
348  */
349  if (sign && dst < dend)
350  *dst++= '-';
351  if (decpt <= 0)
352  {
353  if (dst < dend)
354  *dst++= '0';
355  if (len > 0 && dst < dend)
356  *dst++= '.';
357  for (; decpt < 0 && dst < dend; decpt++)
358  *dst++= '0';
359  }
360 
361  for (i= 1; i <= len && dst < dend; i++)
362  {
363  *dst++= *src++;
364  if (i == decpt && i < len && dst < dend)
365  *dst++= '.';
366  }
367  while (i++ <= decpt && dst < dend)
368  *dst++= '0';
369  }
370  else
371  {
372  /* 'e' format */
373  int decpt_sign= 0;
374 
375  if (--decpt < 0)
376  {
377  decpt= -decpt;
378  width--;
379  decpt_sign= 1;
380  }
381  width-= 1 + exp_len; /* eNNN */
382 
383  if (len > 1)
384  width--;
385 
386  if (width <= 0)
387  {
388  /* Overflow */
389  if (error != NULL)
390  *error= true;
391  width= 0;
392  }
393 
394  /* Do we have to truncate any digits? */
395  if (width < len)
396  {
397  /* Yes, re-convert with a smaller width */
398  free(res);
399  res= dtoa(x, 4, width, &decpt, &sign, &end);
400  src= res;
401  len= end - res;
402  if (--decpt < 0)
403  decpt= -decpt;
404  }
405  /*
406  At this point we are sure we have enough space to put all digits
407  returned by dtoa
408  */
409  if (sign && dst < dend)
410  *dst++= '-';
411  if (dst < dend)
412  *dst++= *src++;
413  if (len > 1 && dst < dend)
414  {
415  *dst++= '.';
416  while (src < end && dst < dend)
417  *dst++= *src++;
418  }
419  if (dst < dend)
420  *dst++= 'e';
421  if (decpt_sign && dst < dend)
422  *dst++= '-';
423 
424  if (decpt >= 100 && dst < dend)
425  {
426  *dst++= decpt / 100 + '0';
427  decpt%= 100;
428  if (dst < dend)
429  *dst++= decpt / 10 + '0';
430  }
431  else if (decpt >= 10 && dst < dend)
432  *dst++= decpt / 10 + '0';
433  if (dst < dend)
434  *dst++= decpt % 10 + '0';
435 
436  }
437 
438 end:
439  free(res);
440  *dst= '\0';
441 
442  return dst - to;
443 }
444 
463 double my_strtod(const char *str, char **end, int *error)
464 {
465  double res;
466  assert(str != NULL && end != NULL && *end != NULL && error != NULL);
467 
468  res= my_strtod_int(str, end, error);
469  return (*error == 0) ? res : (res < 0 ? -DBL_MAX : DBL_MAX);
470 }
471 
472 
473 double my_atof(const char *nptr)
474 {
475  int error;
476  const char *end= nptr+65535; /* Should be enough */
477  return (my_strtod(nptr, (char**) &end, &error));
478 }
479 
480 
481 /****************************************************************
482  *
483  * The author of this software is David M. Gay.
484  *
485  * Copyright (C) 1991, 2000, 2001 by Lucent Technologies.
486  *
487  * Permission to use, copy, modify, and distribute this software for any
488  * purpose without fee is hereby granted, provided that this entire notice
489  * is included in all copies of any software which is or includes a copy
490  * or modification of this software and in all copies of the supporting
491  * documentation for such software.
492  *
493  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
494  * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
495  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
496  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
497  *
498  ***************************************************************/
499 /* Please send bug reports to David M. Gay (dmg at acm dot org,
500  * with " at " changed at "@" and " dot " changed to "."). */
501 
502 /*
503  Original copy of the software is located at http://www.netlib.org/fp/dtoa.c
504  It was adjusted to serve MySQL server needs:
505  * strtod() was modified to not expect a zero-terminated string.
506  It now honors 'se' (end of string) argument as the input parameter,
507  not just as the output one.
508  * in dtoa(), in case of overflow/underflow/NaN result string now contains "0";
509  decpt is set to DTOA_OVERFLOW to indicate overflow.
510  * support for VAX, IBM mainframe and 16-bit hardware removed
511  * we always assume that 64-bit integer type is available
512  * support for Kernigan-Ritchie style headers (pre-ANSI compilers)
513  removed
514  * all gcc warnings ironed out
515  * we always assume multithreaded environment, so we had to change
516  memory allocation procedures to use stack in most cases;
517  malloc is used as the last resort.
518  * pow5mult rewritten to use pre-calculated pow5 list instead of
519  the one generated on the fly.
520 */
521 
522 
523 /*
524  On a machine with IEEE extended-precision registers, it is
525  necessary to specify double-precision (53-bit) rounding precision
526  before invoking strtod or dtoa. If the machine uses (the equivalent
527  of) Intel 80x87 arithmetic, the call
528  _control87(PC_53, MCW_PC);
529  does this with many compilers. Whether this or another call is
530  appropriate depends on the compiler; for this to work, it may be
531  necessary to #include "float.h" or another system-dependent header
532  file.
533 */
534 
535 typedef int32_t Long;
536 typedef uint32_t ULong;
537 typedef int64_t LLong;
538 typedef uint64_t ULLong;
539 
540 typedef union { double d; ULong L[2]; } U;
541 
542 #if defined(WORDS_BIGENDIAN) || (defined(__FLOAT_WORD_ORDER) && \
543  (__FLOAT_WORD_ORDER == __BIG_ENDIAN))
544 #define word0(x) ((U*)&x)->L[0]
545 #define word1(x) ((U*)&x)->L[1]
546 #else
547 #define word0(x) ((U*)&x)->L[1]
548 #define word1(x) ((U*)&x)->L[0]
549 #endif
550 
551 #define dval(x) ((U*)&x)->d
552 
553 /* #define P DBL_MANT_DIG */
554 /* Ten_pmax= floor(P*log(2)/log(5)) */
555 /* Bletch= (highest power of 2 < DBL_MAX_10_EXP) / 16 */
556 /* Quick_max= floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
557 /* Int_max= floor(P*log(FLT_RADIX)/log(10) - 1) */
558 
559 #define Exp_shift 20
560 #define Exp_shift1 20
561 #define Exp_msk1 0x100000
562 #define Exp_mask 0x7ff00000
563 #define P 53
564 #define Bias 1023
565 #define Emin (-1022)
566 #define Exp_1 0x3ff00000
567 #define Exp_11 0x3ff00000
568 #define Ebits 11
569 #define Frac_mask 0xfffff
570 #define Frac_mask1 0xfffff
571 #define Ten_pmax 22
572 #define Bletch 0x10
573 #define Bndry_mask 0xfffff
574 #define Bndry_mask1 0xfffff
575 #define LSB 1
576 #define Sign_bit 0x80000000
577 #define Log2P 1
578 #define Tiny1 1
579 #define Quick_max 14
580 #define Int_max 14
581 
582 #ifndef Flt_Rounds
583 #ifdef FLT_ROUNDS
584 #define Flt_Rounds FLT_ROUNDS
585 #else
586 #define Flt_Rounds 1
587 #endif
588 #endif /*Flt_Rounds*/
589 
590 #define rounded_product(a,b) a*= b
591 #define rounded_quotient(a,b) a/= b
592 
593 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
594 #define Big1 0xffffffff
595 #define FFFFFFFF 0xffffffffUL
596 
597 /* Arbitrary-length integer */
598 
599 typedef struct Bigint
600 {
601  union {
602  ULong *x; /* points right after this Bigint object */
603  } p;
604  int k; /* 2^k = maxwds */
605  int maxwds; /* maximum length in 32-bit words */
606  int sign; /* not zero if number is negative */
607  int wds; /* current length in 32-bit words */
608 } Bigint;
609 
610 static Bigint *Bcopy(Bigint* dst, Bigint* src)
611 {
612  dst->sign= src->sign;
613  dst->wds= src->wds;
614 
615  assert(dst->maxwds >= src->wds);
616 
617  memcpy(dst->p.x, src->p.x, src->wds*sizeof(ULong));
618 
619  return dst;
620 }
621 
622 static Bigint *Balloc(int k)
623 {
624  int x= 1 << k;
625  Bigint* rv= (Bigint*) malloc(sizeof(Bigint));
626 
627  rv->p.x= (ULong*)malloc(x * sizeof(ULong));
628 
629  rv->k= k;
630  rv->maxwds= x;
631 
632  rv->sign= rv->wds= 0;
633 
634  return rv;
635 }
636 
637 
638 /*
639  If object was allocated on stack, try putting it to the free
640  list. Otherwise call free().
641 */
642 
643 static void Bfree(Bigint *v)
644 {
645  if(!v)
646  return;
647  free(v->p.x);
648  free(v);
649 }
650 
651 /* Bigint arithmetic functions */
652 
653 /* Multiply by m and add a */
654 
655 static Bigint *multadd(Bigint *b, int m, int a)
656 {
657  int i, wds;
658  ULong *x;
659  ULLong carry, y;
660  Bigint *b1;
661 
662  wds= b->wds;
663  x= b->p.x;
664  i= 0;
665  carry= a;
666  do
667  {
668  y= *x * (ULLong)m + carry;
669  carry= y >> 32;
670  *x++= (ULong)(y & FFFFFFFF);
671  }
672  while (++i < wds);
673  if (carry)
674  {
675  if (wds >= b->maxwds)
676  {
677  b1= Balloc(b->k+1);
678  Bcopy(b1, b);
679  Bfree(b);
680  b= b1;
681  }
682  b->p.x[wds++]= (ULong) carry;
683  b->wds= wds;
684  }
685  return b;
686 }
687 
688 
689 static Bigint *s2b(const char *s, int nd0, int nd, ULong y9)
690 {
691  Bigint *b;
692  int i, k;
693  Long x, y;
694 
695  x= (nd + 8) / 9;
696  for (k= 0, y= 1; x > y; y <<= 1, k++) ;
697  b= Balloc(k);
698  b->p.x[0]= y9;
699  b->wds= 1;
700 
701  i= 9;
702  if (9 < nd0)
703  {
704  s+= 9;
705  do
706  b= multadd(b, 10, *s++ - '0');
707  while (++i < nd0);
708  s++;
709  }
710  else
711  s+= 10;
712  for(; i < nd; i++)
713  b= multadd(b, 10, *s++ - '0');
714  return b;
715 }
716 
717 
718 static int hi0bits(ULong x)
719 {
720  int k= 0;
721 
722  if (!(x & 0xffff0000))
723  {
724  k= 16;
725  x<<= 16;
726  }
727  if (!(x & 0xff000000))
728  {
729  k+= 8;
730  x<<= 8;
731  }
732  if (!(x & 0xf0000000))
733  {
734  k+= 4;
735  x<<= 4;
736  }
737  if (!(x & 0xc0000000))
738  {
739  k+= 2;
740  x<<= 2;
741  }
742  if (!(x & 0x80000000))
743  {
744  k++;
745  if (!(x & 0x40000000))
746  return 32;
747  }
748  return k;
749 }
750 
751 
752 static int lo0bits(ULong *y)
753 {
754  int k;
755  ULong x= *y;
756 
757  if (x & 7)
758  {
759  if (x & 1)
760  return 0;
761  if (x & 2)
762  {
763  *y= x >> 1;
764  return 1;
765  }
766  *y= x >> 2;
767  return 2;
768  }
769  k= 0;
770  if (!(x & 0xffff))
771  {
772  k= 16;
773  x>>= 16;
774  }
775  if (!(x & 0xff))
776  {
777  k+= 8;
778  x>>= 8;
779  }
780  if (!(x & 0xf))
781  {
782  k+= 4;
783  x>>= 4;
784  }
785  if (!(x & 0x3))
786  {
787  k+= 2;
788  x>>= 2;
789  }
790  if (!(x & 1))
791  {
792  k++;
793  x>>= 1;
794  if (!x)
795  return 32;
796  }
797  *y= x;
798  return k;
799 }
800 
801 
802 /* Convert integer to Bigint number */
803 
804 static Bigint *i2b(int i)
805 {
806  Bigint *b;
807 
808  b= Balloc(1);
809  b->p.x[0]= i;
810  b->wds= 1;
811  return b;
812 }
813 
814 
815 /* Multiply two Bigint numbers */
816 
817 static Bigint *mult(Bigint *a, Bigint *b)
818 {
819  Bigint *c;
820  int k, wa, wb, wc;
821  ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
822  ULong y;
823  ULLong carry, z;
824 
825  if (a->wds < b->wds)
826  {
827  c= a;
828  a= b;
829  b= c;
830  }
831  k= a->k;
832  wa= a->wds;
833  wb= b->wds;
834  wc= wa + wb;
835  if (wc > a->maxwds)
836  k++;
837  c= Balloc(k);
838  for (x= c->p.x, xa= x + wc; x < xa; x++)
839  *x= 0;
840  xa= a->p.x;
841  xae= xa + wa;
842  xb= b->p.x;
843  xbe= xb + wb;
844  xc0= c->p.x;
845  for (; xb < xbe; xc0++)
846  {
847  if ((y= *xb++))
848  {
849  x= xa;
850  xc= xc0;
851  carry= 0;
852  do
853  {
854  z= *x++ * (ULLong)y + *xc + carry;
855  carry= z >> 32;
856  *xc++= (ULong) (z & FFFFFFFF);
857  }
858  while (x < xae);
859  *xc= (ULong) carry;
860  }
861  }
862  for (xc0= c->p.x, xc= xc0 + wc; wc > 0 && !*--xc; --wc) ;
863  c->wds= wc;
864  return c;
865 }
866 
867 
868 /*
869  Precalculated array of powers of 5: tested to be enough for
870  vasting majority of dtoa_r cases.
871 */
872 
873 static ULong powers5[]=
874 {
875  625UL,
876 
877  390625UL,
878 
879  2264035265UL, 35UL,
880 
881  2242703233UL, 762134875UL, 1262UL,
882 
883  3211403009UL, 1849224548UL, 3668416493UL, 3913284084UL, 1593091UL,
884 
885  781532673UL, 64985353UL, 253049085UL, 594863151UL, 3553621484UL,
886  3288652808UL, 3167596762UL, 2788392729UL, 3911132675UL, 590UL,
887 
888  2553183233UL, 3201533787UL, 3638140786UL, 303378311UL, 1809731782UL,
889  3477761648UL, 3583367183UL, 649228654UL, 2915460784UL, 487929380UL,
890  1011012442UL, 1677677582UL, 3428152256UL, 1710878487UL, 1438394610UL,
891  2161952759UL, 4100910556UL, 1608314830UL, 349175UL
892 };
893 
894 
895 static Bigint p5_a[]=
896 {
897  /* { x } - k - maxwds - sign - wds */
898  { { powers5 }, 1, 1, 0, 1 },
899  { { powers5 + 1 }, 1, 1, 0, 1 },
900  { { powers5 + 2 }, 1, 2, 0, 2 },
901  { { powers5 + 4 }, 2, 3, 0, 3 },
902  { { powers5 + 7 }, 3, 5, 0, 5 },
903  { { powers5 + 12 }, 4, 10, 0, 10 },
904  { { powers5 + 22 }, 5, 19, 0, 19 }
905 };
906 
907 #define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1)
908 
909 static Bigint *pow5mult(Bigint *b, int k)
910 {
911  Bigint *b1, *p5, *p51;
912  int i;
913  static int p05[3]= { 5, 25, 125 };
914 
915  if ((i= k & 3))
916  b= multadd(b, p05[i-1], 0);
917 
918  if (!(k>>= 2))
919  return b;
920  p5= p5_a;
921  for (;;)
922  {
923  if (k & 1)
924  {
925  b1= mult(b, p5);
926  Bfree(b);
927  b= b1;
928  }
929  if (!(k>>= 1))
930  break;
931  /* Calculate next power of 5 */
932  if (p5 < p5_a + P5A_MAX)
933  ++p5;
934  else if (p5 == p5_a + P5A_MAX)
935  p5= mult(p5, p5);
936  else
937  {
938  p51= mult(p5, p5);
939  Bfree(p5);
940  p5= p51;
941  }
942  }
943  return b;
944 }
945 
946 
947 static Bigint *lshift(Bigint *b, int k)
948 {
949  int i, k1, n, n1;
950  Bigint *b1;
951  ULong *x, *x1, *xe, z;
952 
953  n= k >> 5;
954  k1= b->k;
955  n1= n + b->wds + 1;
956  for (i= b->maxwds; n1 > i; i<<= 1)
957  k1++;
958  b1= Balloc(k1);
959  x1= b1->p.x;
960  for (i= 0; i < n; i++)
961  *x1++= 0;
962  x= b->p.x;
963  xe= x + b->wds;
964  if (k&= 0x1f)
965  {
966  k1= 32 - k;
967  z= 0;
968  do
969  {
970  *x1++= *x << k | z;
971  z= *x++ >> k1;
972  }
973  while (x < xe);
974  if ((*x1= z))
975  ++n1;
976  }
977  else
978  do
979  *x1++= *x++;
980  while (x < xe);
981  b1->wds= n1 - 1;
982  Bfree(b);
983  return b1;
984 }
985 
986 
987 static int cmp(Bigint *a, Bigint *b)
988 {
989  ULong *xa, *xa0, *xb, *xb0;
990  int i, j;
991 
992  i= a->wds;
993  j= b->wds;
994  if (i-= j)
995  return i;
996  xa0= a->p.x;
997  xa= xa0 + j;
998  xb0= b->p.x;
999  xb= xb0 + j;
1000  for (;;)
1001  {
1002  if (*--xa != *--xb)
1003  return *xa < *xb ? -1 : 1;
1004  if (xa <= xa0)
1005  break;
1006  }
1007  return 0;
1008 }
1009 
1010 
1011 static Bigint *diff(Bigint *a, Bigint *b)
1012 {
1013  Bigint *c;
1014  int i, wa, wb;
1015  ULong *xa, *xae, *xb, *xbe, *xc;
1016  ULLong borrow, y;
1017 
1018  i= cmp(a,b);
1019  if (!i)
1020  {
1021  c= Balloc(0);
1022  c->wds= 1;
1023  c->p.x[0]= 0;
1024  return c;
1025  }
1026  if (i < 0)
1027  {
1028  c= a;
1029  a= b;
1030  b= c;
1031  i= 1;
1032  }
1033  else
1034  i= 0;
1035  c= Balloc(a->k);
1036  c->sign= i;
1037  wa= a->wds;
1038  xa= a->p.x;
1039  xae= xa + wa;
1040  wb= b->wds;
1041  xb= b->p.x;
1042  xbe= xb + wb;
1043  xc= c->p.x;
1044  borrow= 0;
1045  do
1046  {
1047  y= (ULLong)*xa++ - *xb++ - borrow;
1048  borrow= y >> 32 & (ULong)1;
1049  *xc++= (ULong) (y & FFFFFFFF);
1050  }
1051  while (xb < xbe);
1052  while (xa < xae)
1053  {
1054  y= *xa++ - borrow;
1055  borrow= y >> 32 & (ULong)1;
1056  *xc++= (ULong) (y & FFFFFFFF);
1057  }
1058  while (!*--xc)
1059  wa--;
1060  c->wds= wa;
1061  return c;
1062 }
1063 
1064 
1065 static double ulp(double x)
1066 {
1067  Long L;
1068  double a;
1069 
1070  L= (word0(x) & Exp_mask) - (P - 1)*Exp_msk1;
1071  word0(a) = L;
1072  word1(a) = 0;
1073  return dval(a);
1074 }
1075 
1076 
1077 static double b2d(Bigint *a, int *e)
1078 {
1079  ULong *xa, *xa0, w, y, z;
1080  int k;
1081  double d;
1082 #define d0 word0(d)
1083 #define d1 word1(d)
1084 
1085  xa0= a->p.x;
1086  xa= xa0 + a->wds;
1087  y= *--xa;
1088  k= hi0bits(y);
1089  *e= 32 - k;
1090  if (k < Ebits)
1091  {
1092  d0= Exp_1 | y >> (Ebits - k);
1093  w= xa > xa0 ? *--xa : 0;
1094  d1= y << ((32-Ebits) + k) | w >> (Ebits - k);
1095  goto ret_d;
1096  }
1097  z= xa > xa0 ? *--xa : 0;
1098  if (k-= Ebits)
1099  {
1100  d0= Exp_1 | y << k | z >> (32 - k);
1101  y= xa > xa0 ? *--xa : 0;
1102  d1= z << k | y >> (32 - k);
1103  }
1104  else
1105  {
1106  d0= Exp_1 | y;
1107  d1= z;
1108  }
1109  ret_d:
1110 #undef d0
1111 #undef d1
1112  return dval(d);
1113 }
1114 
1115 
1116 static Bigint *d2b(double d, int *e, int *bits)
1117 {
1118  Bigint *b;
1119  int de, k;
1120  ULong *x, y, z;
1121  int i;
1122 #define d0 word0(d)
1123 #define d1 word1(d)
1124 
1125  b= Balloc(1);
1126  x= b->p.x;
1127 
1128  z= d0 & Frac_mask;
1129  d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1130  if ((de= (int)(d0 >> Exp_shift)))
1131  z|= Exp_msk1;
1132  if ((y= d1))
1133  {
1134  if ((k= lo0bits(&y)))
1135  {
1136  x[0]= y | z << (32 - k);
1137  z>>= k;
1138  }
1139  else
1140  x[0]= y;
1141  i= b->wds= (x[1]= z) ? 2 : 1;
1142  }
1143  else
1144  {
1145  k= lo0bits(&z);
1146  x[0]= z;
1147  i= b->wds= 1;
1148  k+= 32;
1149  }
1150  if (de)
1151  {
1152  *e= de - Bias - (P-1) + k;
1153  *bits= P - k;
1154  }
1155  else
1156  {
1157  *e= de - Bias - (P-1) + 1 + k;
1158  *bits= 32*i - hi0bits(x[i-1]);
1159  }
1160  return b;
1161 #undef d0
1162 #undef d1
1163 }
1164 
1165 
1166 static double ratio(Bigint *a, Bigint *b)
1167 {
1168  double da, db;
1169  int k, ka, kb;
1170 
1171  dval(da)= b2d(a, &ka);
1172  dval(db)= b2d(b, &kb);
1173  k= ka - kb + 32*(a->wds - b->wds);
1174  if (k > 0)
1175  word0(da)+= k*Exp_msk1;
1176  else
1177  {
1178  k= -k;
1179  word0(db)+= k*Exp_msk1;
1180  }
1181  return dval(da) / dval(db);
1182 }
1183 
1184 static const double tens[] =
1185 {
1186  1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1187  1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1188  1e20, 1e21, 1e22
1189 };
1190 
1191 static const double bigtens[]= { 1e16, 1e32, 1e64, 1e128, 1e256 };
1192 static const double tinytens[]=
1193 { 1e-16, 1e-32, 1e-64, 1e-128,
1194  9007199254740992.*9007199254740992.e-256 /* = 2^106 * 1e-53 */
1195 };
1196 /*
1197  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow
1198  flag unnecessarily. It leads to a song and dance at the end of strtod.
1199 */
1200 #define Scale_Bit 0x10
1201 #define n_bigtens 5
1202 
1203 /*
1204  strtod for IEEE--arithmetic machines.
1205 
1206  This strtod returns a nearest machine number to the input decimal
1207  string (or sets errno to EOVERFLOW). Ties are broken by the IEEE round-even
1208  rule.
1209 
1210  Inspired loosely by William D. Clinger's paper "How to Read Floating
1211  Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
1212 
1213  Modifications:
1214 
1215  1. We only require IEEE (not IEEE double-extended).
1216  2. We get by with floating-point arithmetic in a case that
1217  Clinger missed -- when we're computing d * 10^n
1218  for a small integer d and the integer n is not too
1219  much larger than 22 (the maximum integer k for which
1220  we can represent 10^k exactly), we may be able to
1221  compute (d*10^k) * 10^(e-k) with just one roundoff.
1222  3. Rather than a bit-at-a-time adjustment of the binary
1223  result in the hard case, we use floating-point
1224  arithmetic to determine the adjustment to within
1225  one bit; only in really hard cases do we need to
1226  compute a second residual.
1227  4. Because of 3., we don't need a large table of powers of 10
1228  for ten-to-e (just some small tables, e.g. of 10^k
1229  for 0 <= k <= 22).
1230 */
1231 
1232 static double my_strtod_int(const char *s00, char **se, int *error)
1233 {
1234  int scale;
1235  int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1236  e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1237  const char *s, *s0, *s1, *end = *se;
1238  double aadj, aadj1, adj, rv, rv0;
1239  Long L;
1240  ULong y, z;
1241  Bigint *bb= NULL, *bb1= NULL, *bd= NULL, *bd0= NULL, *bs= NULL, *delta= NULL;
1242 #ifdef SET_INEXACT
1243  int inexact, oldinexact;
1244 #endif
1245 
1246  c= 0;
1247  *error= 0;
1248 
1249  sign= nz0= nz= 0;
1250  dval(rv)= 0.;
1251  for (s= s00; s < end; s++)
1252  switch (*s) {
1253  case '-':
1254  sign= 1;
1255  /* no break */
1256  case '+':
1257  s++;
1258  goto break2;
1259  case '\t':
1260  case '\n':
1261  case '\v':
1262  case '\f':
1263  case '\r':
1264  case ' ':
1265  continue;
1266  default:
1267  goto break2;
1268  }
1269  break2:
1270  if (s >= end)
1271  goto ret0;
1272 
1273  if (*s == '0')
1274  {
1275  nz0= 1;
1276  while (++s < end && *s == '0') ;
1277  if (s >= end)
1278  goto ret;
1279  }
1280  s0= s;
1281  y= z= 0;
1282  for (nd= nf= 0; s < end && (c= *s) >= '0' && c <= '9'; nd++, s++)
1283  if (nd < 9)
1284  y= 10*y + c - '0';
1285  else if (nd < 16)
1286  z= 10*z + c - '0';
1287  nd0= nd;
1288  if (s < end - 1 && c == '.')
1289  {
1290  c= *++s;
1291  if (!nd)
1292  {
1293  for (; s < end && c == '0'; c= *++s)
1294  nz++;
1295  if (s < end && c > '0' && c <= '9')
1296  {
1297  s0= s;
1298  nf+= nz;
1299  nz= 0;
1300  goto have_dig;
1301  }
1302  goto dig_done;
1303  }
1304  for (; s < end && c >= '0' && c <= '9'; c = *++s)
1305  {
1306  have_dig:
1307  nz++;
1308  if (c-= '0')
1309  {
1310  nf+= nz;
1311  for (i= 1; i < nz; i++)
1312  if (nd++ < 9)
1313  y*= 10;
1314  else if (nd <= DBL_DIG + 1)
1315  z*= 10;
1316  if (nd++ < 9)
1317  y= 10*y + c;
1318  else if (nd <= DBL_DIG + 1)
1319  z= 10*z + c;
1320  nz= 0;
1321  }
1322  }
1323  }
1324  dig_done:
1325  e= 0;
1326  if (s < end && (c == 'e' || c == 'E'))
1327  {
1328  if (!nd && !nz && !nz0)
1329  goto ret0;
1330  s00= s;
1331  esign= 0;
1332  if (++s < end)
1333  switch (c= *s) {
1334  case '-':
1335  esign= 1;
1336  case '+':
1337  c= *++s;
1338  }
1339  if (s < end && c >= '0' && c <= '9')
1340  {
1341  while (s < end && c == '0')
1342  c= *++s;
1343  if (s < end && c > '0' && c <= '9') {
1344  L= c - '0';
1345  s1= s;
1346  while (++s < end && (c= *s) >= '0' && c <= '9')
1347  L= 10*L + c - '0';
1348  if (s - s1 > 8 || L > 19999)
1349  /* Avoid confusion from exponents
1350  * so large that e might overflow.
1351  */
1352  e= 19999; /* safe for 16 bit ints */
1353  else
1354  e= (int)L;
1355  if (esign)
1356  e= -e;
1357  }
1358  else
1359  e= 0;
1360  }
1361  else
1362  s= s00;
1363  }
1364  if (!nd)
1365  {
1366  if (!nz && !nz0)
1367  {
1368  ret0:
1369  s= s00;
1370  sign= 0;
1371  }
1372  goto ret;
1373  }
1374  e1= e -= nf;
1375 
1376  /*
1377  Now we have nd0 digits, starting at s0, followed by a
1378  decimal point, followed by nd-nd0 digits. The number we're
1379  after is the integer represented by those digits times
1380  10**e
1381  */
1382 
1383  if (!nd0)
1384  nd0= nd;
1385  k= nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1386  dval(rv)= y;
1387  if (k > 9)
1388  {
1389 #ifdef SET_INEXACT
1390  if (k > DBL_DIG)
1391  oldinexact = get_inexact();
1392 #endif
1393  dval(rv)= tens[k - 9] * dval(rv) + z;
1394  }
1395  bd0= 0;
1396  if (nd <= DBL_DIG
1397  )
1398  {
1399  if (!e)
1400  goto ret;
1401  if (e > 0)
1402  {
1403  if (e <= Ten_pmax)
1404  {
1405  /* rv = */ rounded_product(dval(rv), tens[e]);
1406  goto ret;
1407  }
1408  i= DBL_DIG - nd;
1409  if (e <= Ten_pmax + i)
1410  {
1411  /*
1412  A fancier test would sometimes let us do
1413  this for larger i values.
1414  */
1415  e-= i;
1416  dval(rv)*= tens[i];
1417  /* rv = */ rounded_product(dval(rv), tens[e]);
1418  goto ret;
1419  }
1420  }
1421 #ifndef Inaccurate_Divide
1422  else if (e >= -Ten_pmax)
1423  {
1424  /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1425  goto ret;
1426  }
1427 #endif
1428  }
1429  e1+= nd - k;
1430 
1431 #ifdef SET_INEXACT
1432  inexact= 1;
1433  if (k <= DBL_DIG)
1434  oldinexact= get_inexact();
1435 #endif
1436  scale= 0;
1437 
1438  /* Get starting approximation = rv * 10**e1 */
1439 
1440  if (e1 > 0)
1441  {
1442  if ((i= e1 & 15))
1443  dval(rv)*= tens[i];
1444  if (e1&= ~15)
1445  {
1446  if (e1 > DBL_MAX_10_EXP)
1447  {
1448  ovfl:
1449  *error= EOVERFLOW;
1450  /* Can't trust HUGE_VAL */
1451  word0(rv)= Exp_mask;
1452  word1(rv)= 0;
1453 #ifdef SET_INEXACT
1454  /* set overflow bit */
1455  dval(rv0)= 1e300;
1456  dval(rv0)*= dval(rv0);
1457 #endif
1458  if (bd0)
1459  goto retfree;
1460  goto ret;
1461  }
1462  e1>>= 4;
1463  for(j= 0; e1 > 1; j++, e1>>= 1)
1464  if (e1 & 1)
1465  dval(rv)*= bigtens[j];
1466  /* The last multiplication could overflow. */
1467  word0(rv)-= P*Exp_msk1;
1468  dval(rv)*= bigtens[j];
1469  if ((z= word0(rv) & Exp_mask) > Exp_msk1 * (DBL_MAX_EXP + Bias - P))
1470  goto ovfl;
1471  if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P))
1472  {
1473  /* set to largest number (Can't trust DBL_MAX) */
1474  word0(rv)= Big0;
1475  word1(rv)= Big1;
1476  }
1477  else
1478  word0(rv)+= P*Exp_msk1;
1479  }
1480  }
1481  else if (e1 < 0)
1482  {
1483  e1= -e1;
1484  if ((i= e1 & 15))
1485  dval(rv)/= tens[i];
1486  if ((e1>>= 4))
1487  {
1488  if (e1 >= 1 << n_bigtens)
1489  goto undfl;
1490  if (e1 & Scale_Bit)
1491  scale= 2 * P;
1492  for(j= 0; e1 > 0; j++, e1>>= 1)
1493  if (e1 & 1)
1494  dval(rv)*= tinytens[j];
1495  if (scale && (j = 2 * P + 1 - ((word0(rv) & Exp_mask) >> Exp_shift)) > 0)
1496  {
1497  /* scaled rv is denormal; zap j low bits */
1498  if (j >= 32)
1499  {
1500  word1(rv)= 0;
1501  if (j >= 53)
1502  word0(rv)= (P + 2) * Exp_msk1;
1503  else
1504  word0(rv)&= 0xffffffff << (j - 32);
1505  }
1506  else
1507  word1(rv)&= 0xffffffff << j;
1508  }
1509  if (!dval(rv))
1510  {
1511  undfl:
1512  dval(rv)= 0.;
1513  if (bd0)
1514  goto retfree;
1515  goto ret;
1516  }
1517  }
1518  }
1519 
1520  /* Now the hard part -- adjusting rv to the correct value.*/
1521 
1522  /* Put digits into bd: true value = bd * 10^e */
1523 
1524  bd0= s2b(s0, nd0, nd, y);
1525 
1526  for(;;)
1527  {
1528  bd= Balloc(bd0->k);
1529  Bcopy(bd, bd0);
1530  bb= d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
1531  bs= i2b(1);
1532 
1533  if (e >= 0)
1534  {
1535  bb2= bb5= 0;
1536  bd2= bd5= e;
1537  }
1538  else
1539  {
1540  bb2= bb5= -e;
1541  bd2= bd5= 0;
1542  }
1543  if (bbe >= 0)
1544  bb2+= bbe;
1545  else
1546  bd2-= bbe;
1547  bs2= bb2;
1548  j= bbe - scale;
1549  i= j + bbbits - 1; /* logb(rv) */
1550  if (i < Emin) /* denormal */
1551  j+= P - Emin;
1552  else
1553  j= P + 1 - bbbits;
1554  bb2+= j;
1555  bd2+= j;
1556  bd2+= scale;
1557  i= bb2 < bd2 ? bb2 : bd2;
1558  if (i > bs2)
1559  i= bs2;
1560  if (i > 0)
1561  {
1562  bb2-= i;
1563  bd2-= i;
1564  bs2-= i;
1565  }
1566  if (bb5 > 0)
1567  {
1568  bs= pow5mult(bs, bb5);
1569  bb1= mult(bs, bb);
1570  Bfree(bb);
1571  bb= bb1;
1572  }
1573  if (bb2 > 0)
1574  bb= lshift(bb, bb2);
1575  if (bd5 > 0)
1576  bd= pow5mult(bd, bd5);
1577  if (bd2 > 0)
1578  bd= lshift(bd, bd2);
1579  if (bs2 > 0)
1580  bs= lshift(bs, bs2);
1581  delta= diff(bb, bd);
1582  dsign= delta->sign;
1583  delta->sign= 0;
1584  i= cmp(delta, bs);
1585 
1586  if (i < 0)
1587  {
1588  /*
1589  Error is less than half an ulp -- check for special case of mantissa
1590  a power of two.
1591  */
1592  if (dsign || word1(rv) || word0(rv) & Bndry_mask ||
1593  (word0(rv) & Exp_mask) <= (2 * P + 1) * Exp_msk1)
1594  {
1595 #ifdef SET_INEXACT
1596  if (!delta->x[0] && delta->wds <= 1)
1597  inexact= 0;
1598 #endif
1599  break;
1600  }
1601  if (!delta->p.x[0] && delta->wds <= 1)
1602  {
1603  /* exact result */
1604 #ifdef SET_INEXACT
1605  inexact= 0;
1606 #endif
1607  break;
1608  }
1609  delta= lshift(delta, Log2P);
1610  if (cmp(delta, bs) > 0)
1611  goto drop_down;
1612  break;
1613  }
1614  if (i == 0)
1615  {
1616  /* exactly half-way between */
1617  if (dsign)
1618  {
1619  if ((word0(rv) & Bndry_mask1) == Bndry_mask1 &&
1620  word1(rv) ==
1621  ((scale && (y = word0(rv) & Exp_mask) <= 2 * P * Exp_msk1) ?
1622  (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
1623  0xffffffff))
1624  {
1625  /*boundary case -- increment exponent*/
1626  word0(rv)= (word0(rv) & Exp_mask) + Exp_msk1;
1627  word1(rv) = 0;
1628  dsign = 0;
1629  break;
1630  }
1631  }
1632  else if (!(word0(rv) & Bndry_mask) && !word1(rv))
1633  {
1634  drop_down:
1635  /* boundary case -- decrement exponent */
1636  if (scale)
1637  {
1638  L= word0(rv) & Exp_mask;
1639  if (L <= (2 *P + 1) * Exp_msk1)
1640  {
1641  if (L > (P + 2) * Exp_msk1)
1642  /* round even ==> accept rv */
1643  break;
1644  /* rv = smallest denormal */
1645  goto undfl;
1646  }
1647  }
1648  L= (word0(rv) & Exp_mask) - Exp_msk1;
1649  word0(rv)= L | Bndry_mask1;
1650  word1(rv)= 0xffffffff;
1651  break;
1652  }
1653  if (!(word1(rv) & LSB))
1654  break;
1655  if (dsign)
1656  dval(rv)+= ulp(dval(rv));
1657  else
1658  {
1659  dval(rv)-= ulp(dval(rv));
1660  if (!dval(rv))
1661  goto undfl;
1662  }
1663  dsign= 1 - dsign;
1664  break;
1665  }
1666  if ((aadj= ratio(delta, bs)) <= 2.)
1667  {
1668  if (dsign)
1669  aadj= aadj1= 1.;
1670  else if (word1(rv) || word0(rv) & Bndry_mask)
1671  {
1672  if (word1(rv) == Tiny1 && !word0(rv))
1673  goto undfl;
1674  aadj= 1.;
1675  aadj1= -1.;
1676  }
1677  else
1678  {
1679  /* special case -- power of FLT_RADIX to be rounded down... */
1680  if (aadj < 2. / FLT_RADIX)
1681  aadj= 1. / FLT_RADIX;
1682  else
1683  aadj*= 0.5;
1684  aadj1= -aadj;
1685  }
1686  }
1687  else
1688  {
1689  aadj*= 0.5;
1690  aadj1= dsign ? aadj : -aadj;
1691  if (Flt_Rounds == 0)
1692  aadj1+= 0.5;
1693  }
1694  y= word0(rv) & Exp_mask;
1695 
1696  /* Check for overflow */
1697 
1698  if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1))
1699  {
1700  dval(rv0)= dval(rv);
1701  word0(rv)-= P * Exp_msk1;
1702  adj= aadj1 * ulp(dval(rv));
1703  dval(rv)+= adj;
1704  if ((word0(rv) & Exp_mask) >= Exp_msk1 * (DBL_MAX_EXP + Bias - P))
1705  {
1706  if (word0(rv0) == Big0 && word1(rv0) == Big1)
1707  goto ovfl;
1708  word0(rv)= Big0;
1709  word1(rv)= Big1;
1710  goto cont;
1711  }
1712  else
1713  word0(rv)+= P * Exp_msk1;
1714  }
1715  else
1716  {
1717  if (scale && y <= 2 * P * Exp_msk1)
1718  {
1719  if (aadj <= 0x7fffffff)
1720  {
1721  if ((z= (ULong) aadj) <= 0)
1722  z= 1;
1723  aadj= z;
1724  aadj1= dsign ? aadj : -aadj;
1725  }
1726  word0(aadj1)+= (2 * P + 1) * Exp_msk1 - y;
1727  }
1728  adj = aadj1 * ulp(dval(rv));
1729  dval(rv) += adj;
1730  }
1731  z= word0(rv) & Exp_mask;
1732 #ifndef SET_INEXACT
1733  if (!scale)
1734  if (y == z)
1735  {
1736  /* Can we stop now? */
1737  L= (Long)aadj;
1738  aadj-= L;
1739  /* The tolerances below are conservative. */
1740  if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1741  {
1742  if (aadj < .4999999 || aadj > .5000001)
1743  break;
1744  }
1745  else if (aadj < .4999999 / FLT_RADIX)
1746  break;
1747  }
1748 #endif
1749  cont:
1750  Bfree(bb);
1751  Bfree(bd);
1752  Bfree(bs);
1753  Bfree(delta);
1754  }
1755 #ifdef SET_INEXACT
1756  if (inexact)
1757  {
1758  if (!oldinexact)
1759  {
1760  word0(rv0)= Exp_1 + (70 << Exp_shift);
1761  word1(rv0)= 0;
1762  dval(rv0)+= 1.;
1763  }
1764  }
1765  else if (!oldinexact)
1766  clear_inexact();
1767 #endif
1768  if (scale)
1769  {
1770  word0(rv0)= Exp_1 - 2 * P * Exp_msk1;
1771  word1(rv0)= 0;
1772  dval(rv)*= dval(rv0);
1773  }
1774 #ifdef SET_INEXACT
1775  if (inexact && !(word0(rv) & Exp_mask))
1776  {
1777  /* set underflow bit */
1778  dval(rv0)= 1e-300;
1779  dval(rv0)*= dval(rv0);
1780  }
1781 #endif
1782  retfree:
1783  Bfree(bb);
1784  Bfree(bd);
1785  Bfree(bs);
1786  Bfree(bd0);
1787  Bfree(delta);
1788  ret:
1789  *se= (char *)s;
1790  return sign ? -dval(rv) : dval(rv);
1791 }
1792 
1793 
1794 static int quorem(Bigint *b, Bigint *S)
1795 {
1796  int n;
1797  ULong *bx, *bxe, q, *sx, *sxe;
1798  ULLong borrow, carry, y, ys;
1799 
1800  n= S->wds;
1801  if (b->wds < n)
1802  return 0;
1803  sx= S->p.x;
1804  sxe= sx + --n;
1805  bx= b->p.x;
1806  bxe= bx + n;
1807  q= *bxe / (*sxe + 1); /* ensure q <= true quotient */
1808  if (q)
1809  {
1810  borrow= 0;
1811  carry= 0;
1812  do
1813  {
1814  ys= *sx++ * (ULLong)q + carry;
1815  carry= ys >> 32;
1816  y= *bx - (ys & FFFFFFFF) - borrow;
1817  borrow= y >> 32 & (ULong)1;
1818  *bx++= (ULong) (y & FFFFFFFF);
1819  }
1820  while (sx <= sxe);
1821  if (!*bxe)
1822  {
1823  bx= b->p.x;
1824  while (--bxe > bx && !*bxe)
1825  --n;
1826  b->wds= n;
1827  }
1828  }
1829  if (cmp(b, S) >= 0)
1830  {
1831  q++;
1832  borrow= 0;
1833  carry= 0;
1834  bx= b->p.x;
1835  sx= S->p.x;
1836  do
1837  {
1838  ys= *sx++ + carry;
1839  carry= ys >> 32;
1840  y= *bx - (ys & FFFFFFFF) - borrow;
1841  borrow= y >> 32 & (ULong)1;
1842  *bx++= (ULong) (y & FFFFFFFF);
1843  }
1844  while (sx <= sxe);
1845  bx= b->p.x;
1846  bxe= bx + n;
1847  if (!*bxe)
1848  {
1849  while (--bxe > bx && !*bxe)
1850  --n;
1851  b->wds= n;
1852  }
1853  }
1854  return q;
1855 }
1856 
1857 
1858 /*
1859  dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1860 
1861  Inspired by "How to Print Floating-Point Numbers Accurately" by
1862  Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
1863 
1864  Modifications:
1865  1. Rather than iterating, we use a simple numeric overestimate
1866  to determine k= floor(log10(d)). We scale relevant
1867  quantities using O(log2(k)) rather than O(k) multiplications.
1868  2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1869  try to generate digits strictly left to right. Instead, we
1870  compute with fewer bits and propagate the carry if necessary
1871  when rounding the final digit up. This is often faster.
1872  3. Under the assumption that input will be rounded nearest,
1873  mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1874  That is, we allow equality in stopping tests when the
1875  round-nearest rule will give the same floating-point value
1876  as would satisfaction of the stopping test with strict
1877  inequality.
1878  4. We remove common factors of powers of 2 from relevant
1879  quantities.
1880  5. When converting floating-point integers less than 1e16,
1881  we use floating-point arithmetic rather than resorting
1882  to multiple-precision integers.
1883  6. When asked to produce fewer than 15 digits, we first try
1884  to get by with floating-point arithmetic; we resort to
1885  multiple-precision integer arithmetic only if we cannot
1886  guarantee that the floating-point calculation has given
1887  the correctly rounded result. For k requested digits and
1888  "uniformly" distributed input, the probability is
1889  something like 10^(k-15) that we must resort to the Long
1890  calculation.
1891  */
1892 
1893 static char *dtoa(double d, int mode, int ndigits, int *decpt, int *sign,
1894  char **rve)
1895 {
1896  /*
1897  Arguments ndigits, decpt, sign are similar to those
1898  of ecvt and fcvt; trailing zeros are suppressed from
1899  the returned string. If not null, *rve is set to point
1900  to the end of the return value. If d is +-Infinity or NaN,
1901  then *decpt is set to DTOA_OVERFLOW.
1902 
1903  mode:
1904  0 ==> shortest string that yields d when read in
1905  and rounded to nearest.
1906 
1907  1 ==> like 0, but with Steele & White stopping rule;
1908  e.g. with IEEE P754 arithmetic , mode 0 gives
1909  1e23 whereas mode 1 gives 9.999999999999999e22.
1910  2 ==> cmax(1,ndigits) significant digits. This gives a
1911  return value similar to that of ecvt, except
1912  that trailing zeros are suppressed.
1913  3 ==> through ndigits past the decimal point. This
1914  gives a return value similar to that from fcvt,
1915  except that trailing zeros are suppressed, and
1916  ndigits can be negative.
1917  4,5 ==> similar to 2 and 3, respectively, but (in
1918  round-nearest mode) with the tests of mode 0 to
1919  possibly return a shorter string that rounds to d.
1920  6-9 ==> Debugging modes similar to mode - 4: don't try
1921  fast floating-point estimate (if applicable).
1922 
1923  Values of mode other than 0-9 are treated as mode 0.
1924 
1925  Sufficient space is allocated to the return value
1926  to hold the suppressed trailing zeros.
1927  */
1928 
1929  int bbits, b2, b5, be, dig, i, ieps, ilim=0, ilim0, ilim1= 0,
1930  j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1931  spec_case, try_quick;
1932  Long L;
1933  int denorm;
1934  ULong x;
1935  Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
1936  double d2, ds, eps;
1937  char *s, *s0;
1938 
1939  if (word0(d) & Sign_bit)
1940  {
1941  /* set sign for everything, including 0's and NaNs */
1942  *sign= 1;
1943  word0(d) &= ~Sign_bit; /* clear sign bit */
1944  }
1945  else
1946  *sign= 0;
1947 
1948  /* If infinity, set decpt to DTOA_OVERFLOW, if 0 set it to 1 */
1949  if ((((word0(d) & Exp_mask) == Exp_mask) && ((*decpt= DTOA_OVERFLOW) != 0)) ||
1950  (!dval(d) && ((*decpt= 1) != 0)))
1951  {
1952  /* Infinity, NaN, 0 */
1953  char *res= (char*) malloc(2);
1954  res[0]= '0';
1955  res[1]= '\0';
1956  if (rve)
1957  *rve= res + 1;
1958  return res;
1959  }
1960 
1961 
1962  b= d2b(dval(d), &be, &bbits);
1963  if ((i= (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
1964  {
1965  dval(d2)= dval(d);
1966  word0(d2) &= Frac_mask1;
1967  word0(d2) |= Exp_11;
1968 
1969  /*
1970  log(x) ~=~ log(1.5) + (x-1.5)/1.5
1971  log10(x) = log(x) / log(10)
1972  ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1973  log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
1974 
1975  This suggests computing an approximation k to log10(d) by
1976 
1977  k= (i - Bias)*0.301029995663981
1978  + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1979 
1980  We want k to be too large rather than too small.
1981  The error in the first-order Taylor series approximation
1982  is in our favor, so we just round up the constant enough
1983  to compensate for any error in the multiplication of
1984  (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1985  and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1986  adding 1e-13 to the constant term more than suffices.
1987  Hence we adjust the constant term to 0.1760912590558.
1988  (We could get a more accurate k by invoking log10,
1989  but this is probably not worthwhile.)
1990  */
1991 
1992  i-= Bias;
1993  denorm= 0;
1994  }
1995  else
1996  {
1997  /* d is denormalized */
1998 
1999  i= bbits + be + (Bias + (P-1) - 1);
2000  x= i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
2001  : word1(d) << (32 - i);
2002  dval(d2)= x;
2003  word0(d2)-= 31*Exp_msk1; /* adjust exponent */
2004  i-= (Bias + (P-1) - 1) + 1;
2005  denorm= 1;
2006  }
2007  ds= (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2008  k= (int)ds;
2009  if (ds < 0. && ds != k)
2010  k--; /* want k= floor(ds) */
2011  k_check= 1;
2012  if (k >= 0 && k <= Ten_pmax)
2013  {
2014  if (dval(d) < tens[k])
2015  k--;
2016  k_check= 0;
2017  }
2018  j= bbits - i - 1;
2019  if (j >= 0)
2020  {
2021  b2= 0;
2022  s2= j;
2023  }
2024  else
2025  {
2026  b2= -j;
2027  s2= 0;
2028  }
2029  if (k >= 0)
2030  {
2031  b5= 0;
2032  s5= k;
2033  s2+= k;
2034  }
2035  else
2036  {
2037  b2-= k;
2038  b5= -k;
2039  s5= 0;
2040  }
2041  if (mode < 0 || mode > 9)
2042  mode= 0;
2043 
2044  try_quick= 1;
2045 
2046  if (mode > 5)
2047  {
2048  mode-= 4;
2049  try_quick= 0;
2050  }
2051  leftright= 1;
2052  switch (mode) {
2053  case 0:
2054  case 1:
2055  ilim= ilim1= -1;
2056  i= 18;
2057  ndigits= 0;
2058  break;
2059  case 2:
2060  leftright= 0;
2061  /* no break */
2062  case 4:
2063  if (ndigits <= 0)
2064  ndigits= 1;
2065  ilim= ilim1= i= ndigits;
2066  break;
2067  case 3:
2068  leftright= 0;
2069  /* no break */
2070  case 5:
2071  i= ndigits + k + 1;
2072  ilim= i;
2073  ilim1= i - 1;
2074  if (i <= 0)
2075  i= 1;
2076  }
2077  s= s0= (char*)malloc(i+1); /* +1 for trailing '\0' appended
2078  at end of function */
2079 
2080  if (ilim >= 0 && ilim <= Quick_max && try_quick)
2081  {
2082  /* Try to get by with floating-point arithmetic. */
2083  i= 0;
2084  dval(d2)= dval(d);
2085  k0= k;
2086  ilim0= ilim;
2087  ieps= 2; /* conservative */
2088  if (k > 0)
2089  {
2090  ds= tens[k&0xf];
2091  j= k >> 4;
2092  if (j & Bletch)
2093  {
2094  /* prevent overflows */
2095  j&= Bletch - 1;
2096  dval(d)/= bigtens[n_bigtens-1];
2097  ieps++;
2098  }
2099  for (; j; j>>= 1, i++)
2100  {
2101  if (j & 1)
2102  {
2103  ieps++;
2104  ds*= bigtens[i];
2105  }
2106  }
2107  dval(d)/= ds;
2108  }
2109  else if ((j1= -k))
2110  {
2111  dval(d)*= tens[j1 & 0xf];
2112  for (j= j1 >> 4; j; j>>= 1, i++)
2113  {
2114  if (j & 1)
2115  {
2116  ieps++;
2117  dval(d)*= bigtens[i];
2118  }
2119  }
2120  }
2121  if (k_check && dval(d) < 1. && ilim > 0)
2122  {
2123  if (ilim1 <= 0)
2124  goto fast_failed;
2125  ilim= ilim1;
2126  k--;
2127  dval(d)*= 10.;
2128  ieps++;
2129  }
2130  dval(eps)= ieps*dval(d) + 7.;
2131  word0(eps)-= (P-1)*Exp_msk1;
2132  if (ilim == 0)
2133  {
2134  S= mhi= 0;
2135  dval(d)-= 5.;
2136  if (dval(d) > dval(eps))
2137  goto one_digit;
2138  if (dval(d) < -dval(eps))
2139  goto no_digits;
2140  goto fast_failed;
2141  }
2142  if (leftright)
2143  {
2144  /* Use Steele & White method of only generating digits needed. */
2145  dval(eps)= 0.5/tens[ilim-1] - dval(eps);
2146  for (i= 0;;)
2147  {
2148  L= (Long) dval(d);
2149  dval(d)-= L;
2150  *s++= '0' + (int)L;
2151  if (dval(d) < dval(eps))
2152  goto ret1;
2153  if (1. - dval(d) < dval(eps))
2154  goto bump_up;
2155  if (++i >= ilim)
2156  break;
2157  dval(eps)*= 10.;
2158  dval(d)*= 10.;
2159  }
2160  }
2161  else
2162  {
2163  /* Generate ilim digits, then fix them up. */
2164  dval(eps)*= tens[ilim-1];
2165  for (i= 1;; i++, dval(d)*= 10.)
2166  {
2167  L= (Long)(dval(d));
2168  if (!(dval(d)-= L))
2169  ilim= i;
2170  *s++= '0' + (int)L;
2171  if (i == ilim)
2172  {
2173  if (dval(d) > 0.5 + dval(eps))
2174  goto bump_up;
2175  else if (dval(d) < 0.5 - dval(eps))
2176  {
2177  while (*--s == '0') {}
2178  s++;
2179  goto ret1;
2180  }
2181  break;
2182  }
2183  }
2184  }
2185  fast_failed:
2186  s= s0;
2187  dval(d)= dval(d2);
2188  k= k0;
2189  ilim= ilim0;
2190  }
2191 
2192  /* Do we have a "small" integer? */
2193 
2194  if (be >= 0 && k <= Int_max)
2195  {
2196  /* Yes. */
2197  ds= tens[k];
2198  if (ndigits < 0 && ilim <= 0)
2199  {
2200  S= mhi= 0;
2201  if (ilim < 0 || dval(d) <= 5*ds)
2202  goto no_digits;
2203  goto one_digit;
2204  }
2205  for (i= 1;; i++, dval(d)*= 10.)
2206  {
2207  L= (Long)(dval(d) / ds);
2208  dval(d)-= L*ds;
2209  *s++= '0' + (int)L;
2210  if (!dval(d))
2211  {
2212  break;
2213  }
2214  if (i == ilim)
2215  {
2216  dval(d)+= dval(d);
2217  if (dval(d) > ds || (dval(d) == ds && L & 1))
2218  {
2219 bump_up:
2220  while (*--s == '9')
2221  if (s == s0)
2222  {
2223  k++;
2224  *s= '0';
2225  break;
2226  }
2227  ++*s++;
2228  }
2229  break;
2230  }
2231  }
2232  goto ret1;
2233  }
2234 
2235  m2= b2;
2236  m5= b5;
2237  mhi= mlo= 0;
2238  if (leftright)
2239  {
2240  i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
2241  b2+= i;
2242  s2+= i;
2243  mhi= i2b(1);
2244  }
2245  if (m2 > 0 && s2 > 0)
2246  {
2247  i= m2 < s2 ? m2 : s2;
2248  b2-= i;
2249  m2-= i;
2250  s2-= i;
2251  }
2252  if (b5 > 0)
2253  {
2254  if (leftright)
2255  {
2256  if (m5 > 0)
2257  {
2258  mhi= pow5mult(mhi, m5);
2259  b1= mult(mhi, b);
2260  Bfree(b);
2261  b= b1;
2262  }
2263  if ((j= b5 - m5))
2264  b= pow5mult(b, j);
2265  }
2266  else
2267  b= pow5mult(b, b5);
2268  }
2269  S= i2b(1);
2270  if (s5 > 0)
2271  S= pow5mult(S, s5);
2272 
2273  /* Check for special case that d is a normalized power of 2. */
2274 
2275  spec_case= 0;
2276  if ((mode < 2 || leftright)
2277  )
2278  {
2279  if (!word1(d) && !(word0(d) & Bndry_mask) &&
2280  word0(d) & (Exp_mask & ~Exp_msk1)
2281  )
2282  {
2283  /* The special case */
2284  b2+= Log2P;
2285  s2+= Log2P;
2286  spec_case= 1;
2287  }
2288  }
2289 
2290  /*
2291  Arrange for convenient computation of quotients:
2292  shift left if necessary so divisor has 4 leading 0 bits.
2293 
2294  Perhaps we should just compute leading 28 bits of S once
2295  a nd for all and pass them and a shift to quorem, so it
2296  can do shifts and ors to compute the numerator for q.
2297  */
2298  if ((i= ((s5 ? 32 - hi0bits(S->p.x[S->wds-1]) : 1) + s2) & 0x1f))
2299  i= 32 - i;
2300  if (i > 4)
2301  {
2302  i-= 4;
2303  b2+= i;
2304  m2+= i;
2305  s2+= i;
2306  }
2307  else if (i < 4)
2308  {
2309  i+= 28;
2310  b2+= i;
2311  m2+= i;
2312  s2+= i;
2313  }
2314  if (b2 > 0)
2315  b= lshift(b, b2);
2316  if (s2 > 0)
2317  S= lshift(S, s2);
2318  if (k_check)
2319  {
2320  if (cmp(b,S) < 0)
2321  {
2322  k--;
2323  /* we botched the k estimate */
2324  b= multadd(b, 10, 0);
2325  if (leftright)
2326  mhi= multadd(mhi, 10, 0);
2327  ilim= ilim1;
2328  }
2329  }
2330  if (ilim <= 0 && (mode == 3 || mode == 5))
2331  {
2332  if (ilim < 0 || cmp(b,S= multadd(S,5,0)) <= 0)
2333  {
2334  /* no digits, fcvt style */
2335 no_digits:
2336  k= -1 - ndigits;
2337  goto ret;
2338  }
2339 one_digit:
2340  *s++= '1';
2341  k++;
2342  goto ret;
2343  }
2344  if (leftright)
2345  {
2346  if (m2 > 0)
2347  mhi= lshift(mhi, m2);
2348 
2349  /*
2350  Compute mlo -- check for special case that d is a normalized power of 2.
2351  */
2352 
2353  mlo= mhi;
2354  if (spec_case)
2355  {
2356  mhi= Balloc(mhi->k);
2357  Bcopy(mhi, mlo);
2358  mhi= lshift(mhi, Log2P);
2359  }
2360 
2361  for (i= 1;;i++)
2362  {
2363  dig= quorem(b,S) + '0';
2364  /* Do we yet have the shortest decimal string that will round to d? */
2365  j= cmp(b, mlo);
2366  delta= diff(S, mhi);
2367  j1= delta->sign ? 1 : cmp(b, delta);
2368  Bfree(delta);
2369  if (j1 == 0 && mode != 1 && !(word1(d) & 1)
2370  )
2371  {
2372  if (dig == '9')
2373  goto round_9_up;
2374  if (j > 0)
2375  dig++;
2376  *s++= dig;
2377  goto ret;
2378  }
2379  if (j < 0 || (j == 0 && mode != 1 && !(word1(d) & 1)))
2380  {
2381  if (!b->p.x[0] && b->wds <= 1)
2382  {
2383  goto accept_dig;
2384  }
2385  if (j1 > 0)
2386  {
2387  b= lshift(b, 1);
2388  j1= cmp(b, S);
2389  if ((j1 > 0 || (j1 == 0 && dig & 1))
2390  && dig++ == '9')
2391  goto round_9_up;
2392  }
2393 accept_dig:
2394  *s++= dig;
2395  goto ret;
2396  }
2397  if (j1 > 0)
2398  {
2399  if (dig == '9')
2400  { /* possible if i == 1 */
2401 round_9_up:
2402  *s++= '9';
2403  goto roundoff;
2404  }
2405  *s++= dig + 1;
2406  goto ret;
2407  }
2408  *s++= dig;
2409  if (i == ilim)
2410  break;
2411  b= multadd(b, 10, 0);
2412  if (mlo == mhi)
2413  mlo= mhi= multadd(mhi, 10, 0);
2414  else
2415  {
2416  mlo= multadd(mlo, 10, 0);
2417  mhi= multadd(mhi, 10, 0);
2418  }
2419  }
2420  }
2421  else
2422  for (i= 1;; i++)
2423  {
2424  *s++= dig= quorem(b,S) + '0';
2425  if (!b->p.x[0] && b->wds <= 1)
2426  {
2427  goto ret;
2428  }
2429  if (i >= ilim)
2430  break;
2431  b= multadd(b, 10, 0);
2432  }
2433 
2434  /* Round off last digit */
2435 
2436  b= lshift(b, 1);
2437  j= cmp(b, S);
2438  if (j > 0 || (j == 0 && dig & 1))
2439  {
2440 roundoff:
2441  while (*--s == '9')
2442  if (s == s0)
2443  {
2444  k++;
2445  *s++= '1';
2446  goto ret;
2447  }
2448  ++*s++;
2449  }
2450  else
2451  {
2452  while (*--s == '0') {}
2453  s++;
2454  }
2455 ret:
2456  Bfree(S);
2457  if (mhi)
2458  {
2459  if (mlo && mlo != mhi)
2460  Bfree(mlo);
2461  Bfree(mhi);
2462  }
2463 ret1:
2464  Bfree(b);
2465  *s= 0;
2466  *decpt= k + 1;
2467  if (rve)
2468  *rve= s;
2469  return s0;
2470 }
2471 
2472 } /* namespace internal */
2473 } /* namespace drizzled */