1 /* mpfr_beta -- beta function 2 3 Copyright 2017-2020 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 nonpositive 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