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