Line data Source code
1 : /* Copyright (C) 2018 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_bern
19 :
20 : /********************************************************************/
21 : /** **/
22 : /** BERNOULLI NUMBERS B_2k **/
23 : /** **/
24 : /********************************************************************/
25 :
26 : /* D = divisorsu(n). Return a/b = \sum_{p-1 | 2n: p prime} 1/p
27 : * B_2k + a/b in Z [Clausen-von Staudt] */
28 : static GEN
29 82523 : fracB2k(GEN D)
30 : {
31 82523 : GEN a = utoipos(5), b = utoipos(6); /* 1/2 + 1/3 */
32 82685 : long i, l = lg(D);
33 622908 : for (i = 2; i < l; i++) /* skip 1 */
34 : {
35 540462 : ulong p = 2*D[i] + 1; /* a/b += 1/p */
36 540462 : if (uisprime(p)) { a = addii(muliu(a,p), b); b = muliu(b,p); }
37 : }
38 82446 : return mkfrac(a,b);
39 : }
40 : /* precision needed to compute B_k for all k <= N */
41 : long
42 1888 : bernbitprec(long N)
43 : { /* 1.612086 ~ log(8Pi) / 2 */
44 1888 : const double log2PI = 1.83787706641;
45 1888 : double logN = log((double)N);
46 1888 : double t = (N + 4) * logN - N*(1 + log2PI) + 1.612086;
47 1888 : return (long)ceil(t / M_LN2) + 10;
48 : }
49 : static long
50 21 : bernprec(long N) { return nbits2prec(bernbitprec(N)); }
51 : /* \sum_{k > M} k^(-n) <= M^(1-n) / (n-1) < 2^-b */
52 : static long
53 1391 : zetamaxpow(long n)
54 : {
55 1391 : long M = (long)ceil(n / (2 * M_PI * M_E));
56 1391 : return M | 1; /* make it odd */
57 : }
58 : /* v * zeta(k) using r precomputed odd powers */
59 : static GEN
60 82513 : bern_zeta(GEN v, long k, GEN pow, long r, long p)
61 : {
62 82513 : GEN z, s = gel(pow, r);
63 : long j;
64 10364907 : for (j = r - 2; j >= 3; j -= 2) s = addii(s, gel(pow,j));
65 81816 : z = s = itor(s, nbits2prec(p));
66 82524 : shiftr_inplace(s, -p); /* zeta(k)(1 - 2^(-k)) - 1*/
67 82520 : s = addrs(s, 1); shiftr_inplace(s, -k);
68 : /* divide by 1 - 2^(-k): s + s/2^k + s/2^(2k) + ... */
69 364724 : for (; k < p; k <<= 1) s = addrr(s, shiftr(s, -k));
70 82483 : return addrr(v, mulrr(v, addrr(z, s)));
71 : }
72 : /* z * j^2 */
73 : static GEN
74 50921690 : muliu2(GEN z, ulong j)
75 50921690 : { return (j | HIGHMASK)? mulii(z, sqru(j)): muliu(z, j*j); }
76 : /* 1 <= m <= n, set y[1] = B_{2m}, ... y[n-m+1] = B_{2n} in Q */
77 : static void
78 304 : bernset(GEN *y, long m, long n)
79 : {
80 304 : long i, j, k, p, prec, r, N = n << 1; /* up to B_N */
81 : GEN u, b, v, t;
82 304 : p = bernbitprec(N); prec = nbits2prec(p);
83 304 : u = sqrr(Pi2n(1, prec)); /* (2Pi)^2 */
84 304 : v = divrr(mpfactr(N, prec), powru(u, n)); shiftr_inplace(v,1);
85 304 : r = zetamaxpow(N);
86 304 : t = cgetg(r+1, t_VEC); b = int2n(p); /* fixed point */
87 5630 : for (j = 3; j <= r; j += 2)
88 : {
89 5326 : GEN z = cgeti(prec);
90 5326 : pari_sp av2 = avma;
91 5326 : affii(divii(b, powuu(j, N)), z);
92 5326 : gel(t,j) = z; set_avma(av2);
93 : }
94 304 : y += n - m;
95 304 : for (i = n, k = N;; i--)
96 82181 : { /* set B_n, k = 2i */
97 82485 : pari_sp av2 = avma;
98 82485 : GEN z = fracB2k(divisorsu(i)), B = bern_zeta(v, k, t, r, p);
99 : long j;
100 : /* B = v * zeta(k), v = 2*k! / (2Pi)^k */
101 82469 : if (!odd(i)) setsigne(B, -1); /* B ~ B_n */
102 82470 : B = roundr(addrr(B, fractor(z,LOWDEFAULTPREC))); /* B - z = B_n */
103 82518 : *y-- = gclone(gsub(B, z));
104 82532 : if (i == m) break;
105 82228 : affrr(divrunextu(mulrr(v,u), k-1), v);
106 10443601 : for (j = r; j >= 3; j -= 2) affii(muliu2(gel(t,j), j), gel(t,j));
107 82214 : set_avma(av2); k -= 2;
108 82187 : if (((N - k) & 0x7f) == 0x7e)
109 : { /* reduce precision if possible */
110 1087 : long p2 = p, prec2 = prec;
111 1087 : p = bernbitprec(k); prec = nbits2prec(p); if (prec2 == prec) continue;
112 1087 : setprec(v, prec); r = zetamaxpow(k);
113 158417 : for (j = 3; j <= r; j += 2) affii(shifti(gel(t,j), p - p2), gel(t,j));
114 1087 : set_avma(av2);
115 : }
116 : }
117 304 : }
118 : /* need B[2..2*nb] as t_INT or t_FRAC */
119 : void
120 234813 : constbern(long nb)
121 : {
122 234813 : const pari_sp av = avma;
123 : long i, l;
124 : GEN B;
125 : pari_timer T;
126 :
127 234813 : l = bernzone? lg(bernzone): 0;
128 234813 : if (l > nb) return;
129 :
130 304 : nb = maxss(nb, l + 127);
131 304 : B = cgetg_block(nb+1, t_VEC);
132 304 : if (bernzone)
133 7040 : { for (i = 1; i < l; i++) gel(B,i) = gel(bernzone,i); }
134 : else
135 : {
136 256 : gel(B,1) = gclone(mkfracss(1,6));
137 256 : gel(B,2) = gclone(mkfracss(-1,30));
138 256 : gel(B,3) = gclone(mkfracss(1,42));
139 256 : gel(B,4) = gclone(mkfracss(-1,30));
140 256 : gel(B,5) = gclone(mkfracss(5,66));
141 256 : gel(B,6) = gclone(mkfracss(-691,2730));
142 256 : gel(B,7) = gclone(mkfracss(7,6));
143 256 : gel(B,8) = gclone(mkfracss(-3617,510));
144 256 : gel(B,9) = gclone(mkfracss(43867,798));
145 256 : gel(B,10)= gclone(mkfracss(-174611,330));
146 256 : gel(B,11)= gclone(mkfracss(854513,138));
147 256 : gel(B,12)= gclone(mkfracss(-236364091,2730));
148 256 : gel(B,13)= gclone(mkfracss(8553103,6)); /* B_26 */
149 256 : l = 14;
150 : }
151 304 : set_avma(av);
152 304 : if (DEBUGLEVEL) {
153 0 : err_printf("caching Bernoulli numbers 2*%ld to 2*%ld\n", l, nb);
154 0 : timer_start(&T);
155 : }
156 304 : bernset((GEN*)B + l, l, nb);
157 304 : if (DEBUGLEVEL) timer_printf(&T, "Bernoulli");
158 304 : swap(B, bernzone); guncloneNULL(B);
159 304 : set_avma(av);
160 : #if 0
161 : if (nb > 200000)
162 : #endif
163 : {
164 304 : const ulong p = 4294967291UL;
165 304 : long n = 2 * nb + 2;
166 304 : GEN t = const_vecsmall(n+1, 1);
167 304 : t[1] = evalvarn(0); t[2] = 0;
168 304 : t = Flx_shift(Flx_invLaplace(t, p), -1); /* t = (exp(x)-1)/x */
169 304 : t = Flx_Laplace(Flxn_inv(t, n, p), p);
170 93162 : for (i = 1; i <= nb; i++)
171 92858 : if (Rg_to_Fl(bernfrac(2*i), p) != uel(t,2*i+2))
172 : {
173 0 : gunclone(bernzone); bernzone = NULL;
174 0 : pari_err_BUG(stack_sprintf("B_{2*%ld}", i));
175 : }
176 304 : set_avma(av);
177 : }
178 : }
179 : /* Obsolete, kept for backward compatibility */
180 : void
181 0 : mpbern(long n, long prec) { (void)prec; constbern(n); }
182 :
183 : /* assume n even > 0, if iz != NULL, assume iz = 1/zeta(n) */
184 : static GEN
185 21 : bernreal_using_zeta(long n, long prec)
186 : {
187 21 : GEN pi2 = Pi2n(1, prec+EXTRAPRECWORD);
188 21 : GEN iz = inv_szeta_euler(n, prec);
189 21 : GEN z = divrr(mpfactr(n, prec), mulrr(powru(pi2, n), iz));
190 21 : shiftr_inplace(z, 1); /* 2 * n! * zeta(n) / (2Pi)^n */
191 21 : if ((n & 3) == 0) setsigne(z, -1);
192 21 : return z;
193 : }
194 : /* assume n even > 0, B = NULL or good approximation to B_n */
195 : static GEN
196 14 : bernfrac_i(long n, GEN B)
197 : {
198 14 : GEN z = fracB2k(divisorsu(n >> 1));
199 14 : if (!B) B = bernreal_using_zeta(n, bernprec(n));
200 14 : B = roundr( gadd(B, fractor(z,LOWDEFAULTPREC)) );
201 14 : return gsub(B, z);
202 : }
203 : GEN
204 9178165 : bernfrac(long n)
205 : {
206 : pari_sp av;
207 : long k;
208 9178165 : if (n <= 1)
209 : {
210 4354 : if (n < 0) pari_err_DOMAIN("bernfrac", "index", "<", gen_0, stoi(n));
211 4347 : return n? mkfrac(gen_m1,gen_2): gen_1;
212 : }
213 9173811 : if (odd(n)) return gen_0;
214 9170034 : k = n >> 1;
215 9170034 : if (!bernzone) constbern(0);
216 9170034 : if (bernzone && k < lg(bernzone)) return gel(bernzone, k);
217 21 : av = avma;
218 21 : return gerepileupto(av, bernfrac_i(n, NULL));
219 : }
220 : GEN
221 77 : bernvec(long n)
222 : {
223 : long i, l;
224 : GEN y;
225 77 : if (n < 0) return cgetg(1, t_VEC);
226 77 : constbern(n);
227 77 : l = n+2; y = cgetg(l, t_VEC); gel(y,1) = gen_1;
228 53536 : for (i = 2; i < l; i++) gel(y,i) = gel(bernzone,i-1);
229 77 : return y;
230 : }
231 :
232 : /* x := pol_x(v); B_k(x) = \sum_{i=0}^k binomial(k, i) B_i x^{k-i} */
233 : static GEN
234 2156 : bernpol_i(long k, long v)
235 : {
236 : GEN B, C;
237 : long i;
238 2156 : if (v < 0) v = 0;
239 2156 : constbern(k >> 1); /* cache B_2, ..., B_2[k/2] */
240 2156 : C = vecbinomial(k);
241 2156 : B = cgetg(k + 3, t_POL);
242 14777 : for (i = 0; i <= k; ++i) gel(B, k-i+2) = gmul(gel(C,i+1), bernfrac(i));
243 2156 : B[1] = evalsigne(1) | evalvarn(v);
244 2156 : return B;
245 : }
246 : GEN
247 2086 : bernpol(long k, long v)
248 : {
249 2086 : pari_sp av = avma;
250 2086 : if (k < 0) pari_err_DOMAIN("bernpol", "index", "<", gen_0, stoi(k));
251 2079 : return gerepileupto(av, bernpol_i(k, v));
252 : }
253 : /* x := pol_x(v); return 1^e + ... + x^e = x^e + (B_{e+1}(x) - B_{e+1})/(e+1) */
254 : static GEN
255 49 : faulhaber(long e, long v)
256 : {
257 : GEN B;
258 49 : if (e == 0) return pol_x(v);
259 35 : B = RgX_integ(bernpol_i(e, v)); /* (B_{e+1}(x) - B_{e+1}) / (e+1) */
260 35 : gel(B,e+2) = gaddgs(gel(B,e+2), 1); /* add x^e, in place */
261 35 : return B;
262 : }
263 : /* sum_v T(v), T a polynomial expression in v */
264 : GEN
265 49 : sumformal(GEN T, long v)
266 : {
267 49 : pari_sp av = avma, av2;
268 : long i, t, d;
269 : GEN R;
270 :
271 49 : T = simplify_shallow(T);
272 49 : t = typ(T);
273 49 : if (is_scalar_t(t))
274 14 : return gerepileupto(av, monomialcopy(T, 1, v < 0? 0: v));
275 35 : if (t != t_POL) pari_err_TYPE("sumformal [not a t_POL]", T);
276 28 : if (v < 0) v = varn(T);
277 28 : av2 = avma;
278 28 : R = gen_0;
279 28 : d = poldegree(T,v);
280 91 : for (i = d; i >= 0; i--)
281 : {
282 63 : GEN c = polcoef_i(T, i, v);
283 63 : if (gequal0(c)) continue;
284 42 : R = gadd(R, gmul(c, faulhaber(i, v)));
285 42 : if (gc_needed(av2,3))
286 : {
287 0 : if(DEBUGMEM>1) pari_warn(warnmem,"sumformal, i = %ld/%ld", i,d);
288 0 : R = gerepileupto(av2, R);
289 : }
290 : }
291 28 : return gerepileupto(av, R);
292 : }
293 :
294 : /* 1/zeta(n) using Euler product. Assume n > 0. */
295 : GEN
296 1194 : inv_szeta_euler(long n, long prec)
297 : {
298 1194 : long bit = prec2nbits(prec);
299 : GEN z, res;
300 : pari_sp av, av2;
301 : double A, D, lba;
302 : ulong p, lim;
303 : forprime_t S;
304 :
305 1194 : if (n > bit) return real_1(prec);
306 :
307 1187 : lba = prec2nbits_mul(prec, M_LN2);
308 1187 : D = exp((lba - log((double)(n-1))) / (n-1));
309 1187 : lim = 1 + (ulong)ceil(D);
310 1187 : if (lim < 3) return subir(gen_1,real2n(-n,prec));
311 1187 : res = cgetr(prec); av = avma; incrprec(prec);
312 :
313 1187 : (void)u_forprime_init(&S, 3, lim);
314 1187 : av2 = avma; A = n / M_LN2; z = subir(gen_1, real2n(-n, prec));
315 86205 : while ((p = u_forprime_next(&S)))
316 : {
317 85018 : long l = bit - (long)floor(A * log((double)p));
318 : GEN h;
319 :
320 85018 : if (l < BITS_IN_LONG) l = BITS_IN_LONG;
321 85018 : l = minss(prec, nbits2prec(l));
322 85018 : h = divrr(z, rpowuu(p, (ulong)n, l));
323 85018 : z = subrr(z, h);
324 85018 : if (gc_needed(av,1))
325 : {
326 0 : if (DEBUGMEM>1) pari_warn(warnmem,"inv_szeta_euler, p = %lu/%lu", p,lim);
327 0 : z = gerepileuptoleaf(av2, z);
328 : }
329 : }
330 1187 : affrr(z, res); set_avma(av); return res;
331 : }
332 :
333 : /* Return B_n */
334 : GEN
335 7917 : bernreal(long n, long prec)
336 : {
337 : pari_sp av;
338 : GEN B;
339 : long p, k;
340 7917 : if (n < 0) pari_err_DOMAIN("bernreal", "index", "<", gen_0, stoi(n));
341 7910 : if (n == 0) return real_1(prec);
342 7910 : if (n == 1) return real_m2n(-1,prec); /*-1/2*/
343 7910 : if (odd(n)) return real_0(prec);
344 :
345 7910 : k = n >> 1;
346 7910 : if (!bernzone) constbern(0);
347 7910 : if (k < lg(bernzone)) return fractor(gel(bernzone,k), prec);
348 7 : p = bernprec(n); av = avma;
349 7 : B = bernreal_using_zeta(n, minss(p, prec));
350 7 : if (p < prec) B = fractor(bernfrac_i(n, B), prec);
351 7 : return gerepileuptoleaf(av, B);
352 : }
353 :
354 : GEN
355 49 : eulerpol(long k, long v)
356 : {
357 49 : pari_sp av = avma;
358 : GEN B, E;
359 49 : if (k < 0) pari_err_DOMAIN("eulerpol", "index", "<", gen_0, stoi(k));
360 42 : k++; B = bernpol_i(k, v);
361 42 : E = RgX_Rg_mul(RgX_sub(B, RgX_rescale(B, gen_2)), uutoQ(2,k));
362 42 : return gerepileupto(av, E);
363 : }
364 :
365 : /*******************************************************************/
366 : /** HARMONIC NUMBERS **/
367 : /*******************************************************************/
368 : /* 1/a + ... + 1/(b-1); a < b <= 2^(BIL-1) */
369 : static GEN
370 4676 : hrec(ulong a, ulong b)
371 : {
372 : ulong m;
373 4676 : switch(b - a)
374 : {
375 4501 : case 1: retmkfrac(gen_1, utoipos(a));
376 133 : case 2: if (a < 65536) retmkfrac(utoipos(2*a + 1), utoipos(a * a + a));
377 0 : retmkfrac(utoipos(2*a + 1), muluu(a, a+1));
378 : }
379 42 : m = (a + b) >> 1;
380 42 : return gadd(hrec(a, m), hrec(m, b));
381 : }
382 : /* exact Harmonic number H_n, n < 2^(BIL-1).
383 : * Could use H_n = sum_k 2^(-k) H^odd_{n \ 2^k} */
384 : GEN
385 4592 : harmonic(ulong n)
386 : {
387 4592 : pari_sp av = avma;
388 4592 : return n? gerepileupto(av, hrec(1, n+1)): gen_0;
389 : }
390 :
391 : /* 1/a^k + ... + 1/(b-1)^k; a < b */
392 : static GEN
393 77 : hreck(ulong a, ulong b, ulong k)
394 : {
395 : ulong m;
396 77 : switch(b - a)
397 : {
398 : GEN x, y;
399 14 : case 1: retmkfrac(gen_1, powuu(a, k));
400 28 : case 2:
401 28 : x = powuu(a, k); y = powuu(a + 1, k);
402 28 : retmkfrac(addii(x, y), mulii(x, y));
403 : }
404 35 : m = (a + b) >> 1;
405 35 : return gadd(hreck(a, m, k), hreck(m, b, k));
406 : }
407 : GEN
408 56 : harmonic0(ulong n, GEN k)
409 : {
410 56 : pari_sp av = avma;
411 : ulong r;
412 56 : if (!n) return gen_0;
413 49 : if (n & HIGHBIT) pari_err_OVERFLOW("harmonic");
414 49 : if (!k) return harmonic(n);
415 21 : if (typ(k) != t_INT) pari_err_TYPE("harmonic", k);
416 14 : if (signe(k) < 0)
417 : {
418 7 : GEN H = poleval(faulhaber(-itos(k), 0), utoipos(n));
419 7 : return gerepileuptoint(av, H);
420 : }
421 7 : r = itou(k);
422 7 : if (!r) return utoipos(n);
423 7 : if (r == 1) return harmonic(n);
424 7 : return gerepileupto(av, hreck(1, n+1, r));
425 : }
426 :
427 :
428 : /**************************************************************/
429 : /* Euler numbers */
430 : /**************************************************************/
431 :
432 : /* precision needed to compute E_k for all k <= N */
433 : static long
434 798 : eulerbitprec(long N)
435 : { /* 1.1605 ~ log(32/Pi) / 2 */
436 798 : const double logPIS2 = 0.4515827;
437 798 : double t = (N + 1) * log((double)N) - N*(1 + logPIS2) + 1.1605;
438 798 : return (long)ceil(t / M_LN2) + 10;
439 : }
440 : static long
441 14 : eulerprec(long N) { return nbits2prec(eulerbitprec(N)); }
442 :
443 : /* sum_{k > M, k odd} (-1)^((k-1)/2)k^(-n) < M^(-n) < 2^-b */
444 : static long
445 784 : lfun4maxpow(long n)
446 : {
447 784 : long M = (long)ceil(2 * n / (M_E * M_PI));
448 784 : return M | 1; /* make it odd */
449 : }
450 :
451 : /* lfun4(k) using r precomputed odd powers */
452 : static GEN
453 49791 : euler_lfun4(GEN v, GEN pow, long r, long p)
454 : {
455 49791 : GEN s = ((r & 3L) == 1)? gel(pow, r): negi(gel(pow, r));
456 : long j;
457 40556894 : for (j = r - 2; j >= 3; j -= 2)
458 40507103 : s = ((j & 3L) == 1)? addii(s, gel(pow,j)): subii(s, gel(pow,j));
459 49791 : s = mulri(v, s); shiftr_inplace(s, -p);
460 49791 : return addrr(v, s);
461 : }
462 :
463 : /* 1 <= m <= n, set y[1] = E_{2m}, ... y[n-m+1] = E_{2n} in Z */
464 : static void
465 21 : eulerset(GEN *y, long m, long n)
466 : {
467 21 : long i, j, k, p, prec, r, N = n << 1, N1 = N + 1; /* up to E_N */
468 : GEN b, u, v, t;
469 21 : p = eulerbitprec(N); prec = nbits2prec(p);
470 21 : u = sqrr(Pi2n(-1, prec)); /* (Pi/2)^2 */
471 21 : v = divrr(mpfactr(N, prec), mulrr(powru(u, n), Pi2n(-2,prec)));
472 21 : r = lfun4maxpow(N1);
473 21 : t = cgetg(r+1, t_VEC); b = int2n(p); /* fixed point */
474 11921 : for (j = 3; j <= r; j += 2)
475 : {
476 11900 : GEN z = cgeti(prec);
477 11900 : pari_sp av2 = avma;
478 11900 : affii(divii(b, powuu(j, N+1)), z);
479 11900 : gel(t,j) = z; set_avma(av2);
480 : }
481 21 : y += n - m;
482 21 : for (i = n, k = N1;; i--)
483 49770 : { /* set E_n, k = 2i + 1 */
484 49791 : pari_sp av2 = avma;
485 49791 : GEN E = euler_lfun4(v, t, r, p);
486 : long j;
487 : /* E = v * lfun4(k), v = (4/Pi)*k! / (Pi/2)^k */
488 49791 : E = roundr(E); if (odd(i)) setsigne(E, -1); /* E ~ E_n */
489 49791 : *y-- = gclone(E);
490 49791 : if (i == m) break;
491 49770 : affrr(divrunextu(mulrr(v,u), k-2), v);
492 40606202 : for (j = r; j >= 3; j -= 2) affii(muliu2(gel(t,j), j), gel(t,j));
493 49770 : set_avma(av2); k -= 2;
494 49770 : if (((N1 - k) & 0x7f) == 0x7e)
495 : { /* reduce precision if possible */
496 763 : long p2 = p, prec2 = prec;
497 763 : p = eulerbitprec(k); prec = nbits2prec(p); if (prec2 == prec) continue;
498 763 : setprec(v, prec); r = lfun4maxpow(k);
499 622923 : for (j = 3; j <= r; j += 2) affii(shifti(gel(t,j), p - p2), gel(t,j));
500 763 : set_avma(av2);
501 : }
502 : }
503 21 : }
504 :
505 : /* need E[2..2*nb] as t_INT */
506 : static void
507 28 : constreuler(long nb)
508 : {
509 28 : const pari_sp av = avma;
510 : long i, l;
511 : GEN E;
512 : pari_timer T;
513 :
514 28 : l = eulerzone? lg(eulerzone): 0;
515 28 : if (l > nb) return;
516 :
517 21 : nb = maxss(nb, l + 127);
518 21 : E = cgetg_block(nb+1, t_VEC);
519 21 : if (eulerzone)
520 896 : { for (i = 1; i < l; i++) gel(E,i) = gel(eulerzone,i); }
521 : else
522 : {
523 14 : gel(E,1) = gclone(stoi(-1));
524 14 : gel(E,2) = gclone(stoi(5));
525 14 : gel(E,3) = gclone(stoi(-61));
526 14 : gel(E,4) = gclone(stoi(1385));
527 14 : gel(E,5) = gclone(stoi(-50521));
528 14 : gel(E,6) = gclone(stoi(2702765));
529 14 : gel(E,7) = gclone(stoi(-199360981));
530 14 : l = 8;
531 : }
532 21 : set_avma(av);
533 21 : if (DEBUGLEVEL) {
534 0 : err_printf("caching Euler numbers 2*%ld to 2*%ld\n", l, nb);
535 0 : timer_start(&T);
536 : }
537 21 : eulerset((GEN*)E + l, l, nb);
538 21 : if (DEBUGLEVEL) timer_printf(&T, "Euler");
539 21 : swap(E, eulerzone); guncloneNULL(E);
540 21 : set_avma(av);
541 : }
542 :
543 : /* 1/lfun(-4,n) using Euler product. Assume n > 0. */
544 : static GEN
545 14 : inv_lfun4(long n, long prec)
546 : {
547 14 : long bit = prec2nbits(prec);
548 : GEN z, res;
549 : pari_sp av, av2;
550 : double A;
551 : ulong p, lim;
552 : forprime_t S;
553 :
554 14 : if (n > bit) return real_1(prec);
555 :
556 14 : lim = 1 + (ulong)ceil(exp2((double)bit / n));
557 14 : res = cgetr(prec); av = avma; incrprec(prec);
558 :
559 14 : (void)u_forprime_init(&S, 3, lim);
560 14 : av2 = avma; A = n / M_LN2; z = real_1(prec);
561 369 : while ((p = u_forprime_next(&S)))
562 : {
563 355 : long l = bit - (long)floor(A * log((double)p));
564 : GEN h;
565 :
566 355 : if (l < BITS_IN_LONG) l = BITS_IN_LONG;
567 355 : l = minss(prec, nbits2prec(l));
568 355 : h = rpowuu(p, (ulong)n, l); if ((p & 3UL) == 1) setsigne(h, -1);
569 355 : z = addrr(z, divrr(z, h)); /* z *= 1 - chi_{-4}(p) / p^n */
570 355 : if (gc_needed(av,1))
571 : {
572 0 : if (DEBUGMEM>1) pari_warn(warnmem,"inv_lfun4, p = %lu/%lu", p,lim);
573 0 : z = gerepileuptoleaf(av2, z);
574 : }
575 : }
576 14 : affrr(z, res); set_avma(av); return res;
577 : }
578 : /* assume n even > 0, E_n = (-1)^(n/2) (4/Pi) n! lfun4(n+1) / (Pi/2)^n */
579 : static GEN
580 14 : eulerreal_using_lfun4(long n, long prec)
581 : {
582 14 : GEN pisur2 = Pi2n(-1, prec+EXTRAPRECWORD);
583 14 : GEN iz = inv_lfun4(n+1, prec);
584 14 : GEN z = divrr(mpfactr(n, prec), mulrr(powru(pisur2, n+1), iz));
585 14 : if ((n & 3L) == 2) setsigne(z, -1);
586 14 : shiftr_inplace(z, 1); return z;
587 : }
588 : /* Euler numbers: 1, 0, -1, 0, 5, 0, -61,... */
589 : GEN
590 217 : eulerfrac(long n)
591 : {
592 : pari_sp av;
593 : long k;
594 : GEN E;
595 217 : if (n <= 0)
596 : {
597 21 : if (n < 0) pari_err_DOMAIN("eulerfrac", "index", "<", gen_0, stoi(n));
598 14 : return gen_1;
599 : }
600 196 : if (odd(n)) return gen_0;
601 189 : k = n >> 1;
602 189 : if (!eulerzone) constreuler(0);
603 189 : if (eulerzone && k < lg(eulerzone)) return gel(eulerzone, k);
604 14 : av = avma; E = eulerreal_using_lfun4(n, eulerprec(n));
605 14 : return gerepileuptoleaf(av, roundr(E));
606 : }
607 : GEN
608 14 : eulervec(long n)
609 : {
610 : long i, l;
611 : GEN y;
612 14 : if (n < 0) return cgetg(1, t_VEC);
613 14 : constreuler(n);
614 14 : l = n+2; y = cgetg(l, t_VEC); gel(y,1) = gen_1;
615 49224 : for (i = 2; i < l; i++) gel(y,i) = gel(eulerzone,i-1);
616 14 : return y;
617 : }
618 :
619 : /* Return E_n */
620 : GEN
621 21 : eulerreal(long n, long prec)
622 : {
623 21 : pari_sp av = avma;
624 : GEN B;
625 : long p, k;
626 21 : if (n < 0) pari_err_DOMAIN("eulerreal", "index", "<", gen_0, stoi(n));
627 14 : if (n == 0) return real_1(prec);
628 14 : if (odd(n)) return real_0(prec);
629 :
630 14 : k = n >> 1;
631 14 : if (!eulerzone) constreuler(0);
632 14 : if (k < lg(eulerzone)) return itor(gel(eulerzone,k), prec);
633 0 : p = eulerprec(n);
634 0 : B = eulerreal_using_lfun4(n, minss(p, prec));
635 0 : if (p < prec) B = itor(roundr(B), prec);
636 0 : return gerepileuptoleaf(av, B);
637 : }
|
