← Back to branch summary

~drizzle-trunk/drizzle/development

« back to all changes in this revision

Viewing changes to mystrings/dtoa.c

Merging Monty's work

Show diffs side-by-side

added added

removed removed

Lines of Context:
mystrings/dtoa.c
528
528
  file.
529
529
*/
530
530
 
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
531
typedef int32_t Long;
541
532
typedef uint32_t ULong;
542
533
typedef int64_t LLong;
592
583
#endif
593
584
#endif /*Flt_Rounds*/
594
585
 
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
586
#define rounded_product(a,b) a*= b
604
587
#define rounded_quotient(a,b) a/= b
605
588
 
1331
1314
#ifdef SET_INEXACT
1332
1315
  int inexact, oldinexact;
1333
1316
#endif
1334
 
#ifdef Honor_FLT_ROUNDS
1335
 
  int rounding;
1336
 
#endif
1337
1317
  Stack_alloc alloc;
1338
1318
 
1339
1319
  c= 0;
1491
1471
  }
1492
1472
  bd0= 0;
1493
1473
  if (nd <= DBL_DIG
1494
 
#ifndef Honor_FLT_ROUNDS
1495
 
    && Flt_Rounds == 1
1496
 
#endif
1497
1474
      )
1498
1475
  {
1499
1476
    if (!e)
1502
1479
    {
1503
1480
      if (e <= Ten_pmax)
1504
1481
      {
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
1482
        /* rv = */ rounded_product(dval(rv), tens[e]);
1514
1483
        goto ret;
1515
1484
      }
1520
1489
          A fancier test would sometimes let us do
1521
1490
          this for larger i values.
1522
1491
        */
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
1492
        e-= i;
1532
1493
        dval(rv)*= tens[i];
1533
1494
        /* rv = */ rounded_product(dval(rv), tens[e]);
1537
1498
#ifndef Inaccurate_Divide
1538
1499
    else if (e >= -Ten_pmax)
1539
1500
    {
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
1501
      /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1549
1502
      goto ret;
1550
1503
    }
1558
1511
    oldinexact= get_inexact();
1559
1512
#endif
1560
1513
  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
1514
 
1572
1515
  /* Get starting approximation = rv * 10**e1 */
1573
1516
 
1582
1525
 ovfl:
1583
1526
        *error= EOVERFLOW;
1584
1527
        /* 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
1528
        word0(rv)= Exp_mask;
1599
1529
        word1(rv)= 0;
1600
 
#endif /*Honor_FLT_ROUNDS*/
1601
1530
#ifdef SET_INEXACT
1602
1531
        /* set overflow bit */
1603
1532
        dval(rv0)= 1e300;
1693
1622
    else
1694
1623
      bd2-= bbe;
1695
1624
    bs2= bb2;
1696
 
#ifdef Honor_FLT_ROUNDS
1697
 
    if (rounding != 1)
1698
 
      bs2++;
1699
 
#endif
1700
1625
    j= bbe - scale;
1701
1626
    i= j + bbbits - 1;  /* logb(rv) */
1702
1627
    if (i < Emin)  /* denormal */
1734
1659
    dsign= delta->sign;
1735
1660
    delta->sign= 0;
1736
1661
    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
1662
 
1804
1663
    if (i < 0)
1805
1664
    {
1906
1765
    {
1907
1766
      aadj*= 0.5;
1908
1767
      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
1768
      if (Flt_Rounds == 0)
1921
1769
        aadj1+= 0.5;
1922
 
#endif /*Check_FLT_ROUNDS*/
1923
1770
    }
1924
1771
    y= word0(rv) & Exp_mask;
1925
1772
 
2146
1993
          4,5 ==> similar to 2 and 3, respectively, but (in
2147
1994
                round-nearest mode) with the tests of mode 0 to
2148
1995
                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
1996
          6-9 ==> Debugging modes similar to mode - 4:  don't try
2153
1997
                fast floating-point estimate (if applicable).
2154
1998
 
2167
2011
  Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
2168
2012
  double d2, ds, eps;
2169
2013
  char *s, *s0;
2170
 
#ifdef Honor_FLT_ROUNDS
2171
 
  int rounding;
2172
 
#endif
2173
2014
  Stack_alloc alloc;
2174
2015
  
2175
2016
  alloc.begin= alloc.free= buf;
2198
2039
    return res;
2199
2040
  }
2200
2041
  
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
2042
 
2212
2043
  b= d2b(dval(d), &be, &bbits, &alloc);
2213
2044
  if ((i= (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))))
2291
2122
  if (mode < 0 || mode > 9)
2292
2123
    mode= 0;
2293
2124
 
2294
 
#ifdef Check_FLT_ROUNDS
2295
 
  try_quick= Rounding == 1;
2296
 
#else
2297
2125
  try_quick= 1;
2298
 
#endif
2299
2126
 
2300
2127
  if (mode > 5)
2301
2128
  {
2330
2157
  }
2331
2158
  s= s0= dtoa_alloc(i, &alloc);
2332
2159
 
2333
 
#ifdef Honor_FLT_ROUNDS
2334
 
  if (mode > 1 && rounding != 1)
2335
 
    leftright= 0;
2336
 
#endif
2337
 
 
2338
2160
  if (ilim >= 0 && ilim <= Quick_max && try_quick)
2339
2161
  {
2340
2162
    /* Try to get by with floating-point arithmetic. */
2464
2286
    {
2465
2287
      L= (Long)(dval(d) / ds);
2466
2288
      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
2289
      *s++= '0' + (int)L;
2476
2290
      if (!dval(d))
2477
2291
      {
2479
2293
      }
2480
2294
      if (i == ilim)
2481
2295
      {
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
2296
        dval(d)+= dval(d);
2492
2297
        if (dval(d) > ds || (dval(d) == ds && L & 1))
2493
2298
        {
2549
2354
 
2550
2355
  spec_case= 0;
2551
2356
  if ((mode < 2 || leftright)
2552
 
#ifdef Honor_FLT_ROUNDS
2553
 
      && rounding == 1
2554
 
#endif
2555
2357
     )
2556
2358
  {
2557
2359
    if (!word1(d) && !(word0(d) & Bndry_mask) &&
2645
2447
      j1= delta->sign ? 1 : cmp(b, delta);
2646
2448
      Bfree(delta, &alloc);
2647
2449
      if (j1 == 0 && mode != 1 && !(word1(d) & 1)
2648
 
#ifdef Honor_FLT_ROUNDS
2649
 
          && rounding >= 1
2650
 
#endif
2651
2450
         )
2652
2451
      {
2653
2452
        if (dig == '9')
2663
2462
        {
2664
2463
          goto accept_dig;
2665
2464
        }
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
2465
        if (j1 > 0)
2674
2466
        {
2675
2467
          b= lshift(b, 1, &alloc);
2684
2476
      }
2685
2477
      if (j1 > 0)
2686
2478
      {
2687
 
#ifdef Honor_FLT_ROUNDS
2688
 
        if (!rounding)
2689
 
          goto accept_dig;
2690
 
#endif
2691
2479
        if (dig == '9')
2692
2480
        { /* possible if i == 1 */
2693
2481
round_9_up:
2697
2485
        *s++= dig + 1;
2698
2486
        goto ret;
2699
2487
      }
2700
 
#ifdef Honor_FLT_ROUNDS
2701
 
keep_dig:
2702
 
#endif
2703
2488
      *s++= dig;
2704
2489
      if (i == ilim)
2705
2490
        break;
2728
2513
 
2729
2514
  /* Round off last digit */
2730
2515
 
2731
 
#ifdef Honor_FLT_ROUNDS
2732
 
  switch (rounding) {
2733
 
  case 0: goto trimzeros;
2734
 
  case 2: goto roundoff;
2735
 
  }
2736
 
#endif
2737
2516
  b= lshift(b, 1, &alloc);
2738
2517
  j= cmp(b, S);
2739
2518
  if (j > 0 || (j == 0 && dig & 1))
2750
2529
  }
2751
2530
  else
2752
2531
  {
2753
 
#ifdef Honor_FLT_ROUNDS
2754
 
trimzeros:
2755
 
#endif
2756
2532
    while (*--s == '0');
2757
2533
    s++;
2758
2534
  }