← Back to branch summary

~drizzle-trunk/drizzle/development

« back to all changes in this revision

Viewing changes to mystrings/dtoa.cc

  • Committer: Eric Day
  • Date: 2009-08-27 07:26:22 UTC
  • mto: This revision was merged to the branch mainline in revision 1131.
  • Revision ID: eday@oddments.org-20090827072622-72te13ua0wdlc2ky
Reworked listen interface to not require binding of TCP ports.

Show diffs side-by-side

added added

removed removed

Lines of Context:
mystrings/dtoa.cc
35
35
 
36
36
 ***************************************************************/
37
37
 
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 *))
 
38
#include <mystrings/m_string.h>  /* for memcpy and NOT_FIXED_DEC */
 
39
#include <stdlib.h>
 
40
 
 
41
#include <algorithm>
 
42
 
 
43
using namespace std;
47
44
 
48
45
/* Magic value returned by dtoa() to indicate overflow */
49
46
#define DTOA_OVERFLOW 9999
50
47
 
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);
 
48
static double my_strtod_int(const char *, char **, int *);
 
49
static char *dtoa(double, int, int, int *, int *, char **);
54
50
 
55
51
/**
56
52
   @brief
87
83
{
88
84
  int decpt, sign, len, i;
89
85
  char *res, *src, *end, *dst= to;
90
 
  char buf[DTOA_BUFF_SIZE];
91
86
  assert(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
92
 
  
93
 
  res= dtoa(x, 5, precision, &decpt, &sign, &end, buf, sizeof(buf));
 
87
 
 
88
  res= dtoa(x, 5, precision, &decpt, &sign, &end);
94
89
 
95
90
  if (decpt == DTOA_OVERFLOW)
96
91
  {
97
 
    dtoa_free(res, buf, sizeof(buf));
 
92
    free(res);
98
93
    *to++= '0';
99
94
    *to= '\0';
100
95
    if (error != NULL)
129
124
  {
130
125
    if (len <= decpt)
131
126
      *dst++= '.';
132
 
    
133
 
    for (i= precision - cmax(0, (len - decpt)); i > 0; i--)
 
127
 
 
128
    for (i= precision - max(0, (len - decpt)); i > 0; i--)
134
129
      *dst++= '0';
135
130
  }
136
 
  
 
131
 
137
132
  *dst= '\0';
138
133
  if (error != NULL)
139
134
    *error= false;
140
135
 
141
 
  dtoa_free(res, buf, sizeof(buf));
 
136
  free(res);
142
137
 
143
138
  return dst - to;
144
139
}
196
191
     my_gcvt(55, ..., 1, ...);
197
192
 
198
193
   We do our best to minimize such cases by:
199
 
   
 
194
 
200
195
   - passing to dtoa() the field width as the number of significant digits
201
 
   
 
196
 
202
197
   - removing the sign of the number early (and decreasing the width before
203
198
     passing it to dtoa())
204
 
   
 
199
 
205
200
   - choosing the proper format to preserve the most number of significant
206
201
     digits.
207
202
*/
211
206
{
212
207
  int decpt, sign, len, exp_len;
213
208
  char *res, *src, *end, *dst= to, *dend= dst + width;
214
 
  char buf[DTOA_BUFF_SIZE];
215
209
  bool have_space, force_e_format;
216
210
  assert(width > 0 && to != NULL);
217
 
  
 
211
 
218
212
  /* We want to remove '-' from equations early */
219
213
  if (x < 0.)
220
214
    width--;
221
215
 
222
 
  res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? width : cmin(width, FLT_DIG),
223
 
            &decpt, &sign, &end, buf, sizeof(buf));
 
216
  res= dtoa(x, 4, type == MY_GCVT_ARG_DOUBLE ? min(width, DBL_DIG) : 
 
217
            min(width, FLT_DIG), &decpt, &sign, &end);
 
218
 
224
219
  if (decpt == DTOA_OVERFLOW)
225
220
  {
226
 
    dtoa_free(res, buf, sizeof(buf));
 
221
    free(res);
227
222
    *to++= '0';
228
223
    *to= '\0';
229
224
    if (error != NULL)
243
238
     to count it here.
244
239
   */
245
240
  exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
246
 
  
 
241
 
247
242
  /*
248
243
     Do we have enough space for all digits in the 'f' format?
249
244
     Let 'len' be the number of significant digits returned by dtoa,
296
291
       ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
297
292
                                            (len > 1 || !force_e_format)))) &&
298
293
         !force_e_format)) &&
299
 
      
 
294
 
300
295
       /*
301
296
         Use the 'e' format in some cases even if we have enough space for the
302
297
         'f' one. See comment for MAX_DECPT_FOR_F_FORMAT.
318
313
          *error= true;
319
314
        width= decpt;
320
315
      }
321
 
      
 
316
 
322
317
      /*
323
318
        We want to truncate (len - width) least significant digits after the
324
319
        decimal point. For this we are calling dtoa with mode=5, passing the
325
320
        number of significant digits = (len-decpt) - (len-width) = width-decpt
326
321
      */
327
 
      dtoa_free(res, buf, sizeof(buf));
328
 
      res= dtoa(x, 5, width - decpt, &decpt, &sign, &end, buf, sizeof(buf));
 
322
      free(res);
 
323
      res= dtoa(x, 5, width - decpt, &decpt, &sign, &end);
329
324
      src= res;
330
325
      len= end - res;
331
326
    }
336
331
      *dst++= '0';
337
332
      goto end;
338
333
    }
339
 
    
 
334
 
340
335
    /*
341
336
      At this point we are sure we have enough space to put all digits
342
337
      returned by dtoa
385
380
        *error= true;
386
381
      width= 0;
387
382
    }
388
 
      
 
383
 
389
384
    /* Do we have to truncate any digits? */
390
385
    if (width < len)
391
386
    {
392
387
      /* Yes, re-convert with a smaller width */
393
 
      dtoa_free(res, buf, sizeof(buf));
394
 
      res= dtoa(x, 4, width, &decpt, &sign, &end, buf, sizeof(buf));
 
388
      free(res);
 
389
      res= dtoa(x, 4, width, &decpt, &sign, &end);
395
390
      src= res;
396
391
      len= end - res;
397
392
      if (--decpt < 0)
431
426
  }
432
427
 
433
428
end:
434
 
  dtoa_free(res, buf, sizeof(buf));
 
429
  free(res);
435
430
  *dst= '\0';
436
431
 
437
432
  return dst - to;
450
445
                  rejected character.
451
446
   @param error   Upon return is set to EOVERFLOW in case of underflow or
452
447
                  overflow.
453
 
   
 
448
 
454
449
   @return        The resulting double value. In case of underflow, 0.0 is
455
450
                  returned. In case overflow, signed DBL_MAX is returned.
456
451
*/
457
452
 
458
453
double my_strtod(const char *str, char **end, int *error)
459
454
{
460
 
  char buf[DTOA_BUFF_SIZE];
461
455
  double res;
462
456
  assert(str != NULL && end != NULL && *end != NULL && error != NULL);
463
457
 
464
 
  res= my_strtod_int(str, end, error, buf, sizeof(buf));
 
458
  res= my_strtod_int(str, end, error);
465
459
  return (*error == 0) ? res : (res < 0 ? -DBL_MAX : DBL_MAX);
466
460
}
467
461
 
528
522
  file.
529
523
*/
530
524
 
531
 
/*
532
 
  #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
533
 
       and dtoa should round accordingly.
534
 
  #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
535
 
       and Honor_FLT_ROUNDS is not #defined.
536
 
 
537
 
  TODO: check if we can get rid of the above two
538
 
*/
539
 
 
540
525
typedef int32_t Long;
541
526
typedef uint32_t ULong;
542
527
typedef int64_t LLong;
592
577
#endif
593
578
#endif /*Flt_Rounds*/
594
579
 
595
 
#ifdef Honor_FLT_ROUNDS
596
 
#define Rounding rounding
597
 
#undef Check_FLT_ROUNDS
598
 
#define Check_FLT_ROUNDS
599
 
#else
600
 
#define Rounding Flt_Rounds
601
 
#endif
602
 
 
603
580
#define rounded_product(a,b) a*= b
604
581
#define rounded_quotient(a,b) a/= b
605
582
 
607
584
#define Big1 0xffffffff
608
585
#define FFFFFFFF 0xffffffffUL
609
586
 
610
 
/* This is tested to be enough for dtoa */
611
 
 
612
 
#define Kmax 15
613
 
 
614
 
#define Bcopy(x,y) memcpy(&x->sign, &y->sign, \
615
 
                          2*sizeof(int) + y->wds*sizeof(ULong))
616
 
 
617
587
/* Arbitrary-length integer */
618
588
 
619
589
typedef struct Bigint
628
598
  int wds;                 /* current length in 32-bit words */
629
599
} Bigint;
630
600
 
631
 
 
632
 
/* A simple stack-memory based allocator for Bigints */
633
 
 
634
 
typedef struct Stack_alloc
 
601
static Bigint *Bcopy(Bigint* dst, Bigint* src)
635
602
{
636
 
  char *begin;
637
 
  char *free;
638
 
  char *end;
639
 
  /*
640
 
    Having list of free blocks lets us reduce maximum required amount
641
 
    of memory from ~4000 bytes to < 1680 (tested on x86).
642
 
  */
643
 
  Bigint *freelist[Kmax+1];
644
 
} Stack_alloc;
645
 
 
646
 
 
647
 
/*
648
 
  Try to allocate object on stack, and resort to malloc if all
649
 
  stack memory is used.
650
 
*/
651
 
 
652
 
static Bigint *Balloc(int k, Stack_alloc *alloc)
 
603
  dst->sign= src->sign;
 
604
  dst->wds= src->wds;
 
605
 
 
606
  assert(dst->maxwds >= src->wds);
 
607
 
 
608
  memcpy(dst->p.x, src->p.x, src->wds*sizeof(ULong));
 
609
 
 
610
  return dst;
 
611
}
 
612
 
 
613
static Bigint *Balloc(int k)
653
614
{
654
615
  Bigint *rv;
655
 
  if (k <= Kmax &&  alloc->freelist[k])
656
 
  {
657
 
    rv= alloc->freelist[k];
658
 
    alloc->freelist[k]= rv->p.next;
659
 
  }
660
 
  else
661
 
  {
662
 
    int x, len;
663
 
 
664
 
    x= 1 << k;
665
 
    len= sizeof(Bigint) + x * sizeof(ULong);
666
 
 
667
 
    if (alloc->free + len <= alloc->end)
668
 
    {
669
 
      rv= (Bigint*) alloc->free;
670
 
      alloc->free+= len;
671
 
    }
672
 
    else
673
 
      rv= (Bigint*) malloc(len);
674
 
 
675
 
    rv->k= k;
676
 
    rv->maxwds= x;
677
 
  }
 
616
 
 
617
  /* TODO: some malloc failure checking */
 
618
 
 
619
  int x= 1 << k;
 
620
  rv= (Bigint*) malloc(sizeof(Bigint));
 
621
 
 
622
  rv->p.x= (ULong*)malloc(x * sizeof(ULong));
 
623
 
 
624
  rv->k= k;
 
625
  rv->maxwds= x;
 
626
 
678
627
  rv->sign= rv->wds= 0;
679
 
  rv->p.x= (ULong*) (rv + 1);
 
628
 
680
629
  return rv;
681
630
}
682
631
 
686
635
  list. Otherwise call free().
687
636
*/
688
637
 
689
 
static void Bfree(Bigint *v, Stack_alloc *alloc)
690
 
{
691
 
  char *gptr= (char*) v;                       /* generic pointer */
692
 
  if (gptr < alloc->begin || gptr >= alloc->end)
693
 
    free(gptr);
694
 
  else if (v->k <= Kmax)
695
 
  {
696
 
    /*
697
 
      Maintain free lists only for stack objects: this way we don't
698
 
      have to bother with freeing lists in the end of dtoa;
699
 
      heap should not be used normally anyway.
700
 
    */
701
 
    v->p.next= alloc->freelist[v->k];
702
 
    alloc->freelist[v->k]= v;
703
 
  }
704
 
}
705
 
 
706
 
 
707
 
/*
708
 
  This is to place return value of dtoa in: tries to use stack
709
 
  as well, but passes by free lists management and just aligns len by
710
 
  sizeof(ULong).
711
 
*/
712
 
 
713
 
static char *dtoa_alloc(int i, Stack_alloc *alloc)
714
 
{
715
 
  char *rv;
716
 
  int aligned_size= (i + sizeof(ULong) - 1) / sizeof(ULong) * sizeof(ULong);
717
 
  if (alloc->free + aligned_size <= alloc->end)
718
 
  {
719
 
    rv= alloc->free;
720
 
    alloc->free+= aligned_size;
721
 
  }
722
 
  else
723
 
    rv= malloc(i);
724
 
  return rv;
725
 
}
726
 
 
727
 
 
728
 
/*
729
 
  dtoa_free() must be used to free values s returned by dtoa()
730
 
  This is the counterpart of dtoa_alloc()
731
 
*/
732
 
 
733
 
static void dtoa_free(char *gptr, char *buf, size_t buf_size)
734
 
{
735
 
  if (gptr < buf || gptr >= buf + buf_size)
736
 
    free(gptr);
737
 
}
738
 
 
 
638
static void Bfree(Bigint *v)
 
639
{
 
640
  if(!v)
 
641
    return;
 
642
  free(v->p.x);
 
643
  free(v);
 
644
}
739
645
 
740
646
/* Bigint arithmetic functions */
741
647
 
742
648
/* Multiply by m and add a */
743
649
 
744
 
static Bigint *multadd(Bigint *b, int m, int a, Stack_alloc *alloc)
 
650
static Bigint *multadd(Bigint *b, int m, int a)
745
651
{
746
652
  int i, wds;
747
653
  ULong *x;
763
669
  {
764
670
    if (wds >= b->maxwds)
765
671
    {
766
 
      b1= Balloc(b->k+1, alloc);
 
672
      b1= Balloc(b->k+1);
767
673
      Bcopy(b1, b);
768
 
      Bfree(b, alloc);
 
674
      Bfree(b);
769
675
      b= b1;
770
676
    }
771
677
    b->p.x[wds++]= (ULong) carry;
775
681
}
776
682
 
777
683
 
778
 
static Bigint *s2b(const char *s, int nd0, int nd, ULong y9, Stack_alloc *alloc)
 
684
static Bigint *s2b(const char *s, int nd0, int nd, ULong y9)
779
685
{
780
686
  Bigint *b;
781
687
  int i, k;
783
689
 
784
690
  x= (nd + 8) / 9;
785
691
  for (k= 0, y= 1; x > y; y <<= 1, k++) ;
786
 
  b= Balloc(k, alloc);
 
692
  b= Balloc(k);
787
693
  b->p.x[0]= y9;
788
694
  b->wds= 1;
789
 
  
 
695
 
790
696
  i= 9;
791
697
  if (9 < nd0)
792
698
  {
793
699
    s+= 9;
794
700
    do
795
 
      b= multadd(b, 10, *s++ - '0', alloc);
 
701
      b= multadd(b, 10, *s++ - '0');
796
702
    while (++i < nd0);
797
703
    s++;
798
704
  }
799
705
  else
800
706
    s+= 10;
801
707
  for(; i < nd; i++)
802
 
    b= multadd(b, 10, *s++ - '0', alloc);
 
708
    b= multadd(b, 10, *s++ - '0');
803
709
  return b;
804
710
}
805
711
 
890
796
 
891
797
/* Convert integer to Bigint number */
892
798
 
893
 
static Bigint *i2b(int i, Stack_alloc *alloc)
 
799
static Bigint *i2b(int i)
894
800
{
895
801
  Bigint *b;
896
802
 
897
 
  b= Balloc(1, alloc);
 
803
  b= Balloc(1);
898
804
  b->p.x[0]= i;
899
805
  b->wds= 1;
900
806
  return b;
903
809
 
904
810
/* Multiply two Bigint numbers */
905
811
 
906
 
static Bigint *mult(Bigint *a, Bigint *b, Stack_alloc *alloc)
 
812
static Bigint *mult(Bigint *a, Bigint *b)
907
813
{
908
814
  Bigint *c;
909
815
  int k, wa, wb, wc;
923
829
  wc= wa + wb;
924
830
  if (wc > a->maxwds)
925
831
    k++;
926
 
  c= Balloc(k, alloc);
 
832
  c= Balloc(k);
927
833
  for (x= c->p.x, xa= x + wc; x < xa; x++)
928
834
    *x= 0;
929
835
  xa= a->p.x;
995
901
 
996
902
#define P5A_MAX (sizeof(p5_a)/sizeof(*p5_a) - 1)
997
903
 
998
 
static Bigint *pow5mult(Bigint *b, int k, Stack_alloc *alloc)
 
904
static Bigint *pow5mult(Bigint *b, int k)
999
905
{
1000
906
  Bigint *b1, *p5, *p51;
1001
907
  int i;
1002
908
  static int p05[3]= { 5, 25, 125 };
1003
909
 
1004
910
  if ((i= k & 3))
1005
 
    b= multadd(b, p05[i-1], 0, alloc);
 
911
    b= multadd(b, p05[i-1], 0);
1006
912
 
1007
913
  if (!(k>>= 2))
1008
914
    return b;
1011
917
  {
1012
918
    if (k & 1)
1013
919
    {
1014
 
      b1= mult(b, p5, alloc);
1015
 
      Bfree(b, alloc);
 
920
      b1= mult(b, p5);
 
921
      Bfree(b);
1016
922
      b= b1;
1017
923
    }
1018
924
    if (!(k>>= 1))
1021
927
    if (p5 < p5_a + P5A_MAX)
1022
928
      ++p5;
1023
929
    else if (p5 == p5_a + P5A_MAX)
1024
 
      p5= mult(p5, p5, alloc);
 
930
      p5= mult(p5, p5);
1025
931
    else
1026
932
    {
1027
 
      p51= mult(p5, p5, alloc);
1028
 
      Bfree(p5, alloc);
 
933
      p51= mult(p5, p5);
 
934
      Bfree(p5);
1029
935
      p5= p51;
1030
936
    }
1031
937
  }
1033
939
}
1034
940
 
1035
941
 
1036
 
static Bigint *lshift(Bigint *b, int k, Stack_alloc *alloc)
 
942
static Bigint *lshift(Bigint *b, int k)
1037
943
{
1038
944
  int i, k1, n, n1;
1039
945
  Bigint *b1;
1044
950
  n1= n + b->wds + 1;
1045
951
  for (i= b->maxwds; n1 > i; i<<= 1)
1046
952
    k1++;
1047
 
  b1= Balloc(k1, alloc);
 
953
  b1= Balloc(k1);
1048
954
  x1= b1->p.x;
1049
955
  for (i= 0; i < n; i++)
1050
956
    *x1++= 0;
1068
974
      *x1++= *x++;
1069
975
    while (x < xe);
1070
976
  b1->wds= n1 - 1;
1071
 
  Bfree(b, alloc);
 
977
  Bfree(b);
1072
978
  return b1;
1073
979
}
1074
980
 
1097
1003
}
1098
1004
 
1099
1005
 
1100
 
static Bigint *diff(Bigint *a, Bigint *b, Stack_alloc *alloc)
 
1006
static Bigint *diff(Bigint *a, Bigint *b)
1101
1007
{
1102
1008
  Bigint *c;
1103
1009
  int i, wa, wb;
1107
1013
  i= cmp(a,b);
1108
1014
  if (!i)
1109
1015
  {
1110
 
    c= Balloc(0, alloc);
 
1016
    c= Balloc(0);
1111
1017
    c->wds= 1;
1112
1018
    c->p.x[0]= 0;
1113
1019
    return c;
1121
1027
  }
1122
1028
  else
1123
1029
    i= 0;
1124
 
  c= Balloc(a->k, alloc);
 
1030
  c= Balloc(a->k);
1125
1031
  c->sign= i;
1126
1032
  wa= a->wds;
1127
1033
  xa= a->p.x;
1202
1108
}
1203
1109
 
1204
1110
 
1205
 
static Bigint *d2b(double d, int *e, int *bits, Stack_alloc *alloc)
 
1111
static Bigint *d2b(double d, int *e, int *bits)
1206
1112
{
1207
1113
  Bigint *b;
1208
1114
  int de, k;
1211
1117
#define d0 word0(d)
1212
1118
#define d1 word1(d)
1213
1119
 
1214
 
  b= Balloc(1, alloc);
 
1120
  b= Balloc(1);
1215
1121
  x= b->p.x;
1216
1122
 
1217
1123
  z= d0 & Frac_mask;
1283
1189
  9007199254740992.*9007199254740992.e-256 /* = 2^106 * 1e-53 */
1284
1190
};
1285
1191
/*
1286
 
  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow 
 
1192
  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow
1287
1193
  flag unnecessarily.  It leads to a song and dance at the end of strtod.
1288
1194
*/
1289
1195
#define Scale_Bit 0x10
1291
1197
 
1292
1198
/*
1293
1199
  strtod for IEEE--arithmetic machines.
1294
 
 
 
1200
 
1295
1201
  This strtod returns a nearest machine number to the input decimal
1296
1202
  string (or sets errno to EOVERFLOW). Ties are broken by the IEEE round-even
1297
1203
  rule.
1298
 
 
 
1204
 
1299
1205
  Inspired loosely by William D. Clinger's paper "How to Read Floating
1300
1206
  Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
1301
 
 
 
1207
 
1302
1208
  Modifications:
1303
 
 
 
1209
 
1304
1210
   1. We only require IEEE (not IEEE double-extended).
1305
1211
   2. We get by with floating-point arithmetic in a case that
1306
1212
     Clinger missed -- when we're computing d * 10^n
1318
1224
     for 0 <= k <= 22).
1319
1225
*/
1320
1226
 
1321
 
static double my_strtod_int(const char *s00, char **se, int *error, char *buf, size_t buf_size)
 
1227
static double my_strtod_int(const char *s00, char **se, int *error)
1322
1228
{
1323
1229
  int scale;
1324
1230
  int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1331
1237
#ifdef SET_INEXACT
1332
1238
  int inexact, oldinexact;
1333
1239
#endif
1334
 
#ifdef Honor_FLT_ROUNDS
1335
 
  int rounding;
1336
 
#endif
1337
 
  Stack_alloc alloc;
1338
1240
 
1339
1241
  c= 0;
1340
1242
  *error= 0;
1341
1243
 
1342
 
  alloc.begin= alloc.free= buf;
1343
 
  alloc.end= buf + buf_size;
1344
 
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
1345
 
 
1346
1244
  sign= nz0= nz= 0;
1347
1245
  dval(rv)= 0.;
1348
1246
  for (s= s00; s < end; s++)
1366
1264
 break2:
1367
1265
  if (s >= end)
1368
1266
    goto ret0;
1369
 
  
 
1267
 
1370
1268
  if (*s == '0')
1371
1269
  {
1372
1270
    nz0= 1;
1491
1389
  }
1492
1390
  bd0= 0;
1493
1391
  if (nd <= DBL_DIG
1494
 
#ifndef Honor_FLT_ROUNDS
1495
 
    && Flt_Rounds == 1
1496
 
#endif
1497
1392
      )
1498
1393
  {
1499
1394
    if (!e)
1502
1397
    {
1503
1398
      if (e <= Ten_pmax)
1504
1399
      {
1505
 
#ifdef Honor_FLT_ROUNDS
1506
 
        /* round correctly FLT_ROUNDS = 2 or 3 */
1507
 
        if (sign)
1508
 
        {
1509
 
          rv= -rv;
1510
 
          sign= 0;
1511
 
        }
1512
 
#endif
1513
1400
        /* rv = */ rounded_product(dval(rv), tens[e]);
1514
1401
        goto ret;
1515
1402
      }
1520
1407
          A fancier test would sometimes let us do
1521
1408
          this for larger i values.
1522
1409
        */
1523
 
#ifdef Honor_FLT_ROUNDS
1524
 
        /* round correctly FLT_ROUNDS = 2 or 3 */
1525
 
        if (sign)
1526
 
        {
1527
 
          rv= -rv;
1528
 
          sign= 0;
1529
 
        }
1530
 
#endif
1531
1410
        e-= i;
1532
1411
        dval(rv)*= tens[i];
1533
1412
        /* rv = */ rounded_product(dval(rv), tens[e]);
1537
1416
#ifndef Inaccurate_Divide
1538
1417
    else if (e >= -Ten_pmax)
1539
1418
    {
1540
 
#ifdef Honor_FLT_ROUNDS
1541
 
      /* round correctly FLT_ROUNDS = 2 or 3 */
1542
 
      if (sign)
1543
 
      {
1544
 
        rv= -rv;
1545
 
        sign= 0;
1546
 
      }
1547
 
#endif
1548
1419
      /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1549
1420
      goto ret;
1550
1421
    }
1558
1429
    oldinexact= get_inexact();
1559
1430
#endif
1560
1431
  scale= 0;
1561
 
#ifdef Honor_FLT_ROUNDS
1562
 
  if ((rounding= Flt_Rounds) >= 2)
1563
 
  {
1564
 
    if (sign)
1565
 
      rounding= rounding == 2 ? 0 : 2;
1566
 
    else
1567
 
      if (rounding != 2)
1568
 
        rounding= 0;
1569
 
  }
1570
 
#endif
1571
1432
 
1572
1433
  /* Get starting approximation = rv * 10**e1 */
1573
1434
 
1582
1443
 ovfl:
1583
1444
        *error= EOVERFLOW;
1584
1445
        /* Can't trust HUGE_VAL */
1585
 
#ifdef Honor_FLT_ROUNDS
1586
 
        switch (rounding)
1587
 
        {
1588
 
        case 0: /* toward 0 */
1589
 
        case 3: /* toward -infinity */
1590
 
          word0(rv)= Big0;
1591
 
          word1(rv)= Big1;
1592
 
          break;
1593
 
        default:
1594
 
          word0(rv)= Exp_mask;
1595
 
          word1(rv)= 0;
1596
 
        }
1597
 
#else /*Honor_FLT_ROUNDS*/
1598
1446
        word0(rv)= Exp_mask;
1599
1447
        word1(rv)= 0;
1600
 
#endif /*Honor_FLT_ROUNDS*/
1601
1448
#ifdef SET_INEXACT
1602
1449
        /* set overflow bit */
1603
1450
        dval(rv0)= 1e300;
1669
1516
 
1670
1517
  /* Put digits into bd: true value = bd * 10^e */
1671
1518
 
1672
 
  bd0= s2b(s0, nd0, nd, y, &alloc);
 
1519
  bd0= s2b(s0, nd0, nd, y);
1673
1520
 
1674
1521
  for(;;)
1675
1522
  {
1676
 
    bd= Balloc(bd0->k, &alloc);
 
1523
    bd= Balloc(bd0->k);
1677
1524
    Bcopy(bd, bd0);
1678
 
    bb= d2b(dval(rv), &bbe, &bbbits, &alloc);  /* rv = bb * 2^bbe */
1679
 
    bs= i2b(1, &alloc);
 
1525
    bb= d2b(dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
 
1526
    bs= i2b(1);
1680
1527
 
1681
1528
    if (e >= 0)
1682
1529
    {
1693
1540
    else
1694
1541
      bd2-= bbe;
1695
1542
    bs2= bb2;
1696
 
#ifdef Honor_FLT_ROUNDS
1697
 
    if (rounding != 1)
1698
 
      bs2++;
1699
 
#endif
1700
1543
    j= bbe - scale;
1701
1544
    i= j + bbbits - 1;  /* logb(rv) */
1702
1545
    if (i < Emin)  /* denormal */
1717
1560
    }
1718
1561
    if (bb5 > 0)
1719
1562
    {
1720
 
      bs= pow5mult(bs, bb5, &alloc);
1721
 
      bb1= mult(bs, bb, &alloc);
1722
 
      Bfree(bb, &alloc);
 
1563
      bs= pow5mult(bs, bb5);
 
1564
      bb1= mult(bs, bb);
 
1565
      Bfree(bb);
1723
1566
      bb= bb1;
1724
1567
    }
1725
1568
    if (bb2 > 0)
1726
 
      bb= lshift(bb, bb2, &alloc);
 
1569
      bb= lshift(bb, bb2);
1727
1570
    if (bd5 > 0)
1728
 
      bd= pow5mult(bd, bd5, &alloc);
 
1571
      bd= pow5mult(bd, bd5);
1729
1572
    if (bd2 > 0)
1730
 
      bd= lshift(bd, bd2, &alloc);
 
1573
      bd= lshift(bd, bd2);
1731
1574
    if (bs2 > 0)
1732
 
      bs= lshift(bs, bs2, &alloc);
1733
 
    delta= diff(bb, bd, &alloc);
 
1575
      bs= lshift(bs, bs2);
 
1576
    delta= diff(bb, bd);
1734
1577
    dsign= delta->sign;
1735
1578
    delta->sign= 0;
1736
1579
    i= cmp(delta, bs);
1737
 
#ifdef Honor_FLT_ROUNDS
1738
 
    if (rounding != 1)
1739
 
    {
1740
 
      if (i < 0)
1741
 
      {
1742
 
        /* Error is less than an ulp */
1743
 
        if (!delta->x[0] && delta->wds <= 1)
1744
 
        {
1745
 
          /* exact */
1746
 
#ifdef SET_INEXACT
1747
 
          inexact= 0;
1748
 
#endif
1749
 
          break;
1750
 
        }
1751
 
        if (rounding)
1752
 
        {
1753
 
          if (dsign)
1754
 
          {
1755
 
            adj= 1.;
1756
 
            goto apply_adj;
1757
 
          }
1758
 
        }
1759
 
        else if (!dsign)
1760
 
        {
1761
 
          adj= -1.;
1762
 
          if (!word1(rv) && !(word0(rv) & Frac_mask))
1763
 
          {
1764
 
            y= word0(rv) & Exp_mask;
1765
 
            if (!scale || y > 2*P*Exp_msk1)
1766
 
            {
1767
 
              delta= lshift(delta,Log2P);
1768
 
              if (cmp(delta, bs) <= 0)
1769
 
              adj= -0.5;
1770
 
            }
1771
 
          }
1772
 
 apply_adj:
1773
 
          if (scale && (y= word0(rv) & Exp_mask) <= 2 * P * Exp_msk1)
1774
 
            word0(adj)+= (2 * P + 1) * Exp_msk1 - y;
1775
 
          dval(rv)+= adj * ulp(dval(rv));
1776
 
        }
1777
 
        break;
1778
 
      }
1779
 
      adj= ratio(delta, bs);
1780
 
      if (adj < 1.)
1781
 
        adj= 1.;
1782
 
      if (adj <= 0x7ffffffe)
1783
 
      {
1784
 
        /* adj = rounding ? ceil(adj) : floor(adj); */
1785
 
        y= adj;
1786
 
        if (y != adj)
1787
 
        {
1788
 
          if (!((rounding >> 1) ^ dsign))
1789
 
            y++;
1790
 
          adj= y;
1791
 
        }
1792
 
      }
1793
 
      if (scale && (y= word0(rv) & Exp_mask) <= 2 * P * Exp_msk1)
1794
 
        word0(adj)+= (2 * P + 1) * Exp_msk1 - y;
1795
 
      adj*= ulp(dval(rv));
1796
 
      if (dsign)
1797
 
        dval(rv)+= adj;
1798
 
      else
1799
 
        dval(rv)-= adj;
1800
 
      goto cont;
1801
 
    }
1802
 
#endif /*Honor_FLT_ROUNDS*/
1803
1580
 
1804
1581
    if (i < 0)
1805
1582
    {
1824
1601
#endif
1825
1602
        break;
1826
1603
      }
1827
 
      delta= lshift(delta, Log2P, &alloc);
 
1604
      delta= lshift(delta, Log2P);
1828
1605
      if (cmp(delta, bs) > 0)
1829
1606
        goto drop_down;
1830
1607
      break;
1906
1683
    {
1907
1684
      aadj*= 0.5;
1908
1685
      aadj1= dsign ? aadj : -aadj;
1909
 
#ifdef Check_FLT_ROUNDS
1910
 
      switch (Rounding)
1911
 
      {
1912
 
      case 2: /* towards +infinity */
1913
 
        aadj1-= 0.5;
1914
 
        break;
1915
 
      case 0: /* towards 0 */
1916
 
      case 3: /* towards -infinity */
1917
 
        aadj1+= 0.5;
1918
 
      }
1919
 
#else
1920
1686
      if (Flt_Rounds == 0)
1921
1687
        aadj1+= 0.5;
1922
 
#endif /*Check_FLT_ROUNDS*/
1923
1688
    }
1924
1689
    y= word0(rv) & Exp_mask;
1925
1690
 
1977
1742
      }
1978
1743
#endif
1979
1744
 cont:
1980
 
    Bfree(bb, &alloc);
1981
 
    Bfree(bd, &alloc);
1982
 
    Bfree(bs, &alloc);
1983
 
    Bfree(delta, &alloc);
 
1745
    Bfree(bb);
 
1746
    Bfree(bd);
 
1747
    Bfree(bs);
 
1748
    Bfree(delta);
1984
1749
  }
1985
1750
#ifdef SET_INEXACT
1986
1751
  if (inexact)
2010
1775
  }
2011
1776
#endif
2012
1777
 retfree:
2013
 
  Bfree(bb, &alloc);
2014
 
  Bfree(bd, &alloc);
2015
 
  Bfree(bs, &alloc);
2016
 
  Bfree(bd0, &alloc);
2017
 
  Bfree(delta, &alloc);
 
1778
  Bfree(bb);
 
1779
  Bfree(bd);
 
1780
  Bfree(bs);
 
1781
  Bfree(bd0);
 
1782
  Bfree(delta);
2018
1783
 ret:
2019
1784
  *se= (char *)s;
2020
1785
  return sign ? -dval(rv) : dval(rv);
2121
1886
 */
2122
1887
 
2123
1888
static char *dtoa(double d, int mode, int ndigits, int *decpt, int *sign,
2124
 
                  char **rve, char *buf, size_t buf_size)
 
1889
                  char **rve)
2125
1890
{
2126
1891
  /*
2127
1892
    Arguments ndigits, decpt, sign are similar to those
2133
1898
    mode:
2134
1899
          0 ==> shortest string that yields d when read in
2135
1900
                and rounded to nearest.
 
1901
 
2136
1902
          1 ==> like 0, but with Steele & White stopping rule;
2137
1903
                e.g. with IEEE P754 arithmetic , mode 0 gives
2138
1904
                1e23 whereas mode 1 gives 9.999999999999999e22.
2146
1912
          4,5 ==> similar to 2 and 3, respectively, but (in
2147
1913
                round-nearest mode) with the tests of mode 0 to
2148
1914
                possibly return a shorter string that rounds to d.
2149
 
                With IEEE arithmetic and compilation with
2150
 
                -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2151
 
                as modes 2 and 3 when FLT_ROUNDS != 1.
2152
1915
          6-9 ==> Debugging modes similar to mode - 4:  don't try
2153
1916
                fast floating-point estimate (if applicable).
2154
1917
 
2167
1930
  Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
2168
1931
  double d2, ds, eps;
2169
1932
  char *s, *s0;
2170
 
#ifdef Honor_FLT_ROUNDS
2171
 
  int rounding;
2172
 
#endif
2173
 
  Stack_alloc alloc;
2174
 
  
2175
 
  alloc.begin= alloc.free= buf;
2176
 
  alloc.end= buf + buf_size;
2177
 
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
2178
1933
 
2179
1934
  if (word0(d) & Sign_bit)
2180
1935
  {
2190
1945
      (!dval(d) && (*decpt= 1)))
2191
1946
  {
2192
1947
    /* Infinity, NaN, 0 */
2193
 
    char *res= (char*) dtoa_alloc(2, &alloc);
 
1948
    char *res= (char*) malloc(2);
2194
1949
    res[0]= '0';
2195
1950
    res[1]= '\0';
2196
1951
    if (rve)
2197
1952
      *rve= res + 1;
2198
1953
    return res;
2199
1954
  }
2200
 
  
2201
 
#ifdef Honor_FLT_ROUNDS
2202
 
  if ((rounding= Flt_Rounds) >= 2)
2203
 
  {
2204
 
    if (*sign)
2205
 
      rounding= rounding == 2 ? 0 : 2;
2206
 
    else
2207
 
      if (rounding != 2)
2208
 
        rounding= 0;
2209
 
  }
2210
 
#endif
2211
 
 
2212
 
  b= d2b(dval(d), &be, &bbits, &alloc);
 
1955
 
 
1956
 
 
1957
  b= d2b(dval(d), &be, &bbits);
2213
1958
  if ((i= (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
2214
1959
  {
2215
1960
    dval(d2)= dval(d);
2221
1966
      log10(x)      =  log(x) / log(10)
2222
1967
                   ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2223
1968
      log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
2224
 
     
 
1969
 
2225
1970
      This suggests computing an approximation k to log10(d) by
2226
 
     
 
1971
 
2227
1972
      k= (i - Bias)*0.301029995663981
2228
1973
           + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2229
 
     
 
1974
 
2230
1975
      We want k to be too large rather than too small.
2231
1976
      The error in the first-order Taylor series approximation
2232
1977
      is in our favor, so we just round up the constant enough
2291
2036
  if (mode < 0 || mode > 9)
2292
2037
    mode= 0;
2293
2038
 
2294
 
#ifdef Check_FLT_ROUNDS
2295
 
  try_quick= Rounding == 1;
2296
 
#else
2297
2039
  try_quick= 1;
2298
 
#endif
2299
2040
 
2300
2041
  if (mode > 5)
2301
2042
  {
2328
2069
    if (i <= 0)
2329
2070
      i= 1;
2330
2071
  }
2331
 
  s= s0= dtoa_alloc(i, &alloc);
2332
 
 
2333
 
#ifdef Honor_FLT_ROUNDS
2334
 
  if (mode > 1 && rounding != 1)
2335
 
    leftright= 0;
2336
 
#endif
 
2072
  s= s0= (char*)malloc(i+1); /* +1 for trailing '\0' appended
 
2073
                                at end of function */
2337
2074
 
2338
2075
  if (ilim >= 0 && ilim <= Quick_max && try_quick)
2339
2076
  {
2432
2169
            goto bump_up;
2433
2170
          else if (dval(d) < 0.5 - dval(eps))
2434
2171
          {
2435
 
            while (*--s == '0');
 
2172
            while (*--s == '0') {}
2436
2173
            s++;
2437
2174
            goto ret1;
2438
2175
          }
2464
2201
    {
2465
2202
      L= (Long)(dval(d) / ds);
2466
2203
      dval(d)-= L*ds;
2467
 
#ifdef Check_FLT_ROUNDS
2468
 
      /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2469
 
      if (dval(d) < 0)
2470
 
      {
2471
 
        L--;
2472
 
        dval(d)+= ds;
2473
 
      }
2474
 
#endif
2475
2204
      *s++= '0' + (int)L;
2476
2205
      if (!dval(d))
2477
2206
      {
2479
2208
      }
2480
2209
      if (i == ilim)
2481
2210
      {
2482
 
#ifdef Honor_FLT_ROUNDS
2483
 
        if (mode > 1)
2484
 
        {
2485
 
          switch (rounding) {
2486
 
          case 0: goto ret1;
2487
 
          case 2: goto bump_up;
2488
 
          }
2489
 
        }
2490
 
#endif
2491
2211
        dval(d)+= dval(d);
2492
2212
        if (dval(d) > ds || (dval(d) == ds && L & 1))
2493
2213
        {
2515
2235
    i = denorm ? be + (Bias + (P-1) - 1 + 1) : 1 + P - bbits;
2516
2236
    b2+= i;
2517
2237
    s2+= i;
2518
 
    mhi= i2b(1, &alloc);
 
2238
    mhi= i2b(1);
2519
2239
  }
2520
2240
  if (m2 > 0 && s2 > 0)
2521
2241
  {
2530
2250
    {
2531
2251
      if (m5 > 0)
2532
2252
      {
2533
 
        mhi= pow5mult(mhi, m5, &alloc);
2534
 
        b1= mult(mhi, b, &alloc);
2535
 
        Bfree(b, &alloc);
 
2253
        mhi= pow5mult(mhi, m5);
 
2254
        b1= mult(mhi, b);
 
2255
        Bfree(b);
2536
2256
        b= b1;
2537
2257
      }
2538
2258
      if ((j= b5 - m5))
2539
 
        b= pow5mult(b, j, &alloc);
 
2259
        b= pow5mult(b, j);
2540
2260
    }
2541
2261
    else
2542
 
      b= pow5mult(b, b5, &alloc);
 
2262
      b= pow5mult(b, b5);
2543
2263
  }
2544
 
  S= i2b(1, &alloc);
 
2264
  S= i2b(1);
2545
2265
  if (s5 > 0)
2546
 
    S= pow5mult(S, s5, &alloc);
 
2266
    S= pow5mult(S, s5);
2547
2267
 
2548
2268
  /* Check for special case that d is a normalized power of 2. */
2549
2269
 
2550
2270
  spec_case= 0;
2551
2271
  if ((mode < 2 || leftright)
2552
 
#ifdef Honor_FLT_ROUNDS
2553
 
      && rounding == 1
2554
 
#endif
2555
2272
     )
2556
2273
  {
2557
2274
    if (!word1(d) && !(word0(d) & Bndry_mask) &&
2568
2285
  /*
2569
2286
    Arrange for convenient computation of quotients:
2570
2287
    shift left if necessary so divisor has 4 leading 0 bits.
2571
 
    
 
2288
 
2572
2289
    Perhaps we should just compute leading 28 bits of S once
2573
2290
    a nd for all and pass them and a shift to quorem, so it
2574
2291
    can do shifts and ors to compute the numerator for q.
2590
2307
    s2+= i;
2591
2308
  }
2592
2309
  if (b2 > 0)
2593
 
    b= lshift(b, b2, &alloc);
 
2310
    b= lshift(b, b2);
2594
2311
  if (s2 > 0)
2595
 
    S= lshift(S, s2, &alloc);
 
2312
    S= lshift(S, s2);
2596
2313
  if (k_check)
2597
2314
  {
2598
2315
    if (cmp(b,S) < 0)
2599
2316
    {
2600
2317
      k--;
2601
2318
      /* we botched the k estimate */
2602
 
      b= multadd(b, 10, 0, &alloc);
 
2319
      b= multadd(b, 10, 0);
2603
2320
      if (leftright)
2604
 
        mhi= multadd(mhi, 10, 0, &alloc);
 
2321
        mhi= multadd(mhi, 10, 0);
2605
2322
      ilim= ilim1;
2606
2323
    }
2607
2324
  }
2608
2325
  if (ilim <= 0 && (mode == 3 || mode == 5))
2609
2326
  {
2610
 
    if (ilim < 0 || cmp(b,S= multadd(S,5,0, &alloc)) <= 0)
 
2327
    if (ilim < 0 || cmp(b,S= multadd(S,5,0)) <= 0)
2611
2328
    {
2612
2329
      /* no digits, fcvt style */
2613
2330
no_digits:
2622
2339
  if (leftright)
2623
2340
  {
2624
2341
    if (m2 > 0)
2625
 
      mhi= lshift(mhi, m2, &alloc);
 
2342
      mhi= lshift(mhi, m2);
2626
2343
 
2627
2344
    /*
2628
2345
      Compute mlo -- check for special case that d is a normalized power of 2.
2631
2348
    mlo= mhi;
2632
2349
    if (spec_case)
2633
2350
    {
2634
 
      mhi= Balloc(mhi->k, &alloc);
 
2351
      mhi= Balloc(mhi->k);
2635
2352
      Bcopy(mhi, mlo);
2636
 
      mhi= lshift(mhi, Log2P, &alloc);
 
2353
      mhi= lshift(mhi, Log2P);
2637
2354
    }
2638
2355
 
2639
2356
    for (i= 1;;i++)
2641
2358
      dig= quorem(b,S) + '0';
2642
2359
      /* Do we yet have the shortest decimal string that will round to d? */
2643
2360
      j= cmp(b, mlo);
2644
 
      delta= diff(S, mhi, &alloc);
 
2361
      delta= diff(S, mhi);
2645
2362
      j1= delta->sign ? 1 : cmp(b, delta);
2646
 
      Bfree(delta, &alloc);
 
2363
      Bfree(delta);
2647
2364
      if (j1 == 0 && mode != 1 && !(word1(d) & 1)
2648
 
#ifdef Honor_FLT_ROUNDS
2649
 
          && rounding >= 1
2650
 
#endif
2651
2365
         )
2652
2366
      {
2653
2367
        if (dig == '9')
2663
2377
        {
2664
2378
          goto accept_dig;
2665
2379
        }
2666
 
#ifdef Honor_FLT_ROUNDS
2667
 
        if (mode > 1)
2668
 
          switch (rounding) {
2669
 
          case 0: goto accept_dig;
2670
 
          case 2: goto keep_dig;
2671
 
          }
2672
 
#endif /*Honor_FLT_ROUNDS*/
2673
2380
        if (j1 > 0)
2674
2381
        {
2675
 
          b= lshift(b, 1, &alloc);
 
2382
          b= lshift(b, 1);
2676
2383
          j1= cmp(b, S);
2677
2384
          if ((j1 > 0 || (j1 == 0 && dig & 1))
2678
2385
              && dig++ == '9')
2684
2391
      }
2685
2392
      if (j1 > 0)
2686
2393
      {
2687
 
#ifdef Honor_FLT_ROUNDS
2688
 
        if (!rounding)
2689
 
          goto accept_dig;
2690
 
#endif
2691
2394
        if (dig == '9')
2692
2395
        { /* possible if i == 1 */
2693
2396
round_9_up:
2697
2400
        *s++= dig + 1;
2698
2401
        goto ret;
2699
2402
      }
2700
 
#ifdef Honor_FLT_ROUNDS
2701
 
keep_dig:
2702
 
#endif
2703
2403
      *s++= dig;
2704
2404
      if (i == ilim)
2705
2405
        break;
2706
 
      b= multadd(b, 10, 0, &alloc);
 
2406
      b= multadd(b, 10, 0);
2707
2407
      if (mlo == mhi)
2708
 
        mlo= mhi= multadd(mhi, 10, 0, &alloc);
 
2408
        mlo= mhi= multadd(mhi, 10, 0);
2709
2409
      else
2710
2410
      {
2711
 
        mlo= multadd(mlo, 10, 0, &alloc);
2712
 
        mhi= multadd(mhi, 10, 0, &alloc);
 
2411
        mlo= multadd(mlo, 10, 0);
 
2412
        mhi= multadd(mhi, 10, 0);
2713
2413
      }
2714
2414
    }
2715
2415
  }
2723
2423
      }
2724
2424
      if (i >= ilim)
2725
2425
        break;
2726
 
      b= multadd(b, 10, 0, &alloc);
 
2426
      b= multadd(b, 10, 0);
2727
2427
    }
2728
2428
 
2729
2429
  /* Round off last digit */
2730
2430
 
2731
 
#ifdef Honor_FLT_ROUNDS
2732
 
  switch (rounding) {
2733
 
  case 0: goto trimzeros;
2734
 
  case 2: goto roundoff;
2735
 
  }
2736
 
#endif
2737
 
  b= lshift(b, 1, &alloc);
 
2431
  b= lshift(b, 1);
2738
2432
  j= cmp(b, S);
2739
2433
  if (j > 0 || (j == 0 && dig & 1))
2740
2434
  {
2750
2444
  }
2751
2445
  else
2752
2446
  {
2753
 
#ifdef Honor_FLT_ROUNDS
2754
 
trimzeros:
2755
 
#endif
2756
 
    while (*--s == '0');
 
2447
    while (*--s == '0') {}
2757
2448
    s++;
2758
2449
  }
2759
2450
ret:
2760
 
  Bfree(S, &alloc);
 
2451
  Bfree(S);
2761
2452
  if (mhi)
2762
2453
  {
2763
2454
    if (mlo && mlo != mhi)
2764
 
      Bfree(mlo, &alloc);
2765
 
    Bfree(mhi, &alloc);
 
2455
      Bfree(mlo);
 
2456
    Bfree(mhi);
2766
2457
  }
2767
2458
ret1:
2768
 
  Bfree(b, &alloc);
 
2459
  Bfree(b);
2769
2460
  *s= 0;
2770
2461
  *decpt= k + 1;
2771
2462
  if (rve)