From: Javier Goizueta <jgoizueta@...>
Date: 2004-10-22T03:27:46+09:00
Subject: Re: A problem with bigdecimal.c, and proposed solution (long)

> Last minute note:
> I've realized that in the unpatched version of bigdecimal.c:
> 
>   (BigDecimal('1')/3).to_f == 1.0/3
> 
> That means that adding more than necessary "figures" degrades
> the result. This could be solved by being more careful when
> computing fig:
> 
>   fig =(DOUBLE_DIGITS+K + 4-n + BASE_FIG - 1) / BASE_FIG;

I've done a small change in the while loop of  that improves the
accuracy: accumulate div by multiplying, not dividing:

    while(ind_m < mm) {
        div *=(flt_t)((S_INT)BASE);
        *d += ((flt_t) ((S_INT)m->frac[ind_m++])) / div;
    }

Using this change and the previous one properly applied, the VpVtoD
function is working quite well now; Here's a patch from Rev.1.43 with
all the changes:


--------------------------------------------------------------BEGIN OF
PATCH
--- Rev.1.43/bigdecimal.c	2004-10-19 12:24:40.000000000 +0200
+++ bigdecimal.c	2004-10-21 20:21:58.000000000 +0200
@@ -1370,7 +1370,7 @@
                 /* The value of BASE**2 + BASE must be represented */
                 /* within one U_LONG. */
 static U_LONG HALF_BASE = 5000L;/* =BASE/2  */
-static S_LONG DBLE_FIG = 8;    /* figure of double */
+static S_LONG DOUBLE_DIGITS = DBL_DIG;    /* figure of double */
 static U_LONG BASE1 = 1000L;    /* =BASE/10  */

 static Real *VpConstOne;    /* constant 1.0 */
@@ -1515,7 +1515,7 @@
 VP_EXPORT U_LONG
 VpDblFig(void)
 {
-    return DBLE_FIG;
+    return DOUBLE_DIGITS+2;
 }

 VP_EXPORT U_LONG
@@ -1760,7 +1760,7 @@
  *   by one U_LONG word(LONG) in the computer used.
  *
  * [Returns]
- * DBLE_FIG   ... OK
+ * DOUBLE_DIGITS   ... OK
  */
 VP_EXPORT U_LONG
 VpInit(U_LONG BaseVal)
@@ -1799,15 +1799,6 @@
     gnAlloc = 0;
 #endif /* _DEBUG */

-    /* Determine # of digits available in one 'double'. */
-
-    v = 1.0;
-    DBLE_FIG = 0;
-    while(v + 1.0 > 1.0) {
-        ++DBLE_FIG;
-        v /= 10;
-    }
-
 #ifdef _DEBUG
     if(gfDebug) {
         printf("VpInit: BaseVal   = %lu\n", BaseVal);
@@ -1819,7 +1810,7 @@
     }
 #endif /* _DEBUG */

-    return DBLE_FIG;
+    return DOUBLE_DIGITS;
 }

 VP_EXPORT Real *
@@ -3401,7 +3392,7 @@
  * [Output]
  *   *d  ... fraction part of m(d = 0.xxxxxxx). where # of 'x's is fig.
  *   *e  ... U_LONG,exponent of m.
- * DBLE_FIG ... Number of digits in a double variable.
+ * DOUBLE_DIGITS ... Number of digits in a double variable.
  *
  *  m -> d*10**e, 0<d<BASE
  * [Returns]
@@ -3413,7 +3404,8 @@
 VP_EXPORT int
 VpVtoD(double *d, S_LONG *e, Real *m)
 {
-    U_LONG ind_m, mm, fig;
+    int n;
+    U_LONG ind_m, mm, fig, d0;
     double div;
     int    f = 1;

@@ -3448,14 +3440,20 @@
         goto Exit;
     }
     /* Normal number */
-    fig =(DBLE_FIG + BASE_FIG - 1) / BASE_FIG;
+    d0 = m->frac[0];
+    n = 1;
+    while (d0>=10) {
+      n++;
+      d0 /= 10;
+    }
+    fig = (DOUBLE_DIGITS+2 + 4-n + BASE_FIG - 1) / BASE_FIG;
     ind_m = 0;
     mm = Min(fig,(m->Prec));
     *d = 0.0;
     div = 1.;
     while(ind_m < mm) {
-        div /=(double)((S_INT)BASE);
-        *d = *d +((double) ((S_INT)m->frac[ind_m++])) * div;
+        div *=(double)((S_INT)BASE);
+        *d += ((double) ((S_INT)m->frac[ind_m++])) / div;
     }
     *e = m->exponent * ((S_INT)BASE_FIG);
     *d *= VpGetSign(m);
@@ -3465,7 +3463,7 @@
     if(gfDebug) {
         VPrint(stdout, " VpVtoD: m=%\n", m);
         printf("   d=%e * 10 **%ld\n", *d, *e);
-        printf("   DBLE_FIG = %ld\n", DBLE_FIG);
+        printf("   DOUBLE_DIGITS = %ld\n", DOUBLE_DIGITS);
     }
 #endif /*_DEBUG */
     return f;
@@ -3663,7 +3661,7 @@
     }
     VpDtoV(y, sqrt(val));    /* y <- sqrt(val) */
     y->exponent += n;
-    n = (DBLE_FIG + BASE_FIG - 1) / BASE_FIG;
+    n = (DOUBLE_DIGITS + BASE_FIG) / BASE_FIG;
     y->MaxPrec = (U_LONG)Min(n , y_prec);
     f->MaxPrec = y->MaxPrec + 1;
     n = y_prec*((S_LONG)BASE_FIG);
--------------------------------------------------------------END OF
PATCH


Note: in my previous post,
  
  format "%.16f",BigDecimal('1.234567890123456').to_f

meant to be

  format "%.15f",BigDecimal('1.234567890123456').to_f

--Javier