1 /* mpfr_zeta -- compute the Riemann Zeta function
2
3 Copyright 2003, 2004, 2005, 2006, 2007, 2008, 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 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 /*
27 Parameters:
28 s - the input floating-point number
29 n, p - parameters from the algorithm
30 tc - an array of p floating-point numbers tc[1]..tc[p]
31 Output:
32 b is the result, i.e.
33 sum(tc[i]*product((s+2j)*(s+2j-1)/n^2,j=1..i-1), i=1..p)*s*n^(-s-1)
34 */
35 static void
mpfr_zeta_part_b(mpfr_t b,mpfr_srcptr s,int n,int p,mpfr_t * tc)36 mpfr_zeta_part_b (mpfr_t b, mpfr_srcptr s, int n, int p, mpfr_t *tc)
37 {
38 mpfr_t s1, d, u;
39 unsigned long n2;
40 int l, t;
41 MPFR_GROUP_DECL (group);
42
43 if (p == 0)
44 {
45 MPFR_SET_ZERO (b);
46 MPFR_SET_POS (b);
47 return;
48 }
49
50 n2 = n * n;
51 MPFR_GROUP_INIT_3 (group, MPFR_PREC (b), s1, d, u);
52
53 /* t equals 2p-2, 2p-3, ... ; s1 equals s+t */
54 t = 2 * p - 2;
55 mpfr_set (d, tc[p], MPFR_RNDN);
56 for (l = 1; l < p; l++)
57 {
58 mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l) */
59 mpfr_mul (d, d, s1, MPFR_RNDN);
60 t = t - 1;
61 mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l-1) */
62 mpfr_mul (d, d, s1, MPFR_RNDN);
63 t = t - 1;
64 mpfr_div_ui (d, d, n2, MPFR_RNDN);
65 mpfr_add (d, d, tc[p-l], MPFR_RNDN);
66 /* since s is positive and the tc[i] have alternate signs,
67 the following is unlikely */
68 if (MPFR_UNLIKELY (mpfr_cmpabs (d, tc[p-l]) > 0))
69 mpfr_set (d, tc[p-l], MPFR_RNDN);
70 }
71 mpfr_mul (d, d, s, MPFR_RNDN);
72 mpfr_add (s1, s, __gmpfr_one, MPFR_RNDN);
73 mpfr_neg (s1, s1, MPFR_RNDN);
74 mpfr_ui_pow (u, n, s1, MPFR_RNDN);
75 mpfr_mul (b, d, u, MPFR_RNDN);
76
77 MPFR_GROUP_CLEAR (group);
78 }
79
80 /* Input: p - an integer
81 Output: fills tc[1..p], tc[i] = bernoulli(2i)/(2i)!
82 tc[1]=1/12, tc[2]=-1/720, tc[3]=1/30240, ...
83 */
84 static void
mpfr_zeta_c(int p,mpfr_t * tc)85 mpfr_zeta_c (int p, mpfr_t *tc)
86 {
87 mpfr_t d;
88 int k, l;
89
90 if (p > 0)
91 {
92 mpfr_init2 (d, MPFR_PREC (tc[1]));
93 mpfr_div_ui (tc[1], __gmpfr_one, 12, MPFR_RNDN);
94 for (k = 2; k <= p; k++)
95 {
96 mpfr_set_ui (d, k-1, MPFR_RNDN);
97 mpfr_div_ui (d, d, 12*k+6, MPFR_RNDN);
98 for (l=2; l < k; l++)
99 {
100 mpfr_div_ui (d, d, 4*(2*k-2*l+3)*(2*k-2*l+2), MPFR_RNDN);
101 mpfr_add (d, d, tc[l], MPFR_RNDN);
102 }
103 mpfr_div_ui (tc[k], d, 24, MPFR_RNDN);
104 MPFR_CHANGE_SIGN (tc[k]);
105 }
106 mpfr_clear (d);
107 }
108 }
109
110 /* Input: s - a floating-point number
111 n - an integer
112 Output: sum - a floating-point number approximating sum(1/i^s, i=1..n-1) */
113 static void
mpfr_zeta_part_a(mpfr_t sum,mpfr_srcptr s,int n)114 mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n)
115 {
116 mpfr_t u, s1;
117 int i;
118 MPFR_GROUP_DECL (group);
119
120 MPFR_GROUP_INIT_2 (group, MPFR_PREC (sum), u, s1);
121
122 mpfr_neg (s1, s, MPFR_RNDN);
123 mpfr_ui_pow (u, n, s1, MPFR_RNDN);
124 mpfr_div_2ui (u, u, 1, MPFR_RNDN);
125 mpfr_set (sum, u, MPFR_RNDN);
126 for (i=n-1; i>1; i--)
127 {
128 mpfr_ui_pow (u, i, s1, MPFR_RNDN);
129 mpfr_add (sum, sum, u, MPFR_RNDN);
130 }
131 mpfr_add (sum, sum, __gmpfr_one, MPFR_RNDN);
132
133 MPFR_GROUP_CLEAR (group);
134 }
135
136 /* Input: s - a floating-point number >= 1/2.
137 rnd_mode - a rounding mode.
138 Assumes s is neither NaN nor Infinite.
139 Output: z - Zeta(s) rounded to the precision of z with direction rnd_mode
140 */
141 static int
mpfr_zeta_pos(mpfr_t z,mpfr_srcptr s,mpfr_rnd_t rnd_mode)142 mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode)
143 {
144 mpfr_t b, c, z_pre, f, s1;
145 double beta, sd, dnep;
146 mpfr_t *tc1;
147 mpfr_prec_t precz, precs, d, dint;
148 int p, n, l, add;
149 int inex;
150 MPFR_GROUP_DECL (group);
151 MPFR_ZIV_DECL (loop);
152
153 MPFR_ASSERTD (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0);
154
155 precz = MPFR_PREC (z);
156 precs = MPFR_PREC (s);
157
158 /* Zeta(x) = 1+1/2^x+1/3^x+1/4^x+1/5^x+O(1/6^x)
159 so with 2^(EXP(x)-1) <= x < 2^EXP(x)
160 So for x > 2^3, k^x > k^8, so 2/k^x < 2/k^8
161 Zeta(x) = 1 + 1/2^x*(1+(2/3)^x+(2/4)^x+...)
162 = 1 + 1/2^x*(1+sum((2/k)^x,k=3..infinity))
163 <= 1 + 1/2^x*(1+sum((2/k)^8,k=3..infinity))
164 And sum((2/k)^8,k=3..infinity) = -257+128*Pi^8/4725 ~= 0.0438035
165 So Zeta(x) <= 1 + 1/2^x*2 for x >= 8
166 The error is < 2^(-x+1) <= 2^(-2^(EXP(x)-1)+1) */
167 if (MPFR_GET_EXP (s) > 3)
168 {
169 mpfr_exp_t err;
170 err = MPFR_GET_EXP (s) - 1;
171 if (err > (mpfr_exp_t) (sizeof (mpfr_exp_t)*CHAR_BIT-2))
172 err = MPFR_EMAX_MAX;
173 else
174 err = ((mpfr_exp_t)1) << err;
175 err = 1 - (-err+1); /* GET_EXP(one) - (-err+1) = err :) */
176 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, err, 0, 1,
177 rnd_mode, {});
178 }
179
180 d = precz + MPFR_INT_CEIL_LOG2(precz) + 10;
181
182 /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */
183 dint = (mpfr_uexp_t) MPFR_GET_EXP (s);
184 mpfr_init2 (s1, MAX (precs, dint));
185 inex = mpfr_sub (s1, s, __gmpfr_one, MPFR_RNDN);
186 MPFR_ASSERTD (inex == 0);
187
188 /* case s=1 should have already been handled */
189 MPFR_ASSERTD (!MPFR_IS_ZERO (s1));
190
191 MPFR_GROUP_INIT_4 (group, MPFR_PREC_MIN, b, c, z_pre, f);
192
193 MPFR_ZIV_INIT (loop, d);
194 for (;;)
195 {
196 /* Principal loop: we compute, in z_pre,
197 an approximation of Zeta(s), that we send to can_round */
198 if (MPFR_GET_EXP (s1) <= -(mpfr_exp_t) ((mpfr_prec_t) (d-3)/2))
199 /* Branch 1: when s-1 is very small, one
200 uses the approximation Zeta(s)=1/(s-1)+gamma,
201 where gamma is Euler's constant */
202 {
203 dint = MAX (d + 3, precs);
204 MPFR_TRACE (printf ("branch 1\ninternal precision=%lu\n",
205 (unsigned long) dint));
206 MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
207 mpfr_div (z_pre, __gmpfr_one, s1, MPFR_RNDN);
208 mpfr_const_euler (f, MPFR_RNDN);
209 mpfr_add (z_pre, z_pre, f, MPFR_RNDN);
210 }
211 else /* Branch 2 */
212 {
213 size_t size;
214
215 MPFR_TRACE (printf ("branch 2\n"));
216 /* Computation of parameters n, p and working precision */
217 dnep = (double) d * LOG2;
218 sd = mpfr_get_d (s, MPFR_RNDN);
219 /* beta = dnep + 0.61 + sd * log (6.2832 / sd);
220 but a larger value is ok */
221 #define LOG6dot2832 1.83787940484160805532
222 beta = dnep + 0.61 + sd * (LOG6dot2832 - LOG2 *
223 __gmpfr_floor_log2 (sd));
224 if (beta <= 0.0)
225 {
226 p = 0;
227 /* n = 1 + (int) (exp ((dnep - LOG2) / sd)); */
228 n = 1 + (int) __gmpfr_ceil_exp2 ((d - 1.0) / sd);
229 }
230 else
231 {
232 p = 1 + (int) beta / 2;
233 n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832);
234 }
235 MPFR_TRACE (printf ("\nn=%d\np=%d\n",n,p));
236 /* add = 4 + floor(1.5 * log(d) / log (2)).
237 We should have add >= 10, which is always fulfilled since
238 d = precz + 11 >= 12, thus ceil(log2(d)) >= 4 */
239 add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2;
240 MPFR_ASSERTD(add >= 10);
241 dint = d + add;
242 if (dint < precs)
243 dint = precs;
244
245 MPFR_TRACE (printf ("internal precision=%lu\n",
246 (unsigned long) dint));
247
248 size = (p + 1) * sizeof(mpfr_t);
249 tc1 = (mpfr_t*) (*__gmp_allocate_func) (size);
250 for (l=1; l<=p; l++)
251 mpfr_init2 (tc1[l], dint);
252 MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
253
254 MPFR_TRACE (printf ("precision of z = %lu\n",
255 (unsigned long) precz));
256
257 /* Computation of the coefficients c_k */
258 mpfr_zeta_c (p, tc1);
259 /* Computation of the 3 parts of the fonction Zeta. */
260 mpfr_zeta_part_a (z_pre, s, n);
261 mpfr_zeta_part_b (b, s, n, p, tc1);
262 /* s1 = s-1 is already computed above */
263 mpfr_div (c, __gmpfr_one, s1, MPFR_RNDN);
264 mpfr_ui_pow (f, n, s1, MPFR_RNDN);
265 mpfr_div (c, c, f, MPFR_RNDN);
266 MPFR_TRACE (MPFR_DUMP (c));
267 mpfr_add (z_pre, z_pre, c, MPFR_RNDN);
268 mpfr_add (z_pre, z_pre, b, MPFR_RNDN);
269 for (l=1; l<=p; l++)
270 mpfr_clear (tc1[l]);
271 (*__gmp_free_func) (tc1, size);
272 /* End branch 2 */
273 }
274
275 MPFR_TRACE (MPFR_DUMP (z_pre));
276 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, d-3, precz, rnd_mode)))
277 break;
278 MPFR_ZIV_NEXT (loop, d);
279 }
280 MPFR_ZIV_FREE (loop);
281
282 inex = mpfr_set (z, z_pre, rnd_mode);
283
284 MPFR_GROUP_CLEAR (group);
285 mpfr_clear (s1);
286
287 return inex;
288 }
289
290 int
mpfr_zeta(mpfr_t z,mpfr_srcptr s,mpfr_rnd_t rnd_mode)291 mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode)
292 {
293 mpfr_t z_pre, s1, y, p;
294 double sd, eps, m1, c;
295 long add;
296 mpfr_prec_t precz, prec1, precs, precs1;
297 int inex;
298 MPFR_GROUP_DECL (group);
299 MPFR_ZIV_DECL (loop);
300 MPFR_SAVE_EXPO_DECL (expo);
301
302 MPFR_LOG_FUNC (
303 ("s[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (s), mpfr_log_prec, s, rnd_mode),
304 ("z[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (z), mpfr_log_prec, z, inex));
305
306 /* Zero, Nan or Inf ? */
307 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s)))
308 {
309 if (MPFR_IS_NAN (s))
310 {
311 MPFR_SET_NAN (z);
312 MPFR_RET_NAN;
313 }
314 else if (MPFR_IS_INF (s))
315 {
316 if (MPFR_IS_POS (s))
317 return mpfr_set_ui (z, 1, MPFR_RNDN); /* Zeta(+Inf) = 1 */
318 MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */
319 MPFR_RET_NAN;
320 }
321 else /* s iz zero */
322 {
323 MPFR_ASSERTD (MPFR_IS_ZERO (s));
324 return mpfr_set_si_2exp (z, -1, -1, rnd_mode);
325 }
326 }
327
328 /* s is neither Nan, nor Inf, nor Zero */
329
330 /* check tiny s: we have zeta(s) = -1/2 - 1/2 log(2 Pi) s + ... around s=0,
331 and for |s| <= 0.074, we have |zeta(s) + 1/2| <= |s|.
332 Thus if |s| <= 1/4*ulp(1/2), we can deduce the correct rounding
333 (the 1/4 covers the case where |zeta(s)| < 1/2 and rounding to nearest).
334 A sufficient condition is that EXP(s) + 1 < -PREC(z). */
335 if (MPFR_GET_EXP (s) + 1 < - (mpfr_exp_t) MPFR_PREC(z))
336 {
337 int signs = MPFR_SIGN(s);
338
339 MPFR_SAVE_EXPO_MARK (expo);
340 mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */
341 if (rnd_mode == MPFR_RNDA)
342 rnd_mode = MPFR_RNDD; /* the result is around -1/2, thus negative */
343 if ((rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDZ) && signs < 0)
344 {
345 mpfr_nextabove (z); /* z = -1/2 + epsilon */
346 inex = 1;
347 }
348 else if (rnd_mode == MPFR_RNDD && signs > 0)
349 {
350 mpfr_nextbelow (z); /* z = -1/2 - epsilon */
351 inex = -1;
352 }
353 else
354 {
355 if (rnd_mode == MPFR_RNDU) /* s > 0: z = -1/2 */
356 inex = 1;
357 else if (rnd_mode == MPFR_RNDD)
358 inex = -1; /* s < 0: z = -1/2 */
359 else /* (MPFR_RNDZ and s > 0) or MPFR_RNDN: z = -1/2 */
360 inex = (signs > 0) ? 1 : -1;
361 }
362 MPFR_SAVE_EXPO_FREE (expo);
363 return mpfr_check_range (z, inex, rnd_mode);
364 }
365
366 /* Check for case s= -2n */
367 if (MPFR_IS_NEG (s))
368 {
369 mpfr_t tmp;
370 tmp[0] = *s;
371 MPFR_EXP (tmp) = MPFR_GET_EXP (s) - 1;
372 if (mpfr_integer_p (tmp))
373 {
374 MPFR_SET_ZERO (z);
375 MPFR_SET_POS (z);
376 MPFR_RET (0);
377 }
378 }
379
380 /* Check for case s= 1 before changing the exponent range */
381 if (mpfr_cmp (s, __gmpfr_one) ==0)
382 {
383 MPFR_SET_INF (z);
384 MPFR_SET_POS (z);
385 mpfr_set_divby0 ();
386 MPFR_RET (0);
387 }
388
389 MPFR_SAVE_EXPO_MARK (expo);
390
391 /* Compute Zeta */
392 if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */
393 inex = mpfr_zeta_pos (z, s, rnd_mode);
394 else /* use reflection formula
395 zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */
396 {
397 int overflow = 0;
398
399 precz = MPFR_PREC (z);
400 precs = MPFR_PREC (s);
401
402 /* Precision precs1 needed to represent 1 - s, and s + 2,
403 without any truncation */
404 precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s));
405 sd = mpfr_get_d (s, MPFR_RNDN) - 1.0;
406 if (sd < 0.0)
407 sd = -sd; /* now sd = abs(s-1.0) */
408 /* Precision prec1 is the precision on elementary computations;
409 it ensures a final precision prec1 - add for zeta(s) */
410 /* eps = pow (2.0, - (double) precz - 14.0); */
411 eps = __gmpfr_ceil_exp2 (- (double) precz - 14.0);
412 m1 = 1.0 + MAX(1.0 / eps, 2.0 * sd) * (1.0 + eps);
413 c = (1.0 + eps) * (1.0 + eps * MAX(8.0, m1));
414 /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */
415 add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1));
416 prec1 = precz + add;
417 prec1 = MAX (prec1, precs1) + 10;
418
419 MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p);
420 MPFR_ZIV_INIT (loop, prec1);
421 for (;;)
422 {
423 mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN);/* s1 = 1-s */
424 mpfr_zeta_pos (z_pre, s1, MPFR_RNDN); /* zeta(1-s) */
425 mpfr_gamma (y, s1, MPFR_RNDN); /* gamma(1-s) */
426 if (MPFR_IS_INF (y)) /* Zeta(s) < 0 for -4k-2 < s < -4k,
427 Zeta(s) > 0 for -4k < s < -4k+2 */
428 {
429 mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */
430 mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */
431 overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1;
432 break;
433 }
434 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); /* gamma(1-s)*zeta(1-s) */
435 mpfr_const_pi (p, MPFR_RNDD);
436 mpfr_mul (y, s, p, MPFR_RNDN);
437 mpfr_div_2ui (y, y, 1, MPFR_RNDN); /* s*Pi/2 */
438 mpfr_sin (y, y, MPFR_RNDN); /* sin(Pi*s/2) */
439 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
440 mpfr_mul_2ui (y, p, 1, MPFR_RNDN); /* 2*Pi */
441 mpfr_neg (s1, s1, MPFR_RNDN); /* s-1 */
442 mpfr_pow (y, y, s1, MPFR_RNDN); /* (2*Pi)^(s-1) */
443 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
444 mpfr_mul_2ui (z_pre, z_pre, 1, MPFR_RNDN);
445
446 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz,
447 rnd_mode)))
448 break;
449
450 MPFR_ZIV_NEXT (loop, prec1);
451 MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p);
452 }
453 MPFR_ZIV_FREE (loop);
454 if (overflow != 0)
455 {
456 inex = mpfr_overflow (z, rnd_mode, overflow);
457 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
458 }
459 else
460 inex = mpfr_set (z, z_pre, rnd_mode);
461 MPFR_GROUP_CLEAR (group);
462 }
463
464 MPFR_SAVE_EXPO_FREE (expo);
465 return mpfr_check_range (z, inex, rnd_mode);
466 }
467