[src/strtofr.c] Started to review the new error analysis and code

(r12705,12706). Minor improvements.

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@12716 280ebfd0-de03-0410-8827-d642c229c3f4

1 files changed, 25 insertions, 17 deletions

diff --git a/src/strtofr.c b/src/strtofr.c
index 11232fc4f..b120b187e 100644
--- a/src/strtofr.c
+++ b/src/strtofr.c
@@ -678,29 +678,32 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
           /* (z, exp_z) = base^(pstr_size - exp_base) */
           z = result + 2*ysize + 1;
           err = mpfr_mpn_exp (z, &exp_z, pstr->base, exp_z, ysize);
-          /* Now {z, ysize}*2^(exp_z_out - ysize_bits) is an approximation of
-             base^exp_z_in, rounded towards zero, with:
-             * if err=-1, the result is exact
-             * if err=-2, an overflow occurred in the computation of exp_z
+
+          /* Now {z, ysize} * 2^(exp_z_out - ysize_bits) is an approximation
+             to base^exp_z_in (denoted b^e below), rounded toward zero, with:
+             * if err = -1, the result is exact;
+             * if err = -2, an overflow occurred in the computation of exp_z;
              * otherwise the error is bounded by 2^err ulps.
-             Thus the exact value of b^e is between z and z + 2^err (up to a
-             power of 2 multiplier for the exponent). Then the error will be:
+             Thus the exact value of b^e is between z and z + 2^err, where
+             z is {z, ysize} properly scaled by a power of 2. Then the error
+             will be:
                y/b^e - trunc(y/z) = eps1 + eps2
              with
                eps1 = y/b^e - y/z <= 0
                eps2 = y/z - trunc(y/z) >= 0
              thus the errors will (partly) compensate, giving a bound
              max(|eps1|,|eps2|).
-             In addition there is a 3rd error eps3 since y might be the
-             conversion of only part of the characters strings, and/or y
+             In addition, there is a 3rd error eps3 since y might be the
+             conversion of only a part of the character string, and/or y
              might be truncated by the mpn_rshift call above:
                eps3 = y0/b^e - y/b^e >= 0.
           */
-          exact = exact && (err == -1);
           if (err == -2)
             goto underflow; /* FIXME: Sure? */
-          if (err == -1)
-            err = 0;
+          else if (err == -1)
+            err = 0; /* see the note below */
+          else
+            exact = 0;
 
           /* Compute the integer division y/z rounded toward zero.
              The quotient will be put at result + ysize (size: ysize + 1),
@@ -716,9 +719,9 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
           /* The truncation error of the mpn_tdiv_qr call (eps2 above) is at
              most 1 ulp. Idem for the error eps3 which has the same sign,
              thus eps2 + eps3 <= 2 ulps.
-             For the error eps1 coming from the approximation of
-             b^e, we have (still up to a power-of-2 normalization):
-             y/z - y/b^e = y*(b^e-z)/(z*b^e) <= y*2^err/(z*b^e).
+             For the error eps1 coming from the approximation to b^e,
+             we have (still up to a power-of-2 normalization):
+             y/z - y/b^e = y * (b^e-z) / (z * b^e) <= y * 2^err / (z * b^e).
              We have to convert that error in terms of ulp(trunc(y/z)).
              We first have ulp(trunc(y/z)) = ulp(y/z).
              Since both y and z are normalized, the quotient
@@ -726,7 +729,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
              bit:
              * if the quotient has exactly ysize limbs, then 1/2 <= |y/z| < 1
                (up to a power of 2) and since 1/2 <= b^e < 1, the error is at
-               most 2^(err+1) ulps.
+               most 2^(err+1) ulps;
              * if the quotient has one extra bit, then 1 <= |y/z| < 2
                (up to a power of 2) and since 1/2 <= b^e < 1, the error is at
                most 2^(err+2) ulps; but since we will shift the result right
@@ -735,8 +738,13 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
 
              Thus the error is:
              * at most 2^(err+1) ulps for eps1
-             * at most 2 ulps for eps2 + eps3, which is opposite sign
-             and we can bound the error by 2^(err+1) ulps in all cases. */
+             * at most 2 ulps for eps2 + eps3, which is of opposite sign
+             and we can bound the error by 2^(err+1) ulps in all cases.
+
+             Note: If eps1 was 0, the error would be bounded by 2 ulps,
+             thus replacing err = -1 by err = 0 above was the right thing
+             to do, since 2^(0+1) = 2.
+          */
 
           /* exp -= exp_z + ysize_bits with overflow checking
              and check that we can add/subtract 2 to exp without overflow */