1 /* mpfr_gamma -- gamma function
2
3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 #define IS_GAMMA
27 #include "lngamma.c"
28 #undef IS_GAMMA
29
30 /* return a sufficient precision such that 2-x is exact, assuming x < 0 */
31 static mpfr_prec_t
mpfr_gamma_2_minus_x_exact(mpfr_srcptr x)32 mpfr_gamma_2_minus_x_exact (mpfr_srcptr x)
33 {
34 /* Since x < 0, 2-x = 2+y with y := -x.
35 If y < 2, a precision w >= PREC(y) + EXP(2)-EXP(y) = PREC(y) + 2 - EXP(y)
36 is enough, since no overlap occurs in 2+y, so no carry happens.
37 If y >= 2, either ULP(y) <= 2, and we need w >= PREC(y)+1 since a
38 carry can occur, or ULP(y) > 2, and we need w >= EXP(y)-1:
39 (a) if EXP(y) <= 1, w = PREC(y) + 2 - EXP(y)
40 (b) if EXP(y) > 1 and EXP(y)-PREC(y) <= 1, w = PREC(y) + 1
41 (c) if EXP(y) > 1 and EXP(y)-PREC(y) > 1, w = EXP(y) - 1 */
42 return (MPFR_GET_EXP(x) <= 1) ? MPFR_PREC(x) + 2 - MPFR_GET_EXP(x)
43 : ((MPFR_GET_EXP(x) <= MPFR_PREC(x) + 1) ? MPFR_PREC(x) + 1
44 : MPFR_GET_EXP(x) - 1);
45 }
46
47 /* return a sufficient precision such that 1-x is exact, assuming x < 1 */
48 static mpfr_prec_t
mpfr_gamma_1_minus_x_exact(mpfr_srcptr x)49 mpfr_gamma_1_minus_x_exact (mpfr_srcptr x)
50 {
51 if (MPFR_IS_POS(x))
52 return MPFR_PREC(x) - MPFR_GET_EXP(x);
53 else if (MPFR_GET_EXP(x) <= 0)
54 return MPFR_PREC(x) + 1 - MPFR_GET_EXP(x);
55 else if (MPFR_PREC(x) >= MPFR_GET_EXP(x))
56 return MPFR_PREC(x) + 1;
57 else
58 return MPFR_GET_EXP(x);
59 }
60
61 /* returns a lower bound of the number of significant bits of n!
62 (not counting the low zero bits).
63 We know n! >= (n/e)^n*sqrt(2*Pi*n) for n >= 1, and the number of zero bits
64 is floor(n/2) + floor(n/4) + floor(n/8) + ...
65 This approximation is exact for n <= 500000, except for n = 219536, 235928,
66 298981, 355854, 464848, 493725, 498992 where it returns a value 1 too small.
67 */
68 static unsigned long
bits_fac(unsigned long n)69 bits_fac (unsigned long n)
70 {
71 mpfr_t x, y;
72 unsigned long r, k;
73 mpfr_init2 (x, 38);
74 mpfr_init2 (y, 38);
75 mpfr_set_ui (x, n, MPFR_RNDZ);
76 mpfr_set_str_binary (y, "10.101101111110000101010001011000101001"); /* upper bound of e */
77 mpfr_div (x, x, y, MPFR_RNDZ);
78 mpfr_pow_ui (x, x, n, MPFR_RNDZ);
79 mpfr_const_pi (y, MPFR_RNDZ);
80 mpfr_mul_ui (y, y, 2 * n, MPFR_RNDZ);
81 mpfr_sqrt (y, y, MPFR_RNDZ);
82 mpfr_mul (x, x, y, MPFR_RNDZ);
83 mpfr_log2 (x, x, MPFR_RNDZ);
84 r = mpfr_get_ui (x, MPFR_RNDU);
85 for (k = 2; k <= n; k *= 2)
86 r -= n / k;
87 mpfr_clear (x);
88 mpfr_clear (y);
89 return r;
90 }
91
92 /* We use the reflection formula
93 Gamma(1+t) Gamma(1-t) = - Pi t / sin(Pi (1 + t))
94 in order to treat the case x <= 1,
95 i.e. with x = 1-t, then Gamma(x) = -Pi*(1-x)/sin(Pi*(2-x))/GAMMA(2-x)
96 */
97 int
mpfr_gamma(mpfr_ptr gamma,mpfr_srcptr x,mpfr_rnd_t rnd_mode)98 mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
99 {
100 mpfr_t xp, GammaTrial, tmp, tmp2;
101 mpz_t fact;
102 mpfr_prec_t realprec;
103 int compared, is_integer;
104 int inex = 0; /* 0 means: result gamma not set yet */
105 MPFR_GROUP_DECL (group);
106 MPFR_SAVE_EXPO_DECL (expo);
107 MPFR_ZIV_DECL (loop);
108
109 MPFR_LOG_FUNC
110 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
111 ("gamma[%Pu]=%.*Rg inexact=%d",
112 mpfr_get_prec (gamma), mpfr_log_prec, gamma, inex));
113
114 /* Trivial cases */
115 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
116 {
117 if (MPFR_IS_NAN (x))
118 {
119 MPFR_SET_NAN (gamma);
120 MPFR_RET_NAN;
121 }
122 else if (MPFR_IS_INF (x))
123 {
124 if (MPFR_IS_NEG (x))
125 {
126 MPFR_SET_NAN (gamma);
127 MPFR_RET_NAN;
128 }
129 else
130 {
131 MPFR_SET_INF (gamma);
132 MPFR_SET_POS (gamma);
133 MPFR_RET (0); /* exact */
134 }
135 }
136 else /* x is zero */
137 {
138 MPFR_ASSERTD(MPFR_IS_ZERO(x));
139 MPFR_SET_INF(gamma);
140 MPFR_SET_SAME_SIGN(gamma, x);
141 mpfr_set_divby0 ();
142 MPFR_RET (0); /* exact */
143 }
144 }
145
146 /* Check for tiny arguments, where gamma(x) ~ 1/x - euler + ....
147 We know from "Bound on Runs of Zeros and Ones for Algebraic Functions",
148 Proceedings of Arith15, T. Lang and J.-M. Muller, 2001, that the maximal
149 number of consecutive zeroes or ones after the round bit is n-1 for an
150 input of n bits. But we need a more precise lower bound. Assume x has
151 n bits, and 1/x is near a floating-point number y of n+1 bits. We can
152 write x = X*2^e, y = Y/2^f with X, Y integers of n and n+1 bits.
153 Thus X*Y^2^(e-f) is near from 1, i.e., X*Y is near from 2^(f-e).
154 Two cases can happen:
155 (i) either X*Y is exactly 2^(f-e), but this can happen only if X and Y
156 are themselves powers of two, i.e., x is a power of two;
157 (ii) or X*Y is at distance at least one from 2^(f-e), thus
158 |xy-1| >= 2^(e-f), or |y-1/x| >= 2^(e-f)/x = 2^(-f)/X >= 2^(-f-n).
159 Since ufp(y) = 2^(n-f) [ufp = unit in first place], this means
160 that the distance |y-1/x| >= 2^(-2n) ufp(y).
161 Now assuming |gamma(x)-1/x| <= 1, which is true for x <= 1,
162 if 2^(-2n) ufp(y) >= 2, the error is at most 2^(-2n-1) ufp(y),
163 and round(1/x) with precision >= 2n+2 gives the correct result.
164 If x < 2^E, then y > 2^(-E), thus ufp(y) > 2^(-E-1).
165 A sufficient condition is thus EXP(x) + 2 <= -2 MAX(PREC(x),PREC(Y)).
166 */
167 if (MPFR_GET_EXP (x) + 2
168 <= -2 * (mpfr_exp_t) MAX(MPFR_PREC(x), MPFR_PREC(gamma)))
169 {
170 int sign = MPFR_SIGN (x); /* retrieve sign before possible override */
171 int special;
172 MPFR_BLOCK_DECL (flags);
173
174 MPFR_SAVE_EXPO_MARK (expo);
175
176 /* for overflow cases, see below; this needs to be done
177 before x possibly gets overridden. */
178 special =
179 MPFR_GET_EXP (x) == 1 - MPFR_EMAX_MAX &&
180 MPFR_IS_POS_SIGN (sign) &&
181 MPFR_IS_LIKE_RNDD (rnd_mode, sign) &&
182 mpfr_powerof2_raw (x);
183
184 MPFR_BLOCK (flags, inex = mpfr_ui_div (gamma, 1, x, rnd_mode));
185 if (inex == 0) /* x is a power of two */
186 {
187 /* return RND(1/x - euler) = RND(+/- 2^k - eps) with eps > 0 */
188 if (rnd_mode == MPFR_RNDN || MPFR_IS_LIKE_RNDU (rnd_mode, sign))
189 inex = 1;
190 else
191 {
192 mpfr_nextbelow (gamma);
193 inex = -1;
194 }
195 }
196 else if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
197 {
198 /* Overflow in the division 1/x. This is a real overflow, except
199 in RNDZ or RNDD when 1/x = 2^emax, i.e. x = 2^(-emax): due to
200 the "- euler", the rounded value in unbounded exponent range
201 is 0.111...11 * 2^emax (not an overflow). */
202 if (!special)
203 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, flags);
204 }
205 MPFR_SAVE_EXPO_FREE (expo);
206 /* Note: an overflow is possible with an infinite result;
207 in this case, the overflow flag will automatically be
208 restored by mpfr_check_range. */
209 return mpfr_check_range (gamma, inex, rnd_mode);
210 }
211
212 is_integer = mpfr_integer_p (x);
213 /* gamma(x) for x a negative integer gives NaN */
214 if (is_integer && MPFR_IS_NEG(x))
215 {
216 MPFR_SET_NAN (gamma);
217 MPFR_RET_NAN;
218 }
219
220 compared = mpfr_cmp_ui (x, 1);
221 if (compared == 0)
222 return mpfr_set_ui (gamma, 1, rnd_mode);
223
224 /* if x is an integer that fits into an unsigned long, use mpfr_fac_ui
225 if argument is not too large.
226 If precision is p, fac_ui costs O(u*p), whereas gamma costs O(p*M(p)),
227 so for u <= M(p), fac_ui should be faster.
228 We approximate here M(p) by p*log(p)^2, which is not a bad guess.
229 Warning: since the generic code does not handle exact cases,
230 we want all cases where gamma(x) is exact to be treated here.
231 */
232 if (is_integer && mpfr_fits_ulong_p (x, MPFR_RNDN))
233 {
234 unsigned long int u;
235 mpfr_prec_t p = MPFR_PREC(gamma);
236 u = mpfr_get_ui (x, MPFR_RNDN);
237 if (u < 44787929UL && bits_fac (u - 1) <= p + (rnd_mode == MPFR_RNDN))
238 /* bits_fac: lower bound on the number of bits of m,
239 where gamma(x) = (u-1)! = m*2^e with m odd. */
240 return mpfr_fac_ui (gamma, u - 1, rnd_mode);
241 /* if bits_fac(...) > p (resp. p+1 for rounding to nearest),
242 then gamma(x) cannot be exact in precision p (resp. p+1).
243 FIXME: remove the test u < 44787929UL after changing bits_fac
244 to return a mpz_t or mpfr_t. */
245 }
246
247 MPFR_SAVE_EXPO_MARK (expo);
248
249 /* check for overflow: according to (6.1.37) in Abramowitz & Stegun,
250 gamma(x) >= exp(-x) * x^(x-1/2) * sqrt(2*Pi)
251 >= 2 * (x/e)^x / x for x >= 1 */
252 if (compared > 0)
253 {
254 mpfr_t yp;
255 mpfr_exp_t expxp;
256 MPFR_BLOCK_DECL (flags);
257
258 /* 1/e rounded down to 53 bits */
259 #define EXPM1_STR "0.010111100010110101011000110110001011001110111100111"
260 mpfr_init2 (xp, 53);
261 mpfr_init2 (yp, 53);
262 mpfr_set_str_binary (xp, EXPM1_STR);
263 mpfr_mul (xp, x, xp, MPFR_RNDZ);
264 mpfr_sub_ui (yp, x, 2, MPFR_RNDZ);
265 mpfr_pow (xp, xp, yp, MPFR_RNDZ); /* (x/e)^(x-2) */
266 mpfr_set_str_binary (yp, EXPM1_STR);
267 mpfr_mul (xp, xp, yp, MPFR_RNDZ); /* x^(x-2) / e^(x-1) */
268 mpfr_mul (xp, xp, yp, MPFR_RNDZ); /* x^(x-2) / e^x */
269 mpfr_mul (xp, xp, x, MPFR_RNDZ); /* lower bound on x^(x-1) / e^x */
270 MPFR_BLOCK (flags, mpfr_mul_2ui (xp, xp, 1, MPFR_RNDZ));
271 expxp = MPFR_GET_EXP (xp);
272 mpfr_clear (xp);
273 mpfr_clear (yp);
274 MPFR_SAVE_EXPO_FREE (expo);
275 return MPFR_OVERFLOW (flags) || expxp > __gmpfr_emax ?
276 mpfr_overflow (gamma, rnd_mode, 1) :
277 mpfr_gamma_aux (gamma, x, rnd_mode);
278 }
279
280 /* now compared < 0 */
281
282 /* check for underflow: for x < 1,
283 gamma(x) = Pi*(x-1)/sin(Pi*(2-x))/gamma(2-x).
284 Since gamma(2-x) >= 2 * ((2-x)/e)^(2-x) / (2-x), we have
285 |gamma(x)| <= Pi*(1-x)*(2-x)/2/((2-x)/e)^(2-x) / |sin(Pi*(2-x))|
286 <= 12 * ((2-x)/e)^x / |sin(Pi*(2-x))|.
287 To avoid an underflow in ((2-x)/e)^x, we compute the logarithm.
288 */
289 if (MPFR_IS_NEG(x))
290 {
291 int underflow = 0, sgn, ck;
292 mpfr_prec_t w;
293
294 mpfr_init2 (xp, 53);
295 mpfr_init2 (tmp, 53);
296 mpfr_init2 (tmp2, 53);
297 /* we want an upper bound for x * [log(2-x)-1].
298 since x < 0, we need a lower bound on log(2-x) */
299 mpfr_ui_sub (xp, 2, x, MPFR_RNDD);
300 mpfr_log (xp, xp, MPFR_RNDD);
301 mpfr_sub_ui (xp, xp, 1, MPFR_RNDD);
302 mpfr_mul (xp, xp, x, MPFR_RNDU);
303
304 /* we need an upper bound on 1/|sin(Pi*(2-x))|,
305 thus a lower bound on |sin(Pi*(2-x))|.
306 If 2-x is exact, then the error of Pi*(2-x) is (1+u)^2 with u = 2^(-p)
307 thus the error on sin(Pi*(2-x)) is less than 1/2ulp + 3Pi(2-x)u,
308 assuming u <= 1, thus <= u + 3Pi(2-x)u */
309
310 w = mpfr_gamma_2_minus_x_exact (x); /* 2-x is exact for prec >= w */
311 w += 17; /* to get tmp2 small enough */
312 mpfr_set_prec (tmp, w);
313 mpfr_set_prec (tmp2, w);
314 ck = mpfr_ui_sub (tmp, 2, x, MPFR_RNDN);
315 MPFR_ASSERTD (ck == 0); (void) ck; /* use ck to avoid a warning */
316 mpfr_const_pi (tmp2, MPFR_RNDN);
317 mpfr_mul (tmp2, tmp2, tmp, MPFR_RNDN); /* Pi*(2-x) */
318 mpfr_sin (tmp, tmp2, MPFR_RNDN); /* sin(Pi*(2-x)) */
319 sgn = mpfr_sgn (tmp);
320 mpfr_abs (tmp, tmp, MPFR_RNDN);
321 mpfr_mul_ui (tmp2, tmp2, 3, MPFR_RNDU); /* 3Pi(2-x) */
322 mpfr_add_ui (tmp2, tmp2, 1, MPFR_RNDU); /* 3Pi(2-x)+1 */
323 mpfr_div_2ui (tmp2, tmp2, mpfr_get_prec (tmp), MPFR_RNDU);
324 /* if tmp2<|tmp|, we get a lower bound */
325 if (mpfr_cmp (tmp2, tmp) < 0)
326 {
327 mpfr_sub (tmp, tmp, tmp2, MPFR_RNDZ); /* low bnd on |sin(Pi*(2-x))| */
328 mpfr_ui_div (tmp, 12, tmp, MPFR_RNDU); /* upper bound */
329 mpfr_log2 (tmp, tmp, MPFR_RNDU);
330 mpfr_add (xp, tmp, xp, MPFR_RNDU);
331 /* The assert below checks that expo.saved_emin - 2 always
332 fits in a long. FIXME if we want to allow mpfr_exp_t to
333 be a long long, for instance. */
334 MPFR_ASSERTN (MPFR_EMIN_MIN - 2 >= LONG_MIN);
335 underflow = mpfr_cmp_si (xp, expo.saved_emin - 2) <= 0;
336 }
337
338 mpfr_clear (xp);
339 mpfr_clear (tmp);
340 mpfr_clear (tmp2);
341 if (underflow) /* the sign is the opposite of that of sin(Pi*(2-x)) */
342 {
343 MPFR_SAVE_EXPO_FREE (expo);
344 return mpfr_underflow (gamma, (rnd_mode == MPFR_RNDN) ? MPFR_RNDZ : rnd_mode, -sgn);
345 }
346 }
347
348 realprec = MPFR_PREC (gamma);
349 /* we want both 1-x and 2-x to be exact */
350 {
351 mpfr_prec_t w;
352 w = mpfr_gamma_1_minus_x_exact (x);
353 if (realprec < w)
354 realprec = w;
355 w = mpfr_gamma_2_minus_x_exact (x);
356 if (realprec < w)
357 realprec = w;
358 }
359 realprec = realprec + MPFR_INT_CEIL_LOG2 (realprec) + 20;
360 MPFR_ASSERTD(realprec >= 5);
361
362 MPFR_GROUP_INIT_4 (group, realprec + MPFR_INT_CEIL_LOG2 (realprec) + 20,
363 xp, tmp, tmp2, GammaTrial);
364 mpz_init (fact);
365 MPFR_ZIV_INIT (loop, realprec);
366 for (;;)
367 {
368 mpfr_exp_t err_g;
369 int ck;
370 MPFR_GROUP_REPREC_4 (group, realprec, xp, tmp, tmp2, GammaTrial);
371
372 /* reflection formula: gamma(x) = Pi*(x-1)/sin(Pi*(2-x))/gamma(2-x) */
373
374 ck = mpfr_ui_sub (xp, 2, x, MPFR_RNDN); /* 2-x, exact */
375 MPFR_ASSERTD(ck == 0); (void) ck; /* use ck to avoid a warning */
376 mpfr_gamma (tmp, xp, MPFR_RNDN); /* gamma(2-x), error (1+u) */
377 mpfr_const_pi (tmp2, MPFR_RNDN); /* Pi, error (1+u) */
378 mpfr_mul (GammaTrial, tmp2, xp, MPFR_RNDN); /* Pi*(2-x), error (1+u)^2 */
379 err_g = MPFR_GET_EXP(GammaTrial);
380 mpfr_sin (GammaTrial, GammaTrial, MPFR_RNDN); /* sin(Pi*(2-x)) */
381 /* If tmp is +Inf, we compute exp(lngamma(x)). */
382 if (mpfr_inf_p (tmp))
383 {
384 inex = mpfr_explgamma (gamma, x, &expo, tmp, tmp2, rnd_mode);
385 if (inex)
386 goto end;
387 else
388 goto ziv_next;
389 }
390 err_g = err_g + 1 - MPFR_GET_EXP(GammaTrial);
391 /* let g0 the true value of Pi*(2-x), g the computed value.
392 We have g = g0 + h with |h| <= |(1+u^2)-1|*g.
393 Thus sin(g) = sin(g0) + h' with |h'| <= |(1+u^2)-1|*g.
394 The relative error is thus bounded by |(1+u^2)-1|*g/sin(g)
395 <= |(1+u^2)-1|*2^err_g. <= 2.25*u*2^err_g for |u|<=1/4.
396 With the rounding error, this gives (0.5 + 2.25*2^err_g)*u. */
397 ck = mpfr_sub_ui (xp, x, 1, MPFR_RNDN); /* x-1, exact */
398 MPFR_ASSERTD(ck == 0); (void) ck; /* use ck to avoid a warning */
399 mpfr_mul (xp, tmp2, xp, MPFR_RNDN); /* Pi*(x-1), error (1+u)^2 */
400 mpfr_mul (GammaTrial, GammaTrial, tmp, MPFR_RNDN);
401 /* [1 + (0.5 + 2.25*2^err_g)*u]*(1+u)^2 = 1 + (2.5 + 2.25*2^err_g)*u
402 + (0.5 + 2.25*2^err_g)*u*(2u+u^2) + u^2.
403 For err_g <= realprec-2, we have (0.5 + 2.25*2^err_g)*u <=
404 0.5*u + 2.25/4 <= 0.6875 and u^2 <= u/4, thus
405 (0.5 + 2.25*2^err_g)*u*(2u+u^2) + u^2 <= 0.6875*(2u+u/4) + u/4
406 <= 1.8*u, thus the rel. error is bounded by (4.5 + 2.25*2^err_g)*u. */
407 mpfr_div (GammaTrial, xp, GammaTrial, MPFR_RNDN);
408 /* the error is of the form (1+u)^3/[1 + (4.5 + 2.25*2^err_g)*u].
409 For realprec >= 5 and err_g <= realprec-2, [(4.5 + 2.25*2^err_g)*u]^2
410 <= 0.71, and for |y|<=0.71, 1/(1-y) can be written 1+a*y with a<=4.
411 (1+u)^3 * (1+4*(4.5 + 2.25*2^err_g)*u)
412 = 1 + (21 + 9*2^err_g)*u + (57+27*2^err_g)*u^2 + (55+27*2^err_g)*u^3
413 + (18+9*2^err_g)*u^4
414 <= 1 + (21 + 9*2^err_g)*u + (57+27*2^err_g)*u^2 + (56+28*2^err_g)*u^3
415 <= 1 + (21 + 9*2^err_g)*u + (59+28*2^err_g)*u^2
416 <= 1 + (23 + 10*2^err_g)*u.
417 The final error is thus bounded by (23 + 10*2^err_g) ulps,
418 which is <= 2^6 for err_g<=2, and <= 2^(err_g+4) for err_g >= 2. */
419 err_g = (err_g <= 2) ? 6 : err_g + 4;
420
421 if (MPFR_LIKELY (MPFR_CAN_ROUND (GammaTrial, realprec - err_g,
422 MPFR_PREC(gamma), rnd_mode)))
423 break;
424
425 ziv_next:
426 MPFR_ZIV_NEXT (loop, realprec);
427 }
428
429 end:
430 MPFR_ZIV_FREE (loop);
431
432 if (inex == 0)
433 inex = mpfr_set (gamma, GammaTrial, rnd_mode);
434 MPFR_GROUP_CLEAR (group);
435 mpz_clear (fact);
436
437 MPFR_SAVE_EXPO_FREE (expo);
438 return mpfr_check_range (gamma, inex, rnd_mode);
439 }
440