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