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