← Back to branch summary

~drizzle-trunk/drizzle/development

« back to all changes in this revision

Viewing changes to strings/dtoa.c

  • Committer: Brian Aker
  • Date: 2008-07-14 16:09:55 UTC
  • Revision ID: brian@tangent.org-20080714160955-v5nzzyjj5hhv7bz6
Removing a few "additional" ways of saying uint64_t

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
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
 
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
 
  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
 
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
  {
1951
2187
    *sign= 0;
1952
2188
 
1953
2189
  /* 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)))
 
2190
  if (((word0(d) & Exp_mask) == Exp_mask && (*decpt= DTOA_OVERFLOW)) ||
 
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 */