1 /* mpfr_zeta -- compute the Riemann Zeta function
3 Copyright 2003-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
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 3 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.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. */
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
], MPFR_RNDN
);
56 for (l
= 1; l
< p
; l
++)
58 mpfr_add_ui (s1
, s
, t
, MPFR_RNDN
); /* s + (2p-2l) */
59 mpfr_mul (d
, d
, s1
, MPFR_RNDN
);
61 mpfr_add_ui (s1
, s
, t
, MPFR_RNDN
); /* s + (2p-2l-1) */
62 mpfr_mul (d
, d
, s1
, MPFR_RNDN
);
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
);
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
);
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, MPFR_RNDN
);
94 for (k
= 2; k
<= p
; k
++)
96 mpfr_set_ui (d
, k
-1, MPFR_RNDN
);
97 mpfr_div_ui (d
, d
, 12*k
+6, MPFR_RNDN
);
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
);
103 mpfr_div_ui (tc
[k
], d
, 24, MPFR_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
, 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
--)
128 mpfr_ui_pow (u
, i
, s1
, MPFR_RNDN
);
129 mpfr_add (sum
, sum
, u
, MPFR_RNDN
);
131 mpfr_add (sum
, sum
, __gmpfr_one
, MPFR_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
, mpfr_rnd_t rnd_mode
)
144 mpfr_t b
, c
, z_pre
, f
, s1
;
145 double beta
, sd
, dnep
;
147 mpfr_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
> (mpfr_exp_t
) (sizeof (mpfr_exp_t
)*CHAR_BIT
-2))
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,
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
, MPFR_RNDN
);
186 MPFR_ASSERTD (inex
== 0);
188 /* case s=1 should have already been handled */
189 MPFR_ASSERTD (!MPFR_IS_ZERO (s1
));
191 MPFR_GROUP_INIT_4 (group
, MPFR_PREC_MIN
, b
, c
, z_pre
, f
);
193 MPFR_ZIV_INIT (loop
, d
);
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 */
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
);
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
));
227 /* n = 1 + (int) (exp ((dnep - LOG2) / sd)); */
228 n
= 1 + (int) __gmpfr_ceil_exp2 ((d
- 1.0) / sd
);
232 p
= 1 + (int) beta
/ 2;
233 n
= 1 + (int) ((sd
+ 2.0 * (double) p
- 1.0) / 6.2832);
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);
245 MPFR_TRACE (printf ("internal precision=%lu\n",
246 (unsigned long) dint
));
248 size
= (p
+ 1) * sizeof(mpfr_t
);
249 tc1
= (mpfr_t
*) (*__gmp_allocate_func
) (size
);
251 mpfr_init2 (tc1
[l
], dint
);
252 MPFR_GROUP_REPREC_4 (group
, dint
, b
, c
, z_pre
, f
);
254 MPFR_TRACE (printf ("precision of z = %lu\n",
255 (unsigned long) precz
));
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
);
271 (*__gmp_free_func
) (tc1
, size
);
275 MPFR_TRACE (MPFR_DUMP (z_pre
));
276 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre
, d
-3, precz
, rnd_mode
)))
278 MPFR_ZIV_NEXT (loop
, d
);
280 MPFR_ZIV_FREE (loop
);
282 inex
= mpfr_set (z
, z_pre
, rnd_mode
);
284 MPFR_GROUP_CLEAR (group
);
291 mpfr_zeta (mpfr_t z
, mpfr_srcptr s
, mpfr_rnd_t rnd_mode
)
293 mpfr_t z_pre
, s1
, y
, p
;
294 double sd
, eps
, m1
, c
;
296 mpfr_prec_t precz
, prec1
, precs
, precs1
;
298 MPFR_GROUP_DECL (group
);
299 MPFR_ZIV_DECL (loop
);
300 MPFR_SAVE_EXPO_DECL (expo
);
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
));
306 /* Zero, Nan or Inf ? */
307 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s
)))
314 else if (MPFR_IS_INF (s
))
317 return mpfr_set_ui (z
, 1, MPFR_RNDN
); /* Zeta(+Inf) = 1 */
318 MPFR_SET_NAN (z
); /* Zeta(-Inf) = NaN */
323 MPFR_ASSERTD (MPFR_IS_ZERO (s
));
324 return mpfr_set_si_2exp (z
, -1, -1, rnd_mode
);
328 /* s is neither Nan, nor Inf, nor Zero */
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
))
337 int signs
= MPFR_SIGN(s
);
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)
345 mpfr_nextabove (z
); /* z = -1/2 + epsilon */
348 else if (rnd_mode
== MPFR_RNDD
&& signs
> 0)
350 mpfr_nextbelow (z
); /* z = -1/2 - epsilon */
355 if (rnd_mode
== MPFR_RNDU
) /* s > 0: z = -1/2 */
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;
362 MPFR_SAVE_EXPO_FREE (expo
);
363 return mpfr_check_range (z
, inex
, rnd_mode
);
366 /* Check for case s= -2n */
371 MPFR_EXP (tmp
) = MPFR_GET_EXP (s
) - 1;
372 if (mpfr_integer_p (tmp
))
380 /* Check for case s= 1 before changing the exponent range */
381 if (mpfr_cmp (s
, __gmpfr_one
) ==0)
389 MPFR_SAVE_EXPO_MARK (expo
);
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) */
399 precz
= MPFR_PREC (z
);
400 precs
= MPFR_PREC (s
);
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;
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
));
417 prec1
= MAX (prec1
, precs1
) + 10;
419 MPFR_GROUP_INIT_4 (group
, prec1
, z_pre
, s1
, y
, p
);
420 MPFR_ZIV_INIT (loop
, prec1
);
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 */
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;
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
);
446 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre
, prec1
- add
, precz
,
450 MPFR_ZIV_NEXT (loop
, prec1
);
451 MPFR_GROUP_REPREC_4 (group
, prec1
, z_pre
, s1
, y
, p
);
453 MPFR_ZIV_FREE (loop
);
456 inex
= mpfr_overflow (z
, rnd_mode
, overflow
);
457 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, MPFR_FLAGS_OVERFLOW
);
460 inex
= mpfr_set (z
, z_pre
, rnd_mode
);
461 MPFR_GROUP_CLEAR (group
);
464 MPFR_SAVE_EXPO_FREE (expo
);
465 return mpfr_check_range (z
, inex
, rnd_mode
);