1 /* mpfr_digamma -- digamma function of a floating-point number 2 3 Copyright 2009-2018 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 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 #include "mpfr-impl.h" 24 25 /* Put in s an approximation of digamma(x). 26 Assumes x >= 2. 27 Assumes s does not overlap with x. 28 Returns an integer e such that the error is bounded by 2^e ulps 29 of the result s. 30 */ 31 static mpfr_exp_t 32 mpfr_digamma_approx (mpfr_ptr s, mpfr_srcptr x) 33 { 34 mpfr_prec_t p = MPFR_PREC (s); 35 mpfr_t t, u, invxx; 36 mpfr_exp_t e, exps, f, expu; 37 unsigned long n; 38 39 MPFR_ASSERTN(MPFR_IS_POS(x) && (MPFR_EXP(x) >= 2)); 40 41 mpfr_init2 (t, p); 42 mpfr_init2 (u, p); 43 mpfr_init2 (invxx, p); 44 45 mpfr_log (s, x, MPFR_RNDN); /* error <= 1/2 ulp */ 46 mpfr_ui_div (t, 1, x, MPFR_RNDN); /* error <= 1/2 ulp */ 47 mpfr_div_2exp (t, t, 1, MPFR_RNDN); /* exact */ 48 mpfr_sub (s, s, t, MPFR_RNDN); 49 /* error <= 1/2 + 1/2*2^(EXP(olds)-EXP(s)) + 1/2*2^(EXP(t)-EXP(s)). 50 For x >= 2, log(x) >= 2*(1/(2x)), thus olds >= 2t, and olds - t >= olds/2, 51 thus 0 <= EXP(olds)-EXP(s) <= 1, and EXP(t)-EXP(s) <= 0, thus 52 error <= 1/2 + 1/2*2 + 1/2 <= 2 ulps. */ 53 e = 2; /* initial error */ 54 mpfr_mul (invxx, x, x, MPFR_RNDZ); /* invxx = x^2 * (1 + theta) 55 for |theta| <= 2^(-p) */ 56 mpfr_ui_div (invxx, 1, invxx, MPFR_RNDU); /* invxx = 1/x^2 * (1 + theta)^2 */ 57 58 /* in the following we note err=xxx when the ratio between the approximation 59 and the exact result can be written (1 + theta)^xxx for |theta| <= 2^(-p), 60 following Higham's method */ 61 mpfr_set_ui (t, 1, MPFR_RNDN); /* err = 0 */ 62 for (n = 1;; n++) 63 { 64 /* The main term is Bernoulli[2n]/(2n)/x^(2n) = B[n]/(2n+1)!(2n)/x^(2n) 65 = B[n]*t[n]/(2n) where t[n]/t[n-1] = 1/(2n)/(2n+1)/x^2. */ 66 mpfr_mul (t, t, invxx, MPFR_RNDU); /* err = err + 3 */ 67 mpfr_div_ui (t, t, 2 * n, MPFR_RNDU); /* err = err + 1 */ 68 mpfr_div_ui (t, t, 2 * n + 1, MPFR_RNDU); /* err = err + 1 */ 69 /* we thus have err = 5n here */ 70 mpfr_div_ui (u, t, 2 * n, MPFR_RNDU); /* err = 5n+1 */ 71 mpfr_mul_z (u, u, mpfr_bernoulli_cache(n), MPFR_RNDU);/* err = 5n+2, and the 72 absolute error is bounded 73 by 10n+4 ulp(u) [Rule 11] */ 74 /* if the terms 'u' are decreasing by a factor two at least, 75 then the error coming from those is bounded by 76 sum((10n+4)/2^n, n=1..infinity) = 24 */ 77 exps = mpfr_get_exp (s); 78 expu = mpfr_get_exp (u); 79 if (expu < exps - (mpfr_exp_t) p) 80 break; 81 mpfr_sub (s, s, u, MPFR_RNDN); /* error <= 24 + n/2 */ 82 if (mpfr_get_exp (s) < exps) 83 e <<= exps - mpfr_get_exp (s); 84 e ++; /* error in mpfr_sub */ 85 f = 10 * n + 4; 86 while (expu < exps) 87 { 88 f = (1 + f) / 2; 89 expu ++; 90 } 91 e += f; /* total rounding error coming from 'u' term */ 92 } 93 94 mpfr_clear (t); 95 mpfr_clear (u); 96 mpfr_clear (invxx); 97 98 f = 0; 99 while (e > 1) 100 { 101 f++; 102 e = (e + 1) / 2; 103 /* Invariant: 2^f * e does not decrease */ 104 } 105 return f; 106 } 107 108 /* Use the reflection formula Digamma(1-x) = Digamma(x) + Pi * cot(Pi*x), 109 i.e., Digamma(x) = Digamma(1-x) - Pi * cot(Pi*x). 110 Assume x < 1/2. */ 111 static int 112 mpfr_digamma_reflection (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 113 { 114 mpfr_prec_t p = MPFR_PREC(y) + 10, q; 115 mpfr_t t, u, v; 116 mpfr_exp_t e1, expv; 117 int inex; 118 MPFR_ZIV_DECL (loop); 119 120 MPFR_LOG_FUNC 121 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode), 122 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y, inex)); 123 124 /* we want that 1-x is exact with precision q: if 0 < x < 1/2, then 125 q = PREC(x)-EXP(x) is ok, otherwise if -1 <= x < 0, q = PREC(x)-EXP(x) 126 is ok, otherwise for x < -1, PREC(x) is ok if EXP(x) <= PREC(x), 127 otherwise we need EXP(x) */ 128 if (MPFR_EXP(x) < 0) 129 q = MPFR_PREC(x) + 1 - MPFR_EXP(x); 130 else if (MPFR_EXP(x) <= MPFR_PREC(x)) 131 q = MPFR_PREC(x) + 1; 132 else 133 q = MPFR_EXP(x); 134 mpfr_init2 (u, q); 135 MPFR_DBGRES(inex = mpfr_ui_sub (u, 1, x, MPFR_RNDN)); 136 MPFR_ASSERTN(inex == 0); 137 138 /* if x is half an integer, cot(Pi*x) = 0, thus Digamma(x) = Digamma(1-x) */ 139 mpfr_mul_2exp (u, u, 1, MPFR_RNDN); 140 inex = mpfr_integer_p (u); 141 mpfr_div_2exp (u, u, 1, MPFR_RNDN); 142 if (inex) 143 { 144 inex = mpfr_digamma (y, u, rnd_mode); 145 goto end; 146 } 147 148 mpfr_init2 (t, p); 149 mpfr_init2 (v, p); 150 151 MPFR_ZIV_INIT (loop, p); 152 for (;;) 153 { 154 mpfr_const_pi (v, MPFR_RNDN); /* v = Pi*(1+theta) for |theta|<=2^(-p) */ 155 mpfr_mul (t, v, x, MPFR_RNDN); /* (1+theta)^2 */ 156 e1 = MPFR_EXP(t) - (mpfr_exp_t) p + 1; /* bound for t: err(t) <= 2^e1 */ 157 mpfr_cot (t, t, MPFR_RNDN); 158 /* cot(t * (1+h)) = cot(t) - theta * (1 + cot(t)^2) with |theta|<=t*h */ 159 if (MPFR_EXP(t) > 0) 160 e1 = e1 + 2 * MPFR_EXP(t) + 1; 161 else 162 e1 = e1 + 1; 163 /* now theta * (1 + cot(t)^2) <= 2^e1 */ 164 e1 += (mpfr_exp_t) p - MPFR_EXP(t); /* error is now 2^e1 ulps */ 165 mpfr_mul (t, t, v, MPFR_RNDN); 166 e1 ++; 167 mpfr_digamma (v, u, MPFR_RNDN); /* error <= 1/2 ulp */ 168 expv = MPFR_EXP(v); 169 mpfr_sub (v, v, t, MPFR_RNDN); 170 if (MPFR_EXP(v) < MPFR_EXP(t)) 171 e1 += MPFR_EXP(t) - MPFR_EXP(v); /* scale error for t wrt new v */ 172 /* now take into account the 1/2 ulp error for v */ 173 if (expv - MPFR_EXP(v) - 1 > e1) 174 e1 = expv - MPFR_EXP(v) - 1; 175 else 176 e1 ++; 177 e1 ++; /* rounding error for mpfr_sub */ 178 if (MPFR_CAN_ROUND (v, p - e1, MPFR_PREC(y), rnd_mode)) 179 break; 180 MPFR_ZIV_NEXT (loop, p); 181 mpfr_set_prec (t, p); 182 mpfr_set_prec (v, p); 183 } 184 MPFR_ZIV_FREE (loop); 185 186 inex = mpfr_set (y, v, rnd_mode); 187 188 mpfr_clear (t); 189 mpfr_clear (v); 190 end: 191 mpfr_clear (u); 192 193 return inex; 194 } 195 196 /* we have x >= 1/2 here */ 197 static int 198 mpfr_digamma_positive (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 199 { 200 mpfr_prec_t p = MPFR_PREC(y) + 10, q; 201 mpfr_t t, u, x_plus_j; 202 int inex; 203 mpfr_exp_t errt, erru, expt; 204 unsigned long j = 0, min; 205 MPFR_ZIV_DECL (loop); 206 207 MPFR_LOG_FUNC 208 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode), 209 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y, inex)); 210 211 /* compute a precision q such that x+1 is exact */ 212 if (MPFR_PREC(x) < MPFR_EXP(x)) 213 q = MPFR_EXP(x); 214 else 215 q = MPFR_PREC(x) + 1; 216 217 /* for very large x, use |digamma(x) - log(x)| < 1/x < 2^(1-EXP(x)) */ 218 if (MPFR_PREC(y) + 10 < MPFR_EXP(x)) 219 { 220 /* this ensures EXP(x) >= 3, thus x >= 4, thus log(x) > 1 */ 221 mpfr_init2 (t, MPFR_PREC(y) + 10); 222 mpfr_log (t, x, MPFR_RNDZ); 223 if (MPFR_CAN_ROUND (t, MPFR_PREC(y) + 10, MPFR_PREC(y), rnd_mode)) 224 { 225 inex = mpfr_set (y, t, rnd_mode); 226 mpfr_clear (t); 227 return inex; 228 } 229 mpfr_clear (t); 230 } 231 232 mpfr_init2 (x_plus_j, q); 233 234 mpfr_init2 (t, p); 235 mpfr_init2 (u, p); 236 MPFR_ZIV_INIT (loop, p); 237 for(;;) 238 { 239 /* Lower bound for x+j in mpfr_digamma_approx call: since the smallest 240 term of the divergent series for Digamma(x) is about exp(-2*Pi*x), and 241 we want it to be less than 2^(-p), this gives x > p*log(2)/(2*Pi) 242 i.e., x >= 0.1103 p. 243 To be safe, we ensure x >= 0.25 * p. 244 */ 245 min = (p + 3) / 4; 246 if (min < 2) 247 min = 2; 248 249 mpfr_set (x_plus_j, x, MPFR_RNDN); 250 mpfr_set_ui (u, 0, MPFR_RNDN); 251 j = 0; 252 while (mpfr_cmp_ui (x_plus_j, min) < 0) 253 { 254 j ++; 255 mpfr_ui_div (t, 1, x_plus_j, MPFR_RNDN); /* err <= 1/2 ulp */ 256 mpfr_add (u, u, t, MPFR_RNDN); 257 inex = mpfr_add_ui (x_plus_j, x_plus_j, 1, MPFR_RNDZ); 258 if (inex != 0) /* we lost one bit */ 259 { 260 q ++; 261 mpfr_prec_round (x_plus_j, q, MPFR_RNDZ); 262 mpfr_nextabove (x_plus_j); 263 } 264 /* since all terms are positive, the error is bounded by j ulps */ 265 } 266 for (erru = 0; j > 1; erru++, j = (j + 1) / 2); 267 errt = mpfr_digamma_approx (t, x_plus_j); 268 expt = MPFR_EXP(t); 269 mpfr_sub (t, t, u, MPFR_RNDN); 270 if (MPFR_EXP(t) < expt) 271 errt += expt - MPFR_EXP(t); 272 if (MPFR_EXP(t) < MPFR_EXP(u)) 273 erru += MPFR_EXP(u) - MPFR_EXP(t); 274 if (errt > erru) 275 errt = errt + 1; 276 else if (errt == erru) 277 errt = errt + 2; 278 else 279 errt = erru + 1; 280 if (MPFR_CAN_ROUND (t, p - errt, MPFR_PREC(y), rnd_mode)) 281 break; 282 MPFR_ZIV_NEXT (loop, p); 283 mpfr_set_prec (t, p); 284 mpfr_set_prec (u, p); 285 } 286 MPFR_ZIV_FREE (loop); 287 inex = mpfr_set (y, t, rnd_mode); 288 mpfr_clear (t); 289 mpfr_clear (u); 290 mpfr_clear (x_plus_j); 291 return inex; 292 } 293 294 int 295 mpfr_digamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 296 { 297 int inex; 298 MPFR_SAVE_EXPO_DECL (expo); 299 300 MPFR_LOG_FUNC 301 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode), 302 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y, inex)); 303 304 if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x))) 305 { 306 if (MPFR_IS_NAN(x)) 307 { 308 MPFR_SET_NAN(y); 309 MPFR_RET_NAN; 310 } 311 else if (MPFR_IS_INF(x)) 312 { 313 if (MPFR_IS_POS(x)) /* Digamma(+Inf) = +Inf */ 314 { 315 MPFR_SET_SAME_SIGN(y, x); 316 MPFR_SET_INF(y); 317 MPFR_RET(0); 318 } 319 else /* Digamma(-Inf) = NaN */ 320 { 321 MPFR_SET_NAN(y); 322 MPFR_RET_NAN; 323 } 324 } 325 else /* Zero case */ 326 { 327 /* the following works also in case of overlap */ 328 MPFR_SET_INF(y); 329 MPFR_SET_OPPOSITE_SIGN(y, x); 330 MPFR_SET_DIVBY0 (); 331 MPFR_RET(0); 332 } 333 } 334 335 /* Digamma is undefined for negative integers */ 336 if (MPFR_IS_NEG(x) && mpfr_integer_p (x)) 337 { 338 MPFR_SET_NAN(y); 339 MPFR_RET_NAN; 340 } 341 342 /* now x is a normal number */ 343 344 MPFR_SAVE_EXPO_MARK (expo); 345 /* for x very small, we have Digamma(x) = -1/x - gamma + O(x), more precisely 346 -1 < Digamma(x) + 1/x < 0 for -0.2 < x < 0.2, thus: 347 (i) either x is a power of two, then 1/x is exactly representable, and 348 as long as 1/2*ulp(1/x) > 1, we can conclude; 349 (ii) otherwise assume x has <= n bits, and y has <= n+1 bits, then 350 |y + 1/x| >= 2^(-2n) ufp(y), where ufp means unit in first place. 351 Since |Digamma(x) + 1/x| <= 1, if 2^(-2n) ufp(y) >= 2, then 352 |y - Digamma(x)| >= 2^(-2n-1)ufp(y), and rounding -1/x gives the correct result. 353 If x < 2^E, then y > 2^(-E), thus ufp(y) > 2^(-E-1). 354 A sufficient condition is thus EXP(x) <= -2 MAX(PREC(x),PREC(Y)). */ 355 if (MPFR_EXP(x) < -2) 356 { 357 if (MPFR_EXP(x) <= -2 * (mpfr_exp_t) MAX(MPFR_PREC(x), MPFR_PREC(y))) 358 { 359 int signx = MPFR_SIGN(x); 360 inex = mpfr_si_div (y, -1, x, rnd_mode); 361 if (inex == 0) /* x is a power of two */ 362 { /* result always -1/x, except when rounding down */ 363 if (rnd_mode == MPFR_RNDA) 364 rnd_mode = (signx > 0) ? MPFR_RNDD : MPFR_RNDU; 365 if (rnd_mode == MPFR_RNDZ) 366 rnd_mode = (signx > 0) ? MPFR_RNDU : MPFR_RNDD; 367 if (rnd_mode == MPFR_RNDU) 368 inex = 1; 369 else if (rnd_mode == MPFR_RNDD) 370 { 371 mpfr_nextbelow (y); 372 inex = -1; 373 } 374 else /* nearest */ 375 inex = 1; 376 } 377 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); 378 goto end; 379 } 380 } 381 382 if (MPFR_IS_NEG(x)) 383 inex = mpfr_digamma_reflection (y, x, rnd_mode); 384 /* if x < 1/2 we use the reflection formula */ 385 else if (MPFR_EXP(x) < 0) 386 inex = mpfr_digamma_reflection (y, x, rnd_mode); 387 else 388 inex = mpfr_digamma_positive (y, x, rnd_mode); 389 390 end: 391 MPFR_SAVE_EXPO_FREE (expo); 392 return mpfr_check_range (y, inex, rnd_mode); 393 } 394