diff options
author | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2007-07-03 23:35:05 +0000 |
---|---|---|
committer | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2007-07-03 23:35:05 +0000 |
commit | c85e17ec69c20771935ce94a01e15871dfe8e3c4 (patch) | |
tree | ab18eea28b6593571a1b511ec5be34474c37e132 /lngamma.c | |
parent | 6e3671cee18f0eb2a57dc804d6c4a1d447e9f808 (diff) | |
download | mpfr-c85e17ec69c20771935ce94a01e15871dfe8e3c4.tar.gz |
lngamma.c: reformat and replaced mpfr_cmp(...) == 0 by mpfr_equal_p.
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@4628 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'lngamma.c')
-rw-r--r-- | lngamma.c | 62 |
1 files changed, 31 insertions, 31 deletions
@@ -610,37 +610,37 @@ mpfr_lgamma (mpfr_ptr y, int *signp, mpfr_srcptr x, mp_rnd_t rnd) { mpfr_t l, h; int ok, inex2; - mp_prec_t w = MPFR_PREC(y) + 14; - - while (1) - { - mpfr_init2 (l, w); - mpfr_init2 (h, w); - /* we want a lower bound on -log(-x), thus an upper bound on log(-x), - thus an upper bound on -x. */ - mpfr_neg (l, x, GMP_RNDU); /* upper bound on -x */ - mpfr_log (l, l, GMP_RNDU); /* upper bound for log(-x) */ - mpfr_neg (l, l, GMP_RNDD); /* lower bound for -log(-x) */ - mpfr_neg (h, x, GMP_RNDD); /* lower bound on -x */ - mpfr_log (h, h, GMP_RNDD); /* lower bound on log(-x) */ - mpfr_neg (h, h, GMP_RNDU); /* upper bound for -log(-x) */ - mpfr_sub (h, h, x, GMP_RNDU); /* upper bound for -log(-x) - x */ - inex = mpfr_prec_round (l, MPFR_PREC(y), rnd); - inex2 = mpfr_prec_round (h, MPFR_PREC(y), rnd); - /* Caution: we not only need l = h, but both inexact flags should - agree. Indeed, one of the inexact flags might be zero. In that - case if we assume ln|gamma(x)| cannot be exact, the other flag - should be correct. We are conservative here and request that both - inexact flags agree. */ - ok = SAME_SIGN (inex, inex2) && mpfr_cmp (l, h) == 0; - if (ok) - mpfr_set (y, h, rnd); /* exact */ - mpfr_clear (l); - mpfr_clear (h); - if (ok) - return inex; - w += MPFR_INT_CEIL_LOG2 (w) + 3; - } + mp_prec_t w = MPFR_PREC (y) + 14; + + while (1) + { + mpfr_init2 (l, w); + mpfr_init2 (h, w); + /* we want a lower bound on -log(-x), thus an upper bound + on log(-x), thus an upper bound on -x. */ + mpfr_neg (l, x, GMP_RNDU); /* upper bound on -x */ + mpfr_log (l, l, GMP_RNDU); /* upper bound for log(-x) */ + mpfr_neg (l, l, GMP_RNDD); /* lower bound for -log(-x) */ + mpfr_neg (h, x, GMP_RNDD); /* lower bound on -x */ + mpfr_log (h, h, GMP_RNDD); /* lower bound on log(-x) */ + mpfr_neg (h, h, GMP_RNDU); /* upper bound for -log(-x) */ + mpfr_sub (h, h, x, GMP_RNDU); /* upper bound for -log(-x) - x */ + inex = mpfr_prec_round (l, MPFR_PREC (y), rnd); + inex2 = mpfr_prec_round (h, MPFR_PREC (y), rnd); + /* Caution: we not only need l = h, but both inexact flags + should agree. Indeed, one of the inexact flags might be + zero. In that case if we assume ln|gamma(x)| cannot be + exact, the other flag should be correct. We are conservative + here and request that both inexact flags agree. */ + ok = SAME_SIGN (inex, inex2) && mpfr_equal_p (l, h); + if (ok) + mpfr_set (y, h, rnd); /* exact */ + mpfr_clear (l); + mpfr_clear (h); + if (ok) + return inex; + w += MPFR_INT_CEIL_LOG2 (w) + 3; + } } } |