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
mpfr_digamma_approx(mpfr_ptr s,mpfr_srcptr x)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
mpfr_digamma_reflection(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)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
mpfr_digamma_positive(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)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
mpfr_digamma(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)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