1/* mpfr_beta -- beta function 2 3Copyright 2017-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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) */ 28int 29mpfr_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} 346