------- Comment #28 from vincent at vinc17 dot org 2006-10-31 22:15 ------- (In reply to comment #27) > It's likely that I'll end up doing it, so would you please tell me how?
According to the C rationale (I haven't checked), the sign of gamma(x) is -1 if [iff] x < 0 && remainder(floor(x), 2) != 0. But if x is a non-positive integer, the sign of gamma(x) isn't defined. Handle these cases first. The test x < 0 is easy to do. In MPFR, you can compute floor(x) (or trunc(x)) with the precision min(PREC(x),max(EXP(x),MPFR_PREC_MIN)), but then, there's no direct function to decide whether the result is even or odd (I thought we added this, but this isn't the case). The solution can be to divide x by 2 (this is exact, except in case of underflow) and call mpfr_frac directly. If the result is between -0.5 and 0, then gamma(x) is negative. If the result is between -1 and -0.5, then gamma(x) is positive. So, a 2-bit precision for mpfr_frac should be sufficient (as -0.5 is representable in this precision), but choose a directed rounding (not GMP_RNDN) for that. Then you can just do a comparison with -0.5; the case of equality with -0.5 depends on the chosen rounding (if you obtain -0.5, then it is an inexact result since x is not an integer). For instance, if you choose GMP_RNDZ, then a result > -0.5 means that gamma(x) is negative, and a result <= -0.5 means that gamma(x) is positive. -- http://gcc.gnu.org/bugzilla/show_bug.cgi?id=29335
