1 /* mpfr_zeta_ui -- compute the Riemann Zeta function for integer argument.
2
3 Copyright 2005-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 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 int
mpfr_zeta_ui(mpfr_ptr z,unsigned long m,mpfr_rnd_t r)27 mpfr_zeta_ui (mpfr_ptr z, unsigned long m, mpfr_rnd_t r)
28 {
29 MPFR_ZIV_DECL (loop);
30
31 MPFR_LOG_FUNC
32 (("m=%lu rnd=%d prec=%Pd", m, r, mpfr_get_prec (z)),
33 ("z[%Pd]=%.*Rg", mpfr_get_prec (z), mpfr_log_prec, z));
34
35 if (m == 0) /* zeta(0) = -1/2 */
36 return mpfr_set_si_2exp (z, -1, -1, r);
37 else if (m == 1)
38 {
39 MPFR_SET_INF (z);
40 MPFR_SET_POS (z);
41 MPFR_SET_DIVBY0 ();
42 return 0;
43 }
44 else /* m >= 2 */
45 {
46 mpfr_prec_t p = MPFR_PREC(z);
47 unsigned long n, k, err, kbits;
48 mpz_t d, t, s, q;
49 mpfr_t y;
50 int inex;
51 MPFR_SAVE_EXPO_DECL (expo);
52
53 if (r == MPFR_RNDA)
54 r = MPFR_RNDU; /* since the result is always positive */
55
56 MPFR_SAVE_EXPO_MARK (expo);
57
58 if (m >= p) /* 2^(-m) < ulp(1) = 2^(1-p). This means that
59 2^(-m) <= 1/2*ulp(1). We have 3^(-m)+4^(-m)+... < 2^(-m)
60 i.e. zeta(m) < 1+2*2^(-m) for m >= 3 */
61 {
62 if (m == 2) /* necessarily p=2 */
63 inex = mpfr_set_ui_2exp (z, 13, -3, r);
64 else if (r == MPFR_RNDZ || r == MPFR_RNDD ||
65 (r == MPFR_RNDN && m > p))
66 {
67 mpfr_set_ui (z, 1, r);
68 inex = -1;
69 }
70 else
71 {
72 mpfr_set_ui (z, 1, r);
73 mpfr_nextabove (z);
74 inex = 1;
75 }
76 goto end;
77 }
78
79 /* now treat also the case where zeta(m) - (1+1/2^m) < 1/2*ulp(1),
80 and the result is either 1+2^(-m) or 1+2^(-m)+2^(1-p). */
81 mpfr_init2 (y, 31);
82
83 if (m >= p / 2) /* otherwise 4^(-m) > 2^(-p) */
84 {
85 /* the following is a lower bound for log(3)/log(2) */
86 mpfr_set_str_binary (y, "1.100101011100000000011010001110");
87 mpfr_mul_ui (y, y, m, MPFR_RNDZ); /* lower bound for log2(3^m) */
88 if (mpfr_cmp_ui (y, p + 2) >= 0)
89 {
90 mpfr_clear (y);
91 mpfr_set_ui (z, 1, MPFR_RNDZ);
92 mpfr_div_2ui (z, z, m, MPFR_RNDZ);
93 mpfr_add_ui (z, z, 1, MPFR_RNDZ);
94 if (r != MPFR_RNDU)
95 inex = -1;
96 else
97 {
98 mpfr_nextabove (z);
99 inex = 1;
100 }
101 goto end;
102 }
103 }
104
105 mpz_init (s);
106 mpz_init (d);
107 mpz_init (t);
108 mpz_init (q);
109
110 p += MPFR_INT_CEIL_LOG2(p); /* account of the n term in the error */
111
112 p += MPFR_INT_CEIL_LOG2(p) + 15; /* initial value */
113
114 MPFR_ZIV_INIT (loop, p);
115 for(;;)
116 {
117 /* 0.39321985067869744 = log(2)/log(3+sqrt(8)) */
118 n = 1 + (unsigned long) (0.39321985067869744 * (double) p);
119 err = n + 4;
120
121 mpfr_set_prec (y, p);
122
123 /* computation of the d[k] */
124 mpz_set_ui (s, 0);
125 mpz_set_ui (t, 1);
126 mpz_mul_2exp (t, t, 2 * n - 1); /* t[n] */
127 mpz_set (d, t);
128 for (k = n; k > 0; k--)
129 {
130 count_leading_zeros (kbits, k);
131 kbits = GMP_NUMB_BITS - kbits;
132 /* if k^m is too large, use mpz_tdiv_q */
133 if (m * kbits > 2 * GMP_NUMB_BITS)
134 {
135 /* if we know in advance that k^m > d, then floor(d/k^m) will
136 be zero below, so there is no need to compute k^m */
137 kbits = (kbits - 1) * m + 1;
138 /* k^m has at least kbits bits */
139 if (kbits > mpz_sizeinbase (d, 2))
140 mpz_set_ui (q, 0);
141 else
142 {
143 mpz_ui_pow_ui (q, k, m);
144 mpz_tdiv_q (q, d, q);
145 }
146 }
147 else /* use several mpz_tdiv_q_ui calls */
148 {
149 unsigned long km = k, mm = m - 1;
150 while (mm > 0 && km < ULONG_MAX / k)
151 {
152 km *= k;
153 mm --;
154 }
155 mpz_tdiv_q_ui (q, d, km);
156 while (mm > 0)
157 {
158 km = k;
159 mm --;
160 while (mm > 0 && km < ULONG_MAX / k)
161 {
162 km *= k;
163 mm --;
164 }
165 mpz_tdiv_q_ui (q, q, km);
166 }
167 }
168 if (k % 2)
169 mpz_add (s, s, q);
170 else
171 mpz_sub (s, s, q);
172
173 /* we have d[k] = sum(t[i], i=k+1..n)
174 with t[i] = n*(n+i-1)!*4^i/(n-i)!/(2i)!
175 t[k-1]/t[k] = k*(2k-1)/(n-k+1)/(n+k-1)/2 */
176 #if (GMP_NUMB_BITS == 32)
177 #define KMAX 46341 /* max k such that k*(2k-1) < 2^32 */
178 #elif (GMP_NUMB_BITS == 64)
179 #define KMAX 3037000500
180 #endif
181 #ifdef KMAX
182 if (k <= KMAX)
183 mpz_mul_ui (t, t, k * (2 * k - 1));
184 else
185 #endif
186 {
187 mpz_mul_ui (t, t, k);
188 mpz_mul_ui (t, t, 2 * k - 1);
189 }
190 mpz_fdiv_q_2exp (t, t, 1);
191 /* Warning: the test below assumes that an unsigned long
192 has no padding bits. */
193 if (n < 1UL << ((sizeof(unsigned long) * CHAR_BIT) / 2))
194 /* (n - k + 1) * (n + k - 1) < n^2 */
195 mpz_divexact_ui (t, t, (n - k + 1) * (n + k - 1));
196 else
197 {
198 mpz_divexact_ui (t, t, n - k + 1);
199 mpz_divexact_ui (t, t, n + k - 1);
200 }
201 mpz_add (d, d, t);
202 }
203
204 /* multiply by 1/(1-2^(1-m)) = 1 + 2^(1-m) + 2^(2-m) + ... */
205 mpz_fdiv_q_2exp (t, s, m - 1);
206 do
207 {
208 err ++;
209 mpz_add (s, s, t);
210 mpz_fdiv_q_2exp (t, t, m - 1);
211 }
212 while (mpz_cmp_ui (t, 0) > 0);
213
214 /* divide by d[n] */
215 mpz_mul_2exp (s, s, p);
216 mpz_tdiv_q (s, s, d);
217 mpfr_set_z (y, s, MPFR_RNDN);
218 mpfr_div_2ui (y, y, p, MPFR_RNDN);
219
220 err = MPFR_INT_CEIL_LOG2 (err);
221
222 if (MPFR_LIKELY(MPFR_CAN_ROUND (y, p - err, MPFR_PREC(z), r)))
223 break;
224
225 MPFR_ZIV_NEXT (loop, p);
226 }
227 MPFR_ZIV_FREE (loop);
228
229 mpz_clear (d);
230 mpz_clear (t);
231 mpz_clear (q);
232 mpz_clear (s);
233 inex = mpfr_set (z, y, r);
234 mpfr_clear (y);
235
236 end:
237 MPFR_LOG_VAR (z);
238 MPFR_LOG_MSG (("inex = %d before mpfr_check_range\n", inex));
239 MPFR_SAVE_EXPO_FREE (expo);
240 return mpfr_check_range (z, inex, r);
241 }
242 }
243