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