1 /* mpfr_beta -- beta function
2
3 Copyright 2017-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 https://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 /* for MPFR_INT_CEIL_LOG2 */
24 #include "mpfr-impl.h"
25
26 /* use formula (6.2.2) from Abramowitz & Stegun:
27 beta(z,w) = gamma(z)*gamma(w)/gamma(z+w) */
28 int
29 mpfr_beta (mpfr_ptr r, mpfr_srcptr z, mpfr_srcptr w, mpfr_rnd_t rnd_mode)
30 {
31 mpfr_exp_t emin, emax;
32 mpfr_uexp_t pmin;
33 mpfr_prec_t prec;
34 mpfr_t z_plus_w, tmp, tmp2;
35 int inex, w_integer;
36 MPFR_GROUP_DECL (group);
37 MPFR_ZIV_DECL (loop);
38 MPFR_SAVE_EXPO_DECL (expo);
39
40 if (mpfr_less_p (z, w))
41 return mpfr_beta (r, w, z, rnd_mode);
42
43 /* Now, either z and w are unordered (at least one is a NaN), or z >= w. */
44
45 if (MPFR_ARE_SINGULAR (z, w))
46 {
47 /* if z or w is NaN, return NaN */
48 if (MPFR_IS_NAN (z) || MPFR_IS_NAN (w))
49 {
50 MPFR_SET_NAN (r);
51 MPFR_RET_NAN;
52 }
53 else if (MPFR_IS_INF (z) || MPFR_IS_INF (w))
54 {
55 /* Since we have z >= w:
56 if z = +Inf and w > 0, then r = +0 (including w = +Inf);
57 if z = +Inf and w = 0, then r = NaN
58 [beta(z,1/log(z)) tends to +Inf whereas
59 beta(z,1/log(log(z))) tends to +0]
60 if z = +Inf and w < 0:
61 if w is an integer or -Inf: r = NaN
62 if -2k-1 < w < -2k: r = -Inf
63 if -2k-2 < w < -2k-1: r = +Inf
64 if w = -Inf and z is finite and not an integer:
65 beta(z,t) for t going to -Inf oscillates between positive and
66 negative values, with poles around integer values of t, thus
67 beta(z,w) gives NaN;
68 if w = -Inf and z is an integer:
69 beta(z,w) gives +0 for z even > 0, -0 for z odd > 0,
70 NaN for z <= 0;
71 if z = -Inf (then w = -Inf too): r = NaN */
72 if (MPFR_IS_INF (z) && MPFR_IS_POS(z)) /* z = +Inf */
73 {
74 if (mpfr_cmp_ui (w, 0) > 0)
75 {
76 MPFR_SET_ZERO(r);
77 MPFR_SET_POS(r);
78 MPFR_RET(0);
79 }
80 else if (MPFR_IS_ZERO(w) || MPFR_IS_INF(w) || mpfr_integer_p (w))
81 {
82 MPFR_SET_NAN(r);
83 MPFR_RET_NAN;
84 }
85 else
86 {
87 long q;
88 mpfr_t t;
89
90 MPFR_SAVE_EXPO_MARK (expo);
91 mpfr_init2 (t, MPFR_PREC_MIN);
92 mpfr_set_ui (t, 1, MPFR_RNDN);
93 mpfr_fmodquo (t, &q, w, t, MPFR_RNDD);
94 mpfr_clear (t);
95 MPFR_SAVE_EXPO_FREE (expo);
96 /* q contains the low bits of trunc(w) where trunc() rounds
97 toward zero, thus if q is odd, then -2k-2 < w < -2k-1 */
98 MPFR_SET_INF(r);
99 if ((unsigned long) q & 1)
100 MPFR_SET_NEG(r);
101 else
102 MPFR_SET_POS(r);
103 MPFR_RET(0);
104 }
105 }
106 else if (MPFR_IS_INF(w)) /* w = -Inf */
107 {
108 if (mpfr_cmp_ui (z, 0) <= 0 || !mpfr_integer_p (z))
109 {
110 MPFR_SET_NAN(r);
111 MPFR_RET_NAN;
112 }
113 else
114 {
115 MPFR_SET_ZERO(r);
116 if (mpfr_odd_p (z))
117 MPFR_SET_NEG(r);
118 else
119 MPFR_SET_POS(r);
120 MPFR_RET(0);
121 }
122 }
123 }
124 else /* z or w is 0 */
125 {
126 /* If x is not a non-positive integer, Gamma(x) is regular, so that
127 when y -> 0 with either y >= 0 or y <= 0,
128 Beta(x,y) ~ Gamma(x) * Gamma(y) / Gamma(x) = Gamma(y)
129 Gamma(y) tends to an infinity of the same sign as y.
130 Thus Beta(x,y) should be an infinity of the same sign as y.
131 */
132 if (mpfr_cmp_ui (z, 0) != 0) /* then w is +0 or -0 and z > 0 */
133 {
134 /* beta(z,+0) = +Inf, beta(z,-0) = -Inf (see above) */
135 MPFR_SET_INF(r);
136 MPFR_SET_SAME_SIGN(r,w);
137 MPFR_SET_DIVBY0 ();
138 MPFR_RET(0);
139 }
140 else if (mpfr_cmp_ui (w, 0) != 0) /* then z is +0 or -0 and w < 0 */
141 {
142 if (mpfr_integer_p (w))
143 {
144 /* For small u > 0, Beta(2u,w+u) and Beta(2u,w-u) have
145 opposite signs, so that they tend to infinities of
146 opposite signs when u -> 0. Thus the result is NaN. */
147 MPFR_SET_NAN(r);
148 MPFR_RET_NAN;
149 }
150 else
151 {
152 /* beta(+0,w) = +Inf, beta(-0,w) = -Inf (see above) */
153 MPFR_SET_INF(r);
154 MPFR_SET_SAME_SIGN(r,z);
155 MPFR_SET_DIVBY0 ();
156 MPFR_RET(0);
157 }
158 }
159 else /* w = z = 0:
160 beta(+0,+0) = +Inf
161 beta(-0,-0) = -Inf
162 beta(+0,-0) = NaN */
163 {
164 if (MPFR_SIGN(z) == MPFR_SIGN(w))
165 {
166 MPFR_SET_INF(r);
167 MPFR_SET_SAME_SIGN(r,z);
168 MPFR_SET_DIVBY0 ();
169 MPFR_RET(0);
170 }
171 else
172 {
173 MPFR_SET_NAN(r);
174 MPFR_RET_NAN;
175 }
176 }
177 }
178 }
179
180 /* special case when w is a negative integer */
181 w_integer = mpfr_integer_p (w);
182 if (w_integer && MPFR_IS_NEG(w))
183 {
184 /* if z < 0 or z+w > 0, or z is not an integer, return NaN */
185 if (MPFR_IS_NEG(z) || mpfr_cmpabs (z, w) > 0 || !mpfr_integer_p (z))
186 {
187 MPFR_SET_NAN(r);
188 MPFR_RET_NAN;
189 }
190 /* If z+w = 0, the result is 1/z. */
191 if (mpfr_cmpabs (z, w) == 0)
192 return mpfr_ui_div (r, 1, z, rnd_mode);
193 /* Now z is an integer and z+w <= 0: return (-1)^z*beta(z,1-w-z).
194 Since z and w are of opposite signs, |z+w| <= max(|z|,|w|). */
195 emax = MAX (MPFR_EXP(z), MPFR_EXP(w));
196 mpfr_init2 (z_plus_w, (mpfr_prec_t) emax);
197 inex = mpfr_add (z_plus_w, z, w, MPFR_RNDN);
198 MPFR_ASSERTN(inex == 0);
199 inex = mpfr_ui_sub (z_plus_w, 1, z_plus_w, MPFR_RNDN);
200 MPFR_ASSERTN(inex == 0);
201 if (mpfr_odd_p (z))
202 {
203 inex = -mpfr_beta (r, z, z_plus_w, MPFR_INVERT_RND (rnd_mode));
204 MPFR_CHANGE_SIGN(r);
205 }
206 else
207 inex = mpfr_beta (r, z, z_plus_w, rnd_mode);
208 mpfr_clear (z_plus_w);
209 return inex;
210 }
211
212 /* special case when z is a negative integer: here w < z and w is not an
213 integer */
214 if (mpfr_integer_p (z) && MPFR_IS_NEG(z))
215 {
216 MPFR_SET_NAN(r);
217 MPFR_RET_NAN;
218 }
219
220 MPFR_SAVE_EXPO_MARK (expo);
221
222 /* compute the smallest precision such that z + w is exact */
223 emax = MAX (MPFR_EXP(z), MPFR_EXP(w));
224 emin = MIN (MPFR_EXP(z) - MPFR_PREC(z), MPFR_EXP(w) - MPFR_PREC(w));
225 MPFR_ASSERTD (emax >= emin);
226 /* Thus the math value of emax - emin is representable in mpfr_uexp_t. */
227 pmin = (mpfr_uexp_t) emax - emin;
228 /* If z and w have same sign, their sum can have exponent emax + 1. */
229 pmin += 1;
230 if (pmin > MPFR_PREC_MAX) /* FIXME: check if result can differ from NaN. */
231 {
232 MPFR_SAVE_EXPO_FREE (expo);
233 MPFR_SET_NAN(r);
234 MPFR_RET_NAN;
235 }
236 MPFR_ASSERTN (pmin <= MPFR_PREC_MAX); /* detect integer overflow */
237 mpfr_init2 (z_plus_w, (mpfr_prec_t) pmin);
238 inex = mpfr_add (z_plus_w, z, w, MPFR_RNDN);
239 /* if z+w overflows with rounding to nearest, then w must be larger than
240 1/2*ulp(z), thus we have an underflow. */
241 if (MPFR_IS_INF(z_plus_w))
242 {
243 mpfr_clear (z_plus_w);
244 MPFR_SAVE_EXPO_FREE (expo);
245 return mpfr_underflow (r, rnd_mode, 1);
246 }
247 MPFR_ASSERTN(inex == 0);
248
249 /* If z+w is 0 or a negative integer, return +0 when w (and thus z) is not
250 an integer. Indeed, gamma(z) and gamma(w) are regular numbers, and
251 gamma(z+w) is Inf, thus 1/gamma(z+w) is zero. Unless there is a rule
252 to choose the sign of 0, we choose +0. */
253 if (mpfr_cmp_ui (z_plus_w, 0) <= 0 && !w_integer
254 && mpfr_integer_p (z_plus_w))
255 {
256 mpfr_clear (z_plus_w);
257 MPFR_SAVE_EXPO_FREE (expo);
258 MPFR_SET_ZERO(r);
259 MPFR_SET_POS(r);
260 MPFR_RET(0);
261 }
262
263 prec = MPFR_PREC(r);
264 prec += MPFR_INT_CEIL_LOG2 (prec);
265 MPFR_GROUP_INIT_2 (group, prec, tmp, tmp2);
266 MPFR_ZIV_INIT (loop, prec);
267 for (;;)
268 {
269 unsigned int inex2; /* unsigned due to bitwise operations */
270
271 MPFR_GROUP_REPREC_2 (group, prec, tmp, tmp2);
272 inex2 = mpfr_gamma (tmp, z, MPFR_RNDN);
273 /* tmp = gamma(z) * (1 + theta) with |theta| <= 2^-prec */
274 inex2 |= mpfr_gamma (tmp2, w, MPFR_RNDN);
275 /* tmp2 = gamma(w) * (1 + theta2) with |theta2| <= 2^-prec */
276 inex2 |= mpfr_mul (tmp, tmp, tmp2, MPFR_RNDN);
277 /* tmp = gamma(z)*gamma(w) * (1 + theta3)^3 with |theta3| <= 2^-prec */
278 inex2 |= mpfr_gamma (tmp2, z_plus_w, MPFR_RNDN);
279 /* tmp2 = gamma(z+w) * (1 + theta4) with |theta4| <= 2^-prec */
280 inex2 |= mpfr_div (tmp, tmp, tmp2, MPFR_RNDN);
281 /* tmp = gamma(z)*gamma(w)/gamma(z+w) * (1 + theta5)^5
282 with |theta5| <= 2^-prec. For prec >= 3, we have
283 |(1 + theta5)^5 - 1| <= 7 * 2^(-prec), thus the error is bounded
284 by 7 ulps */
285
286 if (MPFR_IS_NAN(tmp)) /* FIXME: most probably gamma(z)*gamma(w) = +-Inf,
287 and gamma(z+w) = +-Inf, can we do better? */
288 {
289 mpfr_clear (z_plus_w);
290 MPFR_ZIV_FREE (loop);
291 MPFR_GROUP_CLEAR (group);
292 MPFR_SAVE_EXPO_FREE (expo);
293 MPFR_SET_NAN(r);
294 MPFR_RET_NAN;
295 }
296
297 MPFR_ASSERTN(mpfr_regular_p (tmp));
298
299 /* if inex2 = 0, then tmp is exactly beta(z,w) */
300 if (inex2 == 0 ||
301 MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - 3, MPFR_PREC(r), rnd_mode)))
302 break;
303
304 /* beta(1,+/-2^(-k)) = +/-2^k is exact, and cannot be detected above
305 since gamma(+/-2^(-k)) is not exact */
306 if (mpfr_cmp_ui (z, 1) == 0)
307 {
308 mpfr_exp_t expw = mpfr_get_exp (w);
309 if (mpfr_cmp_ui_2exp (w, 1, expw - 1) == 0)
310 {
311 /* since z >= w, this will only match w <= 1 */
312 mpfr_set_ui_2exp (tmp, 1, 1 - expw, MPFR_RNDN);
313 break;
314 }
315 else if (mpfr_cmp_si_2exp (w, -1, expw - 1) == 0)
316 {
317 mpfr_set_si_2exp (tmp, -1, 1 - expw, MPFR_RNDN);
318 break;
319 }
320 }
321
322 /* beta(2^k,1) = 1/2^k for k > 0 (k <= 0 was already tested above) */
323 if (mpfr_cmp_ui (w, 1) == 0 &&
324 mpfr_cmp_ui_2exp (z, 1, MPFR_EXP(z) - 1) == 0)
325 {
326 mpfr_set_ui_2exp (tmp, 1, 1 - MPFR_EXP(z), MPFR_RNDN);
327 break;
328 }
329
330 /* beta(2,-0.5) = -4 */
331 if (mpfr_cmp_ui (z, 2) == 0 && mpfr_cmp_si_2exp (w, -1, -1) == 0)
332 {
333 mpfr_set_si_2exp (tmp, -1, 2, MPFR_RNDN);
334 break;
335 }
336
337 MPFR_ZIV_NEXT (loop, prec);
338 }
339 MPFR_ZIV_FREE (loop);
340 inex = mpfr_set (r, tmp, rnd_mode);
341 MPFR_GROUP_CLEAR (group);
342 mpfr_clear (z_plus_w);
343 MPFR_SAVE_EXPO_FREE (expo);
344 return mpfr_check_range (r, inex, rnd_mode);
345 }