14a238c70SJohn Marino /* mpfr_gamma -- gamma function
24a238c70SJohn Marino
3*ab6d115fSJohn Marino Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4*ab6d115fSJohn Marino Contributed by the AriC and Caramel projects, INRIA.
54a238c70SJohn Marino
64a238c70SJohn Marino This file is part of the GNU MPFR Library.
74a238c70SJohn Marino
84a238c70SJohn Marino The GNU MPFR Library is free software; you can redistribute it and/or modify
94a238c70SJohn Marino it under the terms of the GNU Lesser General Public License as published by
104a238c70SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
114a238c70SJohn Marino option) any later version.
124a238c70SJohn Marino
134a238c70SJohn Marino The GNU MPFR Library is distributed in the hope that it will be useful, but
144a238c70SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
154a238c70SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
164a238c70SJohn Marino License for more details.
174a238c70SJohn Marino
184a238c70SJohn Marino You should have received a copy of the GNU Lesser General Public License
194a238c70SJohn Marino along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
204a238c70SJohn Marino http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
214a238c70SJohn Marino 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
224a238c70SJohn Marino
234a238c70SJohn Marino #define MPFR_NEED_LONGLONG_H
244a238c70SJohn Marino #include "mpfr-impl.h"
254a238c70SJohn Marino
264a238c70SJohn Marino #define IS_GAMMA
274a238c70SJohn Marino #include "lngamma.c"
284a238c70SJohn Marino #undef IS_GAMMA
294a238c70SJohn Marino
304a238c70SJohn Marino /* return a sufficient precision such that 2-x is exact, assuming x < 0 */
314a238c70SJohn Marino static mpfr_prec_t
mpfr_gamma_2_minus_x_exact(mpfr_srcptr x)324a238c70SJohn Marino mpfr_gamma_2_minus_x_exact (mpfr_srcptr x)
334a238c70SJohn Marino {
344a238c70SJohn Marino /* Since x < 0, 2-x = 2+y with y := -x.
354a238c70SJohn Marino If y < 2, a precision w >= PREC(y) + EXP(2)-EXP(y) = PREC(y) + 2 - EXP(y)
364a238c70SJohn Marino is enough, since no overlap occurs in 2+y, so no carry happens.
374a238c70SJohn Marino If y >= 2, either ULP(y) <= 2, and we need w >= PREC(y)+1 since a
384a238c70SJohn Marino carry can occur, or ULP(y) > 2, and we need w >= EXP(y)-1:
394a238c70SJohn Marino (a) if EXP(y) <= 1, w = PREC(y) + 2 - EXP(y)
404a238c70SJohn Marino (b) if EXP(y) > 1 and EXP(y)-PREC(y) <= 1, w = PREC(y) + 1
414a238c70SJohn Marino (c) if EXP(y) > 1 and EXP(y)-PREC(y) > 1, w = EXP(y) - 1 */
424a238c70SJohn Marino return (MPFR_GET_EXP(x) <= 1) ? MPFR_PREC(x) + 2 - MPFR_GET_EXP(x)
434a238c70SJohn Marino : ((MPFR_GET_EXP(x) <= MPFR_PREC(x) + 1) ? MPFR_PREC(x) + 1
444a238c70SJohn Marino : MPFR_GET_EXP(x) - 1);
454a238c70SJohn Marino }
464a238c70SJohn Marino
474a238c70SJohn Marino /* return a sufficient precision such that 1-x is exact, assuming x < 1 */
484a238c70SJohn Marino static mpfr_prec_t
mpfr_gamma_1_minus_x_exact(mpfr_srcptr x)494a238c70SJohn Marino mpfr_gamma_1_minus_x_exact (mpfr_srcptr x)
504a238c70SJohn Marino {
514a238c70SJohn Marino if (MPFR_IS_POS(x))
524a238c70SJohn Marino return MPFR_PREC(x) - MPFR_GET_EXP(x);
534a238c70SJohn Marino else if (MPFR_GET_EXP(x) <= 0)
544a238c70SJohn Marino return MPFR_PREC(x) + 1 - MPFR_GET_EXP(x);
554a238c70SJohn Marino else if (MPFR_PREC(x) >= MPFR_GET_EXP(x))
564a238c70SJohn Marino return MPFR_PREC(x) + 1;
574a238c70SJohn Marino else
584a238c70SJohn Marino return MPFR_GET_EXP(x);
594a238c70SJohn Marino }
604a238c70SJohn Marino
614a238c70SJohn Marino /* returns a lower bound of the number of significant bits of n!
624a238c70SJohn Marino (not counting the low zero bits).
634a238c70SJohn Marino We know n! >= (n/e)^n*sqrt(2*Pi*n) for n >= 1, and the number of zero bits
644a238c70SJohn Marino is floor(n/2) + floor(n/4) + floor(n/8) + ...
654a238c70SJohn Marino This approximation is exact for n <= 500000, except for n = 219536, 235928,
664a238c70SJohn Marino 298981, 355854, 464848, 493725, 498992 where it returns a value 1 too small.
674a238c70SJohn Marino */
684a238c70SJohn Marino static unsigned long
bits_fac(unsigned long n)694a238c70SJohn Marino bits_fac (unsigned long n)
704a238c70SJohn Marino {
714a238c70SJohn Marino mpfr_t x, y;
724a238c70SJohn Marino unsigned long r, k;
734a238c70SJohn Marino mpfr_init2 (x, 38);
744a238c70SJohn Marino mpfr_init2 (y, 38);
754a238c70SJohn Marino mpfr_set_ui (x, n, MPFR_RNDZ);
764a238c70SJohn Marino mpfr_set_str_binary (y, "10.101101111110000101010001011000101001"); /* upper bound of e */
774a238c70SJohn Marino mpfr_div (x, x, y, MPFR_RNDZ);
784a238c70SJohn Marino mpfr_pow_ui (x, x, n, MPFR_RNDZ);
794a238c70SJohn Marino mpfr_const_pi (y, MPFR_RNDZ);
804a238c70SJohn Marino mpfr_mul_ui (y, y, 2 * n, MPFR_RNDZ);
814a238c70SJohn Marino mpfr_sqrt (y, y, MPFR_RNDZ);
824a238c70SJohn Marino mpfr_mul (x, x, y, MPFR_RNDZ);
834a238c70SJohn Marino mpfr_log2 (x, x, MPFR_RNDZ);
844a238c70SJohn Marino r = mpfr_get_ui (x, MPFR_RNDU);
854a238c70SJohn Marino for (k = 2; k <= n; k *= 2)
864a238c70SJohn Marino r -= n / k;
874a238c70SJohn Marino mpfr_clear (x);
884a238c70SJohn Marino mpfr_clear (y);
894a238c70SJohn Marino return r;
904a238c70SJohn Marino }
914a238c70SJohn Marino
924a238c70SJohn Marino /* We use the reflection formula
934a238c70SJohn Marino Gamma(1+t) Gamma(1-t) = - Pi t / sin(Pi (1 + t))
944a238c70SJohn Marino in order to treat the case x <= 1,
954a238c70SJohn Marino i.e. with x = 1-t, then Gamma(x) = -Pi*(1-x)/sin(Pi*(2-x))/GAMMA(2-x)
964a238c70SJohn Marino */
974a238c70SJohn Marino int
mpfr_gamma(mpfr_ptr gamma,mpfr_srcptr x,mpfr_rnd_t rnd_mode)984a238c70SJohn Marino mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
994a238c70SJohn Marino {
1004a238c70SJohn Marino mpfr_t xp, GammaTrial, tmp, tmp2;
1014a238c70SJohn Marino mpz_t fact;
1024a238c70SJohn Marino mpfr_prec_t realprec;
103*ab6d115fSJohn Marino int compared, is_integer;
104*ab6d115fSJohn Marino int inex = 0; /* 0 means: result gamma not set yet */
1054a238c70SJohn Marino MPFR_GROUP_DECL (group);
1064a238c70SJohn Marino MPFR_SAVE_EXPO_DECL (expo);
1074a238c70SJohn Marino MPFR_ZIV_DECL (loop);
1084a238c70SJohn Marino
1094a238c70SJohn Marino MPFR_LOG_FUNC
1104a238c70SJohn Marino (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
1114a238c70SJohn Marino ("gamma[%Pu]=%.*Rg inexact=%d",
1124a238c70SJohn Marino mpfr_get_prec (gamma), mpfr_log_prec, gamma, inex));
1134a238c70SJohn Marino
1144a238c70SJohn Marino /* Trivial cases */
1154a238c70SJohn Marino if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
1164a238c70SJohn Marino {
1174a238c70SJohn Marino if (MPFR_IS_NAN (x))
1184a238c70SJohn Marino {
1194a238c70SJohn Marino MPFR_SET_NAN (gamma);
1204a238c70SJohn Marino MPFR_RET_NAN;
1214a238c70SJohn Marino }
1224a238c70SJohn Marino else if (MPFR_IS_INF (x))
1234a238c70SJohn Marino {
1244a238c70SJohn Marino if (MPFR_IS_NEG (x))
1254a238c70SJohn Marino {
1264a238c70SJohn Marino MPFR_SET_NAN (gamma);
1274a238c70SJohn Marino MPFR_RET_NAN;
1284a238c70SJohn Marino }
1294a238c70SJohn Marino else
1304a238c70SJohn Marino {
1314a238c70SJohn Marino MPFR_SET_INF (gamma);
1324a238c70SJohn Marino MPFR_SET_POS (gamma);
1334a238c70SJohn Marino MPFR_RET (0); /* exact */
1344a238c70SJohn Marino }
1354a238c70SJohn Marino }
1364a238c70SJohn Marino else /* x is zero */
1374a238c70SJohn Marino {
1384a238c70SJohn Marino MPFR_ASSERTD(MPFR_IS_ZERO(x));
1394a238c70SJohn Marino MPFR_SET_INF(gamma);
1404a238c70SJohn Marino MPFR_SET_SAME_SIGN(gamma, x);
1414a238c70SJohn Marino mpfr_set_divby0 ();
1424a238c70SJohn Marino MPFR_RET (0); /* exact */
1434a238c70SJohn Marino }
1444a238c70SJohn Marino }
1454a238c70SJohn Marino
1464a238c70SJohn Marino /* Check for tiny arguments, where gamma(x) ~ 1/x - euler + ....
1474a238c70SJohn Marino We know from "Bound on Runs of Zeros and Ones for Algebraic Functions",
1484a238c70SJohn Marino Proceedings of Arith15, T. Lang and J.-M. Muller, 2001, that the maximal
1494a238c70SJohn Marino number of consecutive zeroes or ones after the round bit is n-1 for an
1504a238c70SJohn Marino input of n bits. But we need a more precise lower bound. Assume x has
1514a238c70SJohn Marino n bits, and 1/x is near a floating-point number y of n+1 bits. We can
1524a238c70SJohn Marino write x = X*2^e, y = Y/2^f with X, Y integers of n and n+1 bits.
1534a238c70SJohn Marino Thus X*Y^2^(e-f) is near from 1, i.e., X*Y is near from 2^(f-e).
1544a238c70SJohn Marino Two cases can happen:
1554a238c70SJohn Marino (i) either X*Y is exactly 2^(f-e), but this can happen only if X and Y
1564a238c70SJohn Marino are themselves powers of two, i.e., x is a power of two;
1574a238c70SJohn Marino (ii) or X*Y is at distance at least one from 2^(f-e), thus
1584a238c70SJohn Marino |xy-1| >= 2^(e-f), or |y-1/x| >= 2^(e-f)/x = 2^(-f)/X >= 2^(-f-n).
1594a238c70SJohn Marino Since ufp(y) = 2^(n-f) [ufp = unit in first place], this means
1604a238c70SJohn Marino that the distance |y-1/x| >= 2^(-2n) ufp(y).
1614a238c70SJohn Marino Now assuming |gamma(x)-1/x| <= 1, which is true for x <= 1,
1624a238c70SJohn Marino if 2^(-2n) ufp(y) >= 2, the error is at most 2^(-2n-1) ufp(y),
1634a238c70SJohn Marino and round(1/x) with precision >= 2n+2 gives the correct result.
1644a238c70SJohn Marino If x < 2^E, then y > 2^(-E), thus ufp(y) > 2^(-E-1).
1654a238c70SJohn Marino A sufficient condition is thus EXP(x) + 2 <= -2 MAX(PREC(x),PREC(Y)).
1664a238c70SJohn Marino */
1674a238c70SJohn Marino if (MPFR_GET_EXP (x) + 2
1684a238c70SJohn Marino <= -2 * (mpfr_exp_t) MAX(MPFR_PREC(x), MPFR_PREC(gamma)))
1694a238c70SJohn Marino {
1704a238c70SJohn Marino int sign = MPFR_SIGN (x); /* retrieve sign before possible override */
1714a238c70SJohn Marino int special;
1724a238c70SJohn Marino MPFR_BLOCK_DECL (flags);
1734a238c70SJohn Marino
1744a238c70SJohn Marino MPFR_SAVE_EXPO_MARK (expo);
1754a238c70SJohn Marino
1764a238c70SJohn Marino /* for overflow cases, see below; this needs to be done
1774a238c70SJohn Marino before x possibly gets overridden. */
1784a238c70SJohn Marino special =
1794a238c70SJohn Marino MPFR_GET_EXP (x) == 1 - MPFR_EMAX_MAX &&
1804a238c70SJohn Marino MPFR_IS_POS_SIGN (sign) &&
1814a238c70SJohn Marino MPFR_IS_LIKE_RNDD (rnd_mode, sign) &&
1824a238c70SJohn Marino mpfr_powerof2_raw (x);
1834a238c70SJohn Marino
1844a238c70SJohn Marino MPFR_BLOCK (flags, inex = mpfr_ui_div (gamma, 1, x, rnd_mode));
1854a238c70SJohn Marino if (inex == 0) /* x is a power of two */
1864a238c70SJohn Marino {
1874a238c70SJohn Marino /* return RND(1/x - euler) = RND(+/- 2^k - eps) with eps > 0 */
1884a238c70SJohn Marino if (rnd_mode == MPFR_RNDN || MPFR_IS_LIKE_RNDU (rnd_mode, sign))
1894a238c70SJohn Marino inex = 1;
1904a238c70SJohn Marino else
1914a238c70SJohn Marino {
1924a238c70SJohn Marino mpfr_nextbelow (gamma);
1934a238c70SJohn Marino inex = -1;
1944a238c70SJohn Marino }
1954a238c70SJohn Marino }
1964a238c70SJohn Marino else if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
1974a238c70SJohn Marino {
1984a238c70SJohn Marino /* Overflow in the division 1/x. This is a real overflow, except
1994a238c70SJohn Marino in RNDZ or RNDD when 1/x = 2^emax, i.e. x = 2^(-emax): due to
2004a238c70SJohn Marino the "- euler", the rounded value in unbounded exponent range
2014a238c70SJohn Marino is 0.111...11 * 2^emax (not an overflow). */
2024a238c70SJohn Marino if (!special)
2034a238c70SJohn Marino MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, flags);
2044a238c70SJohn Marino }
2054a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
2064a238c70SJohn Marino /* Note: an overflow is possible with an infinite result;
2074a238c70SJohn Marino in this case, the overflow flag will automatically be
2084a238c70SJohn Marino restored by mpfr_check_range. */
2094a238c70SJohn Marino return mpfr_check_range (gamma, inex, rnd_mode);
2104a238c70SJohn Marino }
2114a238c70SJohn Marino
2124a238c70SJohn Marino is_integer = mpfr_integer_p (x);
2134a238c70SJohn Marino /* gamma(x) for x a negative integer gives NaN */
2144a238c70SJohn Marino if (is_integer && MPFR_IS_NEG(x))
2154a238c70SJohn Marino {
2164a238c70SJohn Marino MPFR_SET_NAN (gamma);
2174a238c70SJohn Marino MPFR_RET_NAN;
2184a238c70SJohn Marino }
2194a238c70SJohn Marino
2204a238c70SJohn Marino compared = mpfr_cmp_ui (x, 1);
2214a238c70SJohn Marino if (compared == 0)
2224a238c70SJohn Marino return mpfr_set_ui (gamma, 1, rnd_mode);
2234a238c70SJohn Marino
2244a238c70SJohn Marino /* if x is an integer that fits into an unsigned long, use mpfr_fac_ui
2254a238c70SJohn Marino if argument is not too large.
2264a238c70SJohn Marino If precision is p, fac_ui costs O(u*p), whereas gamma costs O(p*M(p)),
2274a238c70SJohn Marino so for u <= M(p), fac_ui should be faster.
2284a238c70SJohn Marino We approximate here M(p) by p*log(p)^2, which is not a bad guess.
2294a238c70SJohn Marino Warning: since the generic code does not handle exact cases,
2304a238c70SJohn Marino we want all cases where gamma(x) is exact to be treated here.
2314a238c70SJohn Marino */
2324a238c70SJohn Marino if (is_integer && mpfr_fits_ulong_p (x, MPFR_RNDN))
2334a238c70SJohn Marino {
2344a238c70SJohn Marino unsigned long int u;
2354a238c70SJohn Marino mpfr_prec_t p = MPFR_PREC(gamma);
2364a238c70SJohn Marino u = mpfr_get_ui (x, MPFR_RNDN);
2374a238c70SJohn Marino if (u < 44787929UL && bits_fac (u - 1) <= p + (rnd_mode == MPFR_RNDN))
2384a238c70SJohn Marino /* bits_fac: lower bound on the number of bits of m,
2394a238c70SJohn Marino where gamma(x) = (u-1)! = m*2^e with m odd. */
2404a238c70SJohn Marino return mpfr_fac_ui (gamma, u - 1, rnd_mode);
2414a238c70SJohn Marino /* if bits_fac(...) > p (resp. p+1 for rounding to nearest),
2424a238c70SJohn Marino then gamma(x) cannot be exact in precision p (resp. p+1).
2434a238c70SJohn Marino FIXME: remove the test u < 44787929UL after changing bits_fac
2444a238c70SJohn Marino to return a mpz_t or mpfr_t. */
2454a238c70SJohn Marino }
2464a238c70SJohn Marino
2474a238c70SJohn Marino MPFR_SAVE_EXPO_MARK (expo);
2484a238c70SJohn Marino
2494a238c70SJohn Marino /* check for overflow: according to (6.1.37) in Abramowitz & Stegun,
2504a238c70SJohn Marino gamma(x) >= exp(-x) * x^(x-1/2) * sqrt(2*Pi)
2514a238c70SJohn Marino >= 2 * (x/e)^x / x for x >= 1 */
2524a238c70SJohn Marino if (compared > 0)
2534a238c70SJohn Marino {
2544a238c70SJohn Marino mpfr_t yp;
2554a238c70SJohn Marino mpfr_exp_t expxp;
2564a238c70SJohn Marino MPFR_BLOCK_DECL (flags);
2574a238c70SJohn Marino
2584a238c70SJohn Marino /* 1/e rounded down to 53 bits */
2594a238c70SJohn Marino #define EXPM1_STR "0.010111100010110101011000110110001011001110111100111"
2604a238c70SJohn Marino mpfr_init2 (xp, 53);
2614a238c70SJohn Marino mpfr_init2 (yp, 53);
2624a238c70SJohn Marino mpfr_set_str_binary (xp, EXPM1_STR);
2634a238c70SJohn Marino mpfr_mul (xp, x, xp, MPFR_RNDZ);
2644a238c70SJohn Marino mpfr_sub_ui (yp, x, 2, MPFR_RNDZ);
2654a238c70SJohn Marino mpfr_pow (xp, xp, yp, MPFR_RNDZ); /* (x/e)^(x-2) */
2664a238c70SJohn Marino mpfr_set_str_binary (yp, EXPM1_STR);
2674a238c70SJohn Marino mpfr_mul (xp, xp, yp, MPFR_RNDZ); /* x^(x-2) / e^(x-1) */
2684a238c70SJohn Marino mpfr_mul (xp, xp, yp, MPFR_RNDZ); /* x^(x-2) / e^x */
2694a238c70SJohn Marino mpfr_mul (xp, xp, x, MPFR_RNDZ); /* lower bound on x^(x-1) / e^x */
2704a238c70SJohn Marino MPFR_BLOCK (flags, mpfr_mul_2ui (xp, xp, 1, MPFR_RNDZ));
2714a238c70SJohn Marino expxp = MPFR_GET_EXP (xp);
2724a238c70SJohn Marino mpfr_clear (xp);
2734a238c70SJohn Marino mpfr_clear (yp);
2744a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
2754a238c70SJohn Marino return MPFR_OVERFLOW (flags) || expxp > __gmpfr_emax ?
2764a238c70SJohn Marino mpfr_overflow (gamma, rnd_mode, 1) :
2774a238c70SJohn Marino mpfr_gamma_aux (gamma, x, rnd_mode);
2784a238c70SJohn Marino }
2794a238c70SJohn Marino
2804a238c70SJohn Marino /* now compared < 0 */
2814a238c70SJohn Marino
2824a238c70SJohn Marino /* check for underflow: for x < 1,
2834a238c70SJohn Marino gamma(x) = Pi*(x-1)/sin(Pi*(2-x))/gamma(2-x).
2844a238c70SJohn Marino Since gamma(2-x) >= 2 * ((2-x)/e)^(2-x) / (2-x), we have
2854a238c70SJohn Marino |gamma(x)| <= Pi*(1-x)*(2-x)/2/((2-x)/e)^(2-x) / |sin(Pi*(2-x))|
2864a238c70SJohn Marino <= 12 * ((2-x)/e)^x / |sin(Pi*(2-x))|.
2874a238c70SJohn Marino To avoid an underflow in ((2-x)/e)^x, we compute the logarithm.
2884a238c70SJohn Marino */
2894a238c70SJohn Marino if (MPFR_IS_NEG(x))
2904a238c70SJohn Marino {
2914a238c70SJohn Marino int underflow = 0, sgn, ck;
2924a238c70SJohn Marino mpfr_prec_t w;
2934a238c70SJohn Marino
2944a238c70SJohn Marino mpfr_init2 (xp, 53);
2954a238c70SJohn Marino mpfr_init2 (tmp, 53);
2964a238c70SJohn Marino mpfr_init2 (tmp2, 53);
2974a238c70SJohn Marino /* we want an upper bound for x * [log(2-x)-1].
2984a238c70SJohn Marino since x < 0, we need a lower bound on log(2-x) */
2994a238c70SJohn Marino mpfr_ui_sub (xp, 2, x, MPFR_RNDD);
300*ab6d115fSJohn Marino mpfr_log (xp, xp, MPFR_RNDD);
3014a238c70SJohn Marino mpfr_sub_ui (xp, xp, 1, MPFR_RNDD);
3024a238c70SJohn Marino mpfr_mul (xp, xp, x, MPFR_RNDU);
3034a238c70SJohn Marino
3044a238c70SJohn Marino /* we need an upper bound on 1/|sin(Pi*(2-x))|,
3054a238c70SJohn Marino thus a lower bound on |sin(Pi*(2-x))|.
3064a238c70SJohn Marino If 2-x is exact, then the error of Pi*(2-x) is (1+u)^2 with u = 2^(-p)
3074a238c70SJohn Marino thus the error on sin(Pi*(2-x)) is less than 1/2ulp + 3Pi(2-x)u,
3084a238c70SJohn Marino assuming u <= 1, thus <= u + 3Pi(2-x)u */
3094a238c70SJohn Marino
3104a238c70SJohn Marino w = mpfr_gamma_2_minus_x_exact (x); /* 2-x is exact for prec >= w */
3114a238c70SJohn Marino w += 17; /* to get tmp2 small enough */
3124a238c70SJohn Marino mpfr_set_prec (tmp, w);
3134a238c70SJohn Marino mpfr_set_prec (tmp2, w);
3144a238c70SJohn Marino ck = mpfr_ui_sub (tmp, 2, x, MPFR_RNDN);
3154a238c70SJohn Marino MPFR_ASSERTD (ck == 0); (void) ck; /* use ck to avoid a warning */
3164a238c70SJohn Marino mpfr_const_pi (tmp2, MPFR_RNDN);
3174a238c70SJohn Marino mpfr_mul (tmp2, tmp2, tmp, MPFR_RNDN); /* Pi*(2-x) */
3184a238c70SJohn Marino mpfr_sin (tmp, tmp2, MPFR_RNDN); /* sin(Pi*(2-x)) */
3194a238c70SJohn Marino sgn = mpfr_sgn (tmp);
3204a238c70SJohn Marino mpfr_abs (tmp, tmp, MPFR_RNDN);
3214a238c70SJohn Marino mpfr_mul_ui (tmp2, tmp2, 3, MPFR_RNDU); /* 3Pi(2-x) */
3224a238c70SJohn Marino mpfr_add_ui (tmp2, tmp2, 1, MPFR_RNDU); /* 3Pi(2-x)+1 */
3234a238c70SJohn Marino mpfr_div_2ui (tmp2, tmp2, mpfr_get_prec (tmp), MPFR_RNDU);
3244a238c70SJohn Marino /* if tmp2<|tmp|, we get a lower bound */
3254a238c70SJohn Marino if (mpfr_cmp (tmp2, tmp) < 0)
3264a238c70SJohn Marino {
3274a238c70SJohn Marino mpfr_sub (tmp, tmp, tmp2, MPFR_RNDZ); /* low bnd on |sin(Pi*(2-x))| */
3284a238c70SJohn Marino mpfr_ui_div (tmp, 12, tmp, MPFR_RNDU); /* upper bound */
3294a238c70SJohn Marino mpfr_log2 (tmp, tmp, MPFR_RNDU);
3304a238c70SJohn Marino mpfr_add (xp, tmp, xp, MPFR_RNDU);
3314a238c70SJohn Marino /* The assert below checks that expo.saved_emin - 2 always
3324a238c70SJohn Marino fits in a long. FIXME if we want to allow mpfr_exp_t to
3334a238c70SJohn Marino be a long long, for instance. */
3344a238c70SJohn Marino MPFR_ASSERTN (MPFR_EMIN_MIN - 2 >= LONG_MIN);
3354a238c70SJohn Marino underflow = mpfr_cmp_si (xp, expo.saved_emin - 2) <= 0;
3364a238c70SJohn Marino }
3374a238c70SJohn Marino
3384a238c70SJohn Marino mpfr_clear (xp);
3394a238c70SJohn Marino mpfr_clear (tmp);
3404a238c70SJohn Marino mpfr_clear (tmp2);
3414a238c70SJohn Marino if (underflow) /* the sign is the opposite of that of sin(Pi*(2-x)) */
3424a238c70SJohn Marino {
3434a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
3444a238c70SJohn Marino return mpfr_underflow (gamma, (rnd_mode == MPFR_RNDN) ? MPFR_RNDZ : rnd_mode, -sgn);
3454a238c70SJohn Marino }
3464a238c70SJohn Marino }
3474a238c70SJohn Marino
3484a238c70SJohn Marino realprec = MPFR_PREC (gamma);
3494a238c70SJohn Marino /* we want both 1-x and 2-x to be exact */
3504a238c70SJohn Marino {
3514a238c70SJohn Marino mpfr_prec_t w;
3524a238c70SJohn Marino w = mpfr_gamma_1_minus_x_exact (x);
3534a238c70SJohn Marino if (realprec < w)
3544a238c70SJohn Marino realprec = w;
3554a238c70SJohn Marino w = mpfr_gamma_2_minus_x_exact (x);
3564a238c70SJohn Marino if (realprec < w)
3574a238c70SJohn Marino realprec = w;
3584a238c70SJohn Marino }
3594a238c70SJohn Marino realprec = realprec + MPFR_INT_CEIL_LOG2 (realprec) + 20;
3604a238c70SJohn Marino MPFR_ASSERTD(realprec >= 5);
3614a238c70SJohn Marino
3624a238c70SJohn Marino MPFR_GROUP_INIT_4 (group, realprec + MPFR_INT_CEIL_LOG2 (realprec) + 20,
3634a238c70SJohn Marino xp, tmp, tmp2, GammaTrial);
3644a238c70SJohn Marino mpz_init (fact);
3654a238c70SJohn Marino MPFR_ZIV_INIT (loop, realprec);
3664a238c70SJohn Marino for (;;)
3674a238c70SJohn Marino {
3684a238c70SJohn Marino mpfr_exp_t err_g;
3694a238c70SJohn Marino int ck;
3704a238c70SJohn Marino MPFR_GROUP_REPREC_4 (group, realprec, xp, tmp, tmp2, GammaTrial);
3714a238c70SJohn Marino
3724a238c70SJohn Marino /* reflection formula: gamma(x) = Pi*(x-1)/sin(Pi*(2-x))/gamma(2-x) */
3734a238c70SJohn Marino
3744a238c70SJohn Marino ck = mpfr_ui_sub (xp, 2, x, MPFR_RNDN); /* 2-x, exact */
3754a238c70SJohn Marino MPFR_ASSERTD(ck == 0); (void) ck; /* use ck to avoid a warning */
3764a238c70SJohn Marino mpfr_gamma (tmp, xp, MPFR_RNDN); /* gamma(2-x), error (1+u) */
3774a238c70SJohn Marino mpfr_const_pi (tmp2, MPFR_RNDN); /* Pi, error (1+u) */
3784a238c70SJohn Marino mpfr_mul (GammaTrial, tmp2, xp, MPFR_RNDN); /* Pi*(2-x), error (1+u)^2 */
3794a238c70SJohn Marino err_g = MPFR_GET_EXP(GammaTrial);
3804a238c70SJohn Marino mpfr_sin (GammaTrial, GammaTrial, MPFR_RNDN); /* sin(Pi*(2-x)) */
381*ab6d115fSJohn Marino /* If tmp is +Inf, we compute exp(lngamma(x)). */
382*ab6d115fSJohn Marino if (mpfr_inf_p (tmp))
383*ab6d115fSJohn Marino {
384*ab6d115fSJohn Marino inex = mpfr_explgamma (gamma, x, &expo, tmp, tmp2, rnd_mode);
385*ab6d115fSJohn Marino if (inex)
386*ab6d115fSJohn Marino goto end;
387*ab6d115fSJohn Marino else
388*ab6d115fSJohn Marino goto ziv_next;
389*ab6d115fSJohn Marino }
3904a238c70SJohn Marino err_g = err_g + 1 - MPFR_GET_EXP(GammaTrial);
3914a238c70SJohn Marino /* let g0 the true value of Pi*(2-x), g the computed value.
3924a238c70SJohn Marino We have g = g0 + h with |h| <= |(1+u^2)-1|*g.
3934a238c70SJohn Marino Thus sin(g) = sin(g0) + h' with |h'| <= |(1+u^2)-1|*g.
3944a238c70SJohn Marino The relative error is thus bounded by |(1+u^2)-1|*g/sin(g)
3954a238c70SJohn Marino <= |(1+u^2)-1|*2^err_g. <= 2.25*u*2^err_g for |u|<=1/4.
3964a238c70SJohn Marino With the rounding error, this gives (0.5 + 2.25*2^err_g)*u. */
3974a238c70SJohn Marino ck = mpfr_sub_ui (xp, x, 1, MPFR_RNDN); /* x-1, exact */
3984a238c70SJohn Marino MPFR_ASSERTD(ck == 0); (void) ck; /* use ck to avoid a warning */
3994a238c70SJohn Marino mpfr_mul (xp, tmp2, xp, MPFR_RNDN); /* Pi*(x-1), error (1+u)^2 */
4004a238c70SJohn Marino mpfr_mul (GammaTrial, GammaTrial, tmp, MPFR_RNDN);
4014a238c70SJohn Marino /* [1 + (0.5 + 2.25*2^err_g)*u]*(1+u)^2 = 1 + (2.5 + 2.25*2^err_g)*u
4024a238c70SJohn Marino + (0.5 + 2.25*2^err_g)*u*(2u+u^2) + u^2.
4034a238c70SJohn Marino For err_g <= realprec-2, we have (0.5 + 2.25*2^err_g)*u <=
4044a238c70SJohn Marino 0.5*u + 2.25/4 <= 0.6875 and u^2 <= u/4, thus
4054a238c70SJohn Marino (0.5 + 2.25*2^err_g)*u*(2u+u^2) + u^2 <= 0.6875*(2u+u/4) + u/4
4064a238c70SJohn Marino <= 1.8*u, thus the rel. error is bounded by (4.5 + 2.25*2^err_g)*u. */
4074a238c70SJohn Marino mpfr_div (GammaTrial, xp, GammaTrial, MPFR_RNDN);
4084a238c70SJohn Marino /* the error is of the form (1+u)^3/[1 + (4.5 + 2.25*2^err_g)*u].
4094a238c70SJohn Marino For realprec >= 5 and err_g <= realprec-2, [(4.5 + 2.25*2^err_g)*u]^2
4104a238c70SJohn Marino <= 0.71, and for |y|<=0.71, 1/(1-y) can be written 1+a*y with a<=4.
4114a238c70SJohn Marino (1+u)^3 * (1+4*(4.5 + 2.25*2^err_g)*u)
4124a238c70SJohn Marino = 1 + (21 + 9*2^err_g)*u + (57+27*2^err_g)*u^2 + (55+27*2^err_g)*u^3
4134a238c70SJohn Marino + (18+9*2^err_g)*u^4
4144a238c70SJohn Marino <= 1 + (21 + 9*2^err_g)*u + (57+27*2^err_g)*u^2 + (56+28*2^err_g)*u^3
4154a238c70SJohn Marino <= 1 + (21 + 9*2^err_g)*u + (59+28*2^err_g)*u^2
4164a238c70SJohn Marino <= 1 + (23 + 10*2^err_g)*u.
4174a238c70SJohn Marino The final error is thus bounded by (23 + 10*2^err_g) ulps,
4184a238c70SJohn Marino which is <= 2^6 for err_g<=2, and <= 2^(err_g+4) for err_g >= 2. */
4194a238c70SJohn Marino err_g = (err_g <= 2) ? 6 : err_g + 4;
4204a238c70SJohn Marino
4214a238c70SJohn Marino if (MPFR_LIKELY (MPFR_CAN_ROUND (GammaTrial, realprec - err_g,
4224a238c70SJohn Marino MPFR_PREC(gamma), rnd_mode)))
4234a238c70SJohn Marino break;
424*ab6d115fSJohn Marino
425*ab6d115fSJohn Marino ziv_next:
4264a238c70SJohn Marino MPFR_ZIV_NEXT (loop, realprec);
4274a238c70SJohn Marino }
428*ab6d115fSJohn Marino
429*ab6d115fSJohn Marino end:
4304a238c70SJohn Marino MPFR_ZIV_FREE (loop);
4314a238c70SJohn Marino
432*ab6d115fSJohn Marino if (inex == 0)
4334a238c70SJohn Marino inex = mpfr_set (gamma, GammaTrial, rnd_mode);
4344a238c70SJohn Marino MPFR_GROUP_CLEAR (group);
4354a238c70SJohn Marino mpz_clear (fact);
4364a238c70SJohn Marino
4374a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
4384a238c70SJohn Marino return mpfr_check_range (gamma, inex, rnd_mode);
4394a238c70SJohn Marino }
440