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 82433 : fracB2k(GEN D)
30 : {
31 82433 : GEN a = utoipos(5), b = utoipos(6); /* 1/2 + 1/3 */
32 82450 : long i, l = lg(D);
33 625326 : for (i = 2; i < l; i++) /* skip 1 */
34 : {
35 542911 : ulong p = 2*D[i] + 1; /* a/b += 1/p */
36 542911 : if (uisprime(p)) { a = addii(muliu(a,p), b); b = muliu(b,p); }
37 : }
38 82415 : return mkfrac(a,b);
39 : }
40 : /* precision needed to compute B_k for all k <= N */
41 : long
42 1891 : bernbitprec(long N)
43 : { /* 1.612086 ~ log(8Pi) / 2 */
44 1891 : const double log2PI = 1.83787706641;
45 1891 : double logN = log((double)N);
46 1891 : double t = (N + 4) * logN - N*(1 + log2PI) + 1.612086;
47 1891 : 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 1394 : zetamaxpow(long n)
54 : {
55 1394 : long M = (long)ceil(n / (2 * M_PI * M_E));
56 1394 : return M | 1; /* make it odd */
57 : }
58 : /* v * zeta(k) using r precomputed odd powers */
59 : static GEN
60 82409 : bern_zeta(GEN v, long k, GEN pow, long r, long p)
61 : {
62 82409 : GEN z, s = gel(pow, r);
63 : long j;
64 10379307 : for (j = r - 2; j >= 3; j -= 2) s = addii(s, gel(pow,j));
65 82354 : z = s = itor(s, nbits2prec(p));
66 82410 : shiftr_inplace(s, -p); /* zeta(k)(1 - 2^(-k)) - 1*/
67 82411 : s = addrs(s, 1); shiftr_inplace(s, -k);
68 : /* divide by 1 - 2^(-k): s + s/2^k + s/2^(2k) + ... */
69 364728 : for (; k < p; k <<= 1) s = addrr(s, shiftr(s, -k));
70 82407 : return addrr(v, mulrr(v, addrr(z, s)));
71 : }
72 : /* z * j^2 */
73 : static GEN
74 50934260 : muliu2(GEN z, ulong j)
75 50934260 : { 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 302 : bernset(GEN *y, long m, long n)
79 : {
80 302 : long i, j, k, p, prec, r, N = n << 1; /* up to B_N */
81 : GEN u, b, v, t;
82 302 : p = bernbitprec(N); prec = nbits2prec(p);
83 302 : u = sqrr(Pi2n(1, prec)); /* (2Pi)^2 */
84 302 : v = divrr(mpfactr(N, prec), powru(u, n)); shiftr_inplace(v,1);
85 302 : r = zetamaxpow(N);
86 302 : t = cgetg(r+1, t_VEC); b = int2n(p); /* fixed point */
87 5719 : for (j = 3; j <= r; j += 2)
88 : {
89 5417 : GEN z = cgeti(nbits2lg(p));
90 5417 : pari_sp av2 = avma;
91 5417 : affii(divii(b, powuu(j, N)), z);
92 5417 : gel(t,j) = z; set_avma(av2);
93 : }
94 302 : y += n - m;
95 302 : for (i = n, k = N;; i--)
96 82103 : { /* set B_n, k = 2i */
97 82405 : pari_sp av2 = avma;
98 82405 : 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 82413 : if (!odd(i)) setsigne(B, -1); /* B ~ B_n */
102 82412 : B = roundr(addrr(B, fractor(z,LOWDEFAULTPREC))); /* B - z = B_n */
103 82415 : *y-- = gclone(gsub(B, z));
104 82415 : if (i == m) break;
105 82113 : affrr(divrunextu(mulrr(v,u), k-1), v);
106 10459349 : for (j = r; j >= 3; j -= 2) affii(muliu2(gel(t,j), j), gel(t,j));
107 82100 : set_avma(av2); k -= 2;
108 82102 : if (((N - k) & 0x7f) == 0x7e)
109 : { /* reduce precision if possible */
110 1092 : long p2 = p, prec2 = prec;
111 1092 : p = bernbitprec(k); prec = nbits2prec(p); if (prec2 == prec) continue;
112 1092 : setprec(v, prec); r = zetamaxpow(k);
113 158624 : for (j = 3; j <= r; j += 2) affii(shifti(gel(t,j), p - p2), gel(t,j));
114 1092 : set_avma(av2);
115 : }
116 : }
117 302 : }
118 : /* need B[2..2*nb] as t_INT or t_FRAC */
119 : void
120 7240057 : constbern(long nb)
121 : {
122 7240057 : const pari_sp av = avma;
123 : long i, l;
124 : GEN B;
125 : pari_timer T;
126 :
127 7240057 : l = bernzone? lg(bernzone): 0;
128 7240057 : if (l > nb) return;
129 :
130 302 : nb = maxss(nb, l + 127);
131 302 : B = cgetg_block(nb+1, t_VEC);
132 302 : if (bernzone)
133 8832 : { 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 302 : set_avma(av);
152 302 : if (DEBUGLEVEL) {
153 0 : err_printf("caching Bernoulli numbers 2*%ld to 2*%ld\n", l, nb);
154 0 : timer_start(&T);
155 : }
156 302 : bernset((GEN*)B + l, l, nb);
157 302 : if (DEBUGLEVEL) timer_printf(&T, "Bernoulli");
158 302 : swap(B, bernzone); guncloneNULL(B);
159 302 : set_avma(av);
160 : #if 0
161 : if (nb > 200000)
162 : #endif
163 : {
164 302 : const ulong p = 4294967291UL;
165 302 : long n = 2 * nb + 2;
166 302 : GEN t = const_vecsmall(n+1, 1);
167 302 : t[1] = evalvarn(0); t[2] = 0;
168 302 : t = Flx_shift(Flx_invLaplace(t, p), -1); /* t = (exp(x)-1)/x */
169 302 : t = Flx_Laplace(Flxn_inv(t, n, p), p);
170 94710 : for (i = 1; i <= nb; i++)
171 94408 : 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 302 : 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 51724745 : bernfrac(long n)
205 : {
206 : pari_sp av;
207 : long k;
208 51724745 : if (n <= 1)
209 : {
210 14010318 : if (n < 0) pari_err_DOMAIN("bernfrac", "index", "<", gen_0, stoi(n));
211 14010311 : return n? mkfrac(gen_m1,gen_2): gen_1;
212 : }
213 37714427 : if (odd(n)) return gen_0;
214 23709643 : k = n >> 1;
215 23709643 : if (!bernzone) constbern(0);
216 23709643 : if (bernzone && k < lg(bernzone)) return gel(bernzone, k);
217 25 : av = avma;
218 25 : 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 7002660 : bernpol_i(long k, long v)
235 : {
236 : GEN B, C;
237 : long i;
238 7002660 : if (v < 0) v = 0;
239 7002660 : constbern(k >> 1); /* cache B_2, ..., B_2[k/2] */
240 7002660 : C = vecbinomial(k);
241 7002660 : B = cgetg(k + 3, t_POL);
242 49018592 : for (i = 0; i <= k; ++i) gel(B, k-i+2) = gmul(gel(C,i+1), bernfrac(i));
243 7002660 : B[1] = evalsigne(1) | evalvarn(v);
244 7002660 : return B;
245 : }
246 : GEN
247 2548 : bernpol(long k, long v)
248 : {
249 2548 : pari_sp av = avma;
250 2548 : if (k < 0) pari_err_DOMAIN("bernpol", "index", "<", gen_0, stoi(k));
251 2541 : 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_leaf(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_leaf(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 : /* Compute substpol(x*cotanh(x+O(x^p)),x^2,x) */
379 : static GEN
380 1169 : cotanh_Flxn(ulong p)
381 : {
382 : long i, j;
383 1169 : ulong t = p-1, m = (p+1)>>1;
384 1169 : GEN C = cgetg(m+1, t_VECSMALL);
385 1169 : GEN S = cgetg(m+1, t_VECSMALL);
386 1169 : C[1] = S[1] = 0;
387 1169 : C[2] = S[2] = 1;
388 324793 : for (i = p-1, j = m; j > 2; j--)
389 : {
390 323624 : t = Fl_mul(t, i--, p);
391 323624 : uel(S,j) = t;
392 323624 : t = Fl_mul(t, i--, p);
393 323624 : uel(C,j) = t;
394 : }
395 1169 : return Flxn_div(C, S, m-1 ,p);
396 : }
397 :
398 : long
399 1183 : primeisregular(ulong p)
400 : {
401 1183 : pari_sp av = avma;
402 : GEN T;
403 1183 : long i, m = (p+1)>>1;
404 1183 : if (p == 2 || p==3) return 1;
405 1169 : if (!uisprime(p)) pari_err_PRIME("primeisregular",utoi(p));
406 1169 : T = cotanh_Flxn(p);
407 1169 : if (lg(T)-1 < m) return 0;
408 200592 : for (i = 2; i <= m; i++)
409 199878 : if (uel(T,i)==0) return gc_long(av, 0);
410 714 : return gc_long(av, 1);
411 : }
412 :
413 : /*******************************************************************/
414 : /** HARMONIC NUMBERS **/
415 : /*******************************************************************/
416 : /* 1/a + ... + 1/(b-1); a < b <= 2^(BIL-1) */
417 : static GEN
418 4683 : hrec(ulong a, ulong b)
419 : {
420 : ulong m;
421 4683 : switch(b - a)
422 : {
423 4501 : case 1: if (a==1) return gen_1; else retmkfrac(gen_1, utoipos(a));
424 140 : case 2: if (a < 65536) retmkfrac(utoipos(2*a + 1), utoipos(a * a + a));
425 0 : retmkfrac(utoipos(2*a + 1), muluu(a, a+1));
426 : }
427 42 : m = (a + b) >> 1;
428 42 : return gadd(hrec(a, m), hrec(m, b));
429 : }
430 : /* exact Harmonic number H_n, n < 2^(BIL-1).
431 : * Could use H_n = sum_k 2^(-k) H^odd_{n \ 2^k} */
432 : GEN
433 4599 : harmonic(ulong n)
434 : {
435 4599 : pari_sp av = avma;
436 4599 : return n? gc_upto(av, hrec(1, n+1)): gen_0;
437 : }
438 :
439 : /* 1/a^k + ... + 1/(b-1)^k; a < b */
440 : static GEN
441 77 : hreck(ulong a, ulong b, ulong k)
442 : {
443 : ulong m;
444 77 : switch(b - a)
445 : {
446 : GEN x, y;
447 14 : case 1: retmkfrac(gen_1, powuu(a, k));
448 28 : case 2:
449 28 : x = powuu(a, k); y = powuu(a + 1, k);
450 28 : retmkfrac(addii(x, y), mulii(x, y));
451 : }
452 35 : m = (a + b) >> 1;
453 35 : return gadd(hreck(a, m, k), hreck(m, b, k));
454 : }
455 : GEN
456 56 : harmonic0(ulong n, GEN k)
457 : {
458 56 : pari_sp av = avma;
459 : ulong r;
460 56 : if (!n) return gen_0;
461 49 : if (n & HIGHBIT) pari_err_OVERFLOW("harmonic");
462 49 : if (!k) return harmonic(n);
463 21 : if (typ(k) != t_INT) pari_err_TYPE("harmonic", k);
464 14 : if (signe(k) < 0)
465 : {
466 7 : GEN H = poleval(faulhaber(-itos(k), 0), utoipos(n));
467 7 : return gc_INT(av, H);
468 : }
469 7 : r = itou(k);
470 7 : if (!r) return utoipos(n);
471 7 : if (r == 1) return harmonic(n);
472 7 : return gc_upto(av, hreck(1, n+1, r));
473 : }
474 :
475 :
476 : /**************************************************************/
477 : /* Euler numbers */
478 : /**************************************************************/
479 :
480 : /* precision needed to compute E_k for all k <= N */
481 : static long
482 798 : eulerbitprec(long N)
483 : { /* 1.1605 ~ log(32/Pi) / 2 */
484 798 : const double logPIS2 = 0.4515827;
485 798 : double t = (N + 1) * log((double)N) - N*(1 + logPIS2) + 1.1605;
486 798 : return (long)ceil(t / M_LN2) + 10;
487 : }
488 : static long
489 14 : eulerprec(long N) { return nbits2prec(eulerbitprec(N)); }
490 :
491 : /* sum_{k > M, k odd} (-1)^((k-1)/2)k^(-n) < M^(-n) < 2^-b */
492 : static long
493 784 : lfun4maxpow(long n)
494 : {
495 784 : long M = (long)ceil(2 * n / (M_E * M_PI));
496 784 : return M | 1; /* make it odd */
497 : }
498 :
499 : /* lfun4(k) using r precomputed odd powers */
500 : static GEN
501 49791 : euler_lfun4(GEN v, GEN pow, long r, long p)
502 : {
503 49791 : GEN s = ((r & 3L) == 1)? gel(pow, r): negi(gel(pow, r));
504 : long j;
505 40556894 : for (j = r - 2; j >= 3; j -= 2)
506 40507103 : s = ((j & 3L) == 1)? addii(s, gel(pow,j)): subii(s, gel(pow,j));
507 49791 : s = mulri(v, s); shiftr_inplace(s, -p);
508 49791 : return addrr(v, s);
509 : }
510 :
511 : /* 1 <= m <= n, set y[1] = E_{2m}, ... y[n-m+1] = E_{2n} in Z */
512 : static void
513 21 : eulerset(GEN *y, long m, long n)
514 : {
515 21 : long i, j, k, p, prec, r, N = n << 1, N1 = N + 1; /* up to E_N */
516 : GEN b, u, v, t;
517 21 : p = eulerbitprec(N); prec = nbits2prec(p);
518 21 : u = sqrr(Pi2n(-1, prec)); /* (Pi/2)^2 */
519 21 : v = divrr(mpfactr(N, prec), mulrr(powru(u, n), Pi2n(-2,prec)));
520 21 : r = lfun4maxpow(N1);
521 21 : t = cgetg(r+1, t_VEC); b = int2n(p); /* fixed point */
522 11921 : for (j = 3; j <= r; j += 2)
523 : {
524 11900 : GEN z = cgeti(nbits2lg(p));
525 11900 : pari_sp av2 = avma;
526 11900 : affii(divii(b, powuu(j, N+1)), z);
527 11900 : gel(t,j) = z; set_avma(av2);
528 : }
529 21 : y += n - m;
530 21 : for (i = n, k = N1;; i--)
531 49770 : { /* set E_n, k = 2i + 1 */
532 49791 : pari_sp av2 = avma;
533 49791 : GEN E = euler_lfun4(v, t, r, p);
534 : long j;
535 : /* E = v * lfun4(k), v = (4/Pi)*k! / (Pi/2)^k */
536 49791 : E = roundr(E); if (odd(i)) setsigne(E, -1); /* E ~ E_n */
537 49791 : *y-- = gclone(E);
538 49791 : if (i == m) break;
539 49770 : affrr(divrunextu(mulrr(v,u), k-2), v);
540 40606202 : for (j = r; j >= 3; j -= 2) affii(muliu2(gel(t,j), j), gel(t,j));
541 49770 : set_avma(av2); k -= 2;
542 49770 : if (((N1 - k) & 0x7f) == 0x7e)
543 : { /* reduce precision if possible */
544 763 : long p2 = p, prec2 = prec;
545 763 : p = eulerbitprec(k); prec = nbits2prec(p); if (prec2 == prec) continue;
546 763 : setprec(v, prec); r = lfun4maxpow(k);
547 622923 : for (j = 3; j <= r; j += 2) affii(shifti(gel(t,j), p - p2), gel(t,j));
548 763 : set_avma(av2);
549 : }
550 : }
551 21 : }
552 :
553 : /* need E[2..2*nb] as t_INT */
554 : static void
555 28 : constreuler(long nb)
556 : {
557 28 : const pari_sp av = avma;
558 : long i, l;
559 : GEN E;
560 : pari_timer T;
561 :
562 28 : l = eulerzone? lg(eulerzone): 0;
563 28 : if (l > nb) return;
564 :
565 21 : nb = maxss(nb, l + 127);
566 21 : E = cgetg_block(nb+1, t_VEC);
567 21 : if (eulerzone)
568 896 : { for (i = 1; i < l; i++) gel(E,i) = gel(eulerzone,i); }
569 : else
570 : {
571 14 : gel(E,1) = gclone(stoi(-1));
572 14 : gel(E,2) = gclone(stoi(5));
573 14 : gel(E,3) = gclone(stoi(-61));
574 14 : gel(E,4) = gclone(stoi(1385));
575 14 : gel(E,5) = gclone(stoi(-50521));
576 14 : gel(E,6) = gclone(stoi(2702765));
577 14 : gel(E,7) = gclone(stoi(-199360981));
578 14 : l = 8;
579 : }
580 21 : set_avma(av);
581 21 : if (DEBUGLEVEL) {
582 0 : err_printf("caching Euler numbers 2*%ld to 2*%ld\n", l, nb);
583 0 : timer_start(&T);
584 : }
585 21 : eulerset((GEN*)E + l, l, nb);
586 21 : if (DEBUGLEVEL) timer_printf(&T, "Euler");
587 21 : swap(E, eulerzone); guncloneNULL(E);
588 21 : set_avma(av);
589 : }
590 :
591 : /* 1/lfun(-4,n) using Euler product. Assume n > 0. */
592 : static GEN
593 14 : inv_lfun4(long n, long prec)
594 : {
595 14 : long bit = prec2nbits(prec);
596 : GEN z, res;
597 : pari_sp av, av2;
598 : double A;
599 : ulong p, lim;
600 : forprime_t S;
601 :
602 14 : if (n > bit) return real_1(prec);
603 :
604 14 : lim = 1 + (ulong)ceil(exp2((double)bit / n));
605 14 : res = cgetr(prec); av = avma; incrprec(prec);
606 :
607 14 : (void)u_forprime_init(&S, 3, lim);
608 14 : av2 = avma; A = n / M_LN2; z = real_1(prec);
609 369 : while ((p = u_forprime_next(&S)))
610 : {
611 355 : long l = bit - (long)floor(A * log((double)p));
612 : GEN h;
613 :
614 355 : if (l < BITS_IN_LONG) l = BITS_IN_LONG;
615 355 : l = minss(prec, nbits2prec(l));
616 355 : h = rpowuu(p, (ulong)n, l); if ((p & 3UL) == 1) setsigne(h, -1);
617 355 : z = addrr(z, divrr(z, h)); /* z *= 1 - chi_{-4}(p) / p^n */
618 355 : if (gc_needed(av,1))
619 : {
620 0 : if (DEBUGMEM>1) pari_warn(warnmem,"inv_lfun4, p = %lu/%lu", p,lim);
621 0 : z = gc_leaf(av2, z);
622 : }
623 : }
624 14 : affrr(z, res); set_avma(av); return res;
625 : }
626 : /* assume n even > 0, E_n = (-1)^(n/2) (4/Pi) n! lfun4(n+1) / (Pi/2)^n */
627 : static GEN
628 14 : eulerreal_using_lfun4(long n, long prec)
629 : {
630 14 : GEN pisur2 = Pi2n(-1, prec+EXTRAPREC64);
631 14 : GEN iz = inv_lfun4(n+1, prec);
632 14 : GEN z = divrr(mpfactr(n, prec), mulrr(powru(pisur2, n+1), iz));
633 14 : if ((n & 3L) == 2) setsigne(z, -1);
634 14 : shiftr_inplace(z, 1); return z;
635 : }
636 : /* Euler numbers: 1, 0, -1, 0, 5, 0, -61,... */
637 : GEN
638 217 : eulerfrac(long n)
639 : {
640 : pari_sp av;
641 : long k;
642 : GEN E;
643 217 : if (n <= 0)
644 : {
645 21 : if (n < 0) pari_err_DOMAIN("eulerfrac", "index", "<", gen_0, stoi(n));
646 14 : return gen_1;
647 : }
648 196 : if (odd(n)) return gen_0;
649 189 : k = n >> 1;
650 189 : if (!eulerzone) constreuler(0);
651 189 : if (eulerzone && k < lg(eulerzone)) return gel(eulerzone, k);
652 14 : av = avma; E = eulerreal_using_lfun4(n, eulerprec(n));
653 14 : return gc_leaf(av, roundr(E));
654 : }
655 : GEN
656 14 : eulervec(long n)
657 : {
658 : long i, l;
659 : GEN y;
660 14 : if (n < 0) return cgetg(1, t_VEC);
661 14 : constreuler(n);
662 14 : l = n+2; y = cgetg(l, t_VEC); gel(y,1) = gen_1;
663 49224 : for (i = 2; i < l; i++) gel(y,i) = gel(eulerzone,i-1);
664 14 : return y;
665 : }
666 :
667 : /* Return E_n */
668 : GEN
669 21 : eulerreal(long n, long prec)
670 : {
671 21 : pari_sp av = avma;
672 : GEN B;
673 : long p, k;
674 21 : if (n < 0) pari_err_DOMAIN("eulerreal", "index", "<", gen_0, stoi(n));
675 14 : if (n == 0) return real_1(prec);
676 14 : if (odd(n)) return real_0(prec);
677 :
678 14 : k = n >> 1;
679 14 : if (!eulerzone) constreuler(0);
680 14 : if (k < lg(eulerzone)) return itor(gel(eulerzone,k), prec);
681 0 : p = eulerprec(n);
682 0 : B = eulerreal_using_lfun4(n, minss(p, prec));
683 0 : if (p < prec) B = itor(roundr(B), prec);
684 0 : return gc_leaf(av, B);
685 : }
|
