← Back to branch summary

~drizzle-trunk/drizzle/development

« back to all changes in this revision

Viewing changes to strings/dtoa.c

Renamed more stuff to drizzle.

Show diffs side-by-side

added added

removed removed

Lines of Context:
strings/dtoa.c
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
83
79
                      (including the terminating '\0').
84
80
   @param error       if not NULL, points to a location where the status of
85
81
                      conversion is stored upon return.
86
 
                      false  successful conversion
87
 
                      true   the input number is [-,+]infinity or nan.
 
82
                      FALSE  successful conversion
 
83
                      TRUE   the input number is [-,+]infinity or nan.
88
84
                             The output string in this case is always '0'.
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;
96
 
  assert(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
97
 
 
98
 
  res= dtoa(x, 5, precision, &decpt, &sign, &end);
 
92
  char buf[DTOA_BUFF_SIZE];
 
93
  DBUG_ASSERT(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
 
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)
106
 
      *error= true;
 
103
      *error= TRUE;
107
104
    return 1;
108
105
  }
109
106
 
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
 
    *error= false;
 
141
    *error= FALSE;
145
142
 
146
 
  free(res);
 
143
  dtoa_free(res, buf, sizeof(buf));
147
144
 
148
145
  return dst - to;
149
146
}
184
181
                      'width + 1' bytes.
185
182
   @param error       if not NULL, points to a location where the status of
186
183
                      conversion is stored upon return.
187
 
                      false  successful conversion
188
 
                      true   the input number is [-,+]infinity or nan.
 
184
                      FALSE  successful conversion
 
185
                      TRUE   the input number is [-,+]infinity or nan.
189
186
                             The output string in this case is always '0'.
190
187
   @return            number of written characters (excluding terminating '\0')
191
188
 
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;
220
 
  assert(width > 0 && to != NULL);
221
 
 
 
216
  char buf[DTOA_BUFF_SIZE];
 
217
  my_bool have_space, force_e_format;
 
218
  DBUG_ASSERT(width > 0 && to != NULL);
 
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)
235
 
      *error= true;
 
232
      *error= TRUE;
236
233
    return 1;
237
234
  }
238
235
 
239
236
  if (error != NULL)
240
 
    *error= false;
 
237
    *error= FALSE;
241
238
 
242
239
  src= res;
243
240
  len= end - res;
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.
320
317
      if (width < decpt)
321
318
      {
322
319
        if (error != NULL)
323
 
          *error= true;
 
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
387
384
    {
388
385
      /* Overflow */
389
386
      if (error != NULL)
390
 
        *error= true;
 
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
 
  assert(str != NULL && end != NULL && *end != NULL && error != NULL);
 
464
  DBUG_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
 
532
530
  file.
533
531
*/
534
532
 
535
 
typedef int32_t Long;
536
 
typedef uint32_t ULong;
537
 
typedef int64_t LLong;
538
 
typedef uint64_t ULLong;
 
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
 
 
542
typedef int32 Long;
 
543
typedef uint32 ULong;
 
544
typedef int64 LLong;
 
545
typedef uint64 ULLong;
539
546
 
540
547
typedef union { double d; ULong L[2]; } U;
541
548
 
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
 
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;
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
 
  c= 0;
1252
1341
  *error= 0;
1253
1342
 
 
1343
  alloc.begin= alloc.free= buf;
 
1344
  alloc.end= buf + buf_size;
 
1345
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
 
1346
 
1254
1347
  sign= nz0= nz= 0;
1255
1348
  dval(rv)= 0.;
1256
1349
  for (s= s00; s < end; s++)
1274
1367
 break2:
1275
1368
  if (s >= end)
1276
1369
    goto ret0;
1277
 
 
 
1370
  
1278
1371
  if (*s == '0')
1279
1372
  {
1280
1373
    nz0= 1;
1399
1492
  }
1400
1493
  bd0= 0;
1401
1494
  if (nd <= DBL_DIG
 
1495
#ifndef Honor_FLT_ROUNDS
 
1496
    && Flt_Rounds == 1
 
1497
#endif
1402
1498
      )
1403
1499
  {
1404
1500
    if (!e)
1407
1503
    {
1408
1504
      if (e <= Ten_pmax)
1409
1505
      {
 
1506
#ifdef Honor_FLT_ROUNDS
 
1507
        /* round correctly FLT_ROUNDS = 2 or 3 */
 
1508
        if (sign)
 
1509
        {
 
1510
          rv= -rv;
 
1511
          sign= 0;
 
1512
        }
 
1513
#endif
1410
1514
        /* rv = */ rounded_product(dval(rv), tens[e]);
1411
1515
        goto ret;
1412
1516
      }
1417
1521
          A fancier test would sometimes let us do
1418
1522
          this for larger i values.
1419
1523
        */
 
1524
#ifdef Honor_FLT_ROUNDS
 
1525
        /* round correctly FLT_ROUNDS = 2 or 3 */
 
1526
        if (sign)
 
1527
        {
 
1528
          rv= -rv;
 
1529
          sign= 0;
 
1530
        }
 
1531
#endif
1420
1532
        e-= i;
1421
1533
        dval(rv)*= tens[i];
1422
1534
        /* rv = */ rounded_product(dval(rv), tens[e]);
1426
1538
#ifndef Inaccurate_Divide
1427
1539
    else if (e >= -Ten_pmax)
1428
1540
    {
 
1541
#ifdef Honor_FLT_ROUNDS
 
1542
      /* round correctly FLT_ROUNDS = 2 or 3 */
 
1543
      if (sign)
 
1544
      {
 
1545
        rv= -rv;
 
1546
        sign= 0;
 
1547
      }
 
1548
#endif
1429
1549
      /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1430
1550
      goto ret;
1431
1551
    }
1439
1559
    oldinexact= get_inexact();
1440
1560
#endif
1441
1561
  scale= 0;
 
1562
#ifdef Honor_FLT_ROUNDS
 
1563
  if ((rounding= Flt_Rounds) >= 2)
 
1564
  {
 
1565
    if (sign)
 
1566
      rounding= rounding == 2 ? 0 : 2;
 
1567
    else
 
1568
      if (rounding != 2)
 
1569
        rounding= 0;
 
1570
  }
 
1571
#endif
1442
1572
 
1443
1573
  /* Get starting approximation = rv * 10**e1 */
1444
1574
 
1453
1583
 ovfl:
1454
1584
        *error= EOVERFLOW;
1455
1585
        /* Can't trust HUGE_VAL */
 
1586
#ifdef Honor_FLT_ROUNDS
 
1587
        switch (rounding)
 
1588
        {
 
1589
        case 0: /* toward 0 */
 
1590
        case 3: /* toward -infinity */
 
1591
          word0(rv)= Big0;
 
1592
          word1(rv)= Big1;
 
1593
          break;
 
1594
        default:
 
1595
          word0(rv)= Exp_mask;
 
1596
          word1(rv)= 0;
 
1597
        }
 
1598
#else /*Honor_FLT_ROUNDS*/
1456
1599
        word0(rv)= Exp_mask;
1457
1600
        word1(rv)= 0;
 
1601
#endif /*Honor_FLT_ROUNDS*/
1458
1602
#ifdef SET_INEXACT
1459
1603
        /* set overflow bit */
1460
1604
        dval(rv0)= 1e300;
1526
1670
 
1527
1671
  /* Put digits into bd: true value = bd * 10^e */
1528
1672
 
1529
 
  bd0= s2b(s0, nd0, nd, y);
 
1673
  bd0= s2b(s0, nd0, nd, y, &alloc);
1530
1674
 
1531
1675
  for(;;)
1532
1676
  {
1533
 
    bd= Balloc(bd0->k);
 
1677
    bd= Balloc(bd0->k, &alloc);
1534
1678
    Bcopy(bd, bd0);
1535
 
    bb= d2b(dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
1536
 
    bs= i2b(1);
 
1679
    bb= d2b(dval(rv), &bbe, &bbbits, &alloc);  /* rv = bb * 2^bbe */
 
1680
    bs= i2b(1, &alloc);
1537
1681
 
1538
1682
    if (e >= 0)
1539
1683
    {
1550
1694
    else
1551
1695
      bd2-= bbe;
1552
1696
    bs2= bb2;
 
1697
#ifdef Honor_FLT_ROUNDS
 
1698
    if (rounding != 1)
 
1699
      bs2++;
 
1700
#endif
1553
1701
    j= bbe - scale;
1554
1702
    i= j + bbbits - 1;  /* logb(rv) */
1555
1703
    if (i < Emin)  /* denormal */
1570
1718
    }
1571
1719
    if (bb5 > 0)
1572
1720
    {
1573
 
      bs= pow5mult(bs, bb5);
1574
 
      bb1= mult(bs, bb);
1575
 
      Bfree(bb);
 
1721
      bs= pow5mult(bs, bb5, &alloc);
 
1722
      bb1= mult(bs, bb, &alloc);
 
1723
      Bfree(bb, &alloc);
1576
1724
      bb= bb1;
1577
1725
    }
1578
1726
    if (bb2 > 0)
1579
 
      bb= lshift(bb, bb2);
 
1727
      bb= lshift(bb, bb2, &alloc);
1580
1728
    if (bd5 > 0)
1581
 
      bd= pow5mult(bd, bd5);
 
1729
      bd= pow5mult(bd, bd5, &alloc);
1582
1730
    if (bd2 > 0)
1583
 
      bd= lshift(bd, bd2);
 
1731
      bd= lshift(bd, bd2, &alloc);
1584
1732
    if (bs2 > 0)
1585
 
      bs= lshift(bs, bs2);
1586
 
    delta= diff(bb, bd);
 
1733
      bs= lshift(bs, bs2, &alloc);
 
1734
    delta= diff(bb, bd, &alloc);
1587
1735
    dsign= delta->sign;
1588
1736
    delta->sign= 0;
1589
1737
    i= cmp(delta, bs);
 
1738
#ifdef Honor_FLT_ROUNDS
 
1739
    if (rounding != 1)
 
1740
    {
 
1741
      if (i < 0)
 
1742
      {
 
1743
        /* Error is less than an ulp */
 
1744
        if (!delta->x[0] && delta->wds <= 1)
 
1745
        {
 
1746
          /* exact */
 
1747
#ifdef SET_INEXACT
 
1748
          inexact= 0;
 
1749
#endif
 
1750
          break;
 
1751
        }
 
1752
        if (rounding)
 
1753
        {
 
1754
          if (dsign)
 
1755
          {
 
1756
            adj= 1.;
 
1757
            goto apply_adj;
 
1758
          }
 
1759
        }
 
1760
        else if (!dsign)
 
1761
        {
 
1762
          adj= -1.;
 
1763
          if (!word1(rv) && !(word0(rv) & Frac_mask))
 
1764
          {
 
1765
            y= word0(rv) & Exp_mask;
 
1766
            if (!scale || y > 2*P*Exp_msk1)
 
1767
            {
 
1768
              delta= lshift(delta,Log2P);
 
1769
              if (cmp(delta, bs) <= 0)
 
1770
              adj= -0.5;
 
1771
            }
 
1772
          }
 
1773
 apply_adj:
 
1774
          if (scale && (y= word0(rv) & Exp_mask) <= 2 * P * Exp_msk1)
 
1775
            word0(adj)+= (2 * P + 1) * Exp_msk1 - y;
 
1776
          dval(rv)+= adj * ulp(dval(rv));
 
1777
        }
 
1778
        break;
 
1779
      }
 
1780
      adj= ratio(delta, bs);
 
1781
      if (adj < 1.)
 
1782
        adj= 1.;
 
1783
      if (adj <= 0x7ffffffe)
 
1784
      {
 
1785
        /* adj = rounding ? ceil(adj) : floor(adj); */
 
1786
        y= adj;
 
1787
        if (y != adj)
 
1788
        {
 
1789
          if (!((rounding >> 1) ^ dsign))
 
1790
            y++;
 
1791
          adj= y;
 
1792
        }
 
1793
      }
 
1794
      if (scale && (y= word0(rv) & Exp_mask) <= 2 * P * Exp_msk1)
 
1795
        word0(adj)+= (2 * P + 1) * Exp_msk1 - y;
 
1796
      adj*= ulp(dval(rv));
 
1797
      if (dsign)
 
1798
        dval(rv)+= adj;
 
1799
      else
 
1800
        dval(rv)-= adj;
 
1801
      goto cont;
 
1802
    }
 
1803
#endif /*Honor_FLT_ROUNDS*/
1590
1804
 
1591
1805
    if (i < 0)
1592
1806
    {
1611
1825
#endif
1612
1826
        break;
1613
1827
      }
1614
 
      delta= lshift(delta, Log2P);
 
1828
      delta= lshift(delta, Log2P, &alloc);
1615
1829
      if (cmp(delta, bs) > 0)
1616
1830
        goto drop_down;
1617
1831
      break;
1693
1907
    {
1694
1908
      aadj*= 0.5;
1695
1909
      aadj1= dsign ? aadj : -aadj;
 
1910
#ifdef Check_FLT_ROUNDS
 
1911
      switch (Rounding)
 
1912
      {
 
1913
      case 2: /* towards +infinity */
 
1914
        aadj1-= 0.5;
 
1915
        break;
 
1916
      case 0: /* towards 0 */
 
1917
      case 3: /* towards -infinity */
 
1918
        aadj1+= 0.5;
 
1919
      }
 
1920
#else
1696
1921
      if (Flt_Rounds == 0)
1697
1922
        aadj1+= 0.5;
 
1923
#endif /*Check_FLT_ROUNDS*/
1698
1924
    }
1699
1925
    y= word0(rv) & Exp_mask;
1700
1926
 
1752
1978
      }
1753
1979
#endif
1754
1980
 cont:
1755
 
    Bfree(bb);
1756
 
    Bfree(bd);
1757
 
    Bfree(bs);
1758
 
    Bfree(delta);
 
1981
    Bfree(bb, &alloc);
 
1982
    Bfree(bd, &alloc);
 
1983
    Bfree(bs, &alloc);
 
1984
    Bfree(delta, &alloc);
1759
1985
  }
1760
1986
#ifdef SET_INEXACT
1761
1987
  if (inexact)
1785
2011
  }
1786
2012
#endif
1787
2013
 retfree:
1788
 
  Bfree(bb);
1789
 
  Bfree(bd);
1790
 
  Bfree(bs);
1791
 
  Bfree(bd0);
1792
 
  Bfree(delta);
 
2014
  Bfree(bb, &alloc);
 
2015
  Bfree(bd, &alloc);
 
2016
  Bfree(bs, &alloc);
 
2017
  Bfree(bd0, &alloc);
 
2018
  Bfree(delta, &alloc);
1793
2019
 ret:
1794
2020
  *se= (char *)s;
1795
2021
  return sign ? -dval(rv) : dval(rv);
1896
2122
 */
1897
2123
 
1898
2124
static char *dtoa(double d, int mode, int ndigits, int *decpt, int *sign,
1899
 
                  char **rve)
 
2125
                  char **rve, char *buf, size_t buf_size)
1900
2126
{
1901
2127
  /*
1902
2128
    Arguments ndigits, decpt, sign are similar to those
1908
2134
    mode:
1909
2135
          0 ==> shortest string that yields d when read in
1910
2136
                and rounded to nearest.
1911
 
 
1912
2137
          1 ==> like 0, but with Steele & White stopping rule;
1913
2138
                e.g. with IEEE P754 arithmetic , mode 0 gives
1914
2139
                1e23 whereas mode 1 gives 9.999999999999999e22.
1915
 
          2 ==> cmax(1,ndigits) significant digits.  This gives a
 
2140
          2 ==> max(1,ndigits) significant digits.  This gives a
1916
2141
                return value similar to that of ecvt, except
1917
2142
                that trailing zeros are suppressed.
1918
2143
          3 ==> through ndigits past the decimal point.  This
1922
2147
          4,5 ==> similar to 2 and 3, respectively, but (in
1923
2148
                round-nearest mode) with the tests of mode 0 to
1924
2149
                possibly return a shorter string that rounds to d.
 
2150
                With IEEE arithmetic and compilation with
 
2151
                -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
 
2152
                as modes 2 and 3 when FLT_ROUNDS != 1.
1925
2153
          6-9 ==> Debugging modes similar to mode - 4:  don't try
1926
2154
                fast floating-point estimate (if applicable).
1927
2155
 
1931
2159
    to hold the suppressed trailing zeros.
1932
2160
  */
1933
2161
 
1934
 
  int bbits, b2, b5, be, dig, i, ieps, ilim=0, ilim0, ilim1= 0,
 
2162
  int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1= 0,
1935
2163
    j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1936
2164
    spec_case, try_quick;
1937
2165
  Long L;
1938
2166
  int denorm;
1939
2167
  ULong x;
1940
 
  Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
 
2168
  Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1941
2169
  double d2, ds, eps;
1942
2170
  char *s, *s0;
 
2171
#ifdef Honor_FLT_ROUNDS
 
2172
  int rounding;
 
2173
#endif
 
2174
  Stack_alloc alloc;
 
2175
  
 
2176
  alloc.begin= alloc.free= buf;
 
2177
  alloc.end= buf + buf_size;
 
2178
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
1943
2179
 
1944
2180
  if (word0(d) & Sign_bit)
1945
2181
  {
1955
2191
      (!dval(d) && (*decpt= 1)))
1956
2192
  {
1957
2193
    /* Infinity, NaN, 0 */
1958
 
    char *res= (char*) malloc(2);
 
2194
    char *res= (char*) dtoa_alloc(2, &alloc);
1959
2195
    res[0]= '0';
1960
2196
    res[1]= '\0';
1961
2197
    if (rve)
1962
2198
      *rve= res + 1;
1963
2199
    return res;
1964
2200
  }
1965
 
 
1966
 
 
1967
 
  b= d2b(dval(d), &be, &bbits);
 
2201
  
 
2202
#ifdef Honor_FLT_ROUNDS
 
2203
  if ((rounding= Flt_Rounds) >= 2)
 
2204
  {
 
2205
    if (*sign)
 
2206
      rounding= rounding == 2 ? 0 : 2;
 
2207
    else
 
2208
      if (rounding != 2)
 
2209
        rounding= 0;
 
2210
  }
 
2211
#endif
 
2212
 
 
2213
  b= d2b(dval(d), &be, &bbits, &alloc);
1968
2214
  if ((i= (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
1969
2215
  {
1970
2216
    dval(d2)= dval(d);
1976
2222
      log10(x)      =  log(x) / log(10)
1977
2223
                   ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1978
2224
      log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
1979
 
 
 
2225
     
1980
2226
      This suggests computing an approximation k to log10(d) by
1981
 
 
 
2227
     
1982
2228
      k= (i - Bias)*0.301029995663981
1983
2229
           + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1984
 
 
 
2230
     
1985
2231
      We want k to be too large rather than too small.
1986
2232
      The error in the first-order Taylor series approximation
1987
2233
      is in our favor, so we just round up the constant enough
2046
2292
  if (mode < 0 || mode > 9)
2047
2293
    mode= 0;
2048
2294
 
 
2295
#ifdef Check_FLT_ROUNDS
 
2296
  try_quick= Rounding == 1;
 
2297
#else
2049
2298
  try_quick= 1;
 
2299
#endif
2050
2300
 
2051
2301
  if (mode > 5)
2052
2302
  {
2079
2329
    if (i <= 0)
2080
2330
      i= 1;
2081
2331
  }
2082
 
  s= s0= (char*)malloc(i+1); /* +1 for trailing '\0' appended
2083
 
                                at end of function */
 
2332
  s= s0= dtoa_alloc(i, &alloc);
 
2333
 
 
2334
#ifdef Honor_FLT_ROUNDS
 
2335
  if (mode > 1 && rounding != 1)
 
2336
    leftright= 0;
 
2337
#endif
2084
2338
 
2085
2339
  if (ilim >= 0 && ilim <= Quick_max && try_quick)
2086
2340
  {
2179
2433
            goto bump_up;
2180
2434
          else if (dval(d) < 0.5 - dval(eps))
2181
2435
          {
2182
 
            while (*--s == '0') {}
 
2436
            while (*--s == '0');
2183
2437
            s++;
2184
2438
            goto ret1;
2185
2439
          }
2211
2465
    {
2212
2466
      L= (Long)(dval(d) / ds);
2213
2467
      dval(d)-= L*ds;
 
2468
#ifdef Check_FLT_ROUNDS
 
2469
      /* If FLT_ROUNDS == 2, L will usually be high by 1 */
 
2470
      if (dval(d) < 0)
 
2471
      {
 
2472
        L--;
 
2473
        dval(d)+= ds;
 
2474
      }
 
2475
#endif
2214
2476
      *s++= '0' + (int)L;
2215
2477
      if (!dval(d))
2216
2478
      {
2218
2480
      }
2219
2481
      if (i == ilim)
2220
2482
      {
 
2483
#ifdef Honor_FLT_ROUNDS
 
2484
        if (mode > 1)
 
2485
        {
 
2486
          switch (rounding) {
 
2487
          case 0: goto ret1;
 
2488
          case 2: goto bump_up;
 
2489
          }
 
2490
        }
 
2491
#endif
2221
2492
        dval(d)+= dval(d);
2222
2493
        if (dval(d) > ds || (dval(d) == ds && L & 1))
2223
2494
        {
2245
2516
    i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
2246
2517
    b2+= i;
2247
2518
    s2+= i;
2248
 
    mhi= i2b(1);
 
2519
    mhi= i2b(1, &alloc);
2249
2520
  }
2250
2521
  if (m2 > 0 && s2 > 0)
2251
2522
  {
2260
2531
    {
2261
2532
      if (m5 > 0)
2262
2533
      {
2263
 
        mhi= pow5mult(mhi, m5);
2264
 
        b1= mult(mhi, b);
2265
 
        Bfree(b);
 
2534
        mhi= pow5mult(mhi, m5, &alloc);
 
2535
        b1= mult(mhi, b, &alloc);
 
2536
        Bfree(b, &alloc);
2266
2537
        b= b1;
2267
2538
      }
2268
2539
      if ((j= b5 - m5))
2269
 
        b= pow5mult(b, j);
 
2540
        b= pow5mult(b, j, &alloc);
2270
2541
    }
2271
2542
    else
2272
 
      b= pow5mult(b, b5);
 
2543
      b= pow5mult(b, b5, &alloc);
2273
2544
  }
2274
 
  S= i2b(1);
 
2545
  S= i2b(1, &alloc);
2275
2546
  if (s5 > 0)
2276
 
    S= pow5mult(S, s5);
 
2547
    S= pow5mult(S, s5, &alloc);
2277
2548
 
2278
2549
  /* Check for special case that d is a normalized power of 2. */
2279
2550
 
2280
2551
  spec_case= 0;
2281
2552
  if ((mode < 2 || leftright)
 
2553
#ifdef Honor_FLT_ROUNDS
 
2554
      && rounding == 1
 
2555
#endif
2282
2556
     )
2283
2557
  {
2284
2558
    if (!word1(d) && !(word0(d) & Bndry_mask) &&
2295
2569
  /*
2296
2570
    Arrange for convenient computation of quotients:
2297
2571
    shift left if necessary so divisor has 4 leading 0 bits.
2298
 
 
 
2572
    
2299
2573
    Perhaps we should just compute leading 28 bits of S once
2300
2574
    a nd for all and pass them and a shift to quorem, so it
2301
2575
    can do shifts and ors to compute the numerator for q.
2317
2591
    s2+= i;
2318
2592
  }
2319
2593
  if (b2 > 0)
2320
 
    b= lshift(b, b2);
 
2594
    b= lshift(b, b2, &alloc);
2321
2595
  if (s2 > 0)
2322
 
    S= lshift(S, s2);
 
2596
    S= lshift(S, s2, &alloc);
2323
2597
  if (k_check)
2324
2598
  {
2325
2599
    if (cmp(b,S) < 0)
2326
2600
    {
2327
2601
      k--;
2328
2602
      /* we botched the k estimate */
2329
 
      b= multadd(b, 10, 0);
 
2603
      b= multadd(b, 10, 0, &alloc);
2330
2604
      if (leftright)
2331
 
        mhi= multadd(mhi, 10, 0);
 
2605
        mhi= multadd(mhi, 10, 0, &alloc);
2332
2606
      ilim= ilim1;
2333
2607
    }
2334
2608
  }
2335
2609
  if (ilim <= 0 && (mode == 3 || mode == 5))
2336
2610
  {
2337
 
    if (ilim < 0 || cmp(b,S= multadd(S,5,0)) <= 0)
 
2611
    if (ilim < 0 || cmp(b,S= multadd(S,5,0, &alloc)) <= 0)
2338
2612
    {
2339
2613
      /* no digits, fcvt style */
2340
2614
no_digits:
2349
2623
  if (leftright)
2350
2624
  {
2351
2625
    if (m2 > 0)
2352
 
      mhi= lshift(mhi, m2);
 
2626
      mhi= lshift(mhi, m2, &alloc);
2353
2627
 
2354
2628
    /*
2355
2629
      Compute mlo -- check for special case that d is a normalized power of 2.
2358
2632
    mlo= mhi;
2359
2633
    if (spec_case)
2360
2634
    {
2361
 
      mhi= Balloc(mhi->k);
 
2635
      mhi= Balloc(mhi->k, &alloc);
2362
2636
      Bcopy(mhi, mlo);
2363
 
      mhi= lshift(mhi, Log2P);
 
2637
      mhi= lshift(mhi, Log2P, &alloc);
2364
2638
    }
2365
2639
 
2366
2640
    for (i= 1;;i++)
2368
2642
      dig= quorem(b,S) + '0';
2369
2643
      /* Do we yet have the shortest decimal string that will round to d? */
2370
2644
      j= cmp(b, mlo);
2371
 
      delta= diff(S, mhi);
 
2645
      delta= diff(S, mhi, &alloc);
2372
2646
      j1= delta->sign ? 1 : cmp(b, delta);
2373
 
      Bfree(delta);
 
2647
      Bfree(delta, &alloc);
2374
2648
      if (j1 == 0 && mode != 1 && !(word1(d) & 1)
 
2649
#ifdef Honor_FLT_ROUNDS
 
2650
          && rounding >= 1
 
2651
#endif
2375
2652
         )
2376
2653
      {
2377
2654
        if (dig == '9')
2387
2664
        {
2388
2665
          goto accept_dig;
2389
2666
        }
 
2667
#ifdef Honor_FLT_ROUNDS
 
2668
        if (mode > 1)
 
2669
          switch (rounding) {
 
2670
          case 0: goto accept_dig;
 
2671
          case 2: goto keep_dig;
 
2672
          }
 
2673
#endif /*Honor_FLT_ROUNDS*/
2390
2674
        if (j1 > 0)
2391
2675
        {
2392
 
          b= lshift(b, 1);
 
2676
          b= lshift(b, 1, &alloc);
2393
2677
          j1= cmp(b, S);
2394
2678
          if ((j1 > 0 || (j1 == 0 && dig & 1))
2395
2679
              && dig++ == '9')
2401
2685
      }
2402
2686
      if (j1 > 0)
2403
2687
      {
 
2688
#ifdef Honor_FLT_ROUNDS
 
2689
        if (!rounding)
 
2690
          goto accept_dig;
 
2691
#endif
2404
2692
        if (dig == '9')
2405
2693
        { /* possible if i == 1 */
2406
2694
round_9_up:
2410
2698
        *s++= dig + 1;
2411
2699
        goto ret;
2412
2700
      }
 
2701
#ifdef Honor_FLT_ROUNDS
 
2702
keep_dig:
 
2703
#endif
2413
2704
      *s++= dig;
2414
2705
      if (i == ilim)
2415
2706
        break;
2416
 
      b= multadd(b, 10, 0);
 
2707
      b= multadd(b, 10, 0, &alloc);
2417
2708
      if (mlo == mhi)
2418
 
        mlo= mhi= multadd(mhi, 10, 0);
 
2709
        mlo= mhi= multadd(mhi, 10, 0, &alloc);
2419
2710
      else
2420
2711
      {
2421
 
        mlo= multadd(mlo, 10, 0);
2422
 
        mhi= multadd(mhi, 10, 0);
 
2712
        mlo= multadd(mlo, 10, 0, &alloc);
 
2713
        mhi= multadd(mhi, 10, 0, &alloc);
2423
2714
      }
2424
2715
    }
2425
2716
  }
2433
2724
      }
2434
2725
      if (i >= ilim)
2435
2726
        break;
2436
 
      b= multadd(b, 10, 0);
 
2727
      b= multadd(b, 10, 0, &alloc);
2437
2728
    }
2438
2729
 
2439
2730
  /* Round off last digit */
2440
2731
 
2441
 
  b= lshift(b, 1);
 
2732
#ifdef Honor_FLT_ROUNDS
 
2733
  switch (rounding) {
 
2734
  case 0: goto trimzeros;
 
2735
  case 2: goto roundoff;
 
2736
  }
 
2737
#endif
 
2738
  b= lshift(b, 1, &alloc);
2442
2739
  j= cmp(b, S);
2443
2740
  if (j > 0 || (j == 0 && dig & 1))
2444
2741
  {
2454
2751
  }
2455
2752
  else
2456
2753
  {
2457
 
    while (*--s == '0') {}
 
2754
#ifdef Honor_FLT_ROUNDS
 
2755
trimzeros:
 
2756
#endif
 
2757
    while (*--s == '0');
2458
2758
    s++;
2459
2759
  }
2460
2760
ret:
2461
 
  Bfree(S);
 
2761
  Bfree(S, &alloc);
2462
2762
  if (mhi)
2463
2763
  {
2464
2764
    if (mlo && mlo != mhi)
2465
 
      Bfree(mlo);
2466
 
    Bfree(mhi);
 
2765
      Bfree(mlo, &alloc);
 
2766
    Bfree(mhi, &alloc);
2467
2767
  }
2468
2768
ret1:
2469
 
  Bfree(b);
 
2769
  Bfree(b, &alloc);
2470
2770
  *s= 0;
2471
2771
  *decpt= k + 1;
2472
2772
  if (rve)
2473
2773
    *rve= s;
2474
2774
  return s0;
2475
2775
}
2476
 
 
2477
 
} /* namespace internal */
2478
 
} /* namespace drizzled */