1 /* mpfr_zeta -- compute the Riemann Zeta function
2
3 Copyright 2003-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 https://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 <float.h> /* for DBL_MAX */
24
25 #define MPFR_NEED_LONGLONG_H
26 #include "mpfr-impl.h"
27
28 /*
29 Parameters:
30 s - the input floating-point number
31 n, p - parameters from the algorithm
32 tc - an array of p floating-point numbers tc[1]..tc[p]
33 Output:
34 b is the result, i.e.
35 sum(tc[i]*product((s+2j)*(s+2j-1)/n^2,j=1..i-1), i=1..p)*s*n^(-s-1)
36 */
37 static void
38 mpfr_zeta_part_b (mpfr_ptr b, mpfr_srcptr s, int n, int p, mpfr_t *tc)
39 {
40 mpfr_t s1, d, u;
41 unsigned long n2;
42 int l, t;
43 MPFR_GROUP_DECL (group);
44
45 if (p == 0)
46 {
47 MPFR_SET_ZERO (b);
48 MPFR_SET_POS (b);
49 return;
50 }
51
52 n2 = n * n;
53 MPFR_GROUP_INIT_3 (group, MPFR_PREC (b), s1, d, u);
54
55 /* t equals 2p-2, 2p-3, ... ; s1 equals s+t */
56 t = 2 * p - 2;
57 mpfr_set (d, tc[p], MPFR_RNDN);
58 for (l = 1; l < p; l++)
59 {
60 mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l) */
61 mpfr_mul (d, d, s1, MPFR_RNDN);
62 t = t - 1;
63 mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l-1) */
64 mpfr_mul (d, d, s1, MPFR_RNDN);
65 t = t - 1;
66 mpfr_div_ui (d, d, n2, MPFR_RNDN);
67 mpfr_add (d, d, tc[p-l], MPFR_RNDN);
68 /* since s is positive and the tc[i] have alternate signs,
69 the following is unlikely */
70 if (MPFR_UNLIKELY (mpfr_cmpabs (d, tc[p-l]) > 0))
71 mpfr_set (d, tc[p-l], MPFR_RNDN);
72 }
73 mpfr_mul (d, d, s, MPFR_RNDN);
74 mpfr_add (s1, s, __gmpfr_one, MPFR_RNDN);
75 mpfr_neg (s1, s1, MPFR_RNDN);
76 mpfr_ui_pow (u, n, s1, MPFR_RNDN);
77 mpfr_mul (b, d, u, MPFR_RNDN);
78
79 MPFR_GROUP_CLEAR (group);
80 }
81
82 /* Input: p - an integer
83 Output: fills tc[1..p], tc[i] = bernoulli(2i)/(2i)!
84 tc[1]=1/12, tc[2]=-1/720, tc[3]=1/30240, ...
85 Assumes all the tc[i] have the same precision.
86
87 Uses the recurrence (4.60) from the book "Modern Computer Arithmetic"
88 by Brent and Zimmermann for C_k = bernoulli(2k)/(2k)!:
89 sum(C_k/(2k+1-2j)!/4^(k-j), j=0..k) = 1/(2k)!/4^k
90 If we put together the terms involving C_0 and C_1 we get:
91 sum(D_k/(2k+1-2j)!/4^(k-j), j=1..k) = 0
92 with D_1 = C_0/4/(2k+1)/(2k)+C_1-1/(2k)/4=(k-1)/(12k+6),
93 and D_k = C_k for k >= 2.
94
95 FIXME: we have C_k = (-1)^(k-1) 2/(2pi)^(2k) * zeta(2k),
96 see for example formula (4.65) from the above book,
97 thus since |zeta(2k)-1| < 2^(1-2k) for k >= 2, we have:
98 |C_k - E_k| < E_k * 2^(1-2k) for k >= 2 and E_k := (-1)^(k-1) 2/(2pi)^(2k).
99 Then if 2k-1 >= prec we can evaluate E_k instead, which only requires one
100 multiplication per term, instead of O(k) small divisions.
101 */
102 static void
103 mpfr_zeta_c (int p, mpfr_t *tc)
104 {
105 if (p > 0)
106 {
107 mpfr_t d;
108 int k, l;
109 mpfr_prec_t prec = MPFR_PREC (tc[1]);
110
111 mpfr_init2 (d, prec);
112 mpfr_div_ui (tc[1], __gmpfr_one, 12, MPFR_RNDN);
113 for (k = 2; k <= p; k++)
114 {
115 mpfr_set_ui (d, k-1, MPFR_RNDN);
116 mpfr_div_ui (d, d, 12*k+6, MPFR_RNDN);
117 for (l=2; l < k; l++)
118 {
119 mpfr_div_ui (d, d, 4*(2*k-2*l+3)*(2*k-2*l+2), MPFR_RNDN);
120 mpfr_add (d, d, tc[l], MPFR_RNDN);
121 }
122 mpfr_div_ui (tc[k], d, 24, MPFR_RNDN);
123 MPFR_CHANGE_SIGN (tc[k]);
124 }
125 mpfr_clear (d);
126 }
127 }
128
129 /* Input: s - a floating-point number
130 n - an integer
131 Output: sum - a floating-point number approximating sum(1/i^s, i=1..n-1) */
132 static void
133 mpfr_zeta_part_a (mpfr_ptr sum, mpfr_srcptr s, int n)
134 {
135 mpfr_t u, s1;
136 int i;
137 MPFR_GROUP_DECL (group);
138
139 MPFR_GROUP_INIT_2 (group, MPFR_PREC (sum), u, s1);
140
141 mpfr_neg (s1, s, MPFR_RNDN);
142 mpfr_ui_pow (u, n, s1, MPFR_RNDN);
143 mpfr_div_2ui (u, u, 1, MPFR_RNDN);
144 mpfr_set (sum, u, MPFR_RNDN);
145 for (i=n-1; i>1; i--)
146 {
147 mpfr_ui_pow (u, i, s1, MPFR_RNDN);
148 mpfr_add (sum, sum, u, MPFR_RNDN);
149 }
150 mpfr_add (sum, sum, __gmpfr_one, MPFR_RNDN);
151
152 MPFR_GROUP_CLEAR (group);
153 }
154
155 /* Input: s - a floating-point number >= 1/2.
156 rnd_mode - a rounding mode.
157 Assumes s is neither NaN nor Infinite.
158 Output: z - Zeta(s) rounded to the precision of z with direction rnd_mode
159 */
160 static int
161 mpfr_zeta_pos (mpfr_ptr z, mpfr_srcptr s, mpfr_rnd_t rnd_mode)
162 {
163 mpfr_t b, c, z_pre, f, s1;
164 double beta, sd, dnep;
165 mpfr_t *tc1;
166 mpfr_prec_t precz, precs, d, dint;
167 int p, n, l, add;
168 int inex;
169 MPFR_GROUP_DECL (group);
170 MPFR_ZIV_DECL (loop);
171
172 MPFR_ASSERTD (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0);
173
174 precz = MPFR_PREC (z);
175 precs = MPFR_PREC (s);
176
177 /* Zeta(x) = 1+1/2^x+1/3^x+1/4^x+1/5^x+O(1/6^x)
178 so with 2^(EXP(x)-1) <= x < 2^EXP(x)
179 So for x > 2^3, k^x > k^8, so 2/k^x < 2/k^8
180 Zeta(x) = 1 + 1/2^x*(1+(2/3)^x+(2/4)^x+...)
181 = 1 + 1/2^x*(1+sum((2/k)^x,k=3..infinity))
182 <= 1 + 1/2^x*(1+sum((2/k)^8,k=3..infinity))
183 And sum((2/k)^8,k=3..infinity) = -257+128*Pi^8/4725 ~= 0.0438035
184 So Zeta(x) <= 1 + 1/2^x*2 for x >= 8
185 The error is < 2^(-x+1) <= 2^(-2^(EXP(x)-1)+1) */
186 if (MPFR_GET_EXP (s) > 3)
187 {
188 mpfr_exp_t err;
189 err = MPFR_GET_EXP (s) - 1;
190 if (err > (mpfr_exp_t) (sizeof (mpfr_exp_t)*CHAR_BIT-2))
191 err = MPFR_EMAX_MAX;
192 else
193 err = ((mpfr_exp_t)1) << err;
194 err = 1 - (-err+1); /* GET_EXP(one) - (-err+1) = err :) */
195 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, err, 0, 1,
196 rnd_mode, {});
197 }
198
199 d = precz + MPFR_INT_CEIL_LOG2(precz) + 10;
200
201 /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */
202 dint = (mpfr_uexp_t) MPFR_GET_EXP (s);
203 mpfr_init2 (s1, MAX (precs, dint));
204 inex = mpfr_sub (s1, s, __gmpfr_one, MPFR_RNDN);
205 MPFR_ASSERTD (inex == 0);
206
207 /* case s=1 should have already been handled */
208 MPFR_ASSERTD (!MPFR_IS_ZERO (s1));
209
210 MPFR_GROUP_INIT_4 (group, MPFR_PREC_MIN, b, c, z_pre, f);
211
212 MPFR_ZIV_INIT (loop, d);
213 for (;;)
214 {
215 /* Principal loop: we compute, in z_pre,
216 an approximation of Zeta(s), that we send to can_round */
217 if (MPFR_GET_EXP (s1) <= -(mpfr_exp_t) ((mpfr_prec_t) (d-3)/2))
218 /* Branch 1: when s-1 is very small, one
219 uses the approximation Zeta(s)=1/(s-1)+gamma,
220 where gamma is Euler's constant */
221 {
222 dint = MAX (d + 3, precs);
223 /* branch 1, with internal precision dint */
224 MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
225 mpfr_div (z_pre, __gmpfr_one, s1, MPFR_RNDN);
226 mpfr_const_euler (f, MPFR_RNDN);
227 mpfr_add (z_pre, z_pre, f, MPFR_RNDN);
228 }
229 else /* Branch 2 */
230 {
231 size_t size;
232
233 /* branch 2 */
234 /* Computation of parameters n, p and working precision */
235 dnep = (double) d * LOG2;
236 sd = mpfr_get_d (s, MPFR_RNDN);
237 /* beta = dnep + 0.61 + sd * log (6.2832 / sd);
238 but a larger value is OK */
239 #define LOG6dot2832 1.83787940484160805532
240 beta = dnep + 0.61 + sd * (LOG6dot2832 - LOG2 *
241 __gmpfr_floor_log2 (sd));
242 if (beta <= 0.0)
243 {
244 p = 0;
245 /* n = 1 + (int) (exp ((dnep - LOG2) / sd)); */
246 n = 1 + (int) __gmpfr_ceil_exp2 ((d - 1.0) / sd);
247 }
248 else
249 {
250 p = 1 + (int) beta / 2;
251 n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832);
252 }
253 /* add = 4 + floor(1.5 * log(d) / log (2)).
254 We should have add >= 10, which is always fulfilled since
255 d = precz + 11 >= 12, thus ceil(log2(d)) >= 4 */
256 add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2;
257 MPFR_ASSERTD(add >= 10);
258 dint = d + add;
259 if (dint < precs)
260 dint = precs;
261
262 /* internal precision is dint */
263
264 size = (p + 1) * sizeof(mpfr_t);
265 tc1 = (mpfr_t*) mpfr_allocate_func (size);
266 for (l=1; l<=p; l++)
267 mpfr_init2 (tc1[l], dint);
268 MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
269
270 /* precision of z is precz */
271
272 /* Computation of the coefficients c_k */
273 mpfr_zeta_c (p, tc1);
274 /* Computation of the 3 parts of the function Zeta. */
275 mpfr_zeta_part_a (z_pre, s, n);
276 mpfr_zeta_part_b (b, s, n, p, tc1);
277 /* s1 = s-1 is already computed above */
278 mpfr_div (c, __gmpfr_one, s1, MPFR_RNDN);
279 mpfr_ui_pow (f, n, s1, MPFR_RNDN);
280 mpfr_div (c, c, f, MPFR_RNDN);
281 mpfr_add (z_pre, z_pre, c, MPFR_RNDN);
282 mpfr_add (z_pre, z_pre, b, MPFR_RNDN);
283 for (l=1; l<=p; l++)
284 mpfr_clear (tc1[l]);
285 mpfr_free_func (tc1, size);
286 /* End branch 2 */
287 }
288
289 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, d-3, precz, rnd_mode)))
290 break;
291 MPFR_ZIV_NEXT (loop, d);
292 }
293 MPFR_ZIV_FREE (loop);
294
295 inex = mpfr_set (z, z_pre, rnd_mode);
296
297 MPFR_GROUP_CLEAR (group);
298 mpfr_clear (s1);
299
300 return inex;
301 }
302
303 /* TODO: Check the error analysis. The following (undocumented?) one
304 does not take into account the replacement of sin(Pi*s/2) by sinpi(s/2)
305 in commit fd5d811d81f6d1839d4099cc1bb2cde705981648, which could have
306 reduced the error bound since the multiplication by Pi is now exact. */
307 /* return add = 1 + floor(log(c^3*(13+m1))/log(2))
308 where c = (1+eps)*(1+eps*max(8,m1)),
309 m1 = 1 + max(1/eps,2*sd)*(1+eps),
310 eps = 2^(-precz-14)
311 sd = abs(s-1)
312 */
313 static long
314 compute_add (mpfr_srcptr s, mpfr_prec_t precz)
315 {
316 mpfr_t t, u, m1;
317 long add;
318
319 mpfr_inits2 (64, t, u, m1, (mpfr_ptr) 0);
320 if (mpfr_cmp_ui (s, 1) >= 0)
321 mpfr_sub_ui (t, s, 1, MPFR_RNDU);
322 else
323 mpfr_ui_sub (t, 1, s, MPFR_RNDU);
324 /* now t = sd = abs(s-1), rounded up */
325 mpfr_set_ui_2exp (u, 1, - precz - 14, MPFR_RNDU);
326 /* u = eps */
327 /* since 1/eps = 2^(precz+14), if EXP(sd) >= precz+14, then
328 sd >= 1/2*2^(precz+14) thus 2*sd >= 2^(precz+14) >= 1/eps */
329 if (mpfr_get_exp (t) >= precz + 14)
330 mpfr_mul_2ui (t, t, 1, MPFR_RNDU);
331 else
332 mpfr_set_ui_2exp (t, 1, precz + 14, MPFR_RNDU);
333 /* now t = max(1/eps,2*sd) */
334 mpfr_add_ui (u, u, 1, MPFR_RNDU); /* u = 1+eps, rounded up */
335 mpfr_mul (t, t, u, MPFR_RNDU); /* t = max(1/eps,2*sd)*(1+eps) */
336 mpfr_add_ui (m1, t, 1, MPFR_RNDU);
337 if (mpfr_get_exp (m1) <= 3)
338 mpfr_set_ui (t, 8, MPFR_RNDU);
339 else
340 mpfr_set (t, m1, MPFR_RNDU);
341 /* now t = max(8,m1) */
342 mpfr_div_2ui (t, t, precz + 14, MPFR_RNDU); /* eps*max(8,m1) */
343 mpfr_add_ui (t, t, 1, MPFR_RNDU); /* 1+eps*max(8,m1) */
344 mpfr_mul (t, t, u, MPFR_RNDU); /* t = c */
345 mpfr_add_ui (u, m1, 13, MPFR_RNDU); /* 13+m1 */
346 mpfr_mul (u, u, t, MPFR_RNDU); /* c*(13+m1) */
347 mpfr_sqr (t, t, MPFR_RNDU); /* c^2 */
348 mpfr_mul (u, u, t, MPFR_RNDU); /* c^3*(13+m1) */
349 add = mpfr_get_exp (u);
350 mpfr_clears (t, u, m1, (mpfr_ptr) 0);
351 return add;
352 }
353
354 /* return in z a lower bound (for rnd = RNDD) or upper bound (for rnd = RNDU)
355 of |zeta(s)|/2, using:
356 log(|zeta(s)|/2) = (s-1)*log(2*Pi) + lngamma(1-s)
357 + log(|sin(Pi*s/2)| * zeta(1-s)).
358 Assumes s < 1/2 and s1 = 1-s exactly, thus s1 > 1/2.
359 y and p are temporary variables.
360 At input, p is Pi rounded down.
361 The comments in the code are for rnd = RNDD. */
362 static void
363 mpfr_reflection_overflow (mpfr_ptr z, mpfr_ptr s1, mpfr_srcptr s, mpfr_ptr y,
364 mpfr_ptr p, mpfr_rnd_t rnd)
365 {
366 mpz_t sint;
367
368 MPFR_ASSERTD (rnd == MPFR_RNDD || rnd == MPFR_RNDU);
369
370 /* Since log is increasing, we want lower bounds on |sin(Pi*s/2)| and
371 zeta(1-s). */
372 mpz_init (sint);
373 mpfr_get_z (sint, s, MPFR_RNDD); /* sint = floor(s) */
374 /* We first compute a lower bound of |sin(Pi*s/2)|, which is a periodic
375 function of period 2. Thus:
376 if 2k < s < 2k+1, then |sin(Pi*s/2)| is increasing;
377 if 2k-1 < s < 2k, then |sin(Pi*s/2)| is decreasing.
378 These cases are distinguished by testing bit 0 of floor(s) as if
379 represented in two's complement (or equivalently, as an unsigned
380 integer mod 2):
381 0: sint = 0 mod 2, thus 2k < s < 2k+1 and |sin(Pi*s/2)| is increasing;
382 1: sint = 1 mod 2, thus 2k-1 < s < 2k and |sin(Pi*s/2)| is decreasing.
383 Let's recall that the comments are for rnd = RNDD. */
384 if (mpz_tstbit (sint, 0) == 0) /* |sin(Pi*s/2)| is increasing: round down
385 Pi*s to get a lower bound. */
386 {
387 mpfr_mul (y, p, s, rnd);
388 if (rnd == MPFR_RNDD)
389 mpfr_nextabove (p); /* we will need p rounded above afterwards */
390 }
391 else /* |sin(Pi*s/2)| is decreasing: round up Pi*s to get a lower bound. */
392 {
393 if (rnd == MPFR_RNDD)
394 mpfr_nextabove (p);
395 mpfr_mul (y, p, s, MPFR_INVERT_RND(rnd));
396 }
397 mpfr_div_2ui (y, y, 1, MPFR_RNDN); /* exact, rounding mode doesn't matter */
398 /* The rounding direction of sin depends on its sign. We have:
399 if -4k-2 < s < -4k, then -2k-1 < s/2 < -2k, thus sin(Pi*s/2) < 0;
400 if -4k < s < -4k+2, then -2k < s/2 < -2k+1, thus sin(Pi*s/2) > 0.
401 These cases are distinguished by testing bit 1 of floor(s) as if
402 represented in two's complement (or equivalently, as an unsigned
403 integer mod 4):
404 0: sint = {0,1} mod 4, thus -2k < s/2 < -2k+1 and sin(Pi*s/2) > 0;
405 1: sint = {2,3} mod 4, thus -2k-1 < s/2 < -2k and sin(Pi*s/2) < 0.
406 Let's recall that the comments are for rnd = RNDD. */
407 if (mpz_tstbit (sint, 1) == 0) /* -2k < s/2 < -2k+1; sin(Pi*s/2) > 0 */
408 {
409 /* Round sin down to get a lower bound of |sin(Pi*s/2)|. */
410 mpfr_sin (y, y, rnd);
411 }
412 else /* -2k-1 < s/2 < -2k; sin(Pi*s/2) < 0 */
413 {
414 /* Round sin up to get a lower bound of |sin(Pi*s/2)|. */
415 mpfr_sin (y, y, MPFR_INVERT_RND(rnd));
416 mpfr_abs (y, y, MPFR_RNDN); /* exact, rounding mode doesn't matter */
417 }
418 mpz_clear (sint);
419 /* now y <= |sin(Pi*s/2)| when rnd=RNDD, y >= |sin(Pi*s/2)| when rnd=RNDU */
420 mpfr_zeta_pos (z, s1, rnd); /* zeta(1-s) */
421 mpfr_mul (z, z, y, rnd);
422 /* now z <= |sin(Pi*s/2)|*zeta(1-s) */
423 mpfr_log (z, z, rnd);
424 /* now z <= log(|sin(Pi*s/2)|*zeta(1-s)) */
425 mpfr_lngamma (y, s1, rnd);
426 mpfr_add (z, z, y, rnd);
427 /* z <= lngamma(1-s) + log(|sin(Pi*s/2)|*zeta(1-s)) */
428 /* since s-1 < 0, we want to round log(2*pi) upwards */
429 mpfr_mul_2ui (y, p, 1, MPFR_INVERT_RND(rnd));
430 mpfr_log (y, y, MPFR_INVERT_RND(rnd));
431 mpfr_mul (y, y, s1, MPFR_INVERT_RND(rnd));
432 mpfr_sub (z, z, y, rnd);
433 mpfr_exp (z, z, rnd);
434 if (rnd == MPFR_RNDD)
435 mpfr_nextbelow (p); /* restore original p */
436 }
437
438 int
439 mpfr_zeta (mpfr_ptr z, mpfr_srcptr s, mpfr_rnd_t rnd_mode)
440 {
441 mpfr_t z_pre, s1, y, p;
442 long add;
443 mpfr_prec_t precz, prec1, precs, precs1;
444 int inex;
445 MPFR_GROUP_DECL (group);
446 MPFR_ZIV_DECL (loop);
447 MPFR_SAVE_EXPO_DECL (expo);
448
449 MPFR_LOG_FUNC (
450 ("s[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (s), mpfr_log_prec, s, rnd_mode),
451 ("z[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (z), mpfr_log_prec, z, inex));
452
453 /* Zero, Nan or Inf ? */
454 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s)))
455 {
456 if (MPFR_IS_NAN (s))
457 {
458 MPFR_SET_NAN (z);
459 MPFR_RET_NAN;
460 }
461 else if (MPFR_IS_INF (s))
462 {
463 if (MPFR_IS_POS (s))
464 return mpfr_set_ui (z, 1, MPFR_RNDN); /* Zeta(+Inf) = 1 */
465 MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */
466 MPFR_RET_NAN;
467 }
468 else /* s iz zero */
469 {
470 MPFR_ASSERTD (MPFR_IS_ZERO (s));
471 return mpfr_set_si_2exp (z, -1, -1, rnd_mode);
472 }
473 }
474
475 /* s is neither Nan, nor Inf, nor Zero */
476
477 /* check tiny s: we have zeta(s) = -1/2 - 1/2 log(2 Pi) s + ... around s=0,
478 and for |s| <= 2^(-4), we have |zeta(s) + 1/2| <= |s|.
479 EXP(s) + 1 < -PREC(z) is a sufficient condition to be able to round
480 correctly, for any PREC(z) >= 1 (see algorithms.tex for details). */
481 if (MPFR_GET_EXP (s) + 1 < - (mpfr_exp_t) MPFR_PREC(z))
482 {
483 int signs = MPFR_SIGN(s);
484
485 MPFR_SAVE_EXPO_MARK (expo);
486 mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */
487 if (rnd_mode == MPFR_RNDA)
488 rnd_mode = MPFR_RNDD; /* the result is around -1/2, thus negative */
489 if ((rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDZ) && signs < 0)
490 {
491 mpfr_nextabove (z); /* z = -1/2 + epsilon */
492 inex = 1;
493 }
494 else if (rnd_mode == MPFR_RNDD && signs > 0)
495 {
496 mpfr_nextbelow (z); /* z = -1/2 - epsilon */
497 inex = -1;
498 }
499 else
500 {
501 if (rnd_mode == MPFR_RNDU) /* s > 0: z = -1/2 */
502 inex = 1;
503 else if (rnd_mode == MPFR_RNDD)
504 inex = -1; /* s < 0: z = -1/2 */
505 else /* (MPFR_RNDZ and s > 0) or MPFR_RNDN: z = -1/2 */
506 inex = (signs > 0) ? 1 : -1;
507 }
508 MPFR_SAVE_EXPO_FREE (expo);
509 return mpfr_check_range (z, inex, rnd_mode);
510 }
511
512 /* Check for case s= -2n */
513 if (MPFR_IS_NEG (s))
514 {
515 mpfr_t tmp;
516 tmp[0] = *s;
517 MPFR_EXP (tmp) = MPFR_GET_EXP (s) - 1;
518 if (mpfr_integer_p (tmp))
519 {
520 MPFR_SET_ZERO (z);
521 MPFR_SET_POS (z);
522 MPFR_RET (0);
523 }
524 }
525
526 /* Check for case s=1 before changing the exponent range */
527 if (mpfr_equal_p (s, __gmpfr_one))
528 {
529 MPFR_SET_INF (z);
530 MPFR_SET_POS (z);
531 MPFR_SET_DIVBY0 ();
532 MPFR_RET (0);
533 }
534
535 MPFR_SAVE_EXPO_MARK (expo);
536
537 /* Compute Zeta */
538 if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */
539 inex = mpfr_zeta_pos (z, s, rnd_mode);
540 else /* use reflection formula
541 zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */
542 {
543 int overflow = 0;
544
545 precz = MPFR_PREC (z);
546 precs = MPFR_PREC (s);
547
548 /* Precision precs1 needed to represent 1 - s, and s + 2,
549 without any truncation */
550 precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s));
551 /* Precision prec1 is the precision on elementary computations;
552 it ensures a final precision prec1 - add for zeta(s) */
553 add = compute_add (s, precz);
554 prec1 = precz + add;
555 /* FIXME: To avoid that the working precision (prec1) depends on the
556 input precision, one would need to take into account the error made
557 when s1 is not exactly 1-s when computing zeta(s1) and gamma(s1)
558 below, and also in the case y=Inf (i.e. when gamma(s1) overflows).
559 Make sure that underflows do not occur in intermediate computations.
560 Due to the limited precision, they are probably not possible
561 in practice; add some MPFR_ASSERTN's to be sure that problems
562 do not remain undetected? */
563 prec1 = MAX (prec1, precs1) + 10;
564
565 MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p);
566 MPFR_ZIV_INIT (loop, prec1);
567 for (;;)
568 {
569 mpfr_t z_up;
570
571 mpfr_const_pi (p, MPFR_RNDD); /* p is Pi */
572
573 mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN); /* s1 = 1-s */
574 mpfr_gamma (y, s1, MPFR_RNDN); /* gamma(1-s) */
575 if (MPFR_IS_INF (y)) /* zeta(s) < 0 for -4k-2 < s < -4k,
576 zeta(s) > 0 for -4k < s < -4k+2 */
577 {
578 /* FIXME: An overflow in gamma(s1) does not imply that
579 zeta(s) will overflow. A solution:
580 1. Compute
581 log(|zeta(s)|/2) = (s-1)*log(2*pi) + lngamma(1-s)
582 + log(abs(sin(Pi*s/2)) * zeta(1-s))
583 (possibly sharing computations with the normal case)
584 with a rather good accuracy (see (2)).
585 Memorize the sign of sin(...) for the final sign.
586 2. Take the exponential, ~= |zeta(s)|/2. If there is an
587 overflow, then this means an overflow on the final result
588 (due to the multiplication by 2, which has not been done
589 yet).
590 3. Ziv test.
591 4. Correct the sign from the sign of sin(...).
592 5. Round then multiply by 2. Here, an overflow in either
593 operation means a real overflow. */
594 mpfr_reflection_overflow (z_pre, s1, s, y, p, MPFR_RNDD);
595 /* z_pre is a lower bound of |zeta(s)|/2, thus if it overflows,
596 or has exponent emax, then |zeta(s)| overflows too. */
597 if (MPFR_IS_INF (z_pre) || MPFR_GET_EXP(z_pre) == __gmpfr_emax)
598 { /* determine the sign of overflow */
599 mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */
600 mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */
601 overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1;
602 break;
603 }
604 else /* EXP(z_pre) < __gmpfr_emax */
605 {
606 int ok = 0;
607 mpfr_t z_down;
608 mpfr_init2 (z_up, mpfr_get_prec (z_pre));
609 mpfr_reflection_overflow (z_up, s1, s, y, p, MPFR_RNDU);
610 /* if the lower approximation z_pre does not overflow, but
611 z_up does, we need more precision */
612 if (MPFR_IS_INF (z_up) || MPFR_GET_EXP(z_up) == __gmpfr_emax)
613 goto next_loop;
614 /* check if z_pre and z_up round to the same number */
615 mpfr_init2 (z_down, precz);
616 mpfr_set (z_down, z_pre, rnd_mode);
617 /* Note: it might be that EXP(z_down) = emax here, in that
618 case we will have overflow below when we multiply by 2 */
619 mpfr_prec_round (z_up, precz, rnd_mode);
620 ok = mpfr_equal_p (z_down, z_up);
621 mpfr_clear (z_up);
622 mpfr_clear (z_down);
623 if (ok)
624 {
625 /* get correct sign and multiply by 2 */
626 mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */
627 mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */
628 if (mpfr_cmp_si_2exp (s1, -1, -1) > 0)
629 mpfr_neg (z_pre, z_pre, rnd_mode);
630 mpfr_mul_2ui (z_pre, z_pre, 1, rnd_mode);
631 break;
632 }
633 else
634 goto next_loop;
635 }
636 }
637 mpfr_zeta_pos (z_pre, s1, MPFR_RNDN); /* zeta(1-s) */
638 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); /* gamma(1-s)*zeta(1-s) */
639
640 /* multiply z_pre by 2^s*Pi^(s-1) where p=Pi, s1=1-s */
641 mpfr_mul_2ui (y, p, 1, MPFR_RNDN); /* 2*Pi */
642 mpfr_neg (s1, s1, MPFR_RNDN); /* s-1 */
643 mpfr_pow (y, y, s1, MPFR_RNDN); /* (2*Pi)^(s-1) */
644 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
645 mpfr_mul_2ui (z_pre, z_pre, 1, MPFR_RNDN);
646
647 /* multiply z_pre by sin(Pi*s/2) */
648 mpfr_div_2ui (p, s, 1, MPFR_RNDN); /* p = s/2 */
649 /* Can mpfr_sinpi underflow? While with mpfr_sin, we could not
650 answer in any precision without a theoretical study (though
651 an underflow would have been very unlikely as we have a
652 huge exponent range), with mpfr_sinpi, an underflow could
653 occur only in a huge, unsupported precision. Indeed, if
654 mpfr_sinpi underflows, this means that 0 < |sinpi(s/2)| < m,
655 where m is the minimum representable positive number, and in
656 this case, r being the reduced argument such that |r| <= 1/2,
657 one has |sinpi(r)| > |2r|, so that |2r| < m; this can be
658 possible only if |s/2| > 1/2 (otherwise |s| = |2r| < m and
659 s would not be representable as an MPFR number) and s has
660 non-zero bits of exponent less than the minimum exponent
661 (s/2 - r being an integer), i.e. the precision is at least
662 MPFR_EMAX_MAX + 2. With such a huge precision, there would
663 probably be failures before reaching this point. */
664 mpfr_sinpi (y, p, MPFR_RNDN); /* y = sin(Pi*s/2) */
665 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
666
667 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz,
668 rnd_mode)))
669 break;
670
671 next_loop:
672 MPFR_ZIV_NEXT (loop, prec1);
673 MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p);
674 }
675 MPFR_ZIV_FREE (loop);
676 if (overflow != 0)
677 {
678 inex = mpfr_overflow (z, rnd_mode, overflow);
679 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
680 }
681 else
682 inex = mpfr_set (z, z_pre, rnd_mode);
683 MPFR_GROUP_CLEAR (group);
684 }
685
686 MPFR_SAVE_EXPO_FREE (expo);
687 return mpfr_check_range (z, inex, rnd_mode);
688 }