← Back to branch summary

~drizzle-trunk/drizzle/development

« back to all changes in this revision

Viewing changes to mystrings/dtoa.cc

Moved the last of the libdrizzleclient calls into Protocol.

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 */
 
38
#include <mystrings/m_string.h>  /* for memcpy and NOT_FIXED_DEC */
 
39
#include <stdlib.h>
39
40
 
40
41
/**
41
42
   Appears to suffice to not call malloc() in most cases.
87
88
{
88
89
  int decpt, sign, len, i;
89
90
  char *res, *src, *end, *dst= to;
90
 
  char buf[DTOA_BUFF_SIZE];
 
91
  double alignedbuf[(DTOA_BUFF_SIZE+1)/sizeof(double)];
 
92
  char *buf= (char*)alignedbuf;
91
93
  assert(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
92
 
  
 
94
 
93
95
  res= dtoa(x, 5, precision, &decpt, &sign, &end, buf, sizeof(buf));
94
96
 
95
97
  if (decpt == DTOA_OVERFLOW)
129
131
  {
130
132
    if (len <= decpt)
131
133
      *dst++= '.';
132
 
    
 
134
 
133
135
    for (i= precision - cmax(0, (len - decpt)); i > 0; i--)
134
136
      *dst++= '0';
135
137
  }
136
 
  
 
138
 
137
139
  *dst= '\0';
138
140
  if (error != NULL)
139
141
    *error= false;
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
*/
211
213
{
212
214
  int decpt, sign, len, exp_len;
213
215
  char *res, *src, *end, *dst= to, *dend= dst + width;
214
 
  char buf[DTOA_BUFF_SIZE];
 
216
  double alignedbuf[(DTOA_BUFF_SIZE+1)/sizeof(double)];
 
217
  char *buf= (char*)alignedbuf;
215
218
  bool have_space, force_e_format;
216
219
  assert(width > 0 && to != NULL);
217
 
  
 
220
 
218
221
  /* We want to remove '-' from equations early */
219
222
  if (x < 0.)
220
223
    width--;
243
246
     to count it here.
244
247
   */
245
248
  exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
246
 
  
 
249
 
247
250
  /*
248
251
     Do we have enough space for all digits in the 'f' format?
249
252
     Let 'len' be the number of significant digits returned by dtoa,
296
299
       ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
297
300
                                            (len > 1 || !force_e_format)))) &&
298
301
         !force_e_format)) &&
299
 
      
 
302
 
300
303
       /*
301
304
         Use the 'e' format in some cases even if we have enough space for the
302
305
         'f' one. See comment for MAX_DECPT_FOR_F_FORMAT.
318
321
          *error= true;
319
322
        width= decpt;
320
323
      }
321
 
      
 
324
 
322
325
      /*
323
326
        We want to truncate (len - width) least significant digits after the
324
327
        decimal point. For this we are calling dtoa with mode=5, passing the
336
339
      *dst++= '0';
337
340
      goto end;
338
341
    }
339
 
    
 
342
 
340
343
    /*
341
344
      At this point we are sure we have enough space to put all digits
342
345
      returned by dtoa
385
388
        *error= true;
386
389
      width= 0;
387
390
    }
388
 
      
 
391
 
389
392
    /* Do we have to truncate any digits? */
390
393
    if (width < len)
391
394
    {
450
453
                  rejected character.
451
454
   @param error   Upon return is set to EOVERFLOW in case of underflow or
452
455
                  overflow.
453
 
   
 
456
 
454
457
   @return        The resulting double value. In case of underflow, 0.0 is
455
458
                  returned. In case overflow, signed DBL_MAX is returned.
456
459
*/
457
460
 
458
461
double my_strtod(const char *str, char **end, int *error)
459
462
{
460
 
  char buf[DTOA_BUFF_SIZE];
 
463
  double alignedbuf[(DTOA_BUFF_SIZE+1)/sizeof(double)];
 
464
  char *buf= (char*)alignedbuf;
461
465
  double res;
462
466
  assert(str != NULL && end != NULL && *end != NULL && error != NULL);
463
467
 
528
532
  file.
529
533
*/
530
534
 
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
535
typedef int32_t Long;
541
536
typedef uint32_t ULong;
542
537
typedef int64_t LLong;
592
587
#endif
593
588
#endif /*Flt_Rounds*/
594
589
 
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
590
#define rounded_product(a,b) a*= b
604
591
#define rounded_quotient(a,b) a/= b
605
592
 
662
649
    int x, len;
663
650
 
664
651
    x= 1 << k;
665
 
    len= sizeof(Bigint) + x * sizeof(ULong);
 
652
    len= ALIGN_SIZE(sizeof(Bigint) + x * sizeof(ULong));
666
653
 
667
654
    if (alloc->free + len <= alloc->end)
668
655
    {
713
700
static char *dtoa_alloc(int i, Stack_alloc *alloc)
714
701
{
715
702
  char *rv;
716
 
  int aligned_size= (i + sizeof(ULong) - 1) / sizeof(ULong) * sizeof(ULong);
 
703
  int aligned_size= (i + sizeof(double) - 1) / sizeof(double) * sizeof(double);
717
704
  if (alloc->free + aligned_size <= alloc->end)
718
705
  {
719
706
    rv= alloc->free;
720
707
    alloc->free+= aligned_size;
721
708
  }
722
709
  else
723
 
    rv= malloc(i);
 
710
    rv= (char *)malloc(i);
724
711
  return rv;
725
712
}
726
713
 
786
773
  b= Balloc(k, alloc);
787
774
  b->p.x[0]= y9;
788
775
  b->wds= 1;
789
 
  
 
776
 
790
777
  i= 9;
791
778
  if (9 < nd0)
792
779
  {
1283
1270
  9007199254740992.*9007199254740992.e-256 /* = 2^106 * 1e-53 */
1284
1271
};
1285
1272
/*
1286
 
  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow 
 
1273
  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow
1287
1274
  flag unnecessarily.  It leads to a song and dance at the end of strtod.
1288
1275
*/
1289
1276
#define Scale_Bit 0x10
1291
1278
 
1292
1279
/*
1293
1280
  strtod for IEEE--arithmetic machines.
1294
 
 
 
1281
 
1295
1282
  This strtod returns a nearest machine number to the input decimal
1296
1283
  string (or sets errno to EOVERFLOW). Ties are broken by the IEEE round-even
1297
1284
  rule.
1298
 
 
 
1285
 
1299
1286
  Inspired loosely by William D. Clinger's paper "How to Read Floating
1300
1287
  Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
1301
 
 
 
1288
 
1302
1289
  Modifications:
1303
 
 
 
1290
 
1304
1291
   1. We only require IEEE (not IEEE double-extended).
1305
1292
   2. We get by with floating-point arithmetic in a case that
1306
1293
     Clinger missed -- when we're computing d * 10^n
1331
1318
#ifdef SET_INEXACT
1332
1319
  int inexact, oldinexact;
1333
1320
#endif
1334
 
#ifdef Honor_FLT_ROUNDS
1335
 
  int rounding;
1336
 
#endif
1337
1321
  Stack_alloc alloc;
1338
1322
 
1339
1323
  c= 0;
1366
1350
 break2:
1367
1351
  if (s >= end)
1368
1352
    goto ret0;
1369
 
  
 
1353
 
1370
1354
  if (*s == '0')
1371
1355
  {
1372
1356
    nz0= 1;
1491
1475
  }
1492
1476
  bd0= 0;
1493
1477
  if (nd <= DBL_DIG
1494
 
#ifndef Honor_FLT_ROUNDS
1495
 
    && Flt_Rounds == 1
1496
 
#endif
1497
1478
      )
1498
1479
  {
1499
1480
    if (!e)
1502
1483
    {
1503
1484
      if (e <= Ten_pmax)
1504
1485
      {
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
1486
        /* rv = */ rounded_product(dval(rv), tens[e]);
1514
1487
        goto ret;
1515
1488
      }
1520
1493
          A fancier test would sometimes let us do
1521
1494
          this for larger i values.
1522
1495
        */
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
1496
        e-= i;
1532
1497
        dval(rv)*= tens[i];
1533
1498
        /* rv = */ rounded_product(dval(rv), tens[e]);
1537
1502
#ifndef Inaccurate_Divide
1538
1503
    else if (e >= -Ten_pmax)
1539
1504
    {
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
1505
      /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1549
1506
      goto ret;
1550
1507
    }
1558
1515
    oldinexact= get_inexact();
1559
1516
#endif
1560
1517
  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
1518
 
1572
1519
  /* Get starting approximation = rv * 10**e1 */
1573
1520
 
1582
1529
 ovfl:
1583
1530
        *error= EOVERFLOW;
1584
1531
        /* 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
1532
        word0(rv)= Exp_mask;
1599
1533
        word1(rv)= 0;
1600
 
#endif /*Honor_FLT_ROUNDS*/
1601
1534
#ifdef SET_INEXACT
1602
1535
        /* set overflow bit */
1603
1536
        dval(rv0)= 1e300;
1693
1626
    else
1694
1627
      bd2-= bbe;
1695
1628
    bs2= bb2;
1696
 
#ifdef Honor_FLT_ROUNDS
1697
 
    if (rounding != 1)
1698
 
      bs2++;
1699
 
#endif
1700
1629
    j= bbe - scale;
1701
1630
    i= j + bbbits - 1;  /* logb(rv) */
1702
1631
    if (i < Emin)  /* denormal */
1734
1663
    dsign= delta->sign;
1735
1664
    delta->sign= 0;
1736
1665
    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
1666
 
1804
1667
    if (i < 0)
1805
1668
    {
1906
1769
    {
1907
1770
      aadj*= 0.5;
1908
1771
      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
1772
      if (Flt_Rounds == 0)
1921
1773
        aadj1+= 0.5;
1922
 
#endif /*Check_FLT_ROUNDS*/
1923
1774
    }
1924
1775
    y= word0(rv) & Exp_mask;
1925
1776
 
2146
1997
          4,5 ==> similar to 2 and 3, respectively, but (in
2147
1998
                round-nearest mode) with the tests of mode 0 to
2148
1999
                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
2000
          6-9 ==> Debugging modes similar to mode - 4:  don't try
2153
2001
                fast floating-point estimate (if applicable).
2154
2002
 
2167
2015
  Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
2168
2016
  double d2, ds, eps;
2169
2017
  char *s, *s0;
2170
 
#ifdef Honor_FLT_ROUNDS
2171
 
  int rounding;
2172
 
#endif
2173
2018
  Stack_alloc alloc;
2174
 
  
 
2019
 
2175
2020
  alloc.begin= alloc.free= buf;
2176
2021
  alloc.end= buf + buf_size;
2177
2022
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
2197
2042
      *rve= res + 1;
2198
2043
    return res;
2199
2044
  }
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
 
2045
 
2211
2046
 
2212
2047
  b= d2b(dval(d), &be, &bbits, &alloc);
2213
2048
  if ((i= (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
2221
2056
      log10(x)      =  log(x) / log(10)
2222
2057
                   ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2223
2058
      log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
2224
 
     
 
2059
 
2225
2060
      This suggests computing an approximation k to log10(d) by
2226
 
     
 
2061
 
2227
2062
      k= (i - Bias)*0.301029995663981
2228
2063
           + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2229
 
     
 
2064
 
2230
2065
      We want k to be too large rather than too small.
2231
2066
      The error in the first-order Taylor series approximation
2232
2067
      is in our favor, so we just round up the constant enough
2291
2126
  if (mode < 0 || mode > 9)
2292
2127
    mode= 0;
2293
2128
 
2294
 
#ifdef Check_FLT_ROUNDS
2295
 
  try_quick= Rounding == 1;
2296
 
#else
2297
2129
  try_quick= 1;
2298
 
#endif
2299
2130
 
2300
2131
  if (mode > 5)
2301
2132
  {
2330
2161
  }
2331
2162
  s= s0= dtoa_alloc(i, &alloc);
2332
2163
 
2333
 
#ifdef Honor_FLT_ROUNDS
2334
 
  if (mode > 1 && rounding != 1)
2335
 
    leftright= 0;
2336
 
#endif
2337
 
 
2338
2164
  if (ilim >= 0 && ilim <= Quick_max && try_quick)
2339
2165
  {
2340
2166
    /* Try to get by with floating-point arithmetic. */
2432
2258
            goto bump_up;
2433
2259
          else if (dval(d) < 0.5 - dval(eps))
2434
2260
          {
2435
 
            while (*--s == '0');
 
2261
            while (*--s == '0') {}
2436
2262
            s++;
2437
2263
            goto ret1;
2438
2264
          }
2464
2290
    {
2465
2291
      L= (Long)(dval(d) / ds);
2466
2292
      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
2293
      *s++= '0' + (int)L;
2476
2294
      if (!dval(d))
2477
2295
      {
2479
2297
      }
2480
2298
      if (i == ilim)
2481
2299
      {
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
2300
        dval(d)+= dval(d);
2492
2301
        if (dval(d) > ds || (dval(d) == ds && L & 1))
2493
2302
        {
2549
2358
 
2550
2359
  spec_case= 0;
2551
2360
  if ((mode < 2 || leftright)
2552
 
#ifdef Honor_FLT_ROUNDS
2553
 
      && rounding == 1
2554
 
#endif
2555
2361
     )
2556
2362
  {
2557
2363
    if (!word1(d) && !(word0(d) & Bndry_mask) &&
2568
2374
  /*
2569
2375
    Arrange for convenient computation of quotients:
2570
2376
    shift left if necessary so divisor has 4 leading 0 bits.
2571
 
    
 
2377
 
2572
2378
    Perhaps we should just compute leading 28 bits of S once
2573
2379
    a nd for all and pass them and a shift to quorem, so it
2574
2380
    can do shifts and ors to compute the numerator for q.
2645
2451
      j1= delta->sign ? 1 : cmp(b, delta);
2646
2452
      Bfree(delta, &alloc);
2647
2453
      if (j1 == 0 && mode != 1 && !(word1(d) & 1)
2648
 
#ifdef Honor_FLT_ROUNDS
2649
 
          && rounding >= 1
2650
 
#endif
2651
2454
         )
2652
2455
      {
2653
2456
        if (dig == '9')
2663
2466
        {
2664
2467
          goto accept_dig;
2665
2468
        }
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
2469
        if (j1 > 0)
2674
2470
        {
2675
2471
          b= lshift(b, 1, &alloc);
2684
2480
      }
2685
2481
      if (j1 > 0)
2686
2482
      {
2687
 
#ifdef Honor_FLT_ROUNDS
2688
 
        if (!rounding)
2689
 
          goto accept_dig;
2690
 
#endif
2691
2483
        if (dig == '9')
2692
2484
        { /* possible if i == 1 */
2693
2485
round_9_up:
2697
2489
        *s++= dig + 1;
2698
2490
        goto ret;
2699
2491
      }
2700
 
#ifdef Honor_FLT_ROUNDS
2701
 
keep_dig:
2702
 
#endif
2703
2492
      *s++= dig;
2704
2493
      if (i == ilim)
2705
2494
        break;
2728
2517
 
2729
2518
  /* Round off last digit */
2730
2519
 
2731
 
#ifdef Honor_FLT_ROUNDS
2732
 
  switch (rounding) {
2733
 
  case 0: goto trimzeros;
2734
 
  case 2: goto roundoff;
2735
 
  }
2736
 
#endif
2737
2520
  b= lshift(b, 1, &alloc);
2738
2521
  j= cmp(b, S);
2739
2522
  if (j > 0 || (j == 0 && dig & 1))
2750
2533
  }
2751
2534
  else
2752
2535
  {
2753
 
#ifdef Honor_FLT_ROUNDS
2754
 
trimzeros:
2755
 
#endif
2756
 
    while (*--s == '0');
 
2536
    while (*--s == '0') {}
2757
2537
    s++;
2758
2538
  }
2759
2539
ret: