← Back to branch summary

~drizzle-trunk/drizzle/development

« back to all changes in this revision

Viewing changes to mystrings/dtoa.cc

  • Committer: Monty Taylor
  • Date: 2008-11-16 06:29:53 UTC
  • mto: (584.1.9 devel)
  • mto: This revision was merged to the branch mainline in revision 589.
  • Revision ID: monty@inaugust.com-20081116062953-ivdltjmfe009b5fr
Moved stuff into item/

Show diffs side-by-side

added added

removed removed

Lines of Context:
mystrings/dtoa.cc
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 <m_string.h>  /* for memcpy and NOT_FIXED_DEC */
 
39
 
 
40
/**
 
41
   Appears to suffice to not call malloc() in most cases.
 
42
   @todo
 
43
     see if it is possible to get rid of malloc().
 
44
     this constant is sufficient to avoid malloc() on all inputs I have tried.
 
45
*/
 
46
#define DTOA_BUFF_SIZE (420 * sizeof(void *))
54
47
 
55
48
/* Magic value returned by dtoa() to indicate overflow */
56
49
#define DTOA_OVERFLOW 9999
57
50
 
58
 
static double my_strtod_int(const char *, char **, int *);
59
 
static char *dtoa(double, int, int, int *, int *, char **);
 
51
static double my_strtod_int(const char *, char **, int *, char *, size_t);
 
52
static char *dtoa(double, int, int, int *, int *, char **, char *, size_t);
 
53
static void dtoa_free(char *, char *, size_t);
60
54
 
61
55
/**
62
56
   @brief
93
87
{
94
88
  int decpt, sign, len, i;
95
89
  char *res, *src, *end, *dst= to;
 
90
  char buf[DTOA_BUFF_SIZE];
96
91
  assert(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
97
 
 
98
 
  res= dtoa(x, 5, precision, &decpt, &sign, &end);
 
92
  
 
93
  res= dtoa(x, 5, precision, &decpt, &sign, &end, buf, sizeof(buf));
99
94
 
100
95
  if (decpt == DTOA_OVERFLOW)
101
96
  {
102
 
    free(res);
 
97
    dtoa_free(res, buf, sizeof(buf));
103
98
    *to++= '0';
104
99
    *to= '\0';
105
100
    if (error != NULL)
134
129
  {
135
130
    if (len <= decpt)
136
131
      *dst++= '.';
137
 
 
138
 
    for (i= precision - max(0, (len - decpt)); i > 0; i--)
 
132
    
 
133
    for (i= precision - cmax(0, (len - decpt)); i > 0; i--)
139
134
      *dst++= '0';
140
135
  }
141
 
 
 
136
  
142
137
  *dst= '\0';
143
138
  if (error != NULL)
144
139
    *error= false;
145
140
 
146
 
  free(res);
 
141
  dtoa_free(res, buf, sizeof(buf));
147
142
 
148
143
  return dst - to;
149
144
}
201
196
     my_gcvt(55, ..., 1, ...);
202
197
 
203
198
   We do our best to minimize such cases by:
204
 
 
 
199
   
205
200
   - passing to dtoa() the field width as the number of significant digits
206
 
 
 
201
   
207
202
   - removing the sign of the number early (and decreasing the width before
208
203
     passing it to dtoa())
209
 
 
 
204
   
210
205
   - choosing the proper format to preserve the most number of significant
211
206
     digits.
212
207
*/
216
211
{
217
212
  int decpt, sign, len, exp_len;
218
213
  char *res, *src, *end, *dst= to, *dend= dst + width;
 
214
  char buf[DTOA_BUFF_SIZE];
219
215
  bool have_space, force_e_format;
220
216
  assert(width > 0 && to != NULL);
221
 
 
 
217
  
222
218
  /* We want to remove '-' from equations early */
223
219
  if (x < 0.)
224
220
    width--;
225
221
 
226
 
  res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? min(width, DBL_DIG) : 
227
 
            min(width, FLT_DIG), &decpt, &sign, &end);
228
 
 
 
222
  res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? width : cmin(width, FLT_DIG),
 
223
            &decpt, &sign, &end, buf, sizeof(buf));
229
224
  if (decpt == DTOA_OVERFLOW)
230
225
  {
231
 
    free(res);
 
226
    dtoa_free(res, buf, sizeof(buf));
232
227
    *to++= '0';
233
228
    *to= '\0';
234
229
    if (error != NULL)
248
243
     to count it here.
249
244
   */
250
245
  exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
251
 
 
 
246
  
252
247
  /*
253
248
     Do we have enough space for all digits in the 'f' format?
254
249
     Let 'len' be the number of significant digits returned by dtoa,
301
296
       ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
302
297
                                            (len > 1 || !force_e_format)))) &&
303
298
         !force_e_format)) &&
304
 
 
 
299
      
305
300
       /*
306
301
         Use the 'e' format in some cases even if we have enough space for the
307
302
         'f' one. See comment for MAX_DECPT_FOR_F_FORMAT.
323
318
          *error= true;
324
319
        width= decpt;
325
320
      }
326
 
 
 
321
      
327
322
      /*
328
323
        We want to truncate (len - width) least significant digits after the
329
324
        decimal point. For this we are calling dtoa with mode=5, passing the
330
325
        number of significant digits = (len-decpt) - (len-width) = width-decpt
331
326
      */
332
 
      free(res);
333
 
      res= dtoa(x, 5, width - decpt, &decpt, &sign, &end);
 
327
      dtoa_free(res, buf, sizeof(buf));
 
328
      res= dtoa(x, 5, width - decpt, &decpt, &sign, &end, buf, sizeof(buf));
334
329
      src= res;
335
330
      len= end - res;
336
331
    }
341
336
      *dst++= '0';
342
337
      goto end;
343
338
    }
344
 
 
 
339
    
345
340
    /*
346
341
      At this point we are sure we have enough space to put all digits
347
342
      returned by dtoa
390
385
        *error= true;
391
386
      width= 0;
392
387
    }
393
 
 
 
388
      
394
389
    /* Do we have to truncate any digits? */
395
390
    if (width < len)
396
391
    {
397
392
      /* Yes, re-convert with a smaller width */
398
 
      free(res);
399
 
      res= dtoa(x, 4, width, &decpt, &sign, &end);
 
393
      dtoa_free(res, buf, sizeof(buf));
 
394
      res= dtoa(x, 4, width, &decpt, &sign, &end, buf, sizeof(buf));
400
395
      src= res;
401
396
      len= end - res;
402
397
      if (--decpt < 0)
436
431
  }
437
432
 
438
433
end:
439
 
  free(res);
 
434
  dtoa_free(res, buf, sizeof(buf));
440
435
  *dst= '\0';
441
436
 
442
437
  return dst - to;
455
450
                  rejected character.
456
451
   @param error   Upon return is set to EOVERFLOW in case of underflow or
457
452
                  overflow.
458
 
 
 
453
   
459
454
   @return        The resulting double value. In case of underflow, 0.0 is
460
455
                  returned. In case overflow, signed DBL_MAX is returned.
461
456
*/
462
457
 
463
458
double my_strtod(const char *str, char **end, int *error)
464
459
{
 
460
  char buf[DTOA_BUFF_SIZE];
465
461
  double res;
466
462
  assert(str != NULL && end != NULL && *end != NULL && error != NULL);
467
463
 
468
 
  res= my_strtod_int(str, end, error);
 
464
  res= my_strtod_int(str, end, error, buf, sizeof(buf));
469
465
  return (*error == 0) ? res : (res < 0 ? -DBL_MAX : DBL_MAX);
470
466
}
471
467
 
594
590
#define Big1 0xffffffff
595
591
#define FFFFFFFF 0xffffffffUL
596
592
 
 
593
/* This is tested to be enough for dtoa */
 
594
 
 
595
#define Kmax 15
 
596
 
 
597
#define Bcopy(x,y) memcpy(&x->sign, &y->sign, \
 
598
                          2*sizeof(int) + y->wds*sizeof(ULong))
 
599
 
597
600
/* Arbitrary-length integer */
598
601
 
599
602
typedef struct Bigint
608
611
  int wds;                 /* current length in 32-bit words */
609
612
} Bigint;
610
613
 
611
 
static Bigint *Bcopy(Bigint* dst, Bigint* src)
 
614
 
 
615
/* A simple stack-memory based allocator for Bigints */
 
616
 
 
617
typedef struct Stack_alloc
612
618
{
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)
 
619
  char *begin;
 
620
  char *free;
 
621
  char *end;
 
622
  /*
 
623
    Having list of free blocks lets us reduce maximum required amount
 
624
    of memory from ~4000 bytes to < 1680 (tested on x86).
 
625
  */
 
626
  Bigint *freelist[Kmax+1];
 
627
} Stack_alloc;
 
628
 
 
629
 
 
630
/*
 
631
  Try to allocate object on stack, and resort to malloc if all
 
632
  stack memory is used.
 
633
*/
 
634
 
 
635
static Bigint *Balloc(int k, Stack_alloc *alloc)
624
636
{
625
637
  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
 
 
 
638
  if (k <= Kmax &&  alloc->freelist[k])
 
639
  {
 
640
    rv= alloc->freelist[k];
 
641
    alloc->freelist[k]= rv->p.next;
 
642
  }
 
643
  else
 
644
  {
 
645
    int x, len;
 
646
 
 
647
    x= 1 << k;
 
648
    len= sizeof(Bigint) + x * sizeof(ULong);
 
649
 
 
650
    if (alloc->free + len <= alloc->end)
 
651
    {
 
652
      rv= (Bigint*) alloc->free;
 
653
      alloc->free+= len;
 
654
    }
 
655
    else
 
656
      rv= (Bigint*) malloc(len);
 
657
 
 
658
    rv->k= k;
 
659
    rv->maxwds= x;
 
660
  }
637
661
  rv->sign= rv->wds= 0;
638
 
 
 
662
  rv->p.x= (ULong*) (rv + 1);
639
663
  return rv;
640
664
}
641
665
 
645
669
  list. Otherwise call free().
646
670
*/
647
671
 
648
 
static void Bfree(Bigint *v)
649
 
{
650
 
  if(!v)
651
 
    return;
652
 
  free(v->p.x);
653
 
  free(v);
654
 
}
 
672
static void Bfree(Bigint *v, Stack_alloc *alloc)
 
673
{
 
674
  char *gptr= (char*) v;                       /* generic pointer */
 
675
  if (gptr < alloc->begin || gptr >= alloc->end)
 
676
    free(gptr);
 
677
  else if (v->k <= Kmax)
 
678
  {
 
679
    /*
 
680
      Maintain free lists only for stack objects: this way we don't
 
681
      have to bother with freeing lists in the end of dtoa;
 
682
      heap should not be used normally anyway.
 
683
    */
 
684
    v->p.next= alloc->freelist[v->k];
 
685
    alloc->freelist[v->k]= v;
 
686
  }
 
687
}
 
688
 
 
689
 
 
690
/*
 
691
  This is to place return value of dtoa in: tries to use stack
 
692
  as well, but passes by free lists management and just aligns len by
 
693
  sizeof(ULong).
 
694
*/
 
695
 
 
696
static char *dtoa_alloc(int i, Stack_alloc *alloc)
 
697
{
 
698
  char *rv;
 
699
  int aligned_size= (i + sizeof(ULong) - 1) / sizeof(ULong) * sizeof(ULong);
 
700
  if (alloc->free + aligned_size <= alloc->end)
 
701
  {
 
702
    rv= alloc->free;
 
703
    alloc->free+= aligned_size;
 
704
  }
 
705
  else
 
706
    rv= (char *)malloc(i);
 
707
  return rv;
 
708
}
 
709
 
 
710
 
 
711
/*
 
712
  dtoa_free() must be used to free values s returned by dtoa()
 
713
  This is the counterpart of dtoa_alloc()
 
714
*/
 
715
 
 
716
static void dtoa_free(char *gptr, char *buf, size_t buf_size)
 
717
{
 
718
  if (gptr < buf || gptr >= buf + buf_size)
 
719
    free(gptr);
 
720
}
 
721
 
655
722
 
656
723
/* Bigint arithmetic functions */
657
724
 
658
725
/* Multiply by m and add a */
659
726
 
660
 
static Bigint *multadd(Bigint *b, int m, int a)
 
727
static Bigint *multadd(Bigint *b, int m, int a, Stack_alloc *alloc)
661
728
{
662
729
  int i, wds;
663
730
  ULong *x;
679
746
  {
680
747
    if (wds >= b->maxwds)
681
748
    {
682
 
      b1= Balloc(b->k+1);
 
749
      b1= Balloc(b->k+1, alloc);
683
750
      Bcopy(b1, b);
684
 
      Bfree(b);
 
751
      Bfree(b, alloc);
685
752
      b= b1;
686
753
    }
687
754
    b->p.x[wds++]= (ULong) carry;
691
758
}
692
759
 
693
760
 
694
 
static Bigint *s2b(const char *s, int nd0, int nd, ULong y9)
 
761
static Bigint *s2b(const char *s, int nd0, int nd, ULong y9, Stack_alloc *alloc)
695
762
{
696
763
  Bigint *b;
697
764
  int i, k;
699
766
 
700
767
  x= (nd + 8) / 9;
701
768
  for (k= 0, y= 1; x > y; y <<= 1, k++) ;
702
 
  b= Balloc(k);
 
769
  b= Balloc(k, alloc);
703
770
  b->p.x[0]= y9;
704
771
  b->wds= 1;
705
 
 
 
772
  
706
773
  i= 9;
707
774
  if (9 < nd0)
708
775
  {
709
776
    s+= 9;
710
777
    do
711
 
      b= multadd(b, 10, *s++ - '0');
 
778
      b= multadd(b, 10, *s++ - '0', alloc);
712
779
    while (++i < nd0);
713
780
    s++;
714
781
  }
715
782
  else
716
783
    s+= 10;
717
784
  for(; i < nd; i++)
718
 
    b= multadd(b, 10, *s++ - '0');
 
785
    b= multadd(b, 10, *s++ - '0', alloc);
719
786
  return b;
720
787
}
721
788
 
806
873
 
807
874
/* Convert integer to Bigint number */
808
875
 
809
 
static Bigint *i2b(int i)
 
876
static Bigint *i2b(int i, Stack_alloc *alloc)
810
877
{
811
878
  Bigint *b;
812
879
 
813
 
  b= Balloc(1);
 
880
  b= Balloc(1, alloc);
814
881
  b->p.x[0]= i;
815
882
  b->wds= 1;
816
883
  return b;
819
886
 
820
887
/* Multiply two Bigint numbers */
821
888
 
822
 
static Bigint *mult(Bigint *a, Bigint *b)
 
889
static Bigint *mult(Bigint *a, Bigint *b, Stack_alloc *alloc)
823
890
{
824
891
  Bigint *c;
825
892
  int k, wa, wb, wc;
839
906
  wc= wa + wb;
840
907
  if (wc > a->maxwds)
841
908
    k++;
842
 
  c= Balloc(k);
 
909
  c= Balloc(k, alloc);
843
910
  for (x= c->p.x, xa= x + wc; x < xa; x++)
844
911
    *x= 0;
845
912
  xa= a->p.x;
911
978
 
912
979
#define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1)
913
980
 
914
 
static Bigint *pow5mult(Bigint *b, int k)
 
981
static Bigint *pow5mult(Bigint *b, int k, Stack_alloc *alloc)
915
982
{
916
983
  Bigint *b1, *p5, *p51;
917
984
  int i;
918
985
  static int p05[3]= { 5, 25, 125 };
919
986
 
920
987
  if ((i= k & 3))
921
 
    b= multadd(b, p05[i-1], 0);
 
988
    b= multadd(b, p05[i-1], 0, alloc);
922
989
 
923
990
  if (!(k>>= 2))
924
991
    return b;
927
994
  {
928
995
    if (k & 1)
929
996
    {
930
 
      b1= mult(b, p5);
931
 
      Bfree(b);
 
997
      b1= mult(b, p5, alloc);
 
998
      Bfree(b, alloc);
932
999
      b= b1;
933
1000
    }
934
1001
    if (!(k>>= 1))
937
1004
    if (p5 < p5_a + P5A_MAX)
938
1005
      ++p5;
939
1006
    else if (p5 == p5_a + P5A_MAX)
940
 
      p5= mult(p5, p5);
 
1007
      p5= mult(p5, p5, alloc);
941
1008
    else
942
1009
    {
943
 
      p51= mult(p5, p5);
944
 
      Bfree(p5);
 
1010
      p51= mult(p5, p5, alloc);
 
1011
      Bfree(p5, alloc);
945
1012
      p5= p51;
946
1013
    }
947
1014
  }
949
1016
}
950
1017
 
951
1018
 
952
 
static Bigint *lshift(Bigint *b, int k)
 
1019
static Bigint *lshift(Bigint *b, int k, Stack_alloc *alloc)
953
1020
{
954
1021
  int i, k1, n, n1;
955
1022
  Bigint *b1;
960
1027
  n1= n + b->wds + 1;
961
1028
  for (i= b->maxwds; n1 > i; i<<= 1)
962
1029
    k1++;
963
 
  b1= Balloc(k1);
 
1030
  b1= Balloc(k1, alloc);
964
1031
  x1= b1->p.x;
965
1032
  for (i= 0; i < n; i++)
966
1033
    *x1++= 0;
984
1051
      *x1++= *x++;
985
1052
    while (x < xe);
986
1053
  b1->wds= n1 - 1;
987
 
  Bfree(b);
 
1054
  Bfree(b, alloc);
988
1055
  return b1;
989
1056
}
990
1057
 
1013
1080
}
1014
1081
 
1015
1082
 
1016
 
static Bigint *diff(Bigint *a, Bigint *b)
 
1083
static Bigint *diff(Bigint *a, Bigint *b, Stack_alloc *alloc)
1017
1084
{
1018
1085
  Bigint *c;
1019
1086
  int i, wa, wb;
1023
1090
  i= cmp(a,b);
1024
1091
  if (!i)
1025
1092
  {
1026
 
    c= Balloc(0);
 
1093
    c= Balloc(0, alloc);
1027
1094
    c->wds= 1;
1028
1095
    c->p.x[0]= 0;
1029
1096
    return c;
1037
1104
  }
1038
1105
  else
1039
1106
    i= 0;
1040
 
  c= Balloc(a->k);
 
1107
  c= Balloc(a->k, alloc);
1041
1108
  c->sign= i;
1042
1109
  wa= a->wds;
1043
1110
  xa= a->p.x;
1118
1185
}
1119
1186
 
1120
1187
 
1121
 
static Bigint *d2b(double d, int *e, int *bits)
 
1188
static Bigint *d2b(double d, int *e, int *bits, Stack_alloc *alloc)
1122
1189
{
1123
1190
  Bigint *b;
1124
1191
  int de, k;
1127
1194
#define d0 word0(d)
1128
1195
#define d1 word1(d)
1129
1196
 
1130
 
  b= Balloc(1);
 
1197
  b= Balloc(1, alloc);
1131
1198
  x= b->p.x;
1132
1199
 
1133
1200
  z= d0 & Frac_mask;
1199
1266
  9007199254740992.*9007199254740992.e-256 /* = 2^106 * 1e-53 */
1200
1267
};
1201
1268
/*
1202
 
  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow
 
1269
  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow 
1203
1270
  flag unnecessarily.  It leads to a song and dance at the end of strtod.
1204
1271
*/
1205
1272
#define Scale_Bit 0x10
1207
1274
 
1208
1275
/*
1209
1276
  strtod for IEEE--arithmetic machines.
1210
 
 
 
1277
 
1211
1278
  This strtod returns a nearest machine number to the input decimal
1212
1279
  string (or sets errno to EOVERFLOW). Ties are broken by the IEEE round-even
1213
1280
  rule.
1214
 
 
 
1281
 
1215
1282
  Inspired loosely by William D. Clinger's paper "How to Read Floating
1216
1283
  Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
1217
 
 
 
1284
 
1218
1285
  Modifications:
1219
 
 
 
1286
 
1220
1287
   1. We only require IEEE (not IEEE double-extended).
1221
1288
   2. We get by with floating-point arithmetic in a case that
1222
1289
     Clinger missed -- when we're computing d * 10^n
1234
1301
     for 0 <= k <= 22).
1235
1302
*/
1236
1303
 
1237
 
static double my_strtod_int(const char *s00, char **se, int *error)
 
1304
static double my_strtod_int(const char *s00, char **se, int *error, char *buf, size_t buf_size)
1238
1305
{
1239
1306
  int scale;
1240
1307
  int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1247
1314
#ifdef SET_INEXACT
1248
1315
  int inexact, oldinexact;
1249
1316
#endif
 
1317
  Stack_alloc alloc;
1250
1318
 
1251
1319
  c= 0;
1252
1320
  *error= 0;
1253
1321
 
 
1322
  alloc.begin= alloc.free= buf;
 
1323
  alloc.end= buf + buf_size;
 
1324
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
 
1325
 
1254
1326
  sign= nz0= nz= 0;
1255
1327
  dval(rv)= 0.;
1256
1328
  for (s= s00; s < end; s++)
1274
1346
 break2:
1275
1347
  if (s >= end)
1276
1348
    goto ret0;
1277
 
 
 
1349
  
1278
1350
  if (*s == '0')
1279
1351
  {
1280
1352
    nz0= 1;
1526
1598
 
1527
1599
  /* Put digits into bd: true value = bd * 10^e */
1528
1600
 
1529
 
  bd0= s2b(s0, nd0, nd, y);
 
1601
  bd0= s2b(s0, nd0, nd, y, &alloc);
1530
1602
 
1531
1603
  for(;;)
1532
1604
  {
1533
 
    bd= Balloc(bd0->k);
 
1605
    bd= Balloc(bd0->k, &alloc);
1534
1606
    Bcopy(bd, bd0);
1535
 
    bb= d2b(dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
1536
 
    bs= i2b(1);
 
1607
    bb= d2b(dval(rv), &bbe, &bbbits, &alloc);  /* rv = bb * 2^bbe */
 
1608
    bs= i2b(1, &alloc);
1537
1609
 
1538
1610
    if (e >= 0)
1539
1611
    {
1570
1642
    }
1571
1643
    if (bb5 > 0)
1572
1644
    {
1573
 
      bs= pow5mult(bs, bb5);
1574
 
      bb1= mult(bs, bb);
1575
 
      Bfree(bb);
 
1645
      bs= pow5mult(bs, bb5, &alloc);
 
1646
      bb1= mult(bs, bb, &alloc);
 
1647
      Bfree(bb, &alloc);
1576
1648
      bb= bb1;
1577
1649
    }
1578
1650
    if (bb2 > 0)
1579
 
      bb= lshift(bb, bb2);
 
1651
      bb= lshift(bb, bb2, &alloc);
1580
1652
    if (bd5 > 0)
1581
 
      bd= pow5mult(bd, bd5);
 
1653
      bd= pow5mult(bd, bd5, &alloc);
1582
1654
    if (bd2 > 0)
1583
 
      bd= lshift(bd, bd2);
 
1655
      bd= lshift(bd, bd2, &alloc);
1584
1656
    if (bs2 > 0)
1585
 
      bs= lshift(bs, bs2);
1586
 
    delta= diff(bb, bd);
 
1657
      bs= lshift(bs, bs2, &alloc);
 
1658
    delta= diff(bb, bd, &alloc);
1587
1659
    dsign= delta->sign;
1588
1660
    delta->sign= 0;
1589
1661
    i= cmp(delta, bs);
1611
1683
#endif
1612
1684
        break;
1613
1685
      }
1614
 
      delta= lshift(delta, Log2P);
 
1686
      delta= lshift(delta, Log2P, &alloc);
1615
1687
      if (cmp(delta, bs) > 0)
1616
1688
        goto drop_down;
1617
1689
      break;
1752
1824
      }
1753
1825
#endif
1754
1826
 cont:
1755
 
    Bfree(bb);
1756
 
    Bfree(bd);
1757
 
    Bfree(bs);
1758
 
    Bfree(delta);
 
1827
    Bfree(bb, &alloc);
 
1828
    Bfree(bd, &alloc);
 
1829
    Bfree(bs, &alloc);
 
1830
    Bfree(delta, &alloc);
1759
1831
  }
1760
1832
#ifdef SET_INEXACT
1761
1833
  if (inexact)
1785
1857
  }
1786
1858
#endif
1787
1859
 retfree:
1788
 
  Bfree(bb);
1789
 
  Bfree(bd);
1790
 
  Bfree(bs);
1791
 
  Bfree(bd0);
1792
 
  Bfree(delta);
 
1860
  Bfree(bb, &alloc);
 
1861
  Bfree(bd, &alloc);
 
1862
  Bfree(bs, &alloc);
 
1863
  Bfree(bd0, &alloc);
 
1864
  Bfree(delta, &alloc);
1793
1865
 ret:
1794
1866
  *se= (char *)s;
1795
1867
  return sign ? -dval(rv) : dval(rv);
1896
1968
 */
1897
1969
 
1898
1970
static char *dtoa(double d, int mode, int ndigits, int *decpt, int *sign,
1899
 
                  char **rve)
 
1971
                  char **rve, char *buf, size_t buf_size)
1900
1972
{
1901
1973
  /*
1902
1974
    Arguments ndigits, decpt, sign are similar to those
1908
1980
    mode:
1909
1981
          0 ==> shortest string that yields d when read in
1910
1982
                and rounded to nearest.
1911
 
 
1912
1983
          1 ==> like 0, but with Steele & White stopping rule;
1913
1984
                e.g. with IEEE P754 arithmetic , mode 0 gives
1914
1985
                1e23 whereas mode 1 gives 9.999999999999999e22.
1940
2011
  Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
1941
2012
  double d2, ds, eps;
1942
2013
  char *s, *s0;
 
2014
  Stack_alloc alloc;
 
2015
  
 
2016
  alloc.begin= alloc.free= buf;
 
2017
  alloc.end= buf + buf_size;
 
2018
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
1943
2019
 
1944
2020
  if (word0(d) & Sign_bit)
1945
2021
  {
1955
2031
      (!dval(d) && (*decpt= 1)))
1956
2032
  {
1957
2033
    /* Infinity, NaN, 0 */
1958
 
    char *res= (char*) malloc(2);
 
2034
    char *res= (char*) dtoa_alloc(2, &alloc);
1959
2035
    res[0]= '0';
1960
2036
    res[1]= '\0';
1961
2037
    if (rve)
1962
2038
      *rve= res + 1;
1963
2039
    return res;
1964
2040
  }
1965
 
 
1966
 
 
1967
 
  b= d2b(dval(d), &be, &bbits);
 
2041
  
 
2042
 
 
2043
  b= d2b(dval(d), &be, &bbits, &alloc);
1968
2044
  if ((i= (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
1969
2045
  {
1970
2046
    dval(d2)= dval(d);
1976
2052
      log10(x)      =  log(x) / log(10)
1977
2053
                   ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1978
2054
      log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
1979
 
 
 
2055
     
1980
2056
      This suggests computing an approximation k to log10(d) by
1981
 
 
 
2057
     
1982
2058
      k= (i - Bias)*0.301029995663981
1983
2059
           + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1984
 
 
 
2060
     
1985
2061
      We want k to be too large rather than too small.
1986
2062
      The error in the first-order Taylor series approximation
1987
2063
      is in our favor, so we just round up the constant enough
2079
2155
    if (i <= 0)
2080
2156
      i= 1;
2081
2157
  }
2082
 
  s= s0= (char*)malloc(i+1); /* +1 for trailing '\0' appended
2083
 
                                at end of function */
 
2158
  s= s0= dtoa_alloc(i, &alloc);
2084
2159
 
2085
2160
  if (ilim >= 0 && ilim <= Quick_max && try_quick)
2086
2161
  {
2245
2320
    i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
2246
2321
    b2+= i;
2247
2322
    s2+= i;
2248
 
    mhi= i2b(1);
 
2323
    mhi= i2b(1, &alloc);
2249
2324
  }
2250
2325
  if (m2 > 0 && s2 > 0)
2251
2326
  {
2260
2335
    {
2261
2336
      if (m5 > 0)
2262
2337
      {
2263
 
        mhi= pow5mult(mhi, m5);
2264
 
        b1= mult(mhi, b);
2265
 
        Bfree(b);
 
2338
        mhi= pow5mult(mhi, m5, &alloc);
 
2339
        b1= mult(mhi, b, &alloc);
 
2340
        Bfree(b, &alloc);
2266
2341
        b= b1;
2267
2342
      }
2268
2343
      if ((j= b5 - m5))
2269
 
        b= pow5mult(b, j);
 
2344
        b= pow5mult(b, j, &alloc);
2270
2345
    }
2271
2346
    else
2272
 
      b= pow5mult(b, b5);
 
2347
      b= pow5mult(b, b5, &alloc);
2273
2348
  }
2274
 
  S= i2b(1);
 
2349
  S= i2b(1, &alloc);
2275
2350
  if (s5 > 0)
2276
 
    S= pow5mult(S, s5);
 
2351
    S= pow5mult(S, s5, &alloc);
2277
2352
 
2278
2353
  /* Check for special case that d is a normalized power of 2. */
2279
2354
 
2295
2370
  /*
2296
2371
    Arrange for convenient computation of quotients:
2297
2372
    shift left if necessary so divisor has 4 leading 0 bits.
2298
 
 
 
2373
    
2299
2374
    Perhaps we should just compute leading 28 bits of S once
2300
2375
    a nd for all and pass them and a shift to quorem, so it
2301
2376
    can do shifts and ors to compute the numerator for q.
2317
2392
    s2+= i;
2318
2393
  }
2319
2394
  if (b2 > 0)
2320
 
    b= lshift(b, b2);
 
2395
    b= lshift(b, b2, &alloc);
2321
2396
  if (s2 > 0)
2322
 
    S= lshift(S, s2);
 
2397
    S= lshift(S, s2, &alloc);
2323
2398
  if (k_check)
2324
2399
  {
2325
2400
    if (cmp(b,S) < 0)
2326
2401
    {
2327
2402
      k--;
2328
2403
      /* we botched the k estimate */
2329
 
      b= multadd(b, 10, 0);
 
2404
      b= multadd(b, 10, 0, &alloc);
2330
2405
      if (leftright)
2331
 
        mhi= multadd(mhi, 10, 0);
 
2406
        mhi= multadd(mhi, 10, 0, &alloc);
2332
2407
      ilim= ilim1;
2333
2408
    }
2334
2409
  }
2335
2410
  if (ilim <= 0 && (mode == 3 || mode == 5))
2336
2411
  {
2337
 
    if (ilim < 0 || cmp(b,S= multadd(S,5,0)) <= 0)
 
2412
    if (ilim < 0 || cmp(b,S= multadd(S,5,0, &alloc)) <= 0)
2338
2413
    {
2339
2414
      /* no digits, fcvt style */
2340
2415
no_digits:
2349
2424
  if (leftright)
2350
2425
  {
2351
2426
    if (m2 > 0)
2352
 
      mhi= lshift(mhi, m2);
 
2427
      mhi= lshift(mhi, m2, &alloc);
2353
2428
 
2354
2429
    /*
2355
2430
      Compute mlo -- check for special case that d is a normalized power of 2.
2358
2433
    mlo= mhi;
2359
2434
    if (spec_case)
2360
2435
    {
2361
 
      mhi= Balloc(mhi->k);
 
2436
      mhi= Balloc(mhi->k, &alloc);
2362
2437
      Bcopy(mhi, mlo);
2363
 
      mhi= lshift(mhi, Log2P);
 
2438
      mhi= lshift(mhi, Log2P, &alloc);
2364
2439
    }
2365
2440
 
2366
2441
    for (i= 1;;i++)
2368
2443
      dig= quorem(b,S) + '0';
2369
2444
      /* Do we yet have the shortest decimal string that will round to d? */
2370
2445
      j= cmp(b, mlo);
2371
 
      delta= diff(S, mhi);
 
2446
      delta= diff(S, mhi, &alloc);
2372
2447
      j1= delta->sign ? 1 : cmp(b, delta);
2373
 
      Bfree(delta);
 
2448
      Bfree(delta, &alloc);
2374
2449
      if (j1 == 0 && mode != 1 && !(word1(d) & 1)
2375
2450
         )
2376
2451
      {
2389
2464
        }
2390
2465
        if (j1 > 0)
2391
2466
        {
2392
 
          b= lshift(b, 1);
 
2467
          b= lshift(b, 1, &alloc);
2393
2468
          j1= cmp(b, S);
2394
2469
          if ((j1 > 0 || (j1 == 0 && dig & 1))
2395
2470
              && dig++ == '9')
2413
2488
      *s++= dig;
2414
2489
      if (i == ilim)
2415
2490
        break;
2416
 
      b= multadd(b, 10, 0);
 
2491
      b= multadd(b, 10, 0, &alloc);
2417
2492
      if (mlo == mhi)
2418
 
        mlo= mhi= multadd(mhi, 10, 0);
 
2493
        mlo= mhi= multadd(mhi, 10, 0, &alloc);
2419
2494
      else
2420
2495
      {
2421
 
        mlo= multadd(mlo, 10, 0);
2422
 
        mhi= multadd(mhi, 10, 0);
 
2496
        mlo= multadd(mlo, 10, 0, &alloc);
 
2497
        mhi= multadd(mhi, 10, 0, &alloc);
2423
2498
      }
2424
2499
    }
2425
2500
  }
2433
2508
      }
2434
2509
      if (i >= ilim)
2435
2510
        break;
2436
 
      b= multadd(b, 10, 0);
 
2511
      b= multadd(b, 10, 0, &alloc);
2437
2512
    }
2438
2513
 
2439
2514
  /* Round off last digit */
2440
2515
 
2441
 
  b= lshift(b, 1);
 
2516
  b= lshift(b, 1, &alloc);
2442
2517
  j= cmp(b, S);
2443
2518
  if (j > 0 || (j == 0 && dig & 1))
2444
2519
  {
2458
2533
    s++;
2459
2534
  }
2460
2535
ret:
2461
 
  Bfree(S);
 
2536
  Bfree(S, &alloc);
2462
2537
  if (mhi)
2463
2538
  {
2464
2539
    if (mlo && mlo != mhi)
2465
 
      Bfree(mlo);
2466
 
    Bfree(mhi);
 
2540
      Bfree(mlo, &alloc);
 
2541
    Bfree(mhi, &alloc);
2467
2542
  }
2468
2543
ret1:
2469
 
  Bfree(b);
 
2544
  Bfree(b, &alloc);
2470
2545
  *s= 0;
2471
2546
  *decpt= k + 1;
2472
2547
  if (rve)
2473
2548
    *rve= s;
2474
2549
  return s0;
2475
2550
}
2476
 
 
2477
 
} /* namespace internal */
2478
 
} /* namespace drizzled */