← Back to branch summary

~drizzle-trunk/drizzle/development

« back to all changes in this revision

Viewing changes to mystrings/dtoa.cc

  • Committer: Brian Aker
  • Date: 2009-02-05 10:38:55 UTC
  • Revision ID: brian@tangent.org-20090205103855-wajzccrbu7zbvmh4
Reworked some classes out of session.h
Also updated ignore file.

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.
89
90
  char *res, *src, *end, *dst= to;
90
91
  char buf[DTOA_BUFF_SIZE];
91
92
  assert(precision >= 0 && precision < NOT_FIXED_DEC && to != NULL);
92
 
  
 
93
 
93
94
  res= dtoa(x, 5, precision, &decpt, &sign, &end, buf, sizeof(buf));
94
95
 
95
96
  if (decpt == DTOA_OVERFLOW)
129
130
  {
130
131
    if (len <= decpt)
131
132
      *dst++= '.';
132
 
    
 
133
 
133
134
    for (i= precision - cmax(0, (len - decpt)); i > 0; i--)
134
135
      *dst++= '0';
135
136
  }
136
 
  
 
137
 
137
138
  *dst= '\0';
138
139
  if (error != NULL)
139
140
    *error= false;
196
197
     my_gcvt(55, ..., 1, ...);
197
198
 
198
199
   We do our best to minimize such cases by:
199
 
   
 
200
 
200
201
   - passing to dtoa() the field width as the number of significant digits
201
 
   
 
202
 
202
203
   - removing the sign of the number early (and decreasing the width before
203
204
     passing it to dtoa())
204
 
   
 
205
 
205
206
   - choosing the proper format to preserve the most number of significant
206
207
     digits.
207
208
*/
214
215
  char buf[DTOA_BUFF_SIZE];
215
216
  bool have_space, force_e_format;
216
217
  assert(width > 0 && to != NULL);
217
 
  
 
218
 
218
219
  /* We want to remove '-' from equations early */
219
220
  if (x < 0.)
220
221
    width--;
243
244
     to count it here.
244
245
   */
245
246
  exp_len= 1 + (decpt >= 101 || decpt <= -99) + (decpt >= 11 || decpt <= -9);
246
 
  
 
247
 
247
248
  /*
248
249
     Do we have enough space for all digits in the 'f' format?
249
250
     Let 'len' be the number of significant digits returned by dtoa,
296
297
       ((decpt <= width && (decpt >= -1 || (decpt == -2 &&
297
298
                                            (len > 1 || !force_e_format)))) &&
298
299
         !force_e_format)) &&
299
 
      
 
300
 
300
301
       /*
301
302
         Use the 'e' format in some cases even if we have enough space for the
302
303
         'f' one. See comment for MAX_DECPT_FOR_F_FORMAT.
318
319
          *error= true;
319
320
        width= decpt;
320
321
      }
321
 
      
 
322
 
322
323
      /*
323
324
        We want to truncate (len - width) least significant digits after the
324
325
        decimal point. For this we are calling dtoa with mode=5, passing the
336
337
      *dst++= '0';
337
338
      goto end;
338
339
    }
339
 
    
 
340
 
340
341
    /*
341
342
      At this point we are sure we have enough space to put all digits
342
343
      returned by dtoa
385
386
        *error= true;
386
387
      width= 0;
387
388
    }
388
 
      
 
389
 
389
390
    /* Do we have to truncate any digits? */
390
391
    if (width < len)
391
392
    {
450
451
                  rejected character.
451
452
   @param error   Upon return is set to EOVERFLOW in case of underflow or
452
453
                  overflow.
453
 
   
 
454
 
454
455
   @return        The resulting double value. In case of underflow, 0.0 is
455
456
                  returned. In case overflow, signed DBL_MAX is returned.
456
457
*/
528
529
  file.
529
530
*/
530
531
 
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
532
typedef int32_t Long;
541
533
typedef uint32_t ULong;
542
534
typedef int64_t LLong;
592
584
#endif
593
585
#endif /*Flt_Rounds*/
594
586
 
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
587
#define rounded_product(a,b) a*= b
604
588
#define rounded_quotient(a,b) a/= b
605
589
 
720
704
    alloc->free+= aligned_size;
721
705
  }
722
706
  else
723
 
    rv= malloc(i);
 
707
    rv= (char *)malloc(i);
724
708
  return rv;
725
709
}
726
710
 
786
770
  b= Balloc(k, alloc);
787
771
  b->p.x[0]= y9;
788
772
  b->wds= 1;
789
 
  
 
773
 
790
774
  i= 9;
791
775
  if (9 < nd0)
792
776
  {
1283
1267
  9007199254740992.*9007199254740992.e-256 /* = 2^106 * 1e-53 */
1284
1268
};
1285
1269
/*
1286
 
  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow 
 
1270
  The factor of 2^53 in tinytens[4] helps us avoid setting the underflow
1287
1271
  flag unnecessarily.  It leads to a song and dance at the end of strtod.
1288
1272
*/
1289
1273
#define Scale_Bit 0x10
1291
1275
 
1292
1276
/*
1293
1277
  strtod for IEEE--arithmetic machines.
1294
 
 
 
1278
 
1295
1279
  This strtod returns a nearest machine number to the input decimal
1296
1280
  string (or sets errno to EOVERFLOW). Ties are broken by the IEEE round-even
1297
1281
  rule.
1298
 
 
 
1282
 
1299
1283
  Inspired loosely by William D. Clinger's paper "How to Read Floating
1300
1284
  Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
1301
 
 
 
1285
 
1302
1286
  Modifications:
1303
 
 
 
1287
 
1304
1288
   1. We only require IEEE (not IEEE double-extended).
1305
1289
   2. We get by with floating-point arithmetic in a case that
1306
1290
     Clinger missed -- when we're computing d * 10^n
1331
1315
#ifdef SET_INEXACT
1332
1316
  int inexact, oldinexact;
1333
1317
#endif
1334
 
#ifdef Honor_FLT_ROUNDS
1335
 
  int rounding;
1336
 
#endif
1337
1318
  Stack_alloc alloc;
1338
1319
 
1339
1320
  c= 0;
1366
1347
 break2:
1367
1348
  if (s >= end)
1368
1349
    goto ret0;
1369
 
  
 
1350
 
1370
1351
  if (*s == '0')
1371
1352
  {
1372
1353
    nz0= 1;
1491
1472
  }
1492
1473
  bd0= 0;
1493
1474
  if (nd <= DBL_DIG
1494
 
#ifndef Honor_FLT_ROUNDS
1495
 
    && Flt_Rounds == 1
1496
 
#endif
1497
1475
      )
1498
1476
  {
1499
1477
    if (!e)
1502
1480
    {
1503
1481
      if (e <= Ten_pmax)
1504
1482
      {
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
1483
        /* rv = */ rounded_product(dval(rv), tens[e]);
1514
1484
        goto ret;
1515
1485
      }
1520
1490
          A fancier test would sometimes let us do
1521
1491
          this for larger i values.
1522
1492
        */
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
1493
        e-= i;
1532
1494
        dval(rv)*= tens[i];
1533
1495
        /* rv = */ rounded_product(dval(rv), tens[e]);
1537
1499
#ifndef Inaccurate_Divide
1538
1500
    else if (e >= -Ten_pmax)
1539
1501
    {
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
1502
      /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1549
1503
      goto ret;
1550
1504
    }
1558
1512
    oldinexact= get_inexact();
1559
1513
#endif
1560
1514
  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
1515
 
1572
1516
  /* Get starting approximation = rv * 10**e1 */
1573
1517
 
1582
1526
 ovfl:
1583
1527
        *error= EOVERFLOW;
1584
1528
        /* 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
1529
        word0(rv)= Exp_mask;
1599
1530
        word1(rv)= 0;
1600
 
#endif /*Honor_FLT_ROUNDS*/
1601
1531
#ifdef SET_INEXACT
1602
1532
        /* set overflow bit */
1603
1533
        dval(rv0)= 1e300;
1693
1623
    else
1694
1624
      bd2-= bbe;
1695
1625
    bs2= bb2;
1696
 
#ifdef Honor_FLT_ROUNDS
1697
 
    if (rounding != 1)
1698
 
      bs2++;
1699
 
#endif
1700
1626
    j= bbe - scale;
1701
1627
    i= j + bbbits - 1;  /* logb(rv) */
1702
1628
    if (i < Emin)  /* denormal */
1734
1660
    dsign= delta->sign;
1735
1661
    delta->sign= 0;
1736
1662
    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
1663
 
1804
1664
    if (i < 0)
1805
1665
    {
1906
1766
    {
1907
1767
      aadj*= 0.5;
1908
1768
      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
1769
      if (Flt_Rounds == 0)
1921
1770
        aadj1+= 0.5;
1922
 
#endif /*Check_FLT_ROUNDS*/
1923
1771
    }
1924
1772
    y= word0(rv) & Exp_mask;
1925
1773
 
2146
1994
          4,5 ==> similar to 2 and 3, respectively, but (in
2147
1995
                round-nearest mode) with the tests of mode 0 to
2148
1996
                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
1997
          6-9 ==> Debugging modes similar to mode - 4:  don't try
2153
1998
                fast floating-point estimate (if applicable).
2154
1999
 
2167
2012
  Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
2168
2013
  double d2, ds, eps;
2169
2014
  char *s, *s0;
2170
 
#ifdef Honor_FLT_ROUNDS
2171
 
  int rounding;
2172
 
#endif
2173
2015
  Stack_alloc alloc;
2174
 
  
 
2016
 
2175
2017
  alloc.begin= alloc.free= buf;
2176
2018
  alloc.end= buf + buf_size;
2177
2019
  memset(alloc.freelist, 0, sizeof(alloc.freelist));
2197
2039
      *rve= res + 1;
2198
2040
    return res;
2199
2041
  }
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
 
2042
 
2211
2043
 
2212
2044
  b= d2b(dval(d), &be, &bbits, &alloc);
2213
2045
  if ((i= (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
2221
2053
      log10(x)      =  log(x) / log(10)
2222
2054
                   ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2223
2055
      log10(d)= (i-Bias)*log(2)/log(10) + log10(d2)
2224
 
     
 
2056
 
2225
2057
      This suggests computing an approximation k to log10(d) by
2226
 
     
 
2058
 
2227
2059
      k= (i - Bias)*0.301029995663981
2228
2060
           + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2229
 
     
 
2061
 
2230
2062
      We want k to be too large rather than too small.
2231
2063
      The error in the first-order Taylor series approximation
2232
2064
      is in our favor, so we just round up the constant enough
2291
2123
  if (mode < 0 || mode > 9)
2292
2124
    mode= 0;
2293
2125
 
2294
 
#ifdef Check_FLT_ROUNDS
2295
 
  try_quick= Rounding == 1;
2296
 
#else
2297
2126
  try_quick= 1;
2298
 
#endif
2299
2127
 
2300
2128
  if (mode > 5)
2301
2129
  {
2330
2158
  }
2331
2159
  s= s0= dtoa_alloc(i, &alloc);
2332
2160
 
2333
 
#ifdef Honor_FLT_ROUNDS
2334
 
  if (mode > 1 && rounding != 1)
2335
 
    leftright= 0;
2336
 
#endif
2337
 
 
2338
2161
  if (ilim >= 0 && ilim <= Quick_max && try_quick)
2339
2162
  {
2340
2163
    /* Try to get by with floating-point arithmetic. */
2432
2255
            goto bump_up;
2433
2256
          else if (dval(d) < 0.5 - dval(eps))
2434
2257
          {
2435
 
            while (*--s == '0');
 
2258
            while (*--s == '0') {}
2436
2259
            s++;
2437
2260
            goto ret1;
2438
2261
          }
2464
2287
    {
2465
2288
      L= (Long)(dval(d) / ds);
2466
2289
      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
2290
      *s++= '0' + (int)L;
2476
2291
      if (!dval(d))
2477
2292
      {
2479
2294
      }
2480
2295
      if (i == ilim)
2481
2296
      {
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
2297
        dval(d)+= dval(d);
2492
2298
        if (dval(d) > ds || (dval(d) == ds && L & 1))
2493
2299
        {
2549
2355
 
2550
2356
  spec_case= 0;
2551
2357
  if ((mode < 2 || leftright)
2552
 
#ifdef Honor_FLT_ROUNDS
2553
 
      && rounding == 1
2554
 
#endif
2555
2358
     )
2556
2359
  {
2557
2360
    if (!word1(d) && !(word0(d) & Bndry_mask) &&
2568
2371
  /*
2569
2372
    Arrange for convenient computation of quotients:
2570
2373
    shift left if necessary so divisor has 4 leading 0 bits.
2571
 
    
 
2374
 
2572
2375
    Perhaps we should just compute leading 28 bits of S once
2573
2376
    a nd for all and pass them and a shift to quorem, so it
2574
2377
    can do shifts and ors to compute the numerator for q.
2645
2448
      j1= delta->sign ? 1 : cmp(b, delta);
2646
2449
      Bfree(delta, &alloc);
2647
2450
      if (j1 == 0 && mode != 1 && !(word1(d) & 1)
2648
 
#ifdef Honor_FLT_ROUNDS
2649
 
          && rounding >= 1
2650
 
#endif
2651
2451
         )
2652
2452
      {
2653
2453
        if (dig == '9')
2663
2463
        {
2664
2464
          goto accept_dig;
2665
2465
        }
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
2466
        if (j1 > 0)
2674
2467
        {
2675
2468
          b= lshift(b, 1, &alloc);
2684
2477
      }
2685
2478
      if (j1 > 0)
2686
2479
      {
2687
 
#ifdef Honor_FLT_ROUNDS
2688
 
        if (!rounding)
2689
 
          goto accept_dig;
2690
 
#endif
2691
2480
        if (dig == '9')
2692
2481
        { /* possible if i == 1 */
2693
2482
round_9_up:
2697
2486
        *s++= dig + 1;
2698
2487
        goto ret;
2699
2488
      }
2700
 
#ifdef Honor_FLT_ROUNDS
2701
 
keep_dig:
2702
 
#endif
2703
2489
      *s++= dig;
2704
2490
      if (i == ilim)
2705
2491
        break;
2728
2514
 
2729
2515
  /* Round off last digit */
2730
2516
 
2731
 
#ifdef Honor_FLT_ROUNDS
2732
 
  switch (rounding) {
2733
 
  case 0: goto trimzeros;
2734
 
  case 2: goto roundoff;
2735
 
  }
2736
 
#endif
2737
2517
  b= lshift(b, 1, &alloc);
2738
2518
  j= cmp(b, S);
2739
2519
  if (j > 0 || (j == 0 && dig & 1))
2750
2530
  }
2751
2531
  else
2752
2532
  {
2753
 
#ifdef Honor_FLT_ROUNDS
2754
 
trimzeros:
2755
 
#endif
2756
 
    while (*--s == '0');
 
2533
    while (*--s == '0') {}
2757
2534
    s++;
2758
2535
  }
2759
2536
ret: