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 81537 : fracB2k(GEN D)
30 : {
31 81537 : GEN a = utoipos(5), b = utoipos(6); /* 1/2 + 1/3 */
32 81599 : long i, l = lg(D);
33 619117 : for (i = 2; i < l; i++) /* skip 1 */
34 : {
35 537599 : ulong p = 2*D[i] + 1; /* a/b += 1/p */
36 537599 : if (uisprime(p)) { a = addii(muliu(a,p), b); b = muliu(b,p); }
37 : }
38 81518 : return mkfrac(a,b);
39 : }
40 : /* precision needed to compute B_k for all k <= N */
41 : long
42 1870 : bernbitprec(long N)
43 : { /* 1.612086 ~ log(8Pi) / 2 */
44 1870 : const double log2PI = 1.83787706641;
45 1870 : double logN = log((double)N);
46 1870 : double t = (N + 4) * logN - N*(1 + log2PI) + 1.612086;
47 1870 : 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 1373 : zetamaxpow(long n)
54 : {
55 1373 : long M = (long)ceil(n / (2 * M_PI * M_E));
56 1373 : return M | 1; /* make it odd */
57 : }
58 : /* v * zeta(k) using r precomputed odd powers */
59 : static GEN
60 81515 : bern_zeta(GEN v, long k, GEN pow, long r, long p)
61 : {
62 81515 : GEN z, s = gel(pow, r);
63 : long j;
64 10360988 : for (j = r - 2; j >= 3; j -= 2) s = addii(s, gel(pow,j));
65 81430 : z = s = itor(s, nbits2prec(p));
66 81511 : shiftr_inplace(s, -p); /* zeta(k)(1 - 2^(-k)) - 1*/
67 81510 : s = addrs(s, 1); shiftr_inplace(s, -k);
68 : /* divide by 1 - 2^(-k): s + s/2^k + s/2^(2k) + ... */
69 361125 : for (; k < p; k <<= 1) s = addrr(s, shiftr(s, -k));
70 81512 : return addrr(v, mulrr(v, addrr(z, s)));
71 : }
72 : /* z * j^2 */
73 : static GEN
74 50916113 : muliu2(GEN z, ulong j)
75 50916113 : { 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 295 : bernset(GEN *y, long m, long n)
79 : {
80 295 : long i, j, k, p, prec, r, N = n << 1; /* up to B_N */
81 : GEN u, b, v, t;
82 295 : p = bernbitprec(N); prec = nbits2prec(p);
83 295 : u = sqrr(Pi2n(1, prec)); /* (2Pi)^2 */
84 295 : v = divrr(mpfactr(N, prec), powru(u, n)); shiftr_inplace(v,1);
85 295 : r = zetamaxpow(N);
86 295 : t = cgetg(r+1, t_VEC); b = int2n(p); /* fixed point */
87 5558 : for (j = 3; j <= r; j += 2)
88 : {
89 5263 : GEN z = cgeti(nbits2lg(p));
90 5263 : pari_sp av2 = avma;
91 5263 : affii(divii(b, powuu(j, N)), z);
92 5263 : gel(t,j) = z; set_avma(av2);
93 : }
94 295 : y += n - m;
95 295 : for (i = n, k = N;; i--)
96 81218 : { /* set B_n, k = 2i */
97 81513 : pari_sp av2 = avma;
98 81513 : 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 81510 : if (!odd(i)) setsigne(B, -1); /* B ~ B_n */
102 81509 : B = roundr(addrr(B, fractor(z,LOWDEFAULTPREC))); /* B - z = B_n */
103 81518 : *y-- = gclone(gsub(B, z));
104 81518 : if (i == m) break;
105 81223 : affrr(divrunextu(mulrr(v,u), k-1), v);
106 10440222 : for (j = r; j >= 3; j -= 2) affii(muliu2(gel(t,j), j), gel(t,j));
107 81208 : set_avma(av2); k -= 2;
108 81213 : if (((N - k) & 0x7f) == 0x7e)
109 : { /* reduce precision if possible */
110 1078 : long p2 = p, prec2 = prec;
111 1078 : p = bernbitprec(k); prec = nbits2prec(p); if (prec2 == prec) continue;
112 1078 : setprec(v, prec); r = zetamaxpow(k);
113 158372 : for (j = 3; j <= r; j += 2) affii(shifti(gel(t,j), p - p2), gel(t,j));
114 1078 : set_avma(av2);
115 : }
116 : }
117 295 : }
118 : /* need B[2..2*nb] as t_INT or t_FRAC */
119 : void
120 229613 : constbern(long nb)
121 : {
122 229613 : const pari_sp av = avma;
123 : long i, l;
124 : GEN B;
125 : pari_timer T;
126 :
127 229613 : l = bernzone? lg(bernzone): 0;
128 229613 : if (l > nb) return;
129 :
130 295 : nb = maxss(nb, l + 127);
131 295 : B = cgetg_block(nb+1, t_VEC);
132 295 : if (bernzone)
133 7040 : { for (i = 1; i < l; i++) gel(B,i) = gel(bernzone,i); }
134 : else
135 : {
136 247 : gel(B,1) = gclone(mkfracss(1,6));
137 247 : gel(B,2) = gclone(mkfracss(-1,30));
138 247 : gel(B,3) = gclone(mkfracss(1,42));
139 247 : gel(B,4) = gclone(mkfracss(-1,30));
140 247 : gel(B,5) = gclone(mkfracss(5,66));
141 247 : gel(B,6) = gclone(mkfracss(-691,2730));
142 247 : gel(B,7) = gclone(mkfracss(7,6));
143 247 : gel(B,8) = gclone(mkfracss(-3617,510));
144 247 : gel(B,9) = gclone(mkfracss(43867,798));
145 247 : gel(B,10)= gclone(mkfracss(-174611,330));
146 247 : gel(B,11)= gclone(mkfracss(854513,138));
147 247 : gel(B,12)= gclone(mkfracss(-236364091,2730));
148 247 : gel(B,13)= gclone(mkfracss(8553103,6)); /* B_26 */
149 247 : l = 14;
150 : }
151 295 : set_avma(av);
152 295 : if (DEBUGLEVEL) {
153 0 : err_printf("caching Bernoulli numbers 2*%ld to 2*%ld\n", l, nb);
154 0 : timer_start(&T);
155 : }
156 295 : bernset((GEN*)B + l, l, nb);
157 295 : if (DEBUGLEVEL) timer_printf(&T, "Bernoulli");
158 295 : swap(B, bernzone); guncloneNULL(B);
159 295 : set_avma(av);
160 : #if 0
161 : if (nb > 200000)
162 : #endif
163 : {
164 295 : const ulong p = 4294967291UL;
165 295 : long n = 2 * nb + 2;
166 295 : GEN t = const_vecsmall(n+1, 1);
167 295 : t[1] = evalvarn(0); t[2] = 0;
168 295 : t = Flx_shift(Flx_invLaplace(t, p), -1); /* t = (exp(x)-1)/x */
169 295 : t = Flx_Laplace(Flxn_inv(t, n, p), p);
170 92022 : for (i = 1; i <= nb; i++)
171 91727 : 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 295 : 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+EXTRAPREC64);
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 8936939 : bernfrac(long n)
205 : {
206 : pari_sp av;
207 : long k;
208 8936939 : if (n <= 1)
209 : {
210 4487 : if (n < 0) pari_err_DOMAIN("bernfrac", "index", "<", gen_0, stoi(n));
211 4480 : return n? mkfrac(gen_m1,gen_2): gen_1;
212 : }
213 8932452 : if (odd(n)) return gen_0;
214 8928542 : k = n >> 1;
215 8928542 : if (!bernzone) constbern(0);
216 8928542 : 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 329 : bernvec(long n)
222 : {
223 : long i, l;
224 : GEN y;
225 329 : if (n < 0) return cgetg(1, t_VEC);
226 329 : constbern(n);
227 329 : l = n+2; y = cgetg(l, t_VEC); gel(y,1) = gen_1;
228 56679 : for (i = 2; i < l; i++) gel(y,i) = gel(bernzone,i-1);
229 329 : 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 2226 : bernpol_i(long k, long v)
235 : {
236 : GEN B, C;
237 : long i;
238 2226 : if (v < 0) v = 0;
239 2226 : constbern(k >> 1); /* cache B_2, ..., B_2[k/2] */
240 2226 : C = vecbinomial(k);
241 2226 : B = cgetg(k + 3, t_POL);
242 15099 : for (i = 0; i <= k; ++i) gel(B, k-i+2) = gmul(gel(C,i+1), bernfrac(i));
243 2226 : B[1] = evalsigne(1) | evalvarn(v);
244 2226 : return B;
245 : }
246 : GEN
247 2114 : bernpol(long k, long v)
248 : {
249 2114 : pari_sp av = avma;
250 2114 : if (k < 0) pari_err_DOMAIN("bernpol", "index", "<", gen_0, stoi(k));
251 2107 : return gerepileupto(av, bernpol_i(k, v));
252 : }
253 : GEN
254 98 : bernpol_eval(long k, GEN x)
255 : {
256 98 : pari_sp av = avma;
257 : GEN B;
258 98 : if (!x) return bernpol(k, 0);
259 49 : if (gequalX(x)) return bernpol(k, varn(x));
260 49 : if (k < 0) pari_err_DOMAIN("bernpol", "index", "<", gen_0, stoi(k));
261 42 : B = bernpol_i(k, fetch_var_higher());
262 42 : return gerepileupto(av, poleval(B, x));
263 : }
264 :
265 : /* x := pol_x(v); return 1^e + ... + x^e = x^e + (B_{e+1}(x) - B_{e+1})/(e+1) */
266 : static GEN
267 49 : faulhaber(long e, long v)
268 : {
269 : GEN B;
270 49 : if (e == 0) return pol_x(v);
271 35 : B = RgX_integ(bernpol_i(e, v)); /* (B_{e+1}(x) - B_{e+1}) / (e+1) */
272 35 : gel(B,e+2) = gaddgs(gel(B,e+2), 1); /* add x^e, in place */
273 35 : return B;
274 : }
275 : /* sum_v T(v), T a polynomial expression in v */
276 : GEN
277 49 : sumformal(GEN T, long v)
278 : {
279 49 : pari_sp av = avma, av2;
280 : long i, t, d;
281 : GEN R;
282 :
283 49 : T = simplify_shallow(T);
284 49 : t = typ(T);
285 49 : if (is_scalar_t(t))
286 14 : return gerepileupto(av, monomialcopy(T, 1, v < 0? 0: v));
287 35 : if (t != t_POL) pari_err_TYPE("sumformal [not a t_POL]", T);
288 28 : if (v < 0) v = varn(T);
289 28 : av2 = avma;
290 28 : R = gen_0;
291 28 : d = poldegree(T,v);
292 91 : for (i = d; i >= 0; i--)
293 : {
294 63 : GEN c = polcoef_i(T, i, v);
295 63 : if (gequal0(c)) continue;
296 42 : R = gadd(R, gmul(c, faulhaber(i, v)));
297 42 : if (gc_needed(av2,3))
298 : {
299 0 : if(DEBUGMEM>1) pari_warn(warnmem,"sumformal, i = %ld/%ld", i,d);
300 0 : R = gerepileupto(av2, R);
301 : }
302 : }
303 28 : return gerepileupto(av, R);
304 : }
305 :
306 : /* 1/zeta(n) using Euler product. Assume n > 0. */
307 : GEN
308 1194 : inv_szeta_euler(long n, long prec)
309 : {
310 1194 : long bit = prec2nbits(prec);
311 : GEN z, res;
312 : pari_sp av, av2;
313 : double A, D, lba;
314 : ulong p, lim;
315 : forprime_t S;
316 :
317 1194 : if (n > bit) return real_1(prec);
318 :
319 1187 : lba = prec2nbits_mul(prec, M_LN2);
320 1187 : D = exp((lba - log((double)(n-1))) / (n-1));
321 1187 : lim = 1 + (ulong)ceil(D);
322 1187 : if (lim < 3) return subir(gen_1,real2n(-n,prec));
323 1187 : res = cgetr(prec); av = avma; incrprec(prec);
324 :
325 1187 : (void)u_forprime_init(&S, 3, lim);
326 1187 : av2 = avma; A = n / M_LN2; z = subir(gen_1, real2n(-n, prec));
327 86205 : while ((p = u_forprime_next(&S)))
328 : {
329 85018 : long l = bit - (long)floor(A * log((double)p));
330 : GEN h;
331 :
332 85018 : if (l < BITS_IN_LONG) l = BITS_IN_LONG;
333 85018 : l = minss(prec, nbits2prec(l));
334 85018 : h = divrr(z, rpowuu(p, (ulong)n, l));
335 85018 : z = subrr(z, h);
336 85018 : if (gc_needed(av,1))
337 : {
338 0 : if (DEBUGMEM>1) pari_warn(warnmem,"inv_szeta_euler, p = %lu/%lu", p,lim);
339 0 : z = gerepileuptoleaf(av2, z);
340 : }
341 : }
342 1187 : affrr(z, res); set_avma(av); return res;
343 : }
344 :
345 : /* Return B_n */
346 : GEN
347 15414 : bernreal(long n, long prec)
348 : {
349 : pari_sp av;
350 : GEN B;
351 : long p, k;
352 15414 : if (n < 0) pari_err_DOMAIN("bernreal", "index", "<", gen_0, stoi(n));
353 15407 : if (n == 0) return real_1(prec);
354 15407 : if (n == 1) return real_m2n(-1,prec); /*-1/2*/
355 15407 : if (odd(n)) return real_0(prec);
356 :
357 15407 : k = n >> 1;
358 15407 : if (!bernzone) constbern(0);
359 15407 : if (k < lg(bernzone)) return fractor(gel(bernzone,k), prec);
360 7 : p = bernprec(n); av = avma;
361 7 : B = bernreal_using_zeta(n, minss(p, prec));
362 7 : if (p < prec) B = fractor(bernfrac_i(n, B), prec);
363 7 : return gerepileuptoleaf(av, B);
364 : }
365 :
366 : GEN
367 49 : eulerpol(long k, long v)
368 : {
369 49 : pari_sp av = avma;
370 : GEN B, E;
371 49 : if (k < 0) pari_err_DOMAIN("eulerpol", "index", "<", gen_0, stoi(k));
372 42 : k++; B = bernpol_i(k, v);
373 42 : E = RgX_Rg_mul(RgX_sub(B, RgX_rescale(B, gen_2)), uutoQ(2,k));
374 42 : return gerepileupto(av, E);
375 : }
376 :
377 : /*******************************************************************/
378 : /** HARMONIC NUMBERS **/
379 : /*******************************************************************/
380 : /* 1/a + ... + 1/(b-1); a < b <= 2^(BIL-1) */
381 : static GEN
382 4683 : hrec(ulong a, ulong b)
383 : {
384 : ulong m;
385 4683 : switch(b - a)
386 : {
387 4501 : case 1: retmkfrac(gen_1, utoipos(a));
388 140 : case 2: if (a < 65536) retmkfrac(utoipos(2*a + 1), utoipos(a * a + a));
389 0 : retmkfrac(utoipos(2*a + 1), muluu(a, a+1));
390 : }
391 42 : m = (a + b) >> 1;
392 42 : return gadd(hrec(a, m), hrec(m, b));
393 : }
394 : /* exact Harmonic number H_n, n < 2^(BIL-1).
395 : * Could use H_n = sum_k 2^(-k) H^odd_{n \ 2^k} */
396 : GEN
397 4599 : harmonic(ulong n)
398 : {
399 4599 : pari_sp av = avma;
400 4599 : return n? gerepileupto(av, hrec(1, n+1)): gen_0;
401 : }
402 :
403 : /* 1/a^k + ... + 1/(b-1)^k; a < b */
404 : static GEN
405 77 : hreck(ulong a, ulong b, ulong k)
406 : {
407 : ulong m;
408 77 : switch(b - a)
409 : {
410 : GEN x, y;
411 14 : case 1: retmkfrac(gen_1, powuu(a, k));
412 28 : case 2:
413 28 : x = powuu(a, k); y = powuu(a + 1, k);
414 28 : retmkfrac(addii(x, y), mulii(x, y));
415 : }
416 35 : m = (a + b) >> 1;
417 35 : return gadd(hreck(a, m, k), hreck(m, b, k));
418 : }
419 : GEN
420 56 : harmonic0(ulong n, GEN k)
421 : {
422 56 : pari_sp av = avma;
423 : ulong r;
424 56 : if (!n) return gen_0;
425 49 : if (n & HIGHBIT) pari_err_OVERFLOW("harmonic");
426 49 : if (!k) return harmonic(n);
427 21 : if (typ(k) != t_INT) pari_err_TYPE("harmonic", k);
428 14 : if (signe(k) < 0)
429 : {
430 7 : GEN H = poleval(faulhaber(-itos(k), 0), utoipos(n));
431 7 : return gerepileuptoint(av, H);
432 : }
433 7 : r = itou(k);
434 7 : if (!r) return utoipos(n);
435 7 : if (r == 1) return harmonic(n);
436 7 : return gerepileupto(av, hreck(1, n+1, r));
437 : }
438 :
439 :
440 : /**************************************************************/
441 : /* Euler numbers */
442 : /**************************************************************/
443 :
444 : /* precision needed to compute E_k for all k <= N */
445 : static long
446 798 : eulerbitprec(long N)
447 : { /* 1.1605 ~ log(32/Pi) / 2 */
448 798 : const double logPIS2 = 0.4515827;
449 798 : double t = (N + 1) * log((double)N) - N*(1 + logPIS2) + 1.1605;
450 798 : return (long)ceil(t / M_LN2) + 10;
451 : }
452 : static long
453 14 : eulerprec(long N) { return nbits2prec(eulerbitprec(N)); }
454 :
455 : /* sum_{k > M, k odd} (-1)^((k-1)/2)k^(-n) < M^(-n) < 2^-b */
456 : static long
457 784 : lfun4maxpow(long n)
458 : {
459 784 : long M = (long)ceil(2 * n / (M_E * M_PI));
460 784 : return M | 1; /* make it odd */
461 : }
462 :
463 : /* lfun4(k) using r precomputed odd powers */
464 : static GEN
465 49791 : euler_lfun4(GEN v, GEN pow, long r, long p)
466 : {
467 49791 : GEN s = ((r & 3L) == 1)? gel(pow, r): negi(gel(pow, r));
468 : long j;
469 40556894 : for (j = r - 2; j >= 3; j -= 2)
470 40507103 : s = ((j & 3L) == 1)? addii(s, gel(pow,j)): subii(s, gel(pow,j));
471 49791 : s = mulri(v, s); shiftr_inplace(s, -p);
472 49791 : return addrr(v, s);
473 : }
474 :
475 : /* 1 <= m <= n, set y[1] = E_{2m}, ... y[n-m+1] = E_{2n} in Z */
476 : static void
477 21 : eulerset(GEN *y, long m, long n)
478 : {
479 21 : long i, j, k, p, prec, r, N = n << 1, N1 = N + 1; /* up to E_N */
480 : GEN b, u, v, t;
481 21 : p = eulerbitprec(N); prec = nbits2prec(p);
482 21 : u = sqrr(Pi2n(-1, prec)); /* (Pi/2)^2 */
483 21 : v = divrr(mpfactr(N, prec), mulrr(powru(u, n), Pi2n(-2,prec)));
484 21 : r = lfun4maxpow(N1);
485 21 : t = cgetg(r+1, t_VEC); b = int2n(p); /* fixed point */
486 11921 : for (j = 3; j <= r; j += 2)
487 : {
488 11900 : GEN z = cgeti(nbits2lg(p));
489 11900 : pari_sp av2 = avma;
490 11900 : affii(divii(b, powuu(j, N+1)), z);
491 11900 : gel(t,j) = z; set_avma(av2);
492 : }
493 21 : y += n - m;
494 21 : for (i = n, k = N1;; i--)
495 49770 : { /* set E_n, k = 2i + 1 */
496 49791 : pari_sp av2 = avma;
497 49791 : GEN E = euler_lfun4(v, t, r, p);
498 : long j;
499 : /* E = v * lfun4(k), v = (4/Pi)*k! / (Pi/2)^k */
500 49791 : E = roundr(E); if (odd(i)) setsigne(E, -1); /* E ~ E_n */
501 49791 : *y-- = gclone(E);
502 49791 : if (i == m) break;
503 49770 : affrr(divrunextu(mulrr(v,u), k-2), v);
504 40606202 : for (j = r; j >= 3; j -= 2) affii(muliu2(gel(t,j), j), gel(t,j));
505 49770 : set_avma(av2); k -= 2;
506 49770 : if (((N1 - k) & 0x7f) == 0x7e)
507 : { /* reduce precision if possible */
508 763 : long p2 = p, prec2 = prec;
509 763 : p = eulerbitprec(k); prec = nbits2prec(p); if (prec2 == prec) continue;
510 763 : setprec(v, prec); r = lfun4maxpow(k);
511 622923 : for (j = 3; j <= r; j += 2) affii(shifti(gel(t,j), p - p2), gel(t,j));
512 763 : set_avma(av2);
513 : }
514 : }
515 21 : }
516 :
517 : /* need E[2..2*nb] as t_INT */
518 : static void
519 28 : constreuler(long nb)
520 : {
521 28 : const pari_sp av = avma;
522 : long i, l;
523 : GEN E;
524 : pari_timer T;
525 :
526 28 : l = eulerzone? lg(eulerzone): 0;
527 28 : if (l > nb) return;
528 :
529 21 : nb = maxss(nb, l + 127);
530 21 : E = cgetg_block(nb+1, t_VEC);
531 21 : if (eulerzone)
532 896 : { for (i = 1; i < l; i++) gel(E,i) = gel(eulerzone,i); }
533 : else
534 : {
535 14 : gel(E,1) = gclone(stoi(-1));
536 14 : gel(E,2) = gclone(stoi(5));
537 14 : gel(E,3) = gclone(stoi(-61));
538 14 : gel(E,4) = gclone(stoi(1385));
539 14 : gel(E,5) = gclone(stoi(-50521));
540 14 : gel(E,6) = gclone(stoi(2702765));
541 14 : gel(E,7) = gclone(stoi(-199360981));
542 14 : l = 8;
543 : }
544 21 : set_avma(av);
545 21 : if (DEBUGLEVEL) {
546 0 : err_printf("caching Euler numbers 2*%ld to 2*%ld\n", l, nb);
547 0 : timer_start(&T);
548 : }
549 21 : eulerset((GEN*)E + l, l, nb);
550 21 : if (DEBUGLEVEL) timer_printf(&T, "Euler");
551 21 : swap(E, eulerzone); guncloneNULL(E);
552 21 : set_avma(av);
553 : }
554 :
555 : /* 1/lfun(-4,n) using Euler product. Assume n > 0. */
556 : static GEN
557 14 : inv_lfun4(long n, long prec)
558 : {
559 14 : long bit = prec2nbits(prec);
560 : GEN z, res;
561 : pari_sp av, av2;
562 : double A;
563 : ulong p, lim;
564 : forprime_t S;
565 :
566 14 : if (n > bit) return real_1(prec);
567 :
568 14 : lim = 1 + (ulong)ceil(exp2((double)bit / n));
569 14 : res = cgetr(prec); av = avma; incrprec(prec);
570 :
571 14 : (void)u_forprime_init(&S, 3, lim);
572 14 : av2 = avma; A = n / M_LN2; z = real_1(prec);
573 369 : while ((p = u_forprime_next(&S)))
574 : {
575 355 : long l = bit - (long)floor(A * log((double)p));
576 : GEN h;
577 :
578 355 : if (l < BITS_IN_LONG) l = BITS_IN_LONG;
579 355 : l = minss(prec, nbits2prec(l));
580 355 : h = rpowuu(p, (ulong)n, l); if ((p & 3UL) == 1) setsigne(h, -1);
581 355 : z = addrr(z, divrr(z, h)); /* z *= 1 - chi_{-4}(p) / p^n */
582 355 : if (gc_needed(av,1))
583 : {
584 0 : if (DEBUGMEM>1) pari_warn(warnmem,"inv_lfun4, p = %lu/%lu", p,lim);
585 0 : z = gerepileuptoleaf(av2, z);
586 : }
587 : }
588 14 : affrr(z, res); set_avma(av); return res;
589 : }
590 : /* assume n even > 0, E_n = (-1)^(n/2) (4/Pi) n! lfun4(n+1) / (Pi/2)^n */
591 : static GEN
592 14 : eulerreal_using_lfun4(long n, long prec)
593 : {
594 14 : GEN pisur2 = Pi2n(-1, prec+EXTRAPREC64);
595 14 : GEN iz = inv_lfun4(n+1, prec);
596 14 : GEN z = divrr(mpfactr(n, prec), mulrr(powru(pisur2, n+1), iz));
597 14 : if ((n & 3L) == 2) setsigne(z, -1);
598 14 : shiftr_inplace(z, 1); return z;
599 : }
600 : /* Euler numbers: 1, 0, -1, 0, 5, 0, -61,... */
601 : GEN
602 217 : eulerfrac(long n)
603 : {
604 : pari_sp av;
605 : long k;
606 : GEN E;
607 217 : if (n <= 0)
608 : {
609 21 : if (n < 0) pari_err_DOMAIN("eulerfrac", "index", "<", gen_0, stoi(n));
610 14 : return gen_1;
611 : }
612 196 : if (odd(n)) return gen_0;
613 189 : k = n >> 1;
614 189 : if (!eulerzone) constreuler(0);
615 189 : if (eulerzone && k < lg(eulerzone)) return gel(eulerzone, k);
616 14 : av = avma; E = eulerreal_using_lfun4(n, eulerprec(n));
617 14 : return gerepileuptoleaf(av, roundr(E));
618 : }
619 : GEN
620 14 : eulervec(long n)
621 : {
622 : long i, l;
623 : GEN y;
624 14 : if (n < 0) return cgetg(1, t_VEC);
625 14 : constreuler(n);
626 14 : l = n+2; y = cgetg(l, t_VEC); gel(y,1) = gen_1;
627 49224 : for (i = 2; i < l; i++) gel(y,i) = gel(eulerzone,i-1);
628 14 : return y;
629 : }
630 :
631 : /* Return E_n */
632 : GEN
633 21 : eulerreal(long n, long prec)
634 : {
635 21 : pari_sp av = avma;
636 : GEN B;
637 : long p, k;
638 21 : if (n < 0) pari_err_DOMAIN("eulerreal", "index", "<", gen_0, stoi(n));
639 14 : if (n == 0) return real_1(prec);
640 14 : if (odd(n)) return real_0(prec);
641 :
642 14 : k = n >> 1;
643 14 : if (!eulerzone) constreuler(0);
644 14 : if (k < lg(eulerzone)) return itor(gel(eulerzone,k), prec);
645 0 : p = eulerprec(n);
646 0 : B = eulerreal_using_lfun4(n, minss(p, prec));
647 0 : if (p < prec) B = itor(roundr(B), prec);
648 0 : return gerepileuptoleaf(av, B);
649 : }
|
