1 /* mpfr_zeta -- compute the Riemann Zeta function
3 Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by Jean-Luc Re'my and the Spaces project, INRIA Lorraine.
6 This file is part of the GNU MPFR Library.
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 2.1 of the License, or (at your
11 option) any later version.
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.
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.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
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]
33 sum(tc[i]*product((s+2j)*(s+2j-1)/n^2,j=1..i-1), i=1..p)*s*n^(-s-1)
36 mpfr_zeta_part_b (mpfr_t b, mpfr_srcptr s, int n, int p, mpfr_t *tc)
41 MPFR_GROUP_DECL (group);
51 MPFR_GROUP_INIT_3 (group, MPFR_PREC (b), s1, d, u);
53 /* t equals 2p-2, 2p-3, ... ; s1 equals s+t */
55 mpfr_set (d, tc[p], GMP_RNDN);
56 for (l = 1; l < p; l++)
58 mpfr_add_ui (s1, s, t, GMP_RNDN); /* s + (2p-2l) */
59 mpfr_mul (d, d, s1, GMP_RNDN);
61 mpfr_add_ui (s1, s, t, GMP_RNDN); /* s + (2p-2l-1) */
62 mpfr_mul (d, d, s1, GMP_RNDN);
64 mpfr_div_ui (d, d, n2, GMP_RNDN);
65 mpfr_add (d, d, tc[p-l], GMP_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], GMP_RNDN);
71 mpfr_mul (d, d, s, GMP_RNDN);
72 mpfr_add (s1, s, __gmpfr_one, GMP_RNDN);
73 mpfr_neg (s1, s1, GMP_RNDN);
74 mpfr_ui_pow (u, n, s1, GMP_RNDN);
75 mpfr_mul (b, d, u, GMP_RNDN);
77 MPFR_GROUP_CLEAR (group);
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, ...
85 mpfr_zeta_c (int p, mpfr_t *tc)
92 mpfr_init2 (d, MPFR_PREC (tc[1]));
93 mpfr_div_ui (tc[1], __gmpfr_one, 12, GMP_RNDN);
94 for (k = 2; k <= p; k++)
96 mpfr_set_ui (d, k-1, GMP_RNDN);
97 mpfr_div_ui (d, d, 12*k+6, GMP_RNDN);
100 mpfr_div_ui (d, d, 4*(2*k-2*l+3)*(2*k-2*l+2), GMP_RNDN);
101 mpfr_add (d, d, tc[l], GMP_RNDN);
103 mpfr_div_ui (tc[k], d, 24, GMP_RNDN);
104 MPFR_CHANGE_SIGN (tc[k]);
110 /* Input: s - a floating-point number
112 Output: sum - a floating-point number approximating sum(1/i^s, i=1..n-1) */
114 mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n)
118 MPFR_GROUP_DECL (group);
120 MPFR_GROUP_INIT_2 (group, MPFR_PREC (sum), u, s1);
122 mpfr_neg (s1, s, GMP_RNDN);
123 mpfr_ui_pow (u, n, s1, GMP_RNDN);
124 mpfr_div_2ui (u, u, 1, GMP_RNDN);
125 mpfr_set (sum, u, GMP_RNDN);
126 for (i=n-1; i>1; i--)
128 mpfr_ui_pow (u, i, s1, GMP_RNDN);
129 mpfr_add (sum, sum, u, GMP_RNDN);
131 mpfr_add (sum, sum, __gmpfr_one, GMP_RNDN);
133 MPFR_GROUP_CLEAR (group);
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
142 mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
144 mpfr_t b, c, z_pre, f, s1;
145 double beta, sd, dnep;
147 mp_prec_t precz, precs, d, dint;
150 MPFR_GROUP_DECL (group);
151 MPFR_ZIV_DECL (loop);
153 MPFR_ASSERTD (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0);
155 precz = MPFR_PREC (z);
156 precs = MPFR_PREC (s);
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)
170 err = MPFR_GET_EXP (s) - 1;
171 if (err > (mp_exp_t) (sizeof (mp_exp_t)*CHAR_BIT-2))
174 err = ((mp_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,
180 d = precz + MPFR_INT_CEIL_LOG2(precz) + 10;
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, GMP_RNDN);
186 MPFR_ASSERTD (inex == 0);
189 if (MPFR_IS_ZERO (s1))
193 MPFR_ASSERTD (inex == 0);
194 goto clear_and_return;
197 MPFR_GROUP_INIT_4 (group, MPFR_PREC_MIN, b, c, z_pre, f);
199 MPFR_ZIV_INIT (loop, d);
202 /* Principal loop: we compute, in z_pre,
203 an approximation of Zeta(s), that we send to can_round */
204 if (MPFR_GET_EXP (s1) <= -(mp_exp_t) ((mpfr_prec_t) (d-3)/2))
205 /* Branch 1: when s-1 is very small, one
206 uses the approximation Zeta(s)=1/(s-1)+gamma,
207 where gamma is Euler's constant */
209 dint = MAX (d + 3, precs);
210 MPFR_TRACE (printf ("branch 1\ninternal precision=%lu\n",
211 (unsigned long) dint));
212 MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
213 mpfr_div (z_pre, __gmpfr_one, s1, GMP_RNDN);
214 mpfr_const_euler (f, GMP_RNDN);
215 mpfr_add (z_pre, z_pre, f, GMP_RNDN);
221 MPFR_TRACE (printf ("branch 2\n"));
222 /* Computation of parameters n, p and working precision */
223 dnep = (double) d * LOG2;
224 sd = mpfr_get_d (s, GMP_RNDN);
225 /* beta = dnep + 0.61 + sd * log (6.2832 / sd);
226 but a larger value is ok */
227 #define LOG6dot2832 1.83787940484160805532
228 beta = dnep + 0.61 + sd * (LOG6dot2832 - LOG2 *
229 __gmpfr_floor_log2 (sd));
233 /* n = 1 + (int) (exp ((dnep - LOG2) / sd)); */
234 n = 1 + (int) __gmpfr_ceil_exp2 ((d - 1.0) / sd);
238 p = 1 + (int) beta / 2;
239 n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832);
241 MPFR_TRACE (printf ("\nn=%d\np=%d\n",n,p));
242 /* add = 4 + floor(1.5 * log(d) / log (2)).
243 We should have add >= 10, which is always fulfilled since
244 d = precz + 11 >= 12, thus ceil(log2(d)) >= 4 */
245 add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2;
246 MPFR_ASSERTD(add >= 10);
251 MPFR_TRACE (printf ("internal precision=%lu\n",
252 (unsigned long) dint));
254 size = (p + 1) * sizeof(mpfr_t);
255 tc1 = (mpfr_t*) (*__gmp_allocate_func) (size);
257 mpfr_init2 (tc1[l], dint);
258 MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
260 MPFR_TRACE (printf ("precision of z = %lu\n",
261 (unsigned long) precz));
263 /* Computation of the coefficients c_k */
264 mpfr_zeta_c (p, tc1);
265 /* Computation of the 3 parts of the fonction Zeta. */
266 mpfr_zeta_part_a (z_pre, s, n);
267 mpfr_zeta_part_b (b, s, n, p, tc1);
268 /* s1 = s-1 is already computed above */
269 mpfr_div (c, __gmpfr_one, s1, GMP_RNDN);
270 mpfr_ui_pow (f, n, s1, GMP_RNDN);
271 mpfr_div (c, c, f, GMP_RNDN);
272 MPFR_TRACE (MPFR_DUMP (c));
273 mpfr_add (z_pre, z_pre, c, GMP_RNDN);
274 mpfr_add (z_pre, z_pre, b, GMP_RNDN);
277 (*__gmp_free_func) (tc1, size);
281 MPFR_TRACE (MPFR_DUMP (z_pre));
282 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, d-3, precz, rnd_mode)))
284 MPFR_ZIV_NEXT (loop, d);
286 MPFR_ZIV_FREE (loop);
288 inex = mpfr_set (z, z_pre, rnd_mode);
290 MPFR_GROUP_CLEAR (group);
298 mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
300 mpfr_t z_pre, s1, y, p;
301 double sd, eps, m1, c;
303 mp_prec_t precz, prec1, precs, precs1;
305 MPFR_GROUP_DECL (group);
306 MPFR_ZIV_DECL (loop);
307 MPFR_SAVE_EXPO_DECL (expo);
309 MPFR_LOG_FUNC (("s[%#R]=%R rnd=%d", s, s, rnd_mode),
310 ("z[%#R]=%R inexact=%d", z, z, inex));
312 /* Zero, Nan or Inf ? */
313 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s)))
320 else if (MPFR_IS_INF (s))
323 return mpfr_set_ui (z, 1, GMP_RNDN); /* Zeta(+Inf) = 1 */
324 MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */
329 MPFR_ASSERTD (MPFR_IS_ZERO (s));
330 mpfr_set_ui (z, 1, rnd_mode);
331 mpfr_div_2ui (z, z, 1, rnd_mode);
332 MPFR_CHANGE_SIGN (z);
337 /* s is neither Nan, nor Inf, nor Zero */
339 /* check tiny s: we have zeta(s) = -1/2 - 1/2 log(2 Pi) s + ... around s=0,
340 and for |s| <= 0.074, we have |zeta(s) + 1/2| <= |s|.
341 Thus if |s| <= 1/4*ulp(1/2), we can deduce the correct rounding
342 (the 1/4 covers the case where |zeta(s)| < 1/2 and rounding to nearest).
343 A sufficient condition is that EXP(s) + 1 < -PREC(z). */
344 if (MPFR_EXP(s) + 1 < - (mp_exp_t) MPFR_PREC(z))
346 int signs = MPFR_SIGN(s);
347 mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */
348 if ((rnd_mode == GMP_RNDU || rnd_mode == GMP_RNDZ) && signs < 0)
350 mpfr_nextabove (z); /* z = -1/2 + epsilon */
353 else if (rnd_mode == GMP_RNDD && signs > 0)
355 mpfr_nextbelow (z); /* z = -1/2 - epsilon */
360 if (rnd_mode == GMP_RNDU) /* s > 0: z = -1/2 */
362 else if (rnd_mode == GMP_RNDD)
363 inex = -1; /* s < 0: z = -1/2 */
364 else /* (GMP_RNDZ and s > 0) or GMP_RNDN: z = -1/2 */
365 inex = (signs > 0) ? 1 : -1;
367 return mpfr_check_range (z, inex, rnd_mode);
370 /* Check for case s= -2n */
375 MPFR_EXP (tmp) = MPFR_EXP (s) - 1;
376 if (mpfr_integer_p (tmp))
384 MPFR_SAVE_EXPO_MARK (expo);
387 if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */
388 inex = mpfr_zeta_pos (z, s, rnd_mode);
389 else /* use reflection formula
390 zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */
394 precz = MPFR_PREC (z);
395 precs = MPFR_PREC (s);
397 /* Precision precs1 needed to represent 1 - s, and s + 2,
398 without any truncation */
399 precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s));
400 sd = mpfr_get_d (s, GMP_RNDN) - 1.0;
402 sd = -sd; /* now sd = abs(s-1.0) */
403 /* Precision prec1 is the precision on elementary computations;
404 it ensures a final precision prec1 - add for zeta(s) */
405 /* eps = pow (2.0, - (double) precz - 14.0); */
406 eps = __gmpfr_ceil_exp2 (- (double) precz - 14.0);
407 m1 = 1.0 + MAX(1.0 / eps, 2.0 * sd) * (1.0 + eps);
408 c = (1.0 + eps) * (1.0 + eps * MAX(8.0, m1));
409 /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */
410 add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1));
412 prec1 = MAX (prec1, precs1) + 10;
414 MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p);
415 MPFR_ZIV_INIT (loop, prec1);
418 mpfr_sub (s1, __gmpfr_one, s, GMP_RNDN);/* s1 = 1-s */
419 mpfr_zeta_pos (z_pre, s1, GMP_RNDN); /* zeta(1-s) */
420 mpfr_gamma (y, s1, GMP_RNDN); /* gamma(1-s) */
421 if (MPFR_IS_INF (y)) /* Zeta(s) < 0 for -4k-2 < s < -4k,
422 Zeta(s) > 0 for -4k < s < -4k+2 */
424 mpfr_div_2ui (s1, s, 2, GMP_RNDN); /* s/4, exact */
425 mpfr_frac (s1, s1, GMP_RNDN); /* exact, -1 < s1 < 0 */
426 overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1;
429 mpfr_mul (z_pre, z_pre, y, GMP_RNDN); /* gamma(1-s)*zeta(1-s) */
430 mpfr_const_pi (p, GMP_RNDD);
431 mpfr_mul (y, s, p, GMP_RNDN);
432 mpfr_div_2ui (y, y, 1, GMP_RNDN); /* s*Pi/2 */
433 mpfr_sin (y, y, GMP_RNDN); /* sin(Pi*s/2) */
434 mpfr_mul (z_pre, z_pre, y, GMP_RNDN);
435 mpfr_mul_2ui (y, p, 1, GMP_RNDN); /* 2*Pi */
436 mpfr_neg (s1, s1, GMP_RNDN); /* s-1 */
437 mpfr_pow (y, y, s1, GMP_RNDN); /* (2*Pi)^(s-1) */
438 mpfr_mul (z_pre, z_pre, y, GMP_RNDN);
439 mpfr_mul_2ui (z_pre, z_pre, 1, GMP_RNDN);
441 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz,
445 MPFR_ZIV_NEXT (loop, prec1);
446 MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p);
448 MPFR_ZIV_FREE (loop);
451 inex = mpfr_overflow (z, rnd_mode, overflow);
452 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
455 inex = mpfr_set (z, z_pre, rnd_mode);
456 MPFR_GROUP_CLEAR (group);
459 MPFR_SAVE_EXPO_FREE (expo);
460 return mpfr_check_range (z, inex, rnd_mode);