← Back to branch summary

~drizzle-trunk/drizzle/development

« back to all changes in this revision

Viewing changes to strings/dtoa.c

Merge/fix in FAQ.

Show diffs side-by-side

added added

removed removed

Lines of Context:
strings/dtoa.c
11
11
 
12
12
   You should have received a copy of the GNU General Public License
13
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 */
 
14
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA */
15
15
 
16
16
/****************************************************************
17
17
 
20
20
 
21
21
  The author of this software is David M. Gay.
22
22
 
23
 
  Copyright (C) 1991, 2000, 2001 by Lucent Technologies.
 
23
  Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
24
24
 
25
25
  Permission to use, copy, modify, and distribute this software for any
26
26
  purpose without fee is hereby granted, provided that this entire notice
35
35
 
36
36
 ***************************************************************/
37
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
 
{
 
38
#include <my_base.h> /* for EOVERFLOW on Windows */
 
39
#include <my_global.h>
 
40
#include <m_string.h>  /* for memcpy and NOT_FIXED_DEC */
 
41
 
 
42
/**
 
43
   Appears to suffice to not call malloc() in most cases.
 
44
   @todo
 
45
     see if it is possible to get rid of malloc().
 
46
     this constant is sufficient to avoid malloc() on all inputs I have tried.
 
47
*/
 
48
#define DTOA_BUFF_SIZE (420 * sizeof(void *))
54
49
 
55
50
/* Magic value returned by dtoa() to indicate overflow */
56
51
#define DTOA_OVERFLOW 9999
57
52
 
58
 
static double my_strtod_int(const char *, char **, int *);
59
 
static char *dtoa(double, int, int, int *, int *, char **);
 
53
static double my_strtod_int(const char *, char **, int *, char *, size_t);
 
54
static char *dtoa(double, int, int, int *, int *, char **, char *, size_t);
 
55
static void dtoa_free(char *, char *, size_t);
60
56
 
61
57
/**
62
58
   @brief
89
85
   @return            number of written characters (excluding terminating '\0')
90
86
*/
91
87
 
92
 
size_t my_fcvt(double x, int precision, char *to, bool *error)
 
88
size_t my_fcvt(double x, int precision, char *to, my_bool *error)
93
89
{
94
90
  int decpt, sign, len, i;
95
91
  char *res, *src, *end, *dst= to;
 
92
  char buf[DTOA_BUFF_SIZE];
96
93
  assert(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
97
 
 
98
 
  res= dtoa(x, 5, precision, &decpt, &sign, &end);
 
94
  
 
95
  res= dtoa(x, 5, precision, &decpt, &sign, &end, buf, sizeof(buf));
99
96
 
100
97
  if (decpt == DTOA_OVERFLOW)
101
98
  {
102
 
    free(res);
 
99
    dtoa_free(res, buf, sizeof(buf));
103
100
    *to++= '0';
104
101
    *to= '\0';
105
102
    if (error != NULL)
134
131
  {
135
132
    if (len <= decpt)
136
133
      *dst++= '.';
137
 
 
 
134
    
138
135
    for (i= precision - max(0, (len - decpt)); i > 0; i--)
139
136
      *dst++= '0';
140
137
  }
141
 
 
 
138
  
142
139
  *dst= '\0';
143
140
  if (error != NULL)
144
141
    *error= false;
145
142
 
146
 
  free(res);
 
143
  dtoa_free(res, buf, sizeof(buf));
147
144
 
148
145
  return dst - to;
149
146
}
201
198
     my_gcvt(55, ..., 1, ...);
202
199
 
203
200
   We do our best to minimize such cases by:
204
 
 
 
201
   
205
202
   - passing to dtoa() the field width as the number of significant digits
206
 
 
 
203
   
207
204
   - removing the sign of the number early (and decreasing the width before
208
205
     passing it to dtoa())
209
 
 
 
206
   
210
207
   - choosing the proper format to preserve the most number of significant
211
208
     digits.
212
209
*/
213
210
 
214
211
size_t my_gcvt(double x, my_gcvt_arg_type type, int width, char *to,
215
 
               bool *error)
 
212
               my_bool *error)
216
213
{
217
214
  int decpt, sign, len, exp_len;
218
215
  char *res, *src, *end, *dst= to, *dend= dst + width;
219
 
  bool have_space, force_e_format;
 
216
  char buf[DTOA_BUFF_SIZE];
 
217
  my_bool have_space, force_e_format;
220
218
  assert(width > 0 && to != NULL);
221
 
 
 
219
  
222
220
  /* We want to remove '-' from equations early */
223
221
  if (x < 0.)
224
222
    width--;
225
223
 
226
 
  res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? min(width, DBL_DIG) : 
227
 
            min(width, FLT_DIG), &decpt, &sign, &end);
228
 
 
 
224
  res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? width : min(width, FLT_DIG),
 
225
            &decpt, &sign, &end, buf, sizeof(buf));
229
226
  if (decpt == DTOA_OVERFLOW)
230
227
  {
231
 
    free(res);
 
228
    dtoa_free(res, buf, sizeof(buf));
232
229
    *to++= '0';
233
230
    *to= '\0';
234
231
    if (error != NULL)
248
245
     to count it here.
249
246
   */
250
247
  exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
251
 
 
 
248
  
252
249
  /*
253
250
     Do we have enough space for all digits in the 'f' format?
254
251
     Let 'len' be the number of significant digits returned by dtoa,
301
298
       ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
302
299
                                            (len > 1 || !force_e_format)))) &&
303
300
         !force_e_format)) &&
304
 
 
 
301
      
305
302
       /*
306
303
         Use the 'e' format in some cases even if we have enough space for the
307
304
         'f' one. See comment for MAX_DECPT_FOR_F_FORMAT.
323
320
          *error= true;
324
321
        width= decpt;
325
322
      }
326
 
 
 
323
      
327
324
      /*
328
325
        We want to truncate (len - width) least significant digits after the
329
326
        decimal point. For this we are calling dtoa with mode=5, passing the
330
327
        number of significant digits = (len-decpt) - (len-width) = width-decpt
331
328
      */
332
 
      free(res);
333
 
      res= dtoa(x, 5, width - decpt, &decpt, &sign, &end);
 
329
      dtoa_free(res, buf, sizeof(buf));
 
330
      res= dtoa(x, 5, width - decpt, &decpt, &sign, &end, buf, sizeof(buf));
334
331
      src= res;
335
332
      len= end - res;
336
333
    }
341
338
      *dst++= '0';
342
339
      goto end;
343
340
    }
344
 
 
 
341
    
345
342
    /*
346
343
      At this point we are sure we have enough space to put all digits
347
344
      returned by dtoa
390
387
        *error= true;
391
388
      width= 0;
392
389
    }
393
 
 
 
390
      
394
391
    /* Do we have to truncate any digits? */
395
392
    if (width < len)
396
393
    {
397
394
      /* Yes, re-convert with a smaller width */
398
 
      free(res);
399
 
      res= dtoa(x, 4, width, &decpt, &sign, &end);
 
395
      dtoa_free(res, buf, sizeof(buf));
 
396
      res= dtoa(x, 4, width, &decpt, &sign, &end, buf, sizeof(buf));
400
397
      src= res;
401
398
      len= end - res;
402
399
      if (--decpt < 0)
436
433
  }
437
434
 
438
435
end:
439
 
  free(res);
 
436
  dtoa_free(res, buf, sizeof(buf));
440
437
  *dst= '\0';
441
438
 
442
439
  return dst - to;
455
452
                  rejected character.
456
453
   @param error   Upon return is set to EOVERFLOW in case of underflow or
457
454
                  overflow.
458
 
 
 
455
   
459
456
   @return        The resulting double value. In case of underflow, 0.0 is
460
457
                  returned. In case overflow, signed DBL_MAX is returned.
461
458
*/
462
459
 
463
460
double my_strtod(const char *str, char **end, int *error)
464
461
{
 
462
  char buf[DTOA_BUFF_SIZE];
465
463
  double res;
466
464
  assert(str != NULL && end != NULL && *end != NULL && error != NULL);
467
465
 
468
 
  res= my_strtod_int(str, end, error);
 
466
  res= my_strtod_int(str, end, error, buf, sizeof(buf));
469
467
  return (*error == 0) ? res : (res < 0 ? -DBL_MAX : DBL_MAX);
470
468
}
471
469
 
482
480
 *
483
481
 * The author of this software is David M. Gay.
484
482
 *
485
 
 * Copyright (C) 1991, 2000, 2001 by Lucent Technologies.
 
483
 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
486
484
 *
487
485
 * Permission to use, copy, modify, and distribute this software for any
488
486
 * purpose without fee is hereby granted, provided that this entire notice
532
530
  file.
533
531
*/
534
532
 
 
533
/*
 
534
  #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
 
535
       and dtoa should round accordingly.
 
536
  #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
 
537
       and Honor_FLT_ROUNDS is not #defined.
 
538
 
 
539
  TODO: check if we can get rid of the above two
 
540
*/
 
541
 
535
542
typedef int32_t Long;
536
543
typedef uint32_t ULong;
537
544
typedef int64_t LLong;
587
594
#endif
588
595
#endif /*Flt_Rounds*/
589
596
 
 
597
#ifdef Honor_FLT_ROUNDS
 
598
#define Rounding rounding
 
599
#undef Check_FLT_ROUNDS
 
600
#define Check_FLT_ROUNDS
 
601
#else
 
602
#define Rounding Flt_Rounds
 
603
#endif
 
604
 
590
605
#define rounded_product(a,b) a*= b
591
606
#define rounded_quotient(a,b) a/= b
592
607
 
594
609
#define Big1 0xffffffff
595
610
#define FFFFFFFF 0xffffffffUL
596
611
 
 
612
/* This is tested to be enough for dtoa */
 
613
 
 
614
#define Kmax 15
 
615
 
 
616
#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign,   \
 
617
                          2*sizeof(int) + y->wds*sizeof(ULong))
 
618
 
597
619
/* Arbitrary-length integer */
598
620
 
599
621
typedef struct Bigint
608
630
  int wds;                 /* current length in 32-bit words */
609
631
} Bigint;
610
632
 
611
 
static Bigint *Bcopy(Bigint* dst, Bigint* src)
 
633
 
 
634
/* A simple stack-memory based allocator for Bigints */
 
635
 
 
636
typedef struct Stack_alloc
612
637
{
613
 
  dst->sign= src->sign;
614
 
  dst->wds= src->wds;
615
 
 
616
 
  assert(dst->maxwds >= src->wds);
617
 
 
618
 
  memcpy(dst->p.x, src->p.x, src->wds*sizeof(ULong));
619
 
 
620
 
  return dst;
621
 
}
622
 
 
623
 
static Bigint *Balloc(int k)
 
638
  char *begin;
 
639
  char *free;
 
640
  char *end;
 
641
  /*
 
642
    Having list of free blocks lets us reduce maximum required amount
 
643
    of memory from ~4000 bytes to < 1680 (tested on x86).
 
644
  */
 
645
  Bigint *freelist[Kmax+1];
 
646
} Stack_alloc;
 
647
 
 
648
 
 
649
/*
 
650
  Try to allocate object on stack, and resort to malloc if all
 
651
  stack memory is used.
 
652
*/
 
653
 
 
654
static Bigint *Balloc(int k, Stack_alloc *alloc)
624
655
{
625
656
  Bigint *rv;
626
 
 
627
 
  /* TODO: some malloc failure checking */
628
 
 
629
 
  int x= 1 << k;
630
 
  rv= (Bigint*) malloc(sizeof(Bigint));
631
 
 
632
 
  rv->p.x= (ULong*)malloc(x * sizeof(ULong));
633
 
 
634
 
  rv->k= k;
635
 
  rv->maxwds= x;
636
 
 
 
657
  if (k <= Kmax &&  alloc->freelist[k])
 
658
  {
 
659
    rv= alloc->freelist[k];
 
660
    alloc->freelist[k]= rv->p.next;
 
661
  }
 
662
  else
 
663
  {
 
664
    int x, len;
 
665
 
 
666
    x= 1 << k;
 
667
    len= sizeof(Bigint) + x * sizeof(ULong);
 
668
 
 
669
    if (alloc->free + len <= alloc->end)
 
670
    {
 
671
      rv= (Bigint*) alloc->free;
 
672
      alloc->free+= len;
 
673
    }
 
674
    else
 
675
      rv= (Bigint*) malloc(len);
 
676
 
 
677
    rv->k= k;
 
678
    rv->maxwds= x;
 
679
  }
637
680
  rv->sign= rv->wds= 0;
638
 
 
 
681
  rv->p.x= (ULong*) (rv + 1);
639
682
  return rv;
640
683
}
641
684
 
645
688
  list. Otherwise call free().
646
689
*/
647
690
 
648
 
static void Bfree(Bigint *v)
649
 
{
650
 
  if(!v)
651
 
    return;
652
 
  free(v->p.x);
653
 
  free(v);
654
 
}
 
691
static void Bfree(Bigint *v, Stack_alloc *alloc)
 
692
{
 
693
  char *gptr= (char*) v;                       /* generic pointer */
 
694
  if (gptr < alloc->begin || gptr >= alloc->end)
 
695
    free(gptr);
 
696
  else if (v->k <= Kmax)
 
697
  {
 
698
    /*
 
699
      Maintain free lists only for stack objects: this way we don't
 
700
      have to bother with freeing lists in the end of dtoa;
 
701
      heap should not be used normally anyway.
 
702
    */
 
703
    v->p.next= alloc->freelist[v->k];
 
704
    alloc->freelist[v->k]= v;
 
705
  }
 
706
}
 
707
 
 
708
 
 
709
/*
 
710
  This is to place return value of dtoa in: tries to use stack
 
711
  as well, but passes by free lists management and just aligns len by
 
712
  sizeof(ULong).
 
713
*/
 
714
 
 
715
static char *dtoa_alloc(int i, Stack_alloc *alloc)
 
716
{
 
717
  char *rv;
 
718
  int aligned_size= (i + sizeof(ULong) - 1) / sizeof(ULong) * sizeof(ULong);
 
719
  if (alloc->free + aligned_size <= alloc->end)
 
720
  {
 
721
    rv= alloc->free;
 
722
    alloc->free+= aligned_size;
 
723
  }
 
724
  else
 
725
    rv= malloc(i);
 
726
  return rv;
 
727
}
 
728
 
 
729
 
 
730
/*
 
731
  dtoa_free() must be used to free values s returned by dtoa()
 
732
  This is the counterpart of dtoa_alloc()
 
733
*/
 
734
 
 
735
static void dtoa_free(char *gptr, char *buf, size_t buf_size)
 
736
{
 
737
  if (gptr < buf || gptr >= buf + buf_size)
 
738
    free(gptr);
 
739
}
 
740
 
655
741
 
656
742
/* Bigint arithmetic functions */
657
743
 
658
744
/* Multiply by m and add a */
659
745
 
660
 
static Bigint *multadd(Bigint *b, int m, int a)
 
746
static Bigint *multadd(Bigint *b, int m, int a, Stack_alloc *alloc)
661
747
{
662
748
  int i, wds;
663
749
  ULong *x;
679
765
  {
680
766
    if (wds >= b->maxwds)
681
767
    {
682
 
      b1= Balloc(b->k+1);
 
768
      b1= Balloc(b->k+1, alloc);
683
769
      Bcopy(b1, b);
684
 
      Bfree(b);
 
770
      Bfree(b, alloc);
685
771
      b= b1;
686
772
    }
687
773
    b->p.x[wds++]= (ULong) carry;
691
777
}
692
778
 
693
779
 
694
 
static Bigint *s2b(const char *s, int nd0, int nd, ULong y9)
 
780
static Bigint *s2b(const char *s, int nd0, int nd, ULong y9, Stack_alloc *alloc)
695
781
{
696
782
  Bigint *b;
697
783
  int i, k;
699
785
 
700
786
  x= (nd + 8) / 9;
701
787
  for (k= 0, y= 1; x > y; y <<= 1, k++) ;
702
 
  b= Balloc(k);
 
788
  b= Balloc(k, alloc);
703
789
  b->p.x[0]= y9;
704
790
  b->wds= 1;
705
 
 
 
791
  
706
792
  i= 9;
707
793
  if (9 < nd0)
708
794
  {
709
795
    s+= 9;
710
796
    do
711
 
      b= multadd(b, 10, *s++ - '0');
 
797
      b= multadd(b, 10, *s++ - '0', alloc);
712
798
    while (++i < nd0);
713
799
    s++;
714
800
  }
715
801
  else
716
802
    s+= 10;
717
803
  for(; i < nd; i++)
718
 
    b= multadd(b, 10, *s++ - '0');
 
804
    b= multadd(b, 10, *s++ - '0', alloc);
719
805
  return b;
720
806
}
721
807
 
722
808
 
723
 
static int hi0bits(ULong x)
 
809
static int hi0bits(register ULong x)
724
810
{
725
 
  int k= 0;
 
811
  register int k= 0;
726
812
 
727
813
  if (!(x & 0xffff0000))
728
814
  {
756
842
 
757
843
static int lo0bits(ULong *y)
758
844
{
759
 
  int k;
760
 
  ULong x= *y;
 
845
  register int k;
 
846
  register ULong x= *y;
761
847
 
762
848
  if (x & 7)
763
849
  {
806
892
 
807
893
/* Convert integer to Bigint number */
808
894
 
809
 
static Bigint *i2b(int i)
 
895
static Bigint *i2b(int i, Stack_alloc *alloc)
810
896
{
811
897
  Bigint *b;
812
898
 
813
 
  b= Balloc(1);
 
899
  b= Balloc(1, alloc);
814
900
  b->p.x[0]= i;
815
901
  b->wds= 1;
816
902
  return b;
819
905
 
820
906
/* Multiply two Bigint numbers */
821
907
 
822
 
static Bigint *mult(Bigint *a, Bigint *b)
 
908
static Bigint *mult(Bigint *a, Bigint *b, Stack_alloc *alloc)
823
909
{
824
910
  Bigint *c;
825
911
  int k, wa, wb, wc;
839
925
  wc= wa + wb;
840
926
  if (wc > a->maxwds)
841
927
    k++;
842
 
  c= Balloc(k);
 
928
  c= Balloc(k, alloc);
843
929
  for (x= c->p.x, xa= x + wc; x < xa; x++)
844
930
    *x= 0;
845
931
  xa= a->p.x;
911
997
 
912
998
#define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1)
913
999
 
914
 
static Bigint *pow5mult(Bigint *b, int k)
 
1000
static Bigint *pow5mult(Bigint *b, int k, Stack_alloc *alloc)
915
1001
{
916
1002
  Bigint *b1, *p5, *p51;
917
1003
  int i;
918
1004
  static int p05[3]= { 5, 25, 125 };
919
1005
 
920
1006
  if ((i= k & 3))
921
 
    b= multadd(b, p05[i-1], 0);
 
1007
    b= multadd(b, p05[i-1], 0, alloc);
922
1008
 
923
1009
  if (!(k>>= 2))
924
1010
    return b;
927
1013
  {
928
1014
    if (k & 1)
929
1015
    {
930
 
      b1= mult(b, p5);
931
 
      Bfree(b);
 
1016
      b1= mult(b, p5, alloc);
 
1017
      Bfree(b, alloc);
932
1018
      b= b1;
933
1019
    }
934
1020
    if (!(k>>= 1))
937
1023
    if (p5 < p5_a + P5A_MAX)
938
1024
      ++p5;
939
1025
    else if (p5 == p5_a + P5A_MAX)
940
 
      p5= mult(p5, p5);
 
1026
      p5= mult(p5, p5, alloc);
941
1027
    else
942
1028
    {
943
 
      p51= mult(p5, p5);
944
 
      Bfree(p5);
 
1029
      p51= mult(p5, p5, alloc);
 
1030
      Bfree(p5, alloc);
945
1031
      p5= p51;
946
1032
    }
947
1033
  }
949
1035
}
950
1036
 
951
1037
 
952
 
static Bigint *lshift(Bigint *b, int k)
 
1038
static Bigint *lshift(Bigint *b, int k, Stack_alloc *alloc)
953
1039
{
954
1040
  int i, k1, n, n1;
955
1041
  Bigint *b1;
960
1046
  n1= n + b->wds + 1;
961
1047
  for (i= b->maxwds; n1 > i; i<<= 1)
962
1048
    k1++;
963
 
  b1= Balloc(k1);
 
1049
  b1= Balloc(k1, alloc);
964
1050
  x1= b1->p.x;
965
1051
  for (i= 0; i < n; i++)
966
1052
    *x1++= 0;
984
1070
      *x1++= *x++;
985
1071
    while (x < xe);
986
1072
  b1->wds= n1 - 1;
987
 
  Bfree(b);
 
1073
  Bfree(b, alloc);
988
1074
  return b1;
989
1075
}
990
1076
 
1013
1099
}
1014
1100
 
1015
1101
 
1016
 
static Bigint *diff(Bigint *a, Bigint *b)
 
1102
static Bigint *diff(Bigint *a, Bigint *b, Stack_alloc *alloc)
1017
1103
{
1018
1104
  Bigint *c;
1019
1105
  int i, wa, wb;
1023
1109
  i= cmp(a,b);
1024
1110
  if (!i)
1025
1111
  {
1026
 
    c= Balloc(0);
 
1112
    c= Balloc(0, alloc);
1027
1113
    c->wds= 1;
1028
1114
    c->p.x[0]= 0;
1029
1115
    return c;
1037
1123
  }
1038
1124
  else
1039
1125
    i= 0;
1040
 
  c= Balloc(a->k);
 
1126
  c= Balloc(a->k, alloc);
1041
1127
  c->sign= i;
1042
1128
  wa= a->wds;
1043
1129
  xa= a->p.x;
1069
1155
 
1070
1156
static double ulp(double x)
1071
1157
{
1072
 
  Long L;
 
1158
  register Long L;
1073
1159
  double a;
1074
1160
 
1075
1161
  L= (word0(x) & Exp_mask) - (P - 1)*Exp_msk1;
1118
1204
}
1119
1205
 
1120
1206
 
1121
 
static Bigint *d2b(double d, int *e, int *bits)
 
1207
static Bigint *d2b(double d, int *e, int *bits, Stack_alloc *alloc)
1122
1208
{
1123
1209
  Bigint *b;
1124
1210
  int de, k;
1127
1213
#define d0 word0(d)
1128
1214
#define d1 word1(d)
1129
1215
 
1130
 
  b= Balloc(1);
 
1216
  b= Balloc(1, alloc);
1131
1217
  x= b->p.x;
1132
1218
 
1133
1219
  z= d0 & Frac_mask;
1199
1285
  9007199254740992.*9007199254740992.e-256 /* = 2^106 * 1e-53 */
1200
1286
};
1201
1287
/*
1202
 
  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow
 
1288
  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow 
1203
1289
  flag unnecessarily.  It leads to a song and dance at the end of strtod.
1204
1290
*/
1205
1291
#define Scale_Bit 0x10
1207
1293
 
1208
1294
/*
1209
1295
  strtod for IEEE--arithmetic machines.
1210
 
 
 
1296
 
1211
1297
  This strtod returns a nearest machine number to the input decimal
1212
1298
  string (or sets errno to EOVERFLOW). Ties are broken by the IEEE round-even
1213
1299
  rule.
1214
 
 
 
1300
 
1215
1301
  Inspired loosely by William D. Clinger's paper "How to Read Floating
1216
1302
  Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
1217
 
 
 
1303
 
1218
1304
  Modifications:
1219
 
 
 
1305
 
1220
1306
   1. We only require IEEE (not IEEE double-extended).
1221
1307
   2. We get by with floating-point arithmetic in a case that
1222
1308
     Clinger missed -- when we're computing d * 10^n
1234
1320
     for 0 <= k <= 22).
1235
1321
*/
1236
1322
 
1237
 
static double my_strtod_int(const char *s00, char **se, int *error)
 
1323
static double my_strtod_int(const char *s00, char **se, int *error, char *buf, size_t buf_size)
1238
1324
{
1239
1325
  int scale;
1240
1326
  int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1243
1329
  double aadj, aadj1, adj, rv, rv0;
1244
1330
  Long L;
1245
1331
  ULong y, z;
1246
 
  Bigint *bb= NULL, *bb1= NULL, *bd= NULL, *bd0= NULL, *bs= NULL, *delta= NULL;
 
1332
  Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1247
1333
#ifdef SET_INEXACT
1248
1334
  int inexact, oldinexact;
1249
1335
#endif
 
1336
#ifdef Honor_FLT_ROUNDS
 
1337
  int rounding;
 
1338
#endif
 
1339
  Stack_alloc alloc;
1250
1340
 
1251
1341
  c= 0;
1252
1342
  *error= 0;
1253
1343
 
 
1344
  alloc.begin= alloc.free= buf;
 
1345
  alloc.end= buf + buf_size;
 
1346
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
 
1347
 
1254
1348
  sign= nz0= nz= 0;
1255
1349
  dval(rv)= 0.;
1256
1350
  for (s= s00; s < end; s++)
1274
1368
 break2:
1275
1369
  if (s >= end)
1276
1370
    goto ret0;
1277
 
 
 
1371
  
1278
1372
  if (*s == '0')
1279
1373
  {
1280
1374
    nz0= 1;
1399
1493
  }
1400
1494
  bd0= 0;
1401
1495
  if (nd <= DBL_DIG
 
1496
#ifndef Honor_FLT_ROUNDS
 
1497
    && Flt_Rounds == 1
 
1498
#endif
1402
1499
      )
1403
1500
  {
1404
1501
    if (!e)
1407
1504
    {
1408
1505
      if (e <= Ten_pmax)
1409
1506
      {
 
1507
#ifdef Honor_FLT_ROUNDS
 
1508
        /* round correctly FLT_ROUNDS = 2 or 3 */
 
1509
        if (sign)
 
1510
        {
 
1511
          rv= -rv;
 
1512
          sign= 0;
 
1513
        }
 
1514
#endif
1410
1515
        /* rv = */ rounded_product(dval(rv), tens[e]);
1411
1516
        goto ret;
1412
1517
      }
1417
1522
          A fancier test would sometimes let us do
1418
1523
          this for larger i values.
1419
1524
        */
 
1525
#ifdef Honor_FLT_ROUNDS
 
1526
        /* round correctly FLT_ROUNDS = 2 or 3 */
 
1527
        if (sign)
 
1528
        {
 
1529
          rv= -rv;
 
1530
          sign= 0;
 
1531
        }
 
1532
#endif
1420
1533
        e-= i;
1421
1534
        dval(rv)*= tens[i];
1422
1535
        /* rv = */ rounded_product(dval(rv), tens[e]);
1426
1539
#ifndef Inaccurate_Divide
1427
1540
    else if (e >= -Ten_pmax)
1428
1541
    {
 
1542
#ifdef Honor_FLT_ROUNDS
 
1543
      /* round correctly FLT_ROUNDS = 2 or 3 */
 
1544
      if (sign)
 
1545
      {
 
1546
        rv= -rv;
 
1547
        sign= 0;
 
1548
      }
 
1549
#endif
1429
1550
      /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1430
1551
      goto ret;
1431
1552
    }
1439
1560
    oldinexact= get_inexact();
1440
1561
#endif
1441
1562
  scale= 0;
 
1563
#ifdef Honor_FLT_ROUNDS
 
1564
  if ((rounding= Flt_Rounds) >= 2)
 
1565
  {
 
1566
    if (sign)
 
1567
      rounding= rounding == 2 ? 0 : 2;
 
1568
    else
 
1569
      if (rounding != 2)
 
1570
        rounding= 0;
 
1571
  }
 
1572
#endif
1442
1573
 
1443
1574
  /* Get starting approximation = rv * 10**e1 */
1444
1575
 
1453
1584
 ovfl:
1454
1585
        *error= EOVERFLOW;
1455
1586
        /* Can't trust HUGE_VAL */
 
1587
#ifdef Honor_FLT_ROUNDS
 
1588
        switch (rounding)
 
1589
        {
 
1590
        case 0: /* toward 0 */
 
1591
        case 3: /* toward -infinity */
 
1592
          word0(rv)= Big0;
 
1593
          word1(rv)= Big1;
 
1594
          break;
 
1595
        default:
 
1596
          word0(rv)= Exp_mask;
 
1597
          word1(rv)= 0;
 
1598
        }
 
1599
#else /*Honor_FLT_ROUNDS*/
1456
1600
        word0(rv)= Exp_mask;
1457
1601
        word1(rv)= 0;
 
1602
#endif /*Honor_FLT_ROUNDS*/
1458
1603
#ifdef SET_INEXACT
1459
1604
        /* set overflow bit */
1460
1605
        dval(rv0)= 1e300;
1526
1671
 
1527
1672
  /* Put digits into bd: true value = bd * 10^e */
1528
1673
 
1529
 
  bd0= s2b(s0, nd0, nd, y);
 
1674
  bd0= s2b(s0, nd0, nd, y, &alloc);
1530
1675
 
1531
1676
  for(;;)
1532
1677
  {
1533
 
    bd= Balloc(bd0->k);
 
1678
    bd= Balloc(bd0->k, &alloc);
1534
1679
    Bcopy(bd, bd0);
1535
 
    bb= d2b(dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
1536
 
    bs= i2b(1);
 
1680
    bb= d2b(dval(rv), &bbe, &bbbits, &alloc);  /* rv = bb * 2^bbe */
 
1681
    bs= i2b(1, &alloc);
1537
1682
 
1538
1683
    if (e >= 0)
1539
1684
    {
1550
1695
    else
1551
1696
      bd2-= bbe;
1552
1697
    bs2= bb2;
 
1698
#ifdef Honor_FLT_ROUNDS
 
1699
    if (rounding != 1)
 
1700
      bs2++;
 
1701
#endif
1553
1702
    j= bbe - scale;
1554
1703
    i= j + bbbits - 1;  /* logb(rv) */
1555
1704
    if (i < Emin)  /* denormal */
1570
1719
    }
1571
1720
    if (bb5 > 0)
1572
1721
    {
1573
 
      bs= pow5mult(bs, bb5);
1574
 
      bb1= mult(bs, bb);
1575
 
      Bfree(bb);
 
1722
      bs= pow5mult(bs, bb5, &alloc);
 
1723
      bb1= mult(bs, bb, &alloc);
 
1724
      Bfree(bb, &alloc);
1576
1725
      bb= bb1;
1577
1726
    }
1578
1727
    if (bb2 > 0)
1579
 
      bb= lshift(bb, bb2);
 
1728
      bb= lshift(bb, bb2, &alloc);
1580
1729
    if (bd5 > 0)
1581
 
      bd= pow5mult(bd, bd5);
 
1730
      bd= pow5mult(bd, bd5, &alloc);
1582
1731
    if (bd2 > 0)
1583
 
      bd= lshift(bd, bd2);
 
1732
      bd= lshift(bd, bd2, &alloc);
1584
1733
    if (bs2 > 0)
1585
 
      bs= lshift(bs, bs2);
1586
 
    delta= diff(bb, bd);
 
1734
      bs= lshift(bs, bs2, &alloc);
 
1735
    delta= diff(bb, bd, &alloc);
1587
1736
    dsign= delta->sign;
1588
1737
    delta->sign= 0;
1589
1738
    i= cmp(delta, bs);
 
1739
#ifdef Honor_FLT_ROUNDS
 
1740
    if (rounding != 1)
 
1741
    {
 
1742
      if (i < 0)
 
1743
      {
 
1744
        /* Error is less than an ulp */
 
1745
        if (!delta->x[0] && delta->wds <= 1)
 
1746
        {
 
1747
          /* exact */
 
1748
#ifdef SET_INEXACT
 
1749
          inexact= 0;
 
1750
#endif
 
1751
          break;
 
1752
        }
 
1753
        if (rounding)
 
1754
        {
 
1755
          if (dsign)
 
1756
          {
 
1757
            adj= 1.;
 
1758
            goto apply_adj;
 
1759
          }
 
1760
        }
 
1761
        else if (!dsign)
 
1762
        {
 
1763
          adj= -1.;
 
1764
          if (!word1(rv) && !(word0(rv) & Frac_mask))
 
1765
          {
 
1766
            y= word0(rv) & Exp_mask;
 
1767
            if (!scale || y > 2*P*Exp_msk1)
 
1768
            {
 
1769
              delta= lshift(delta,Log2P);
 
1770
              if (cmp(delta, bs) <= 0)
 
1771
              adj= -0.5;
 
1772
            }
 
1773
          }
 
1774
 apply_adj:
 
1775
          if (scale && (y= word0(rv) & Exp_mask) <= 2 * P * Exp_msk1)
 
1776
            word0(adj)+= (2 * P + 1) * Exp_msk1 - y;
 
1777
          dval(rv)+= adj * ulp(dval(rv));
 
1778
        }
 
1779
        break;
 
1780
      }
 
1781
      adj= ratio(delta, bs);
 
1782
      if (adj < 1.)
 
1783
        adj= 1.;
 
1784
      if (adj <= 0x7ffffffe)
 
1785
      {
 
1786
        /* adj = rounding ? ceil(adj) : floor(adj); */
 
1787
        y= adj;
 
1788
        if (y != adj)
 
1789
        {
 
1790
          if (!((rounding >> 1) ^ dsign))
 
1791
            y++;
 
1792
          adj= y;
 
1793
        }
 
1794
      }
 
1795
      if (scale && (y= word0(rv) & Exp_mask) <= 2 * P * Exp_msk1)
 
1796
        word0(adj)+= (2 * P + 1) * Exp_msk1 - y;
 
1797
      adj*= ulp(dval(rv));
 
1798
      if (dsign)
 
1799
        dval(rv)+= adj;
 
1800
      else
 
1801
        dval(rv)-= adj;
 
1802
      goto cont;
 
1803
    }
 
1804
#endif /*Honor_FLT_ROUNDS*/
1590
1805
 
1591
1806
    if (i < 0)
1592
1807
    {
1611
1826
#endif
1612
1827
        break;
1613
1828
      }
1614
 
      delta= lshift(delta, Log2P);
 
1829
      delta= lshift(delta, Log2P, &alloc);
1615
1830
      if (cmp(delta, bs) > 0)
1616
1831
        goto drop_down;
1617
1832
      break;
1693
1908
    {
1694
1909
      aadj*= 0.5;
1695
1910
      aadj1= dsign ? aadj : -aadj;
 
1911
#ifdef Check_FLT_ROUNDS
 
1912
      switch (Rounding)
 
1913
      {
 
1914
      case 2: /* towards +infinity */
 
1915
        aadj1-= 0.5;
 
1916
        break;
 
1917
      case 0: /* towards 0 */
 
1918
      case 3: /* towards -infinity */
 
1919
        aadj1+= 0.5;
 
1920
      }
 
1921
#else
1696
1922
      if (Flt_Rounds == 0)
1697
1923
        aadj1+= 0.5;
 
1924
#endif /*Check_FLT_ROUNDS*/
1698
1925
    }
1699
1926
    y= word0(rv) & Exp_mask;
1700
1927
 
1752
1979
      }
1753
1980
#endif
1754
1981
 cont:
1755
 
    Bfree(bb);
1756
 
    Bfree(bd);
1757
 
    Bfree(bs);
1758
 
    Bfree(delta);
 
1982
    Bfree(bb, &alloc);
 
1983
    Bfree(bd, &alloc);
 
1984
    Bfree(bs, &alloc);
 
1985
    Bfree(delta, &alloc);
1759
1986
  }
1760
1987
#ifdef SET_INEXACT
1761
1988
  if (inexact)
1785
2012
  }
1786
2013
#endif
1787
2014
 retfree:
1788
 
  Bfree(bb);
1789
 
  Bfree(bd);
1790
 
  Bfree(bs);
1791
 
  Bfree(bd0);
1792
 
  Bfree(delta);
 
2015
  Bfree(bb, &alloc);
 
2016
  Bfree(bd, &alloc);
 
2017
  Bfree(bs, &alloc);
 
2018
  Bfree(bd0, &alloc);
 
2019
  Bfree(delta, &alloc);
1793
2020
 ret:
1794
2021
  *se= (char *)s;
1795
2022
  return sign ? -dval(rv) : dval(rv);
1896
2123
 */
1897
2124
 
1898
2125
static char *dtoa(double d, int mode, int ndigits, int *decpt, int *sign,
1899
 
                  char **rve)
 
2126
                  char **rve, char *buf, size_t buf_size)
1900
2127
{
1901
2128
  /*
1902
2129
    Arguments ndigits, decpt, sign are similar to those
1908
2135
    mode:
1909
2136
          0 ==> shortest string that yields d when read in
1910
2137
                and rounded to nearest.
1911
 
 
1912
2138
          1 ==> like 0, but with Steele & White stopping rule;
1913
2139
                e.g. with IEEE P754 arithmetic , mode 0 gives
1914
2140
                1e23 whereas mode 1 gives 9.999999999999999e22.
1915
 
          2 ==> cmax(1,ndigits) significant digits.  This gives a
 
2141
          2 ==> max(1,ndigits) significant digits.  This gives a
1916
2142
                return value similar to that of ecvt, except
1917
2143
                that trailing zeros are suppressed.
1918
2144
          3 ==> through ndigits past the decimal point.  This
1922
2148
          4,5 ==> similar to 2 and 3, respectively, but (in
1923
2149
                round-nearest mode) with the tests of mode 0 to
1924
2150
                possibly return a shorter string that rounds to d.
 
2151
                With IEEE arithmetic and compilation with
 
2152
                -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
 
2153
                as modes 2 and 3 when FLT_ROUNDS != 1.
1925
2154
          6-9 ==> Debugging modes similar to mode - 4:  don't try
1926
2155
                fast floating-point estimate (if applicable).
1927
2156
 
1937
2166
  Long L;
1938
2167
  int denorm;
1939
2168
  ULong x;
1940
 
  Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
 
2169
  Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1941
2170
  double d2, ds, eps;
1942
2171
  char *s, *s0;
 
2172
#ifdef Honor_FLT_ROUNDS
 
2173
  int rounding;
 
2174
#endif
 
2175
  Stack_alloc alloc;
 
2176
  
 
2177
  alloc.begin= alloc.free= buf;
 
2178
  alloc.end= buf + buf_size;
 
2179
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
1943
2180
 
1944
2181
  if (word0(d) & Sign_bit)
1945
2182
  {
1951
2188
    *sign= 0;
1952
2189
 
1953
2190
  /* If infinity, set decpt to DTOA_OVERFLOW, if 0 set it to 1 */
1954
 
  if ((((word0(d) & Exp_mask) == Exp_mask) && ((*decpt= DTOA_OVERFLOW) != 0)) ||
1955
 
      (!dval(d) && ((*decpt= 1) != 0)))
 
2191
  if (((word0(d) & Exp_mask) == Exp_mask && (*decpt= DTOA_OVERFLOW)) ||
 
2192
      (!dval(d) && (*decpt= 1)))
1956
2193
  {
1957
2194
    /* Infinity, NaN, 0 */
1958
 
    char *res= (char*) malloc(2);
 
2195
    char *res= (char*) dtoa_alloc(2, &alloc);
1959
2196
    res[0]= '0';
1960
2197
    res[1]= '\0';
1961
2198
    if (rve)
1962
2199
      *rve= res + 1;
1963
2200
    return res;
1964
2201
  }
1965
 
 
1966
 
 
1967
 
  b= d2b(dval(d), &be, &bbits);
 
2202
  
 
2203
#ifdef Honor_FLT_ROUNDS
 
2204
  if ((rounding= Flt_Rounds) >= 2)
 
2205
  {
 
2206
    if (*sign)
 
2207
      rounding= rounding == 2 ? 0 : 2;
 
2208
    else
 
2209
      if (rounding != 2)
 
2210
        rounding= 0;
 
2211
  }
 
2212
#endif
 
2213
 
 
2214
  b= d2b(dval(d), &be, &bbits, &alloc);
1968
2215
  if ((i= (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
1969
2216
  {
1970
2217
    dval(d2)= dval(d);
1976
2223
      log10(x)      =  log(x) / log(10)
1977
2224
                   ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1978
2225
      log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
1979
 
 
 
2226
     
1980
2227
      This suggests computing an approximation k to log10(d) by
1981
 
 
 
2228
     
1982
2229
      k= (i - Bias)*0.301029995663981
1983
2230
           + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1984
 
 
 
2231
     
1985
2232
      We want k to be too large rather than too small.
1986
2233
      The error in the first-order Taylor series approximation
1987
2234
      is in our favor, so we just round up the constant enough
2046
2293
  if (mode < 0 || mode > 9)
2047
2294
    mode= 0;
2048
2295
 
 
2296
#ifdef Check_FLT_ROUNDS
 
2297
  try_quick= Rounding == 1;
 
2298
#else
2049
2299
  try_quick= 1;
 
2300
#endif
2050
2301
 
2051
2302
  if (mode > 5)
2052
2303
  {
2079
2330
    if (i <= 0)
2080
2331
      i= 1;
2081
2332
  }
2082
 
  s= s0= (char*)malloc(i+1); /* +1 for trailing '\0' appended
2083
 
                                at end of function */
 
2333
  s= s0= dtoa_alloc(i, &alloc);
 
2334
 
 
2335
#ifdef Honor_FLT_ROUNDS
 
2336
  if (mode > 1 && rounding != 1)
 
2337
    leftright= 0;
 
2338
#endif
2084
2339
 
2085
2340
  if (ilim >= 0 && ilim <= Quick_max && try_quick)
2086
2341
  {
2179
2434
            goto bump_up;
2180
2435
          else if (dval(d) < 0.5 - dval(eps))
2181
2436
          {
2182
 
            while (*--s == '0') {}
 
2437
            while (*--s == '0');
2183
2438
            s++;
2184
2439
            goto ret1;
2185
2440
          }
2211
2466
    {
2212
2467
      L= (Long)(dval(d) / ds);
2213
2468
      dval(d)-= L*ds;
 
2469
#ifdef Check_FLT_ROUNDS
 
2470
      /* If FLT_ROUNDS == 2, L will usually be high by 1 */
 
2471
      if (dval(d) < 0)
 
2472
      {
 
2473
        L--;
 
2474
        dval(d)+= ds;
 
2475
      }
 
2476
#endif
2214
2477
      *s++= '0' + (int)L;
2215
2478
      if (!dval(d))
2216
2479
      {
2218
2481
      }
2219
2482
      if (i == ilim)
2220
2483
      {
 
2484
#ifdef Honor_FLT_ROUNDS
 
2485
        if (mode > 1)
 
2486
        {
 
2487
          switch (rounding) {
 
2488
          case 0: goto ret1;
 
2489
          case 2: goto bump_up;
 
2490
          }
 
2491
        }
 
2492
#endif
2221
2493
        dval(d)+= dval(d);
2222
2494
        if (dval(d) > ds || (dval(d) == ds && L & 1))
2223
2495
        {
2245
2517
    i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
2246
2518
    b2+= i;
2247
2519
    s2+= i;
2248
 
    mhi= i2b(1);
 
2520
    mhi= i2b(1, &alloc);
2249
2521
  }
2250
2522
  if (m2 > 0 && s2 > 0)
2251
2523
  {
2260
2532
    {
2261
2533
      if (m5 > 0)
2262
2534
      {
2263
 
        mhi= pow5mult(mhi, m5);
2264
 
        b1= mult(mhi, b);
2265
 
        Bfree(b);
 
2535
        mhi= pow5mult(mhi, m5, &alloc);
 
2536
        b1= mult(mhi, b, &alloc);
 
2537
        Bfree(b, &alloc);
2266
2538
        b= b1;
2267
2539
      }
2268
2540
      if ((j= b5 - m5))
2269
 
        b= pow5mult(b, j);
 
2541
        b= pow5mult(b, j, &alloc);
2270
2542
    }
2271
2543
    else
2272
 
      b= pow5mult(b, b5);
 
2544
      b= pow5mult(b, b5, &alloc);
2273
2545
  }
2274
 
  S= i2b(1);
 
2546
  S= i2b(1, &alloc);
2275
2547
  if (s5 > 0)
2276
 
    S= pow5mult(S, s5);
 
2548
    S= pow5mult(S, s5, &alloc);
2277
2549
 
2278
2550
  /* Check for special case that d is a normalized power of 2. */
2279
2551
 
2280
2552
  spec_case= 0;
2281
2553
  if ((mode < 2 || leftright)
 
2554
#ifdef Honor_FLT_ROUNDS
 
2555
      && rounding == 1
 
2556
#endif
2282
2557
     )
2283
2558
  {
2284
2559
    if (!word1(d) && !(word0(d) & Bndry_mask) &&
2295
2570
  /*
2296
2571
    Arrange for convenient computation of quotients:
2297
2572
    shift left if necessary so divisor has 4 leading 0 bits.
2298
 
 
 
2573
    
2299
2574
    Perhaps we should just compute leading 28 bits of S once
2300
2575
    a nd for all and pass them and a shift to quorem, so it
2301
2576
    can do shifts and ors to compute the numerator for q.
2317
2592
    s2+= i;
2318
2593
  }
2319
2594
  if (b2 > 0)
2320
 
    b= lshift(b, b2);
 
2595
    b= lshift(b, b2, &alloc);
2321
2596
  if (s2 > 0)
2322
 
    S= lshift(S, s2);
 
2597
    S= lshift(S, s2, &alloc);
2323
2598
  if (k_check)
2324
2599
  {
2325
2600
    if (cmp(b,S) < 0)
2326
2601
    {
2327
2602
      k--;
2328
2603
      /* we botched the k estimate */
2329
 
      b= multadd(b, 10, 0);
 
2604
      b= multadd(b, 10, 0, &alloc);
2330
2605
      if (leftright)
2331
 
        mhi= multadd(mhi, 10, 0);
 
2606
        mhi= multadd(mhi, 10, 0, &alloc);
2332
2607
      ilim= ilim1;
2333
2608
    }
2334
2609
  }
2335
2610
  if (ilim <= 0 && (mode == 3 || mode == 5))
2336
2611
  {
2337
 
    if (ilim < 0 || cmp(b,S= multadd(S,5,0)) <= 0)
 
2612
    if (ilim < 0 || cmp(b,S= multadd(S,5,0, &alloc)) <= 0)
2338
2613
    {
2339
2614
      /* no digits, fcvt style */
2340
2615
no_digits:
2349
2624
  if (leftright)
2350
2625
  {
2351
2626
    if (m2 > 0)
2352
 
      mhi= lshift(mhi, m2);
 
2627
      mhi= lshift(mhi, m2, &alloc);
2353
2628
 
2354
2629
    /*
2355
2630
      Compute mlo -- check for special case that d is a normalized power of 2.
2358
2633
    mlo= mhi;
2359
2634
    if (spec_case)
2360
2635
    {
2361
 
      mhi= Balloc(mhi->k);
 
2636
      mhi= Balloc(mhi->k, &alloc);
2362
2637
      Bcopy(mhi, mlo);
2363
 
      mhi= lshift(mhi, Log2P);
 
2638
      mhi= lshift(mhi, Log2P, &alloc);
2364
2639
    }
2365
2640
 
2366
2641
    for (i= 1;;i++)
2368
2643
      dig= quorem(b,S) + '0';
2369
2644
      /* Do we yet have the shortest decimal string that will round to d? */
2370
2645
      j= cmp(b, mlo);
2371
 
      delta= diff(S, mhi);
 
2646
      delta= diff(S, mhi, &alloc);
2372
2647
      j1= delta->sign ? 1 : cmp(b, delta);
2373
 
      Bfree(delta);
 
2648
      Bfree(delta, &alloc);
2374
2649
      if (j1 == 0 && mode != 1 && !(word1(d) & 1)
 
2650
#ifdef Honor_FLT_ROUNDS
 
2651
          && rounding >= 1
 
2652
#endif
2375
2653
         )
2376
2654
      {
2377
2655
        if (dig == '9')
2387
2665
        {
2388
2666
          goto accept_dig;
2389
2667
        }
 
2668
#ifdef Honor_FLT_ROUNDS
 
2669
        if (mode > 1)
 
2670
          switch (rounding) {
 
2671
          case 0: goto accept_dig;
 
2672
          case 2: goto keep_dig;
 
2673
          }
 
2674
#endif /*Honor_FLT_ROUNDS*/
2390
2675
        if (j1 > 0)
2391
2676
        {
2392
 
          b= lshift(b, 1);
 
2677
          b= lshift(b, 1, &alloc);
2393
2678
          j1= cmp(b, S);
2394
2679
          if ((j1 > 0 || (j1 == 0 && dig & 1))
2395
2680
              && dig++ == '9')
2401
2686
      }
2402
2687
      if (j1 > 0)
2403
2688
      {
 
2689
#ifdef Honor_FLT_ROUNDS
 
2690
        if (!rounding)
 
2691
          goto accept_dig;
 
2692
#endif
2404
2693
        if (dig == '9')
2405
2694
        { /* possible if i == 1 */
2406
2695
round_9_up:
2410
2699
        *s++= dig + 1;
2411
2700
        goto ret;
2412
2701
      }
 
2702
#ifdef Honor_FLT_ROUNDS
 
2703
keep_dig:
 
2704
#endif
2413
2705
      *s++= dig;
2414
2706
      if (i == ilim)
2415
2707
        break;
2416
 
      b= multadd(b, 10, 0);
 
2708
      b= multadd(b, 10, 0, &alloc);
2417
2709
      if (mlo == mhi)
2418
 
        mlo= mhi= multadd(mhi, 10, 0);
 
2710
        mlo= mhi= multadd(mhi, 10, 0, &alloc);
2419
2711
      else
2420
2712
      {
2421
 
        mlo= multadd(mlo, 10, 0);
2422
 
        mhi= multadd(mhi, 10, 0);
 
2713
        mlo= multadd(mlo, 10, 0, &alloc);
 
2714
        mhi= multadd(mhi, 10, 0, &alloc);
2423
2715
      }
2424
2716
    }
2425
2717
  }
2433
2725
      }
2434
2726
      if (i >= ilim)
2435
2727
        break;
2436
 
      b= multadd(b, 10, 0);
 
2728
      b= multadd(b, 10, 0, &alloc);
2437
2729
    }
2438
2730
 
2439
2731
  /* Round off last digit */
2440
2732
 
2441
 
  b= lshift(b, 1);
 
2733
#ifdef Honor_FLT_ROUNDS
 
2734
  switch (rounding) {
 
2735
  case 0: goto trimzeros;
 
2736
  case 2: goto roundoff;
 
2737
  }
 
2738
#endif
 
2739
  b= lshift(b, 1, &alloc);
2442
2740
  j= cmp(b, S);
2443
2741
  if (j > 0 || (j == 0 && dig & 1))
2444
2742
  {
2454
2752
  }
2455
2753
  else
2456
2754
  {
2457
 
    while (*--s == '0') {}
 
2755
#ifdef Honor_FLT_ROUNDS
 
2756
trimzeros:
 
2757
#endif
 
2758
    while (*--s == '0');
2458
2759
    s++;
2459
2760
  }
2460
2761
ret:
2461
 
  Bfree(S);
 
2762
  Bfree(S, &alloc);
2462
2763
  if (mhi)
2463
2764
  {
2464
2765
    if (mlo && mlo != mhi)
2465
 
      Bfree(mlo);
2466
 
    Bfree(mhi);
 
2766
      Bfree(mlo, &alloc);
 
2767
    Bfree(mhi, &alloc);
2467
2768
  }
2468
2769
ret1:
2469
 
  Bfree(b);
 
2770
  Bfree(b, &alloc);
2470
2771
  *s= 0;
2471
2772
  *decpt= k + 1;
2472
2773
  if (rve)
2473
2774
    *rve= s;
2474
2775
  return s0;
2475
2776
}
2476
 
 
2477
 
} /* namespace internal */
2478
 
} /* namespace drizzled */