← Back to branch summary

~drizzle-trunk/drizzle/development

« back to all changes in this revision

Viewing changes to strings/dtoa.c

Merged in changes. 
Edited a the comment test case so deal with our version bump.

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;
 
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 *))
49
49
 
50
50
/* Magic value returned by dtoa() to indicate overflow */
51
51
#define DTOA_OVERFLOW 9999
52
52
 
53
 
static double my_strtod_int(const char *, char **, int *);
54
 
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);
55
56
 
56
57
/**
57
58
   @brief
78
79
                      (including the terminating '\0').
79
80
   @param error       if not NULL, points to a location where the status of
80
81
                      conversion is stored upon return.
81
 
                      false  successful conversion
82
 
                      true   the input number is [-,+]infinity or nan.
 
82
                      FALSE  successful conversion
 
83
                      TRUE   the input number is [-,+]infinity or nan.
83
84
                             The output string in this case is always '0'.
84
85
   @return            number of written characters (excluding terminating '\0')
85
86
*/
86
87
 
87
 
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)
88
89
{
89
90
  int decpt, sign, len, i;
90
91
  char *res, *src, *end, *dst= to;
91
 
  assert(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
92
 
 
93
 
  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));
94
96
 
95
97
  if (decpt == DTOA_OVERFLOW)
96
98
  {
97
 
    free(res);
 
99
    dtoa_free(res, buf, sizeof(buf));
98
100
    *to++= '0';
99
101
    *to= '\0';
100
102
    if (error != NULL)
101
 
      *error= true;
 
103
      *error= TRUE;
102
104
    return 1;
103
105
  }
104
106
 
129
131
  {
130
132
    if (len <= decpt)
131
133
      *dst++= '.';
132
 
 
 
134
    
133
135
    for (i= precision - max(0, (len - decpt)); i > 0; i--)
134
136
      *dst++= '0';
135
137
  }
136
 
 
 
138
  
137
139
  *dst= '\0';
138
140
  if (error != NULL)
139
 
    *error= false;
 
141
    *error= FALSE;
140
142
 
141
 
  free(res);
 
143
  dtoa_free(res, buf, sizeof(buf));
142
144
 
143
145
  return dst - to;
144
146
}
179
181
                      'width + 1' bytes.
180
182
   @param error       if not NULL, points to a location where the status of
181
183
                      conversion is stored upon return.
182
 
                      false  successful conversion
183
 
                      true   the input number is [-,+]infinity or nan.
 
184
                      FALSE  successful conversion
 
185
                      TRUE   the input number is [-,+]infinity or nan.
184
186
                             The output string in this case is always '0'.
185
187
   @return            number of written characters (excluding terminating '\0')
186
188
 
196
198
     my_gcvt(55, ..., 1, ...);
197
199
 
198
200
   We do our best to minimize such cases by:
199
 
 
 
201
   
200
202
   - passing to dtoa() the field width as the number of significant digits
201
 
 
 
203
   
202
204
   - removing the sign of the number early (and decreasing the width before
203
205
     passing it to dtoa())
204
 
 
 
206
   
205
207
   - choosing the proper format to preserve the most number of significant
206
208
     digits.
207
209
*/
208
210
 
209
211
size_t my_gcvt(double x, my_gcvt_arg_type type, int width, char *to,
210
 
               bool *error)
 
212
               my_bool *error)
211
213
{
212
214
  int decpt, sign, len, exp_len;
213
215
  char *res, *src, *end, *dst= to, *dend= dst + width;
214
 
  bool have_space, force_e_format;
215
 
  assert(width > 0 && to != NULL);
216
 
 
 
216
  char buf[DTOA_BUFF_SIZE];
 
217
  my_bool have_space, force_e_format;
 
218
  DBUG_ASSERT(width > 0 && to != NULL);
 
219
  
217
220
  /* We want to remove '-' from equations early */
218
221
  if (x < 0.)
219
222
    width--;
220
223
 
221
 
  res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? min(width, DBL_DIG) : 
222
 
            min(width, FLT_DIG), &decpt, &sign, &end);
223
 
 
 
224
  res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? width : min(width, FLT_DIG),
 
225
            &decpt, &sign, &end, buf, sizeof(buf));
224
226
  if (decpt == DTOA_OVERFLOW)
225
227
  {
226
 
    free(res);
 
228
    dtoa_free(res, buf, sizeof(buf));
227
229
    *to++= '0';
228
230
    *to= '\0';
229
231
    if (error != NULL)
230
 
      *error= true;
 
232
      *error= TRUE;
231
233
    return 1;
232
234
  }
233
235
 
234
236
  if (error != NULL)
235
 
    *error= false;
 
237
    *error= FALSE;
236
238
 
237
239
  src= res;
238
240
  len= end - res;
243
245
     to count it here.
244
246
   */
245
247
  exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
246
 
 
 
248
  
247
249
  /*
248
250
     Do we have enough space for all digits in the 'f' format?
249
251
     Let 'len' be the number of significant digits returned by dtoa,
296
298
       ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
297
299
                                            (len > 1 || !force_e_format)))) &&
298
300
         !force_e_format)) &&
299
 
 
 
301
      
300
302
       /*
301
303
         Use the 'e' format in some cases even if we have enough space for the
302
304
         'f' one. See comment for MAX_DECPT_FOR_F_FORMAT.
315
317
      if (width < decpt)
316
318
      {
317
319
        if (error != NULL)
318
 
          *error= true;
 
320
          *error= TRUE;
319
321
        width= decpt;
320
322
      }
321
 
 
 
323
      
322
324
      /*
323
325
        We want to truncate (len - width) least significant digits after the
324
326
        decimal point. For this we are calling dtoa with mode=5, passing the
325
327
        number of significant digits = (len-decpt) - (len-width) = width-decpt
326
328
      */
327
 
      free(res);
328
 
      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));
329
331
      src= res;
330
332
      len= end - res;
331
333
    }
336
338
      *dst++= '0';
337
339
      goto end;
338
340
    }
339
 
 
 
341
    
340
342
    /*
341
343
      At this point we are sure we have enough space to put all digits
342
344
      returned by dtoa
382
384
    {
383
385
      /* Overflow */
384
386
      if (error != NULL)
385
 
        *error= true;
 
387
        *error= TRUE;
386
388
      width= 0;
387
389
    }
388
 
 
 
390
      
389
391
    /* Do we have to truncate any digits? */
390
392
    if (width < len)
391
393
    {
392
394
      /* Yes, re-convert with a smaller width */
393
 
      free(res);
394
 
      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));
395
397
      src= res;
396
398
      len= end - res;
397
399
      if (--decpt < 0)
431
433
  }
432
434
 
433
435
end:
434
 
  free(res);
 
436
  dtoa_free(res, buf, sizeof(buf));
435
437
  *dst= '\0';
436
438
 
437
439
  return dst - to;
450
452
                  rejected character.
451
453
   @param error   Upon return is set to EOVERFLOW in case of underflow or
452
454
                  overflow.
453
 
 
 
455
   
454
456
   @return        The resulting double value. In case of underflow, 0.0 is
455
457
                  returned. In case overflow, signed DBL_MAX is returned.
456
458
*/
457
459
 
458
460
double my_strtod(const char *str, char **end, int *error)
459
461
{
 
462
  char buf[DTOA_BUFF_SIZE];
460
463
  double res;
461
 
  assert(str != NULL && end != NULL && *end != NULL && error != NULL);
 
464
  DBUG_ASSERT(str != NULL && end != NULL && *end != NULL && error != NULL);
462
465
 
463
 
  res= my_strtod_int(str, end, error);
 
466
  res= my_strtod_int(str, end, error, buf, sizeof(buf));
464
467
  return (*error == 0) ? res : (res < 0 ? -DBL_MAX : DBL_MAX);
465
468
}
466
469
 
527
530
  file.
528
531
*/
529
532
 
530
 
typedef int32_t Long;
531
 
typedef uint32_t ULong;
532
 
typedef int64_t LLong;
533
 
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;
534
546
 
535
547
typedef union { double d; ULong L[2]; } U;
536
548
 
582
594
#endif
583
595
#endif /*Flt_Rounds*/
584
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
 
585
605
#define rounded_product(a,b) a*= b
586
606
#define rounded_quotient(a,b) a/= b
587
607
 
589
609
#define Big1 0xffffffff
590
610
#define FFFFFFFF 0xffffffffUL
591
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
 
592
619
/* Arbitrary-length integer */
593
620
 
594
621
typedef struct Bigint
603
630
  int wds;                 /* current length in 32-bit words */
604
631
} Bigint;
605
632
 
606
 
static Bigint *Bcopy(Bigint* dst, Bigint* src)
 
633
 
 
634
/* A simple stack-memory based allocator for Bigints */
 
635
 
 
636
typedef struct Stack_alloc
607
637
{
608
 
  dst->sign= src->sign;
609
 
  dst->wds= src->wds;
610
 
 
611
 
  assert(dst->maxwds >= src->wds);
612
 
 
613
 
  memcpy(dst->p.x, src->p.x, src->wds*sizeof(ULong));
614
 
 
615
 
  return dst;
616
 
}
617
 
 
618
 
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)
619
655
{
620
656
  Bigint *rv;
621
 
 
622
 
  /* TODO: some malloc failure checking */
623
 
 
624
 
  int x= 1 << k;
625
 
  rv= (Bigint*) malloc(sizeof(Bigint));
626
 
 
627
 
  rv->p.x= (ULong*)malloc(x * sizeof(ULong));
628
 
 
629
 
  rv->k= k;
630
 
  rv->maxwds= x;
631
 
 
 
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
  }
632
680
  rv->sign= rv->wds= 0;
633
 
 
 
681
  rv->p.x= (ULong*) (rv + 1);
634
682
  return rv;
635
683
}
636
684
 
640
688
  list. Otherwise call free().
641
689
*/
642
690
 
643
 
static void Bfree(Bigint *v)
644
 
{
645
 
  if(!v)
646
 
    return;
647
 
  free(v->p.x);
648
 
  free(v);
649
 
}
 
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
 
650
741
 
651
742
/* Bigint arithmetic functions */
652
743
 
653
744
/* Multiply by m and add a */
654
745
 
655
 
static Bigint *multadd(Bigint *b, int m, int a)
 
746
static Bigint *multadd(Bigint *b, int m, int a, Stack_alloc *alloc)
656
747
{
657
748
  int i, wds;
658
749
  ULong *x;
674
765
  {
675
766
    if (wds >= b->maxwds)
676
767
    {
677
 
      b1= Balloc(b->k+1);
 
768
      b1= Balloc(b->k+1, alloc);
678
769
      Bcopy(b1, b);
679
 
      Bfree(b);
 
770
      Bfree(b, alloc);
680
771
      b= b1;
681
772
    }
682
773
    b->p.x[wds++]= (ULong) carry;
686
777
}
687
778
 
688
779
 
689
 
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)
690
781
{
691
782
  Bigint *b;
692
783
  int i, k;
694
785
 
695
786
  x= (nd + 8) / 9;
696
787
  for (k= 0, y= 1; x > y; y <<= 1, k++) ;
697
 
  b= Balloc(k);
 
788
  b= Balloc(k, alloc);
698
789
  b->p.x[0]= y9;
699
790
  b->wds= 1;
700
 
 
 
791
  
701
792
  i= 9;
702
793
  if (9 < nd0)
703
794
  {
704
795
    s+= 9;
705
796
    do
706
 
      b= multadd(b, 10, *s++ - '0');
 
797
      b= multadd(b, 10, *s++ - '0', alloc);
707
798
    while (++i < nd0);
708
799
    s++;
709
800
  }
710
801
  else
711
802
    s+= 10;
712
803
  for(; i < nd; i++)
713
 
    b= multadd(b, 10, *s++ - '0');
 
804
    b= multadd(b, 10, *s++ - '0', alloc);
714
805
  return b;
715
806
}
716
807
 
801
892
 
802
893
/* Convert integer to Bigint number */
803
894
 
804
 
static Bigint *i2b(int i)
 
895
static Bigint *i2b(int i, Stack_alloc *alloc)
805
896
{
806
897
  Bigint *b;
807
898
 
808
 
  b= Balloc(1);
 
899
  b= Balloc(1, alloc);
809
900
  b->p.x[0]= i;
810
901
  b->wds= 1;
811
902
  return b;
814
905
 
815
906
/* Multiply two Bigint numbers */
816
907
 
817
 
static Bigint *mult(Bigint *a, Bigint *b)
 
908
static Bigint *mult(Bigint *a, Bigint *b, Stack_alloc *alloc)
818
909
{
819
910
  Bigint *c;
820
911
  int k, wa, wb, wc;
834
925
  wc= wa + wb;
835
926
  if (wc > a->maxwds)
836
927
    k++;
837
 
  c= Balloc(k);
 
928
  c= Balloc(k, alloc);
838
929
  for (x= c->p.x, xa= x + wc; x < xa; x++)
839
930
    *x= 0;
840
931
  xa= a->p.x;
906
997
 
907
998
#define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1)
908
999
 
909
 
static Bigint *pow5mult(Bigint *b, int k)
 
1000
static Bigint *pow5mult(Bigint *b, int k, Stack_alloc *alloc)
910
1001
{
911
1002
  Bigint *b1, *p5, *p51;
912
1003
  int i;
913
1004
  static int p05[3]= { 5, 25, 125 };
914
1005
 
915
1006
  if ((i= k & 3))
916
 
    b= multadd(b, p05[i-1], 0);
 
1007
    b= multadd(b, p05[i-1], 0, alloc);
917
1008
 
918
1009
  if (!(k>>= 2))
919
1010
    return b;
922
1013
  {
923
1014
    if (k & 1)
924
1015
    {
925
 
      b1= mult(b, p5);
926
 
      Bfree(b);
 
1016
      b1= mult(b, p5, alloc);
 
1017
      Bfree(b, alloc);
927
1018
      b= b1;
928
1019
    }
929
1020
    if (!(k>>= 1))
932
1023
    if (p5 < p5_a + P5A_MAX)
933
1024
      ++p5;
934
1025
    else if (p5 == p5_a + P5A_MAX)
935
 
      p5= mult(p5, p5);
 
1026
      p5= mult(p5, p5, alloc);
936
1027
    else
937
1028
    {
938
 
      p51= mult(p5, p5);
939
 
      Bfree(p5);
 
1029
      p51= mult(p5, p5, alloc);
 
1030
      Bfree(p5, alloc);
940
1031
      p5= p51;
941
1032
    }
942
1033
  }
944
1035
}
945
1036
 
946
1037
 
947
 
static Bigint *lshift(Bigint *b, int k)
 
1038
static Bigint *lshift(Bigint *b, int k, Stack_alloc *alloc)
948
1039
{
949
1040
  int i, k1, n, n1;
950
1041
  Bigint *b1;
955
1046
  n1= n + b->wds + 1;
956
1047
  for (i= b->maxwds; n1 > i; i<<= 1)
957
1048
    k1++;
958
 
  b1= Balloc(k1);
 
1049
  b1= Balloc(k1, alloc);
959
1050
  x1= b1->p.x;
960
1051
  for (i= 0; i < n; i++)
961
1052
    *x1++= 0;
979
1070
      *x1++= *x++;
980
1071
    while (x < xe);
981
1072
  b1->wds= n1 - 1;
982
 
  Bfree(b);
 
1073
  Bfree(b, alloc);
983
1074
  return b1;
984
1075
}
985
1076
 
1008
1099
}
1009
1100
 
1010
1101
 
1011
 
static Bigint *diff(Bigint *a, Bigint *b)
 
1102
static Bigint *diff(Bigint *a, Bigint *b, Stack_alloc *alloc)
1012
1103
{
1013
1104
  Bigint *c;
1014
1105
  int i, wa, wb;
1018
1109
  i= cmp(a,b);
1019
1110
  if (!i)
1020
1111
  {
1021
 
    c= Balloc(0);
 
1112
    c= Balloc(0, alloc);
1022
1113
    c->wds= 1;
1023
1114
    c->p.x[0]= 0;
1024
1115
    return c;
1032
1123
  }
1033
1124
  else
1034
1125
    i= 0;
1035
 
  c= Balloc(a->k);
 
1126
  c= Balloc(a->k, alloc);
1036
1127
  c->sign= i;
1037
1128
  wa= a->wds;
1038
1129
  xa= a->p.x;
1113
1204
}
1114
1205
 
1115
1206
 
1116
 
static Bigint *d2b(double d, int *e, int *bits)
 
1207
static Bigint *d2b(double d, int *e, int *bits, Stack_alloc *alloc)
1117
1208
{
1118
1209
  Bigint *b;
1119
1210
  int de, k;
1122
1213
#define d0 word0(d)
1123
1214
#define d1 word1(d)
1124
1215
 
1125
 
  b= Balloc(1);
 
1216
  b= Balloc(1, alloc);
1126
1217
  x= b->p.x;
1127
1218
 
1128
1219
  z= d0 & Frac_mask;
1194
1285
  9007199254740992.*9007199254740992.e-256 /* = 2^106 * 1e-53 */
1195
1286
};
1196
1287
/*
1197
 
  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 
1198
1289
  flag unnecessarily.  It leads to a song and dance at the end of strtod.
1199
1290
*/
1200
1291
#define Scale_Bit 0x10
1202
1293
 
1203
1294
/*
1204
1295
  strtod for IEEE--arithmetic machines.
1205
 
 
 
1296
 
1206
1297
  This strtod returns a nearest machine number to the input decimal
1207
1298
  string (or sets errno to EOVERFLOW). Ties are broken by the IEEE round-even
1208
1299
  rule.
1209
 
 
 
1300
 
1210
1301
  Inspired loosely by William D. Clinger's paper "How to Read Floating
1211
1302
  Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
1212
 
 
 
1303
 
1213
1304
  Modifications:
1214
 
 
 
1305
 
1215
1306
   1. We only require IEEE (not IEEE double-extended).
1216
1307
   2. We get by with floating-point arithmetic in a case that
1217
1308
     Clinger missed -- when we're computing d * 10^n
1229
1320
     for 0 <= k <= 22).
1230
1321
*/
1231
1322
 
1232
 
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)
1233
1324
{
1234
1325
  int scale;
1235
1326
  int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1238
1329
  double aadj, aadj1, adj, rv, rv0;
1239
1330
  Long L;
1240
1331
  ULong y, z;
1241
 
  Bigint *bb= NULL, *bb1= NULL, *bd= NULL, *bd0= NULL, *bs= NULL, *delta= NULL;
 
1332
  Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1242
1333
#ifdef SET_INEXACT
1243
1334
  int inexact, oldinexact;
1244
1335
#endif
 
1336
#ifdef Honor_FLT_ROUNDS
 
1337
  int rounding;
 
1338
#endif
 
1339
  Stack_alloc alloc;
1245
1340
 
1246
 
  c= 0;
1247
1341
  *error= 0;
1248
1342
 
 
1343
  alloc.begin= alloc.free= buf;
 
1344
  alloc.end= buf + buf_size;
 
1345
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
 
1346
 
1249
1347
  sign= nz0= nz= 0;
1250
1348
  dval(rv)= 0.;
1251
1349
  for (s= s00; s < end; s++)
1269
1367
 break2:
1270
1368
  if (s >= end)
1271
1369
    goto ret0;
1272
 
 
 
1370
  
1273
1371
  if (*s == '0')
1274
1372
  {
1275
1373
    nz0= 1;
1394
1492
  }
1395
1493
  bd0= 0;
1396
1494
  if (nd <= DBL_DIG
 
1495
#ifndef Honor_FLT_ROUNDS
 
1496
    && Flt_Rounds == 1
 
1497
#endif
1397
1498
      )
1398
1499
  {
1399
1500
    if (!e)
1402
1503
    {
1403
1504
      if (e <= Ten_pmax)
1404
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
1405
1514
        /* rv = */ rounded_product(dval(rv), tens[e]);
1406
1515
        goto ret;
1407
1516
      }
1412
1521
          A fancier test would sometimes let us do
1413
1522
          this for larger i values.
1414
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
1415
1532
        e-= i;
1416
1533
        dval(rv)*= tens[i];
1417
1534
        /* rv = */ rounded_product(dval(rv), tens[e]);
1421
1538
#ifndef Inaccurate_Divide
1422
1539
    else if (e >= -Ten_pmax)
1423
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
1424
1549
      /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1425
1550
      goto ret;
1426
1551
    }
1434
1559
    oldinexact= get_inexact();
1435
1560
#endif
1436
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
1437
1572
 
1438
1573
  /* Get starting approximation = rv * 10**e1 */
1439
1574
 
1448
1583
 ovfl:
1449
1584
        *error= EOVERFLOW;
1450
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*/
1451
1599
        word0(rv)= Exp_mask;
1452
1600
        word1(rv)= 0;
 
1601
#endif /*Honor_FLT_ROUNDS*/
1453
1602
#ifdef SET_INEXACT
1454
1603
        /* set overflow bit */
1455
1604
        dval(rv0)= 1e300;
1521
1670
 
1522
1671
  /* Put digits into bd: true value = bd * 10^e */
1523
1672
 
1524
 
  bd0= s2b(s0, nd0, nd, y);
 
1673
  bd0= s2b(s0, nd0, nd, y, &alloc);
1525
1674
 
1526
1675
  for(;;)
1527
1676
  {
1528
 
    bd= Balloc(bd0->k);
 
1677
    bd= Balloc(bd0->k, &alloc);
1529
1678
    Bcopy(bd, bd0);
1530
 
    bb= d2b(dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
1531
 
    bs= i2b(1);
 
1679
    bb= d2b(dval(rv), &bbe, &bbbits, &alloc);  /* rv = bb * 2^bbe */
 
1680
    bs= i2b(1, &alloc);
1532
1681
 
1533
1682
    if (e >= 0)
1534
1683
    {
1545
1694
    else
1546
1695
      bd2-= bbe;
1547
1696
    bs2= bb2;
 
1697
#ifdef Honor_FLT_ROUNDS
 
1698
    if (rounding != 1)
 
1699
      bs2++;
 
1700
#endif
1548
1701
    j= bbe - scale;
1549
1702
    i= j + bbbits - 1;  /* logb(rv) */
1550
1703
    if (i < Emin)  /* denormal */
1565
1718
    }
1566
1719
    if (bb5 > 0)
1567
1720
    {
1568
 
      bs= pow5mult(bs, bb5);
1569
 
      bb1= mult(bs, bb);
1570
 
      Bfree(bb);
 
1721
      bs= pow5mult(bs, bb5, &alloc);
 
1722
      bb1= mult(bs, bb, &alloc);
 
1723
      Bfree(bb, &alloc);
1571
1724
      bb= bb1;
1572
1725
    }
1573
1726
    if (bb2 > 0)
1574
 
      bb= lshift(bb, bb2);
 
1727
      bb= lshift(bb, bb2, &alloc);
1575
1728
    if (bd5 > 0)
1576
 
      bd= pow5mult(bd, bd5);
 
1729
      bd= pow5mult(bd, bd5, &alloc);
1577
1730
    if (bd2 > 0)
1578
 
      bd= lshift(bd, bd2);
 
1731
      bd= lshift(bd, bd2, &alloc);
1579
1732
    if (bs2 > 0)
1580
 
      bs= lshift(bs, bs2);
1581
 
    delta= diff(bb, bd);
 
1733
      bs= lshift(bs, bs2, &alloc);
 
1734
    delta= diff(bb, bd, &alloc);
1582
1735
    dsign= delta->sign;
1583
1736
    delta->sign= 0;
1584
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*/
1585
1804
 
1586
1805
    if (i < 0)
1587
1806
    {
1606
1825
#endif
1607
1826
        break;
1608
1827
      }
1609
 
      delta= lshift(delta, Log2P);
 
1828
      delta= lshift(delta, Log2P, &alloc);
1610
1829
      if (cmp(delta, bs) > 0)
1611
1830
        goto drop_down;
1612
1831
      break;
1688
1907
    {
1689
1908
      aadj*= 0.5;
1690
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
1691
1921
      if (Flt_Rounds == 0)
1692
1922
        aadj1+= 0.5;
 
1923
#endif /*Check_FLT_ROUNDS*/
1693
1924
    }
1694
1925
    y= word0(rv) & Exp_mask;
1695
1926
 
1747
1978
      }
1748
1979
#endif
1749
1980
 cont:
1750
 
    Bfree(bb);
1751
 
    Bfree(bd);
1752
 
    Bfree(bs);
1753
 
    Bfree(delta);
 
1981
    Bfree(bb, &alloc);
 
1982
    Bfree(bd, &alloc);
 
1983
    Bfree(bs, &alloc);
 
1984
    Bfree(delta, &alloc);
1754
1985
  }
1755
1986
#ifdef SET_INEXACT
1756
1987
  if (inexact)
1780
2011
  }
1781
2012
#endif
1782
2013
 retfree:
1783
 
  Bfree(bb);
1784
 
  Bfree(bd);
1785
 
  Bfree(bs);
1786
 
  Bfree(bd0);
1787
 
  Bfree(delta);
 
2014
  Bfree(bb, &alloc);
 
2015
  Bfree(bd, &alloc);
 
2016
  Bfree(bs, &alloc);
 
2017
  Bfree(bd0, &alloc);
 
2018
  Bfree(delta, &alloc);
1788
2019
 ret:
1789
2020
  *se= (char *)s;
1790
2021
  return sign ? -dval(rv) : dval(rv);
1891
2122
 */
1892
2123
 
1893
2124
static char *dtoa(double d, int mode, int ndigits, int *decpt, int *sign,
1894
 
                  char **rve)
 
2125
                  char **rve, char *buf, size_t buf_size)
1895
2126
{
1896
2127
  /*
1897
2128
    Arguments ndigits, decpt, sign are similar to those
1903
2134
    mode:
1904
2135
          0 ==> shortest string that yields d when read in
1905
2136
                and rounded to nearest.
1906
 
 
1907
2137
          1 ==> like 0, but with Steele & White stopping rule;
1908
2138
                e.g. with IEEE P754 arithmetic , mode 0 gives
1909
2139
                1e23 whereas mode 1 gives 9.999999999999999e22.
1910
 
          2 ==> cmax(1,ndigits) significant digits.  This gives a
 
2140
          2 ==> max(1,ndigits) significant digits.  This gives a
1911
2141
                return value similar to that of ecvt, except
1912
2142
                that trailing zeros are suppressed.
1913
2143
          3 ==> through ndigits past the decimal point.  This
1917
2147
          4,5 ==> similar to 2 and 3, respectively, but (in
1918
2148
                round-nearest mode) with the tests of mode 0 to
1919
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.
1920
2153
          6-9 ==> Debugging modes similar to mode - 4:  don't try
1921
2154
                fast floating-point estimate (if applicable).
1922
2155
 
1926
2159
    to hold the suppressed trailing zeros.
1927
2160
  */
1928
2161
 
1929
 
  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,
1930
2163
    j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1931
2164
    spec_case, try_quick;
1932
2165
  Long L;
1933
2166
  int denorm;
1934
2167
  ULong x;
1935
 
  Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
 
2168
  Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1936
2169
  double d2, ds, eps;
1937
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));
1938
2179
 
1939
2180
  if (word0(d) & Sign_bit)
1940
2181
  {
1950
2191
      (!dval(d) && (*decpt= 1)))
1951
2192
  {
1952
2193
    /* Infinity, NaN, 0 */
1953
 
    char *res= (char*) malloc(2);
 
2194
    char *res= (char*) dtoa_alloc(2, &alloc);
1954
2195
    res[0]= '0';
1955
2196
    res[1]= '\0';
1956
2197
    if (rve)
1957
2198
      *rve= res + 1;
1958
2199
    return res;
1959
2200
  }
1960
 
 
1961
 
 
1962
 
  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);
1963
2214
  if ((i= (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
1964
2215
  {
1965
2216
    dval(d2)= dval(d);
1971
2222
      log10(x)      =  log(x) / log(10)
1972
2223
                   ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1973
2224
      log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
1974
 
 
 
2225
     
1975
2226
      This suggests computing an approximation k to log10(d) by
1976
 
 
 
2227
     
1977
2228
      k= (i - Bias)*0.301029995663981
1978
2229
           + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1979
 
 
 
2230
     
1980
2231
      We want k to be too large rather than too small.
1981
2232
      The error in the first-order Taylor series approximation
1982
2233
      is in our favor, so we just round up the constant enough
2041
2292
  if (mode < 0 || mode > 9)
2042
2293
    mode= 0;
2043
2294
 
 
2295
#ifdef Check_FLT_ROUNDS
 
2296
  try_quick= Rounding == 1;
 
2297
#else
2044
2298
  try_quick= 1;
 
2299
#endif
2045
2300
 
2046
2301
  if (mode > 5)
2047
2302
  {
2074
2329
    if (i <= 0)
2075
2330
      i= 1;
2076
2331
  }
2077
 
  s= s0= (char*)malloc(i+1); /* +1 for trailing '\0' appended
2078
 
                                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
2079
2338
 
2080
2339
  if (ilim >= 0 && ilim <= Quick_max && try_quick)
2081
2340
  {
2174
2433
            goto bump_up;
2175
2434
          else if (dval(d) < 0.5 - dval(eps))
2176
2435
          {
2177
 
            while (*--s == '0') {}
 
2436
            while (*--s == '0');
2178
2437
            s++;
2179
2438
            goto ret1;
2180
2439
          }
2206
2465
    {
2207
2466
      L= (Long)(dval(d) / ds);
2208
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
2209
2476
      *s++= '0' + (int)L;
2210
2477
      if (!dval(d))
2211
2478
      {
2213
2480
      }
2214
2481
      if (i == ilim)
2215
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
2216
2492
        dval(d)+= dval(d);
2217
2493
        if (dval(d) > ds || (dval(d) == ds && L & 1))
2218
2494
        {
2240
2516
    i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
2241
2517
    b2+= i;
2242
2518
    s2+= i;
2243
 
    mhi= i2b(1);
 
2519
    mhi= i2b(1, &alloc);
2244
2520
  }
2245
2521
  if (m2 > 0 && s2 > 0)
2246
2522
  {
2255
2531
    {
2256
2532
      if (m5 > 0)
2257
2533
      {
2258
 
        mhi= pow5mult(mhi, m5);
2259
 
        b1= mult(mhi, b);
2260
 
        Bfree(b);
 
2534
        mhi= pow5mult(mhi, m5, &alloc);
 
2535
        b1= mult(mhi, b, &alloc);
 
2536
        Bfree(b, &alloc);
2261
2537
        b= b1;
2262
2538
      }
2263
2539
      if ((j= b5 - m5))
2264
 
        b= pow5mult(b, j);
 
2540
        b= pow5mult(b, j, &alloc);
2265
2541
    }
2266
2542
    else
2267
 
      b= pow5mult(b, b5);
 
2543
      b= pow5mult(b, b5, &alloc);
2268
2544
  }
2269
 
  S= i2b(1);
 
2545
  S= i2b(1, &alloc);
2270
2546
  if (s5 > 0)
2271
 
    S= pow5mult(S, s5);
 
2547
    S= pow5mult(S, s5, &alloc);
2272
2548
 
2273
2549
  /* Check for special case that d is a normalized power of 2. */
2274
2550
 
2275
2551
  spec_case= 0;
2276
2552
  if ((mode < 2 || leftright)
 
2553
#ifdef Honor_FLT_ROUNDS
 
2554
      && rounding == 1
 
2555
#endif
2277
2556
     )
2278
2557
  {
2279
2558
    if (!word1(d) && !(word0(d) & Bndry_mask) &&
2290
2569
  /*
2291
2570
    Arrange for convenient computation of quotients:
2292
2571
    shift left if necessary so divisor has 4 leading 0 bits.
2293
 
 
 
2572
    
2294
2573
    Perhaps we should just compute leading 28 bits of S once
2295
2574
    a nd for all and pass them and a shift to quorem, so it
2296
2575
    can do shifts and ors to compute the numerator for q.
2312
2591
    s2+= i;
2313
2592
  }
2314
2593
  if (b2 > 0)
2315
 
    b= lshift(b, b2);
 
2594
    b= lshift(b, b2, &alloc);
2316
2595
  if (s2 > 0)
2317
 
    S= lshift(S, s2);
 
2596
    S= lshift(S, s2, &alloc);
2318
2597
  if (k_check)
2319
2598
  {
2320
2599
    if (cmp(b,S) < 0)
2321
2600
    {
2322
2601
      k--;
2323
2602
      /* we botched the k estimate */
2324
 
      b= multadd(b, 10, 0);
 
2603
      b= multadd(b, 10, 0, &alloc);
2325
2604
      if (leftright)
2326
 
        mhi= multadd(mhi, 10, 0);
 
2605
        mhi= multadd(mhi, 10, 0, &alloc);
2327
2606
      ilim= ilim1;
2328
2607
    }
2329
2608
  }
2330
2609
  if (ilim <= 0 && (mode == 3 || mode == 5))
2331
2610
  {
2332
 
    if (ilim < 0 || cmp(b,S= multadd(S,5,0)) <= 0)
 
2611
    if (ilim < 0 || cmp(b,S= multadd(S,5,0, &alloc)) <= 0)
2333
2612
    {
2334
2613
      /* no digits, fcvt style */
2335
2614
no_digits:
2344
2623
  if (leftright)
2345
2624
  {
2346
2625
    if (m2 > 0)
2347
 
      mhi= lshift(mhi, m2);
 
2626
      mhi= lshift(mhi, m2, &alloc);
2348
2627
 
2349
2628
    /*
2350
2629
      Compute mlo -- check for special case that d is a normalized power of 2.
2353
2632
    mlo= mhi;
2354
2633
    if (spec_case)
2355
2634
    {
2356
 
      mhi= Balloc(mhi->k);
 
2635
      mhi= Balloc(mhi->k, &alloc);
2357
2636
      Bcopy(mhi, mlo);
2358
 
      mhi= lshift(mhi, Log2P);
 
2637
      mhi= lshift(mhi, Log2P, &alloc);
2359
2638
    }
2360
2639
 
2361
2640
    for (i= 1;;i++)
2363
2642
      dig= quorem(b,S) + '0';
2364
2643
      /* Do we yet have the shortest decimal string that will round to d? */
2365
2644
      j= cmp(b, mlo);
2366
 
      delta= diff(S, mhi);
 
2645
      delta= diff(S, mhi, &alloc);
2367
2646
      j1= delta->sign ? 1 : cmp(b, delta);
2368
 
      Bfree(delta);
 
2647
      Bfree(delta, &alloc);
2369
2648
      if (j1 == 0 && mode != 1 && !(word1(d) & 1)
 
2649
#ifdef Honor_FLT_ROUNDS
 
2650
          && rounding >= 1
 
2651
#endif
2370
2652
         )
2371
2653
      {
2372
2654
        if (dig == '9')
2382
2664
        {
2383
2665
          goto accept_dig;
2384
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*/
2385
2674
        if (j1 > 0)
2386
2675
        {
2387
 
          b= lshift(b, 1);
 
2676
          b= lshift(b, 1, &alloc);
2388
2677
          j1= cmp(b, S);
2389
2678
          if ((j1 > 0 || (j1 == 0 && dig & 1))
2390
2679
              && dig++ == '9')
2396
2685
      }
2397
2686
      if (j1 > 0)
2398
2687
      {
 
2688
#ifdef Honor_FLT_ROUNDS
 
2689
        if (!rounding)
 
2690
          goto accept_dig;
 
2691
#endif
2399
2692
        if (dig == '9')
2400
2693
        { /* possible if i == 1 */
2401
2694
round_9_up:
2405
2698
        *s++= dig + 1;
2406
2699
        goto ret;
2407
2700
      }
 
2701
#ifdef Honor_FLT_ROUNDS
 
2702
keep_dig:
 
2703
#endif
2408
2704
      *s++= dig;
2409
2705
      if (i == ilim)
2410
2706
        break;
2411
 
      b= multadd(b, 10, 0);
 
2707
      b= multadd(b, 10, 0, &alloc);
2412
2708
      if (mlo == mhi)
2413
 
        mlo= mhi= multadd(mhi, 10, 0);
 
2709
        mlo= mhi= multadd(mhi, 10, 0, &alloc);
2414
2710
      else
2415
2711
      {
2416
 
        mlo= multadd(mlo, 10, 0);
2417
 
        mhi= multadd(mhi, 10, 0);
 
2712
        mlo= multadd(mlo, 10, 0, &alloc);
 
2713
        mhi= multadd(mhi, 10, 0, &alloc);
2418
2714
      }
2419
2715
    }
2420
2716
  }
2428
2724
      }
2429
2725
      if (i >= ilim)
2430
2726
        break;
2431
 
      b= multadd(b, 10, 0);
 
2727
      b= multadd(b, 10, 0, &alloc);
2432
2728
    }
2433
2729
 
2434
2730
  /* Round off last digit */
2435
2731
 
2436
 
  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);
2437
2739
  j= cmp(b, S);
2438
2740
  if (j > 0 || (j == 0 && dig & 1))
2439
2741
  {
2449
2751
  }
2450
2752
  else
2451
2753
  {
2452
 
    while (*--s == '0') {}
 
2754
#ifdef Honor_FLT_ROUNDS
 
2755
trimzeros:
 
2756
#endif
 
2757
    while (*--s == '0');
2453
2758
    s++;
2454
2759
  }
2455
2760
ret:
2456
 
  Bfree(S);
 
2761
  Bfree(S, &alloc);
2457
2762
  if (mhi)
2458
2763
  {
2459
2764
    if (mlo && mlo != mhi)
2460
 
      Bfree(mlo);
2461
 
    Bfree(mhi);
 
2765
      Bfree(mlo, &alloc);
 
2766
    Bfree(mhi, &alloc);
2462
2767
  }
2463
2768
ret1:
2464
 
  Bfree(b);
 
2769
  Bfree(b, &alloc);
2465
2770
  *s= 0;
2466
2771
  *decpt= k + 1;
2467
2772
  if (rve)