14a238c70SJohn Marino /* mpfr_lngamma -- lngamma function
24a238c70SJohn Marino
3*ab6d115fSJohn Marino Copyright 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 /* given a precision p, return alpha, such that the argument reduction
274a238c70SJohn Marino will use k = alpha*p*log(2).
284a238c70SJohn Marino
294a238c70SJohn Marino Warning: we should always have alpha >= log(2)/(2Pi) ~ 0.11,
304a238c70SJohn Marino and the smallest value of alpha multiplied by the smallest working
314a238c70SJohn Marino precision should be >= 4.
324a238c70SJohn Marino */
334a238c70SJohn Marino static void
mpfr_gamma_alpha(mpfr_t s,mpfr_prec_t p)344a238c70SJohn Marino mpfr_gamma_alpha (mpfr_t s, mpfr_prec_t p)
354a238c70SJohn Marino {
364a238c70SJohn Marino if (p <= 100)
374a238c70SJohn Marino mpfr_set_ui_2exp (s, 614, -10, MPFR_RNDN); /* about 0.6 */
384a238c70SJohn Marino else if (p <= 500)
394a238c70SJohn Marino mpfr_set_ui_2exp (s, 819, -10, MPFR_RNDN); /* about 0.8 */
404a238c70SJohn Marino else if (p <= 1000)
414a238c70SJohn Marino mpfr_set_ui_2exp (s, 1331, -10, MPFR_RNDN); /* about 1.3 */
424a238c70SJohn Marino else if (p <= 2000)
434a238c70SJohn Marino mpfr_set_ui_2exp (s, 1741, -10, MPFR_RNDN); /* about 1.7 */
444a238c70SJohn Marino else if (p <= 5000)
454a238c70SJohn Marino mpfr_set_ui_2exp (s, 2253, -10, MPFR_RNDN); /* about 2.2 */
464a238c70SJohn Marino else if (p <= 10000)
474a238c70SJohn Marino mpfr_set_ui_2exp (s, 3482, -10, MPFR_RNDN); /* about 3.4 */
484a238c70SJohn Marino else
494a238c70SJohn Marino mpfr_set_ui_2exp (s, 9, -1, MPFR_RNDN); /* 4.5 */
504a238c70SJohn Marino }
514a238c70SJohn Marino
52*ab6d115fSJohn Marino #ifdef IS_GAMMA
53*ab6d115fSJohn Marino
54*ab6d115fSJohn Marino /* This function is called in case of intermediate overflow/underflow.
55*ab6d115fSJohn Marino The s1 and s2 arguments are temporary MPFR numbers, having the
56*ab6d115fSJohn Marino working precision. If the result could be determined, then the
57*ab6d115fSJohn Marino flags are updated via pexpo, y is set to the result, and the
58*ab6d115fSJohn Marino (non-zero) ternary value is returned. Otherwise 0 is returned
59*ab6d115fSJohn Marino in order to perform the next Ziv iteration. */
604a238c70SJohn Marino static int
mpfr_explgamma(mpfr_ptr y,mpfr_srcptr x,mpfr_save_expo_t * pexpo,mpfr_ptr s1,mpfr_ptr s2,mpfr_rnd_t rnd)61*ab6d115fSJohn Marino mpfr_explgamma (mpfr_ptr y, mpfr_srcptr x, mpfr_save_expo_t *pexpo,
62*ab6d115fSJohn Marino mpfr_ptr s1, mpfr_ptr s2, mpfr_rnd_t rnd)
63*ab6d115fSJohn Marino {
64*ab6d115fSJohn Marino mpfr_t t1, t2;
65*ab6d115fSJohn Marino int inex1, inex2, sign;
66*ab6d115fSJohn Marino MPFR_BLOCK_DECL (flags1);
67*ab6d115fSJohn Marino MPFR_BLOCK_DECL (flags2);
68*ab6d115fSJohn Marino MPFR_GROUP_DECL (group);
69*ab6d115fSJohn Marino
70*ab6d115fSJohn Marino MPFR_BLOCK (flags1, inex1 = mpfr_lgamma (s1, &sign, x, MPFR_RNDD));
71*ab6d115fSJohn Marino MPFR_ASSERTN (inex1 != 0);
72*ab6d115fSJohn Marino /* s1 = RNDD(lngamma(x)), inexact */
73*ab6d115fSJohn Marino if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags1)))
74*ab6d115fSJohn Marino {
75*ab6d115fSJohn Marino if (MPFR_SIGN (s1) > 0)
76*ab6d115fSJohn Marino {
77*ab6d115fSJohn Marino MPFR_SAVE_EXPO_UPDATE_FLAGS (*pexpo, MPFR_FLAGS_OVERFLOW);
78*ab6d115fSJohn Marino return mpfr_overflow (y, rnd, sign);
79*ab6d115fSJohn Marino }
80*ab6d115fSJohn Marino else
81*ab6d115fSJohn Marino {
82*ab6d115fSJohn Marino MPFR_SAVE_EXPO_UPDATE_FLAGS (*pexpo, MPFR_FLAGS_UNDERFLOW);
83*ab6d115fSJohn Marino return mpfr_underflow (y, rnd == MPFR_RNDN ? MPFR_RNDZ : rnd, sign);
84*ab6d115fSJohn Marino }
85*ab6d115fSJohn Marino }
86*ab6d115fSJohn Marino
87*ab6d115fSJohn Marino mpfr_set (s2, s1, MPFR_RNDN); /* exact */
88*ab6d115fSJohn Marino mpfr_nextabove (s2); /* v = RNDU(lngamma(z0)) */
89*ab6d115fSJohn Marino
90*ab6d115fSJohn Marino if (sign < 0)
91*ab6d115fSJohn Marino rnd = MPFR_INVERT_RND (rnd); /* since the result with be negated */
92*ab6d115fSJohn Marino MPFR_GROUP_INIT_2 (group, MPFR_PREC (y), t1, t2);
93*ab6d115fSJohn Marino MPFR_BLOCK (flags1, inex1 = mpfr_exp (t1, s1, rnd));
94*ab6d115fSJohn Marino MPFR_BLOCK (flags2, inex2 = mpfr_exp (t2, s2, rnd));
95*ab6d115fSJohn Marino /* t1 is the rounding with mode 'rnd' of a lower bound on |Gamma(x)|,
96*ab6d115fSJohn Marino t2 is the rounding with mode 'rnd' of an upper bound, thus if both
97*ab6d115fSJohn Marino are equal, so is the wanted result. If t1 and t2 differ or the flags
98*ab6d115fSJohn Marino differ, at some point of Ziv's loop they should agree. */
99*ab6d115fSJohn Marino if (mpfr_equal_p (t1, t2) && flags1 == flags2)
100*ab6d115fSJohn Marino {
101*ab6d115fSJohn Marino MPFR_ASSERTN ((inex1 > 0 && inex2 > 0) || (inex1 < 0 && inex2 < 0));
102*ab6d115fSJohn Marino mpfr_set4 (y, t1, MPFR_RNDN, sign); /* exact */
103*ab6d115fSJohn Marino if (sign < 0)
104*ab6d115fSJohn Marino inex1 = - inex1;
105*ab6d115fSJohn Marino MPFR_SAVE_EXPO_UPDATE_FLAGS (*pexpo, flags1);
106*ab6d115fSJohn Marino }
107*ab6d115fSJohn Marino else
108*ab6d115fSJohn Marino inex1 = 0; /* couldn't determine the result */
109*ab6d115fSJohn Marino MPFR_GROUP_CLEAR (group);
110*ab6d115fSJohn Marino
111*ab6d115fSJohn Marino return inex1;
112*ab6d115fSJohn Marino }
113*ab6d115fSJohn Marino
114*ab6d115fSJohn Marino #else
115*ab6d115fSJohn Marino
116*ab6d115fSJohn Marino static int
unit_bit(mpfr_srcptr x)117*ab6d115fSJohn Marino unit_bit (mpfr_srcptr x)
1184a238c70SJohn Marino {
1194a238c70SJohn Marino mpfr_exp_t expo;
1204a238c70SJohn Marino mpfr_prec_t prec;
1214a238c70SJohn Marino mp_limb_t x0;
1224a238c70SJohn Marino
1234a238c70SJohn Marino expo = MPFR_GET_EXP (x);
1244a238c70SJohn Marino if (expo <= 0)
1254a238c70SJohn Marino return 0; /* |x| < 1 */
1264a238c70SJohn Marino
1274a238c70SJohn Marino prec = MPFR_PREC (x);
1284a238c70SJohn Marino if (expo > prec)
1294a238c70SJohn Marino return 0; /* y is a multiple of 2^(expo-prec), thus an even integer */
1304a238c70SJohn Marino
1314a238c70SJohn Marino /* Now, the unit bit is represented. */
1324a238c70SJohn Marino
133*ab6d115fSJohn Marino prec = MPFR_PREC2LIMBS (prec) * GMP_NUMB_BITS - expo;
1344a238c70SJohn Marino /* number of represented fractional bits (including the trailing 0's) */
1354a238c70SJohn Marino
1364a238c70SJohn Marino x0 = *(MPFR_MANT (x) + prec / GMP_NUMB_BITS);
1374a238c70SJohn Marino /* limb containing the unit bit */
1384a238c70SJohn Marino
1394a238c70SJohn Marino return (x0 >> (prec % GMP_NUMB_BITS)) & 1;
1404a238c70SJohn Marino }
141*ab6d115fSJohn Marino
1424a238c70SJohn Marino #endif
1434a238c70SJohn Marino
1444a238c70SJohn Marino /* lngamma(x) = log(gamma(x)).
1454a238c70SJohn Marino We use formula [6.1.40] from Abramowitz&Stegun:
1464a238c70SJohn Marino lngamma(z) = (z-1/2)*log(z) - z + 1/2*log(2*Pi)
1474a238c70SJohn Marino + sum (Bernoulli[2m]/(2m)/(2m-1)/z^(2m-1),m=1..infinity)
1484a238c70SJohn Marino According to [6.1.42], if the sum is truncated after m=n, the error
1494a238c70SJohn Marino R_n(z) is bounded by |B[2n+2]|*K(z)/(2n+1)/(2n+2)/|z|^(2n+1)
1504a238c70SJohn Marino where K(z) = max (z^2/(u^2+z^2)) for u >= 0.
1514a238c70SJohn Marino For z real, |K(z)| <= 1 thus R_n(z) is bounded by the first neglected term.
1524a238c70SJohn Marino */
1534a238c70SJohn Marino #ifdef IS_GAMMA
1544a238c70SJohn Marino #define GAMMA_FUNC mpfr_gamma_aux
1554a238c70SJohn Marino #else
1564a238c70SJohn Marino #define GAMMA_FUNC mpfr_lngamma_aux
1574a238c70SJohn Marino #endif
1584a238c70SJohn Marino
1594a238c70SJohn Marino static int
GAMMA_FUNC(mpfr_ptr y,mpfr_srcptr z0,mpfr_rnd_t rnd)1604a238c70SJohn Marino GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mpfr_rnd_t rnd)
1614a238c70SJohn Marino {
1624a238c70SJohn Marino mpfr_prec_t precy, w; /* working precision */
1634a238c70SJohn Marino mpfr_t s, t, u, v, z;
1644a238c70SJohn Marino unsigned long m, k, maxm;
1654a238c70SJohn Marino mpz_t *INITIALIZED(B); /* variable B declared as initialized */
166*ab6d115fSJohn Marino int compared;
167*ab6d115fSJohn Marino int inexact = 0; /* 0 means: result y not set yet */
1684a238c70SJohn Marino mpfr_exp_t err_s, err_t;
1694a238c70SJohn Marino unsigned long Bm = 0; /* number of allocated B[] */
1704a238c70SJohn Marino unsigned long oldBm;
1714a238c70SJohn Marino double d;
1724a238c70SJohn Marino MPFR_SAVE_EXPO_DECL (expo);
173*ab6d115fSJohn Marino MPFR_ZIV_DECL (loop);
1744a238c70SJohn Marino
1754a238c70SJohn Marino compared = mpfr_cmp_ui (z0, 1);
1764a238c70SJohn Marino
1774a238c70SJohn Marino MPFR_SAVE_EXPO_MARK (expo);
1784a238c70SJohn Marino
1794a238c70SJohn Marino #ifndef IS_GAMMA /* lngamma or lgamma */
1804a238c70SJohn Marino if (compared == 0 || (compared > 0 && mpfr_cmp_ui (z0, 2) == 0))
1814a238c70SJohn Marino {
1824a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
1834a238c70SJohn Marino return mpfr_set_ui (y, 0, MPFR_RNDN); /* lngamma(1 or 2) = +0 */
1844a238c70SJohn Marino }
1854a238c70SJohn Marino
1864a238c70SJohn Marino /* Deal here with tiny inputs. We have for -0.3 <= x <= 0.3:
1874a238c70SJohn Marino - log|x| - gamma*x <= log|gamma(x)| <= - log|x| - gamma*x + x^2 */
1884a238c70SJohn Marino if (MPFR_EXP(z0) <= - (mpfr_exp_t) MPFR_PREC(y))
1894a238c70SJohn Marino {
1904a238c70SJohn Marino mpfr_t l, h, g;
191*ab6d115fSJohn Marino int ok, inex1, inex2;
1924a238c70SJohn Marino mpfr_prec_t prec = MPFR_PREC(y) + 14;
1934a238c70SJohn Marino MPFR_ZIV_DECL (loop);
1944a238c70SJohn Marino
1954a238c70SJohn Marino MPFR_ZIV_INIT (loop, prec);
1964a238c70SJohn Marino do
1974a238c70SJohn Marino {
1984a238c70SJohn Marino mpfr_init2 (l, prec);
1994a238c70SJohn Marino if (MPFR_IS_POS(z0))
2004a238c70SJohn Marino {
2014a238c70SJohn Marino mpfr_log (l, z0, MPFR_RNDU); /* upper bound for log(z0) */
2024a238c70SJohn Marino mpfr_init2 (h, MPFR_PREC(l));
2034a238c70SJohn Marino }
2044a238c70SJohn Marino else
2054a238c70SJohn Marino {
2064a238c70SJohn Marino mpfr_init2 (h, MPFR_PREC(z0));
2074a238c70SJohn Marino mpfr_neg (h, z0, MPFR_RNDN); /* exact */
2084a238c70SJohn Marino mpfr_log (l, h, MPFR_RNDU); /* upper bound for log(-z0) */
2094a238c70SJohn Marino mpfr_set_prec (h, MPFR_PREC(l));
2104a238c70SJohn Marino }
2114a238c70SJohn Marino mpfr_neg (l, l, MPFR_RNDD); /* lower bound for -log(|z0|) */
2124a238c70SJohn Marino mpfr_set (h, l, MPFR_RNDD); /* exact */
2134a238c70SJohn Marino mpfr_nextabove (h); /* upper bound for -log(|z0|), avoids two calls
2144a238c70SJohn Marino to mpfr_log */
2154a238c70SJohn Marino mpfr_init2 (g, MPFR_PREC(l));
2164a238c70SJohn Marino /* if z0>0, we need an upper approximation of Euler's constant
2174a238c70SJohn Marino for the left bound */
2184a238c70SJohn Marino mpfr_const_euler (g, MPFR_IS_POS(z0) ? MPFR_RNDU : MPFR_RNDD);
2194a238c70SJohn Marino mpfr_mul (g, g, z0, MPFR_RNDD);
2204a238c70SJohn Marino mpfr_sub (l, l, g, MPFR_RNDD);
2214a238c70SJohn Marino mpfr_const_euler (g, MPFR_IS_POS(z0) ? MPFR_RNDD : MPFR_RNDU); /* cached */
2224a238c70SJohn Marino mpfr_mul (g, g, z0, MPFR_RNDU);
2234a238c70SJohn Marino mpfr_sub (h, h, g, MPFR_RNDD);
2244a238c70SJohn Marino mpfr_mul (g, z0, z0, MPFR_RNDU);
2254a238c70SJohn Marino mpfr_add (h, h, g, MPFR_RNDU);
226*ab6d115fSJohn Marino inex1 = mpfr_prec_round (l, MPFR_PREC(y), rnd);
2274a238c70SJohn Marino inex2 = mpfr_prec_round (h, MPFR_PREC(y), rnd);
2284a238c70SJohn Marino /* Caution: we not only need l = h, but both inexact flags should
2294a238c70SJohn Marino agree. Indeed, one of the inexact flags might be zero. In that
2304a238c70SJohn Marino case if we assume lngamma(z0) cannot be exact, the other flag
2314a238c70SJohn Marino should be correct. We are conservative here and request that both
2324a238c70SJohn Marino inexact flags agree. */
233*ab6d115fSJohn Marino ok = SAME_SIGN (inex1, inex2) && mpfr_cmp (l, h) == 0;
2344a238c70SJohn Marino if (ok)
2354a238c70SJohn Marino mpfr_set (y, h, rnd); /* exact */
2364a238c70SJohn Marino mpfr_clear (l);
2374a238c70SJohn Marino mpfr_clear (h);
2384a238c70SJohn Marino mpfr_clear (g);
2394a238c70SJohn Marino if (ok)
2404a238c70SJohn Marino {
241*ab6d115fSJohn Marino MPFR_ZIV_FREE (loop);
2424a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
243*ab6d115fSJohn Marino return mpfr_check_range (y, inex1, rnd);
2444a238c70SJohn Marino }
2454a238c70SJohn Marino /* since we have log|gamma(x)| = - log|x| - gamma*x + O(x^2),
2464a238c70SJohn Marino if x ~ 2^(-n), then we have a n-bit approximation, thus
2474a238c70SJohn Marino we can try again with a working precision of n bits,
2484a238c70SJohn Marino especially when n >> PREC(y).
2494a238c70SJohn Marino Otherwise we would use the reflection formula evaluating x-1,
2504a238c70SJohn Marino which would need precision n. */
2514a238c70SJohn Marino MPFR_ZIV_NEXT (loop, prec);
2524a238c70SJohn Marino }
2534a238c70SJohn Marino while (prec <= -MPFR_EXP(z0));
2544a238c70SJohn Marino MPFR_ZIV_FREE (loop);
2554a238c70SJohn Marino }
2564a238c70SJohn Marino #endif
2574a238c70SJohn Marino
2584a238c70SJohn Marino precy = MPFR_PREC(y);
2594a238c70SJohn Marino
2604a238c70SJohn Marino mpfr_init2 (s, MPFR_PREC_MIN);
2614a238c70SJohn Marino mpfr_init2 (t, MPFR_PREC_MIN);
2624a238c70SJohn Marino mpfr_init2 (u, MPFR_PREC_MIN);
2634a238c70SJohn Marino mpfr_init2 (v, MPFR_PREC_MIN);
2644a238c70SJohn Marino mpfr_init2 (z, MPFR_PREC_MIN);
2654a238c70SJohn Marino
2664a238c70SJohn Marino if (compared < 0)
2674a238c70SJohn Marino {
2684a238c70SJohn Marino mpfr_exp_t err_u;
2694a238c70SJohn Marino
2704a238c70SJohn Marino /* use reflection formula:
2714a238c70SJohn Marino gamma(x) = Pi*(x-1)/sin(Pi*(2-x))/gamma(2-x)
2724a238c70SJohn Marino thus lngamma(x) = log(Pi*(x-1)/sin(Pi*(2-x))) - lngamma(2-x) */
2734a238c70SJohn Marino
2744a238c70SJohn Marino w = precy + MPFR_INT_CEIL_LOG2 (precy);
275*ab6d115fSJohn Marino w += MPFR_INT_CEIL_LOG2 (w) + 14;
276*ab6d115fSJohn Marino MPFR_ZIV_INIT (loop, w);
2774a238c70SJohn Marino while (1)
2784a238c70SJohn Marino {
2794a238c70SJohn Marino MPFR_ASSERTD(w >= 3);
2804a238c70SJohn Marino mpfr_set_prec (s, w);
2814a238c70SJohn Marino mpfr_set_prec (t, w);
2824a238c70SJohn Marino mpfr_set_prec (u, w);
2834a238c70SJohn Marino mpfr_set_prec (v, w);
2844a238c70SJohn Marino /* In the following, we write r for a real of absolute value
2854a238c70SJohn Marino at most 2^(-w). Different instances of r may represent different
2864a238c70SJohn Marino values. */
2874a238c70SJohn Marino mpfr_ui_sub (s, 2, z0, MPFR_RNDD); /* s = (2-z0) * (1+2r) >= 1 */
2884a238c70SJohn Marino mpfr_const_pi (t, MPFR_RNDN); /* t = Pi * (1+r) */
2894a238c70SJohn Marino mpfr_lngamma (u, s, MPFR_RNDN); /* lngamma(2-x) */
2904a238c70SJohn Marino /* Let s = (2-z0) + h. By construction, -(2-z0)*2^(1-w) <= h <= 0.
2914a238c70SJohn Marino We have lngamma(s) = lngamma(2-z0) + h*Psi(z), z in [2-z0+h,2-z0].
2924a238c70SJohn Marino Since 2-z0+h = s >= 1 and |Psi(x)| <= max(1,log(x)) for x >= 1,
2934a238c70SJohn Marino the error on u is bounded by
2944a238c70SJohn Marino ulp(u)/2 + (2-z0)*max(1,log(2-z0))*2^(1-w)
2954a238c70SJohn Marino = (1/2 + (2-z0)*max(1,log(2-z0))*2^(1-E(u))) ulp(u) */
2964a238c70SJohn Marino d = (double) MPFR_GET_EXP(s) * 0.694; /* upper bound for log(2-z0) */
2974a238c70SJohn Marino err_u = MPFR_GET_EXP(s) + __gmpfr_ceil_log2 (d) + 1 - MPFR_GET_EXP(u);
2984a238c70SJohn Marino err_u = (err_u >= 0) ? err_u + 1 : 0;
2994a238c70SJohn Marino /* now the error on u is bounded by 2^err_u ulps */
3004a238c70SJohn Marino
3014a238c70SJohn Marino mpfr_mul (s, s, t, MPFR_RNDN); /* Pi*(2-x) * (1+r)^4 */
3024a238c70SJohn Marino err_s = MPFR_GET_EXP(s); /* 2-x <= 2^err_s */
3034a238c70SJohn Marino mpfr_sin (s, s, MPFR_RNDN); /* sin(Pi*(2-x)) */
3044a238c70SJohn Marino /* the error on s is bounded by 1/2*ulp(s) + [(1+2^(-w))^4-1]*(2-x)
3054a238c70SJohn Marino <= 1/2*ulp(s) + 5*2^(-w)*(2-x) for w >= 3
3064a238c70SJohn Marino <= (1/2 + 5 * 2^(-E(s)) * (2-x)) ulp(s) */
3074a238c70SJohn Marino err_s += 3 - MPFR_GET_EXP(s);
3084a238c70SJohn Marino err_s = (err_s >= 0) ? err_s + 1 : 0;
3094a238c70SJohn Marino /* the error on s is bounded by 2^err_s ulp(s), thus by
3104a238c70SJohn Marino 2^(err_s+1)*2^(-w)*|s| since ulp(s) <= 2^(1-w)*|s|.
3114a238c70SJohn Marino Now n*2^(-w) can always be written |(1+r)^n-1| for some
3124a238c70SJohn Marino |r|<=2^(-w), thus taking n=2^(err_s+1) we see that
3134a238c70SJohn Marino |S - s| <= |(1+r)^(2^(err_s+1))-1| * |s|, where S is the
3144a238c70SJohn Marino true value.
3154a238c70SJohn Marino In fact if ulp(s) <= ulp(S) the same inequality holds for
3164a238c70SJohn Marino |S| instead of |s| in the right hand side, i.e., we can
3174a238c70SJohn Marino write s = (1+r)^(2^(err_s+1)) * S.
3184a238c70SJohn Marino But if ulp(S) < ulp(s), we need to add one ``bit'' to the error,
3194a238c70SJohn Marino to get s = (1+r)^(2^(err_s+2)) * S. This is true since with
3204a238c70SJohn Marino E = n*2^(-w) we have |s - S| <= E * |s|, thus
3214a238c70SJohn Marino |s - S| <= E/(1-E) * |S|.
3224a238c70SJohn Marino Now E/(1-E) is bounded by 2E as long as E<=1/2,
3234a238c70SJohn Marino and 2E can be written (1+r)^(2n)-1 as above.
3244a238c70SJohn Marino */
3254a238c70SJohn Marino err_s += 2; /* exponent of relative error */
3264a238c70SJohn Marino
3274a238c70SJohn Marino mpfr_sub_ui (v, z0, 1, MPFR_RNDN); /* v = (x-1) * (1+r) */
3284a238c70SJohn Marino mpfr_mul (v, v, t, MPFR_RNDN); /* v = Pi*(x-1) * (1+r)^3 */
3294a238c70SJohn Marino mpfr_div (v, v, s, MPFR_RNDN); /* Pi*(x-1)/sin(Pi*(2-x)) */
3304a238c70SJohn Marino mpfr_abs (v, v, MPFR_RNDN);
3314a238c70SJohn Marino /* (1+r)^(3+2^err_s+1) */
3324a238c70SJohn Marino err_s = (err_s <= 1) ? 3 : err_s + 1;
3334a238c70SJohn Marino /* now (1+r)^M with M <= 2^err_s */
3344a238c70SJohn Marino mpfr_log (v, v, MPFR_RNDN);
3354a238c70SJohn Marino /* log(v*(1+e)) = log(v)+log(1+e) where |e| <= 2^(err_s-w).
3364a238c70SJohn Marino Since |log(1+e)| <= 2*e for |e| <= 1/4, the error on v is
3374a238c70SJohn Marino bounded by ulp(v)/2 + 2^(err_s+1-w). */
3384a238c70SJohn Marino if (err_s + 2 > w)
3394a238c70SJohn Marino {
3404a238c70SJohn Marino w += err_s + 2;
3414a238c70SJohn Marino }
3424a238c70SJohn Marino else
3434a238c70SJohn Marino {
3444a238c70SJohn Marino err_s += 1 - MPFR_GET_EXP(v);
3454a238c70SJohn Marino err_s = (err_s >= 0) ? err_s + 1 : 0;
3464a238c70SJohn Marino /* the error on v is bounded by 2^err_s ulps */
3474a238c70SJohn Marino err_u += MPFR_GET_EXP(u); /* absolute error on u */
3484a238c70SJohn Marino err_s += MPFR_GET_EXP(v); /* absolute error on v */
3494a238c70SJohn Marino mpfr_sub (s, v, u, MPFR_RNDN);
3504a238c70SJohn Marino /* the total error on s is bounded by ulp(s)/2 + 2^(err_u-w)
3514a238c70SJohn Marino + 2^(err_s-w) <= ulp(s)/2 + 2^(max(err_u,err_s)+1-w) */
3524a238c70SJohn Marino err_s = (err_s >= err_u) ? err_s : err_u;
3534a238c70SJohn Marino err_s += 1 - MPFR_GET_EXP(s); /* error is 2^err_s ulp(s) */
3544a238c70SJohn Marino err_s = (err_s >= 0) ? err_s + 1 : 0;
3554a238c70SJohn Marino if (mpfr_can_round (s, w - err_s, MPFR_RNDN, MPFR_RNDZ, precy
3564a238c70SJohn Marino + (rnd == MPFR_RNDN)))
3574a238c70SJohn Marino goto end;
3584a238c70SJohn Marino }
359*ab6d115fSJohn Marino MPFR_ZIV_NEXT (loop, w);
3604a238c70SJohn Marino }
361*ab6d115fSJohn Marino MPFR_ZIV_FREE (loop);
3624a238c70SJohn Marino }
3634a238c70SJohn Marino
3644a238c70SJohn Marino /* now z0 > 1 */
3654a238c70SJohn Marino
3664a238c70SJohn Marino MPFR_ASSERTD (compared > 0);
3674a238c70SJohn Marino
3684a238c70SJohn Marino /* since k is O(w), the value of log(z0*...*(z0+k-1)) is about w*log(w),
3694a238c70SJohn Marino so there is a cancellation of ~log(w) in the argument reconstruction */
3704a238c70SJohn Marino w = precy + MPFR_INT_CEIL_LOG2 (precy);
3714a238c70SJohn Marino w += MPFR_INT_CEIL_LOG2 (w) + 13;
372*ab6d115fSJohn Marino MPFR_ZIV_INIT (loop, w);
373*ab6d115fSJohn Marino while (1)
374*ab6d115fSJohn Marino {
3754a238c70SJohn Marino MPFR_ASSERTD (w >= 3);
3764a238c70SJohn Marino
3774a238c70SJohn Marino /* argument reduction: we compute gamma(z0 + k), where the series
3784a238c70SJohn Marino has error term B_{2n}/(z0+k)^(2n) ~ (n/(Pi*e*(z0+k)))^(2n)
3794a238c70SJohn Marino and we need k steps of argument reconstruction. Assuming k is large
3804a238c70SJohn Marino with respect to z0, and k = n, we get 1/(Pi*e)^(2n) ~ 2^(-w), i.e.,
3814a238c70SJohn Marino k ~ w*log(2)/2/log(Pi*e) ~ 0.1616 * w.
3824a238c70SJohn Marino However, since the series is more expensive to compute, the optimal
3834a238c70SJohn Marino value seems to be k ~ 4.5 * w experimentally. */
3844a238c70SJohn Marino mpfr_set_prec (s, 53);
3854a238c70SJohn Marino mpfr_gamma_alpha (s, w);
3864a238c70SJohn Marino mpfr_set_ui_2exp (s, 9, -1, MPFR_RNDU);
3874a238c70SJohn Marino mpfr_mul_ui (s, s, w, MPFR_RNDU);
3884a238c70SJohn Marino if (mpfr_cmp (z0, s) < 0)
3894a238c70SJohn Marino {
3904a238c70SJohn Marino mpfr_sub (s, s, z0, MPFR_RNDU);
3914a238c70SJohn Marino k = mpfr_get_ui (s, MPFR_RNDU);
3924a238c70SJohn Marino if (k < 3)
3934a238c70SJohn Marino k = 3;
3944a238c70SJohn Marino }
3954a238c70SJohn Marino else
3964a238c70SJohn Marino k = 3;
3974a238c70SJohn Marino
3984a238c70SJohn Marino mpfr_set_prec (s, w);
3994a238c70SJohn Marino mpfr_set_prec (t, w);
4004a238c70SJohn Marino mpfr_set_prec (u, w);
4014a238c70SJohn Marino mpfr_set_prec (v, w);
4024a238c70SJohn Marino mpfr_set_prec (z, w);
4034a238c70SJohn Marino
4044a238c70SJohn Marino mpfr_add_ui (z, z0, k, MPFR_RNDN);
4054a238c70SJohn Marino /* z = (z0+k)*(1+t1) with |t1| <= 2^(-w) */
4064a238c70SJohn Marino
4074a238c70SJohn Marino /* z >= 4 ensures the relative error on log(z) is small,
4084a238c70SJohn Marino and also (z-1/2)*log(z)-z >= 0 */
4094a238c70SJohn Marino MPFR_ASSERTD (mpfr_cmp_ui (z, 4) >= 0);
4104a238c70SJohn Marino
4114a238c70SJohn Marino mpfr_log (s, z, MPFR_RNDN); /* log(z) */
4124a238c70SJohn Marino /* we have s = log((z0+k)*(1+t1))*(1+t2) with |t1|, |t2| <= 2^(-w).
4134a238c70SJohn Marino Since w >= 2 and z0+k >= 4, we can write log((z0+k)*(1+t1))
4144a238c70SJohn Marino = log(z0+k) * (1+t3) with |t3| <= 2^(-w), thus we have
4154a238c70SJohn Marino s = log(z0+k) * (1+t4)^2 with |t4| <= 2^(-w) */
4164a238c70SJohn Marino mpfr_mul_2ui (t, z, 1, MPFR_RNDN); /* t = 2z * (1+t5) */
4174a238c70SJohn Marino mpfr_sub_ui (t, t, 1, MPFR_RNDN); /* t = 2z-1 * (1+t6)^3 */
4184a238c70SJohn Marino /* since we can write 2z*(1+t5) = (2z-1)*(1+t5') with
4194a238c70SJohn Marino t5' = 2z/(2z-1) * t5, thus |t5'| <= 8/7 * t5 */
4204a238c70SJohn Marino mpfr_mul (s, s, t, MPFR_RNDN); /* (2z-1)*log(z) * (1+t7)^6 */
4214a238c70SJohn Marino mpfr_div_2ui (s, s, 1, MPFR_RNDN); /* (z-1/2)*log(z) * (1+t7)^6 */
4224a238c70SJohn Marino mpfr_sub (s, s, z, MPFR_RNDN); /* (z-1/2)*log(z)-z */
4234a238c70SJohn Marino /* s = [(z-1/2)*log(z)-z]*(1+u)^14, s >= 1/2 */
4244a238c70SJohn Marino
4254a238c70SJohn Marino mpfr_ui_div (u, 1, z, MPFR_RNDN); /* 1/z * (1+u), u <= 1/4 since z >= 4 */
4264a238c70SJohn Marino
4274a238c70SJohn Marino /* the first term is B[2]/2/z = 1/12/z: t=1/12/z, C[2]=1 */
4284a238c70SJohn Marino mpfr_div_ui (t, u, 12, MPFR_RNDN); /* 1/(12z) * (1+u)^2, t <= 3/128 */
4294a238c70SJohn Marino mpfr_set (v, t, MPFR_RNDN); /* (1+u)^2, v < 2^(-5) */
4304a238c70SJohn Marino mpfr_add (s, s, v, MPFR_RNDN); /* (1+u)^15 */
4314a238c70SJohn Marino
4324a238c70SJohn Marino mpfr_mul (u, u, u, MPFR_RNDN); /* 1/z^2 * (1+u)^3 */
4334a238c70SJohn Marino
4344a238c70SJohn Marino if (Bm == 0)
4354a238c70SJohn Marino {
4364a238c70SJohn Marino B = mpfr_bernoulli_internal ((mpz_t *) 0, 0);
4374a238c70SJohn Marino B = mpfr_bernoulli_internal (B, 1);
4384a238c70SJohn Marino Bm = 2;
4394a238c70SJohn Marino }
4404a238c70SJohn Marino
4414a238c70SJohn Marino /* m <= maxm ensures that 2*m*(2*m+1) <= ULONG_MAX */
4424a238c70SJohn Marino maxm = 1UL << (GMP_NUMB_BITS / 2 - 1);
4434a238c70SJohn Marino
4444a238c70SJohn Marino /* s:(1+u)^15, t:(1+u)^2, t <= 3/128 */
4454a238c70SJohn Marino
4464a238c70SJohn Marino for (m = 2; MPFR_GET_EXP(v) + (mpfr_exp_t) w >= MPFR_GET_EXP(s); m++)
4474a238c70SJohn Marino {
4484a238c70SJohn Marino mpfr_mul (t, t, u, MPFR_RNDN); /* (1+u)^(10m-14) */
4494a238c70SJohn Marino if (m <= maxm)
4504a238c70SJohn Marino {
4514a238c70SJohn Marino mpfr_mul_ui (t, t, 2*(m-1)*(2*m-3), MPFR_RNDN);
4524a238c70SJohn Marino mpfr_div_ui (t, t, 2*m*(2*m-1), MPFR_RNDN);
4534a238c70SJohn Marino mpfr_div_ui (t, t, 2*m*(2*m+1), MPFR_RNDN);
4544a238c70SJohn Marino }
4554a238c70SJohn Marino else
4564a238c70SJohn Marino {
4574a238c70SJohn Marino mpfr_mul_ui (t, t, 2*(m-1), MPFR_RNDN);
4584a238c70SJohn Marino mpfr_mul_ui (t, t, 2*m-3, MPFR_RNDN);
4594a238c70SJohn Marino mpfr_div_ui (t, t, 2*m, MPFR_RNDN);
4604a238c70SJohn Marino mpfr_div_ui (t, t, 2*m-1, MPFR_RNDN);
4614a238c70SJohn Marino mpfr_div_ui (t, t, 2*m, MPFR_RNDN);
4624a238c70SJohn Marino mpfr_div_ui (t, t, 2*m+1, MPFR_RNDN);
4634a238c70SJohn Marino }
4644a238c70SJohn Marino /* (1+u)^(10m-8) */
4654a238c70SJohn Marino /* invariant: t=1/(2m)/(2m-1)/z^(2m-1)/(2m+1)! */
4664a238c70SJohn Marino if (Bm <= m)
4674a238c70SJohn Marino {
4684a238c70SJohn Marino B = mpfr_bernoulli_internal (B, m); /* B[2m]*(2m+1)!, exact */
4694a238c70SJohn Marino Bm ++;
4704a238c70SJohn Marino }
4714a238c70SJohn Marino mpfr_mul_z (v, t, B[m], MPFR_RNDN); /* (1+u)^(10m-7) */
4724a238c70SJohn Marino MPFR_ASSERTD(MPFR_GET_EXP(v) <= - (2 * m + 3));
4734a238c70SJohn Marino mpfr_add (s, s, v, MPFR_RNDN);
4744a238c70SJohn Marino }
4754a238c70SJohn Marino /* m <= 1/2*Pi*e*z ensures that |v[m]| < 1/2^(2m+3) */
4764a238c70SJohn Marino MPFR_ASSERTD ((double) m <= 4.26 * mpfr_get_d (z, MPFR_RNDZ));
4774a238c70SJohn Marino
4784a238c70SJohn Marino /* We have sum([(1+u)^(10m-7)-1]*1/2^(2m+3), m=2..infinity)
4794a238c70SJohn Marino <= 1.46*u for u <= 2^(-3).
4804a238c70SJohn Marino We have 0 < lngamma(z) - [(z - 1/2) ln(z) - z + 1/2 ln(2 Pi)] < 0.021
4814a238c70SJohn Marino for z >= 4, thus since the initial s >= 0.85, the different values of
4824a238c70SJohn Marino s differ by at most one binade, and the total rounding error on s
4834a238c70SJohn Marino in the for-loop is bounded by 2*(m-1)*ulp(final_s).
4844a238c70SJohn Marino The error coming from the v's is bounded by
4854a238c70SJohn Marino 1.46*2^(-w) <= 2*ulp(final_s).
4864a238c70SJohn Marino Thus the total error so far is bounded by [(1+u)^15-1]*s+2m*ulp(s)
4874a238c70SJohn Marino <= (2m+47)*ulp(s).
4884a238c70SJohn Marino Taking into account the truncation error (which is bounded by the last
4894a238c70SJohn Marino term v[] according to 6.1.42 in A&S), the bound is (2m+48)*ulp(s).
4904a238c70SJohn Marino */
4914a238c70SJohn Marino
4924a238c70SJohn Marino /* add 1/2*log(2*Pi) and subtract log(z0*(z0+1)*...*(z0+k-1)) */
4934a238c70SJohn Marino mpfr_const_pi (v, MPFR_RNDN); /* v = Pi*(1+u) */
4944a238c70SJohn Marino mpfr_mul_2ui (v, v, 1, MPFR_RNDN); /* v = 2*Pi * (1+u) */
4954a238c70SJohn Marino if (k)
4964a238c70SJohn Marino {
4974a238c70SJohn Marino unsigned long l;
4984a238c70SJohn Marino mpfr_set (t, z0, MPFR_RNDN); /* t = z0*(1+u) */
4994a238c70SJohn Marino for (l = 1; l < k; l++)
5004a238c70SJohn Marino {
5014a238c70SJohn Marino mpfr_add_ui (u, z0, l, MPFR_RNDN); /* u = (z0+l)*(1+u) */
5024a238c70SJohn Marino mpfr_mul (t, t, u, MPFR_RNDN); /* (1+u)^(2l+1) */
5034a238c70SJohn Marino }
5044a238c70SJohn Marino /* now t: (1+u)^(2k-1) */
5054a238c70SJohn Marino /* instead of computing log(sqrt(2*Pi)/t), we compute
5064a238c70SJohn Marino 1/2*log(2*Pi/t^2), which trades a square root for a square */
5074a238c70SJohn Marino mpfr_mul (t, t, t, MPFR_RNDN); /* (z0*...*(z0+k-1))^2, (1+u)^(4k-1) */
5084a238c70SJohn Marino mpfr_div (v, v, t, MPFR_RNDN);
5094a238c70SJohn Marino /* 2*Pi/(z0*...*(z0+k-1))^2 (1+u)^(4k+1) */
5104a238c70SJohn Marino }
5114a238c70SJohn Marino #ifdef IS_GAMMA
5124a238c70SJohn Marino err_s = MPFR_GET_EXP(s);
5134a238c70SJohn Marino mpfr_exp (s, s, MPFR_RNDN);
514*ab6d115fSJohn Marino /* If s is +Inf, we compute exp(lngamma(z0)). */
515*ab6d115fSJohn Marino if (mpfr_inf_p (s))
516*ab6d115fSJohn Marino {
517*ab6d115fSJohn Marino inexact = mpfr_explgamma (y, z0, &expo, s, t, rnd);
518*ab6d115fSJohn Marino if (inexact)
519*ab6d115fSJohn Marino goto end0;
520*ab6d115fSJohn Marino else
521*ab6d115fSJohn Marino goto ziv_next;
522*ab6d115fSJohn Marino }
5234a238c70SJohn Marino /* before the exponential, we have s = s0 + h where
5244a238c70SJohn Marino |h| <= (2m+48)*ulp(s), thus exp(s0) = exp(s) * exp(-h).
5254a238c70SJohn Marino For |h| <= 1/4, we have |exp(h)-1| <= 1.2*|h| thus
5264a238c70SJohn Marino |exp(s) - exp(s0)| <= 1.2 * exp(s) * (2m+48)* 2^(EXP(s)-w). */
5274a238c70SJohn Marino d = 1.2 * (2.0 * (double) m + 48.0);
5284a238c70SJohn Marino /* the error on s is bounded by d*2^err_s * 2^(-w) */
5294a238c70SJohn Marino mpfr_sqrt (t, v, MPFR_RNDN);
5304a238c70SJohn Marino /* let v0 be the exact value of v. We have v = v0*(1+u)^(4k+1),
5314a238c70SJohn Marino thus t = sqrt(v0)*(1+u)^(2k+3/2). */
5324a238c70SJohn Marino mpfr_mul (s, s, t, MPFR_RNDN);
5334a238c70SJohn Marino /* the error on input s is bounded by (1+u)^(d*2^err_s),
5344a238c70SJohn Marino and that on t is (1+u)^(2k+3/2), thus the
5354a238c70SJohn Marino total error is (1+u)^(d*2^err_s+2k+5/2) */
5364a238c70SJohn Marino err_s += __gmpfr_ceil_log2 (d);
5374a238c70SJohn Marino err_t = __gmpfr_ceil_log2 (2.0 * (double) k + 2.5);
5384a238c70SJohn Marino err_s = (err_s >= err_t) ? err_s + 1 : err_t + 1;
5394a238c70SJohn Marino #else
5404a238c70SJohn Marino mpfr_log (t, v, MPFR_RNDN);
5414a238c70SJohn Marino /* let v0 be the exact value of v. We have v = v0*(1+u)^(4k+1),
5424a238c70SJohn Marino thus log(v) = log(v0) + (4k+1)*log(1+u). Since |log(1+u)/u| <= 1.07
5434a238c70SJohn Marino for |u| <= 2^(-3), the absolute error on log(v) is bounded by
5444a238c70SJohn Marino 1.07*(4k+1)*u, and the rounding error by ulp(t). */
5454a238c70SJohn Marino mpfr_div_2ui (t, t, 1, MPFR_RNDN);
5464a238c70SJohn Marino /* the error on t is now bounded by ulp(t) + 0.54*(4k+1)*2^(-w).
5474a238c70SJohn Marino We have sqrt(2*Pi)/(z0*(z0+1)*...*(z0+k-1)) <= sqrt(2*Pi)/k! <= 0.5
5484a238c70SJohn Marino since k>=3, thus t <= -0.5 and ulp(t) >= 2^(-w).
5494a238c70SJohn Marino Thus the error on t is bounded by (2.16*k+1.54)*ulp(t). */
5504a238c70SJohn Marino err_t = MPFR_GET_EXP(t) + (mpfr_exp_t)
5514a238c70SJohn Marino __gmpfr_ceil_log2 (2.2 * (double) k + 1.6);
5524a238c70SJohn Marino err_s = MPFR_GET_EXP(s) + (mpfr_exp_t)
5534a238c70SJohn Marino __gmpfr_ceil_log2 (2.0 * (double) m + 48.0);
5544a238c70SJohn Marino mpfr_add (s, s, t, MPFR_RNDN); /* this is a subtraction in fact */
5554a238c70SJohn Marino /* the final error in ulp(s) is
5564a238c70SJohn Marino <= 1 + 2^(err_t-EXP(s)) + 2^(err_s-EXP(s))
5574a238c70SJohn Marino <= 2^(1+max(err_t,err_s)-EXP(s)) if err_t <> err_s
5584a238c70SJohn Marino <= 2^(2+max(err_t,err_s)-EXP(s)) if err_t = err_s */
5594a238c70SJohn Marino err_s = (err_t == err_s) ? 1 + err_s : ((err_t > err_s) ? err_t : err_s);
5604a238c70SJohn Marino err_s += 1 - MPFR_GET_EXP(s);
5614a238c70SJohn Marino #endif
562*ab6d115fSJohn Marino if (MPFR_LIKELY (MPFR_CAN_ROUND (s, w - err_s, precy, rnd)))
563*ab6d115fSJohn Marino break;
564*ab6d115fSJohn Marino #ifdef IS_GAMMA
565*ab6d115fSJohn Marino ziv_next:
566*ab6d115fSJohn Marino #endif
567*ab6d115fSJohn Marino MPFR_ZIV_NEXT (loop, w);
5684a238c70SJohn Marino }
5694a238c70SJohn Marino
570*ab6d115fSJohn Marino #ifdef IS_GAMMA
571*ab6d115fSJohn Marino end0:
572*ab6d115fSJohn Marino #endif
5734a238c70SJohn Marino oldBm = Bm;
5744a238c70SJohn Marino while (Bm--)
5754a238c70SJohn Marino mpz_clear (B[Bm]);
5764a238c70SJohn Marino (*__gmp_free_func) (B, oldBm * sizeof (mpz_t));
5774a238c70SJohn Marino
5784a238c70SJohn Marino end:
579*ab6d115fSJohn Marino if (inexact == 0)
5804a238c70SJohn Marino inexact = mpfr_set (y, s, rnd);
581*ab6d115fSJohn Marino MPFR_ZIV_FREE (loop);
5824a238c70SJohn Marino
5834a238c70SJohn Marino mpfr_clear (s);
5844a238c70SJohn Marino mpfr_clear (t);
5854a238c70SJohn Marino mpfr_clear (u);
5864a238c70SJohn Marino mpfr_clear (v);
5874a238c70SJohn Marino mpfr_clear (z);
5884a238c70SJohn Marino
5894a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
5904a238c70SJohn Marino return mpfr_check_range (y, inexact, rnd);
5914a238c70SJohn Marino }
5924a238c70SJohn Marino
5934a238c70SJohn Marino #ifndef IS_GAMMA
5944a238c70SJohn Marino
5954a238c70SJohn Marino int
mpfr_lngamma(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd)5964a238c70SJohn Marino mpfr_lngamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
5974a238c70SJohn Marino {
5984a238c70SJohn Marino int inex;
5994a238c70SJohn Marino
6004a238c70SJohn Marino MPFR_LOG_FUNC
6014a238c70SJohn Marino (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
6024a238c70SJohn Marino ("y[%Pu]=%.*Rg inexact=%d",
6034a238c70SJohn Marino mpfr_get_prec (y), mpfr_log_prec, y, inex));
6044a238c70SJohn Marino
6054a238c70SJohn Marino /* special cases */
6064a238c70SJohn Marino if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
6074a238c70SJohn Marino {
6084a238c70SJohn Marino if (MPFR_IS_NAN (x) || MPFR_IS_NEG (x))
6094a238c70SJohn Marino {
6104a238c70SJohn Marino MPFR_SET_NAN (y);
6114a238c70SJohn Marino MPFR_RET_NAN;
6124a238c70SJohn Marino }
6134a238c70SJohn Marino else /* lngamma(+Inf) = lngamma(+0) = +Inf */
6144a238c70SJohn Marino {
6154a238c70SJohn Marino if (MPFR_IS_ZERO (x))
6164a238c70SJohn Marino mpfr_set_divby0 ();
6174a238c70SJohn Marino MPFR_SET_INF (y);
6184a238c70SJohn Marino MPFR_SET_POS (y);
6194a238c70SJohn Marino MPFR_RET (0); /* exact */
6204a238c70SJohn Marino }
6214a238c70SJohn Marino }
6224a238c70SJohn Marino
6234a238c70SJohn Marino /* if x < 0 and -2k-1 <= x <= -2k, then lngamma(x) = NaN */
6244a238c70SJohn Marino if (MPFR_IS_NEG (x) && (unit_bit (x) == 0 || mpfr_integer_p (x)))
6254a238c70SJohn Marino {
6264a238c70SJohn Marino MPFR_SET_NAN (y);
6274a238c70SJohn Marino MPFR_RET_NAN;
6284a238c70SJohn Marino }
6294a238c70SJohn Marino
6304a238c70SJohn Marino inex = mpfr_lngamma_aux (y, x, rnd);
6314a238c70SJohn Marino return inex;
6324a238c70SJohn Marino }
6334a238c70SJohn Marino
6344a238c70SJohn Marino int
mpfr_lgamma(mpfr_ptr y,int * signp,mpfr_srcptr x,mpfr_rnd_t rnd)6354a238c70SJohn Marino mpfr_lgamma (mpfr_ptr y, int *signp, mpfr_srcptr x, mpfr_rnd_t rnd)
6364a238c70SJohn Marino {
6374a238c70SJohn Marino int inex;
6384a238c70SJohn Marino
6394a238c70SJohn Marino MPFR_LOG_FUNC
6404a238c70SJohn Marino (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
6414a238c70SJohn Marino ("y[%Pu]=%.*Rg signp=%d inexact=%d",
6424a238c70SJohn Marino mpfr_get_prec (y), mpfr_log_prec, y, *signp, inex));
6434a238c70SJohn Marino
6444a238c70SJohn Marino *signp = 1; /* most common case */
6454a238c70SJohn Marino
6464a238c70SJohn Marino if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
6474a238c70SJohn Marino {
6484a238c70SJohn Marino if (MPFR_IS_NAN (x))
6494a238c70SJohn Marino {
6504a238c70SJohn Marino MPFR_SET_NAN (y);
6514a238c70SJohn Marino MPFR_RET_NAN;
6524a238c70SJohn Marino }
6534a238c70SJohn Marino else
6544a238c70SJohn Marino {
6554a238c70SJohn Marino if (MPFR_IS_ZERO (x))
6564a238c70SJohn Marino mpfr_set_divby0 ();
6574a238c70SJohn Marino *signp = MPFR_INT_SIGN (x);
6584a238c70SJohn Marino MPFR_SET_INF (y);
6594a238c70SJohn Marino MPFR_SET_POS (y);
6604a238c70SJohn Marino MPFR_RET (0);
6614a238c70SJohn Marino }
6624a238c70SJohn Marino }
6634a238c70SJohn Marino
6644a238c70SJohn Marino if (MPFR_IS_NEG (x))
6654a238c70SJohn Marino {
6664a238c70SJohn Marino if (mpfr_integer_p (x))
6674a238c70SJohn Marino {
6684a238c70SJohn Marino MPFR_SET_INF (y);
6694a238c70SJohn Marino MPFR_SET_POS (y);
6704a238c70SJohn Marino mpfr_set_divby0 ();
6714a238c70SJohn Marino MPFR_RET (0);
6724a238c70SJohn Marino }
6734a238c70SJohn Marino
6744a238c70SJohn Marino if (unit_bit (x) == 0)
6754a238c70SJohn Marino *signp = -1;
6764a238c70SJohn Marino
6774a238c70SJohn Marino /* For tiny negative x, we have gamma(x) = 1/x - euler + O(x),
6784a238c70SJohn Marino thus |gamma(x)| = -1/x + euler + O(x), and
6794a238c70SJohn Marino log |gamma(x)| = -log(-x) - euler*x + O(x^2).
6804a238c70SJohn Marino More precisely we have for -0.4 <= x < 0:
6814a238c70SJohn Marino -log(-x) <= log |gamma(x)| <= -log(-x) - x.
6824a238c70SJohn Marino Since log(x) is not representable, we may have an instance of the
6834a238c70SJohn Marino Table Maker Dilemma. The only way to ensure correct rounding is to
6844a238c70SJohn Marino compute an interval [l,h] such that l <= -log(-x) and
6854a238c70SJohn Marino -log(-x) - x <= h, and check whether l and h round to the same number
6864a238c70SJohn Marino for the target precision and rounding modes. */
6874a238c70SJohn Marino if (MPFR_EXP(x) + 1 <= - (mpfr_exp_t) MPFR_PREC(y))
6884a238c70SJohn Marino /* since PREC(y) >= 1, this ensures EXP(x) <= -2,
6894a238c70SJohn Marino thus |x| <= 0.25 < 0.4 */
6904a238c70SJohn Marino {
6914a238c70SJohn Marino mpfr_t l, h;
6924a238c70SJohn Marino int ok, inex2;
6934a238c70SJohn Marino mpfr_prec_t w = MPFR_PREC (y) + 14;
6944a238c70SJohn Marino mpfr_exp_t expl;
6954a238c70SJohn Marino
6964a238c70SJohn Marino while (1)
6974a238c70SJohn Marino {
6984a238c70SJohn Marino mpfr_init2 (l, w);
6994a238c70SJohn Marino mpfr_init2 (h, w);
7004a238c70SJohn Marino /* we want a lower bound on -log(-x), thus an upper bound
7014a238c70SJohn Marino on log(-x), thus an upper bound on -x. */
7024a238c70SJohn Marino mpfr_neg (l, x, MPFR_RNDU); /* upper bound on -x */
7034a238c70SJohn Marino mpfr_log (l, l, MPFR_RNDU); /* upper bound for log(-x) */
7044a238c70SJohn Marino mpfr_neg (l, l, MPFR_RNDD); /* lower bound for -log(-x) */
7054a238c70SJohn Marino mpfr_neg (h, x, MPFR_RNDD); /* lower bound on -x */
7064a238c70SJohn Marino mpfr_log (h, h, MPFR_RNDD); /* lower bound on log(-x) */
7074a238c70SJohn Marino mpfr_neg (h, h, MPFR_RNDU); /* upper bound for -log(-x) */
7084a238c70SJohn Marino mpfr_sub (h, h, x, MPFR_RNDU); /* upper bound for -log(-x) - x */
7094a238c70SJohn Marino inex = mpfr_prec_round (l, MPFR_PREC (y), rnd);
7104a238c70SJohn Marino inex2 = mpfr_prec_round (h, MPFR_PREC (y), rnd);
7114a238c70SJohn Marino /* Caution: we not only need l = h, but both inexact flags
7124a238c70SJohn Marino should agree. Indeed, one of the inexact flags might be
7134a238c70SJohn Marino zero. In that case if we assume ln|gamma(x)| cannot be
7144a238c70SJohn Marino exact, the other flag should be correct. We are conservative
7154a238c70SJohn Marino here and request that both inexact flags agree. */
7164a238c70SJohn Marino ok = SAME_SIGN (inex, inex2) && mpfr_equal_p (l, h);
7174a238c70SJohn Marino if (ok)
7184a238c70SJohn Marino mpfr_set (y, h, rnd); /* exact */
7194a238c70SJohn Marino else
7204a238c70SJohn Marino expl = MPFR_EXP (l);
7214a238c70SJohn Marino mpfr_clear (l);
7224a238c70SJohn Marino mpfr_clear (h);
7234a238c70SJohn Marino if (ok)
7244a238c70SJohn Marino return inex;
7254a238c70SJohn Marino /* if ulp(log(-x)) <= |x| there is no reason to loop,
7264a238c70SJohn Marino since the width of [l, h] will be at least |x| */
7274a238c70SJohn Marino if (expl < MPFR_EXP(x) + (mpfr_exp_t) w)
7284a238c70SJohn Marino break;
7294a238c70SJohn Marino w += MPFR_INT_CEIL_LOG2(w) + 3;
7304a238c70SJohn Marino }
7314a238c70SJohn Marino }
7324a238c70SJohn Marino }
7334a238c70SJohn Marino
7344a238c70SJohn Marino inex = mpfr_lngamma_aux (y, x, rnd);
7354a238c70SJohn Marino return inex;
7364a238c70SJohn Marino }
7374a238c70SJohn Marino
7384a238c70SJohn Marino #endif
739