diff options
| -rw-r--r-- | gamma.c | 12 | ||||
| -rw-r--r-- | tests/tgamma.c | 14 |
2 files changed, 24 insertions, 2 deletions
@@ -173,15 +173,23 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) if (compared > 0) { int overflow; + mpfr_t yp; + mpfr_clear_overflow (); mpfr_init2 (xp, 53); + mpfr_init2 (yp, 53); mpfr_set_d (xp, EXPM1, GMP_RNDZ); /* 1/e rounded down */ mpfr_mul (xp, x, xp, GMP_RNDZ); - mpfr_pow (xp, xp, x, GMP_RNDZ); - mpfr_clear_overflow (); + mpfr_sub_ui (yp, x, 2, GMP_RNDZ); + mpfr_pow (xp, xp, yp, GMP_RNDZ); /* (x/e)^(x-2) */ + mpfr_set_d (yp, EXPM1, GMP_RNDZ); + mpfr_mul (xp, xp, yp, GMP_RNDZ); /* x^(x-2) / e^(x-1) */ + mpfr_mul (xp, xp, yp, GMP_RNDZ); /* x^(x-2) / e^x */ + mpfr_mul (xp, xp, x, GMP_RNDZ); /* x^(x-1) / e^x */ mpfr_mul_2exp (xp, xp, 1, GMP_RNDZ); overflow = mpfr_overflow_p (); mpfr_clear (xp); + mpfr_clear (yp); return (overflow) ? mpfr_overflow (gamma, rnd_mode, 1) : mpfr_gamma_aux (gamma, x, rnd_mode); } diff --git a/tests/tgamma.c b/tests/tgamma.c index df0471eb8..9e8ca7746 100644 --- a/tests/tgamma.c +++ b/tests/tgamma.c @@ -374,6 +374,20 @@ special_overflow (void) exit (1); } + mpfr_set_emax (1024); + mpfr_set_prec (x, 53); + mpfr_set_prec (y, 53); + mpfr_set_str_binary (x, "101010110100110011111010000110001000111100000110101E-43"); + mpfr_gamma (x, x, GMP_RNDU); + mpfr_set_str_binary (y, "110000011110001000111110110101011110000100001111111E971"); + if (mpfr_cmp (x, y) != 0) + { + printf ("Error for gamma(4)\n"); + printf ("expected "); mpfr_dump (y); + printf ("got "); mpfr_dump (x); + exit (1); + } + mpfr_clear (y); mpfr_clear (x); set_emin (MPFR_EMIN_MIN); |