[src/gamma.c] Corrections in a proof for tiny arguments. Added a TODO.

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

1 files changed, 18 insertions, 12 deletions

diff --git a/src/gamma.c b/src/gamma.c
index 7f42509e7..9a192eed0 100644
--- a/src/gamma.c
+++ b/src/gamma.c
@@ -158,14 +158,16 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
         }
     }
 
-  /* Check for tiny arguments, where gamma(x) ~ 1/x - euler + ....
+  /* Check for tiny arguments, where gamma(x) ~ 1/x - euler + ... can be
+     approximated by 1/x, with some error term ~= - euler.
      We know from "Bound on Runs of Zeros and Ones for Algebraic Functions",
      Proceedings of Arith15, T. Lang and J.-M. Muller, 2001, that the maximal
-     number of consecutive zeroes or ones after the round bit is n-1 for an
-     input of n bits. But we need a more precise lower bound. Assume x has
-     n bits, and 1/x is near a floating-point number y of n+1 bits. We can
-     write x = X*2^e, y = Y/2^f with X, Y integers of n and n+1 bits.
-     Thus X*Y^2^(e-f) is near from 1, i.e., X*Y is near from 2^(f-e).
+     number of consecutive zeroes or ones after the round bit for 1/x is n-1
+     for an input x of n bits [this is an actually much older result!].
+     But we need a more precise lower bound. Assume x has n bits, and 1/x
+     is near a floating-point number y of n+1 bits. We can write x = X*2^e,
+     y = Y/2^f with X, Y integers of n and n+1 bits respectively.
+     Thus X*Y^2^(e-f) is near 1, i.e., X*Y is near the integer 2^(f-e).
      Two cases can happen:
      (i) either X*Y is exactly 2^(f-e), but this can happen only if X and Y
          are themselves powers of two, i.e., x is a power of two;
@@ -173,13 +175,17 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
           |xy-1| >= 2^(e-f), or |y-1/x| >= 2^(e-f)/x = 2^(-f)/X >= 2^(-f-n).
           Since ufp(y) = 2^(n-f) [ufp = unit in first place], this means
           that the distance |y-1/x| >= 2^(-2n) ufp(y).
-          Now assuming |gamma(x)-1/x| <= 1, which is true for x <= 1,
-          if 2^(-2n) ufp(y) >= 2, the error is at most 2^(-2n-1) ufp(y),
-          and round(1/x) with precision >= 2n+2 gives the correct result.
-          If x < 2^E, then y > 2^(-E), thus ufp(y) > 2^(-E-1).
-          A sufficient condition is thus EXP(x) + 2 <= -2 MAX(PREC(x),PREC(Y)).
+          Now, assuming |gamma(x)-1/x| < 1, which is true for 0 < x <= 1,
+          if 2^(-2n) ufp(y) >= 1, then gamma(x) and 1/x round in the same
+          way, so that rounding 1/x gives the correct result and correct
+          (nonzero) ternary value.
+          If x < 2^E, then y > 2^(-E), thus ufp(y) >= 2^(-E).
+          A sufficient condition is thus EXP(x) <= -2 MAX(PREC(x),PREC(Y)).
   */
-  if (MPFR_GET_EXP (x) + 2
+  /* TODO: Handle the very small arguments that do not satisfy the condition,
+     by using the approximation 1/x - euler and a Ziv loop. Otherwise, after
+     some tests, even Gamma(1+x)/x would be faster than the generic code. */
+  if (MPFR_GET_EXP (x)
       <= -2 * (mpfr_exp_t) MAX(MPFR_PREC(x), MPFR_PREC(gamma)))
     {
       int sign = MPFR_SIGN (x); /* retrieve sign before possible override */