Line data Source code
1 : /* Copyright (C) 2015 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. It is distributed in the hope that it will be useful, but WITHOUT
8 : ANY WARRANTY WHATSOEVER.
9 :
10 : Check the License for details. You should have received a copy of it, along
11 : with the package; see the file 'COPYING'. If not, write to the Free Software
12 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
13 :
14 : /********************************************************************/
15 : /** **/
16 : /** Dirichlet series through Euler product **/
17 : /** **/
18 : /********************************************************************/
19 : #include "pari.h"
20 : #include "paripriv.h"
21 :
22 : static void
23 28 : err_direuler(GEN x)
24 28 : { pari_err_DOMAIN("direuler","constant term","!=", gen_1,x); }
25 :
26 : /* s = t_POL (tolerate t_SER of valuation 0) of constant term = 1
27 : * d = minimal such that p^d > X
28 : * V indexed by 1..X will contain the a_n
29 : * v[1..n] contains the indices nj such that V[nj] != 0 */
30 : static long
31 20468 : dirmuleuler_small(GEN V, GEN v, long n, ulong p, GEN s, long d)
32 : {
33 20468 : long i, j, m = n, D = minss(d+2, lg(s));
34 20468 : ulong q = 1, X = lg(V)-1;
35 :
36 70245 : for (i = 3, q = p; i < D; i++, q *= p) /* q*p does not overflow */
37 : {
38 49777 : GEN aq = gel(s,i);
39 49777 : if (gequal0(aq)) continue;
40 : /* j = 1 */
41 42084 : gel(V,q) = aq;
42 42084 : v[++n] = q;
43 1766674 : for (j = 2; j <= m; j++)
44 : {
45 1724590 : ulong nj = umuluu_le(uel(v,j), q, X);
46 1724590 : if (!nj) continue;
47 133910 : gel(V,nj) = gmul(aq, gel(V,v[j]));
48 133910 : v[++n] = nj;
49 : }
50 : }
51 20468 : return n;
52 : }
53 :
54 : /* ap != 0 for efficiency, p > sqrt(X) */
55 : static void
56 150010 : dirmuleuler_large(GEN V, ulong p, GEN ap)
57 : {
58 150010 : long j, jp, X = lg(V)-1;
59 150010 : gel(V,p) = ap;
60 599879 : for (j = 2, jp = 2*p; jp <= X; j++, jp += p) gel(V,jp) = gmul(ap, gel(V,j));
61 150010 : }
62 :
63 : static ulong
64 9107 : direulertou(GEN a, GEN fl(GEN))
65 : {
66 9107 : if (typ(a) != t_INT)
67 : {
68 49 : a = fl(a);
69 28 : if (typ(a) != t_INT) pari_err_TYPE("direuler", a);
70 : }
71 9086 : return signe(a)<=0 ? 0: itou(a);
72 : }
73 :
74 : static GEN
75 3556 : direuler_Sbad(GEN V, GEN v, GEN Sbad, ulong *n)
76 : {
77 3556 : long i, l = lg(Sbad);
78 3556 : ulong X = lg(V)-1;
79 3556 : GEN pbad = gen_1;
80 9275 : for (i = 1; i < l; i++)
81 : {
82 5754 : GEN ai = gel(Sbad,i);
83 : ulong q;
84 5754 : if (typ(ai) != t_VEC || lg(ai) != 3)
85 14 : pari_err_TYPE("direuler [bad primes]",ai);
86 5740 : q = gtou(gel(ai,1));
87 5733 : if (q <= X)
88 : {
89 4627 : long d = ulogint(X, q) + 1;
90 4627 : GEN s = direuler_factor(gel(ai,2), d);
91 4613 : *n = dirmuleuler_small(V, v, *n, q, s, d);
92 4613 : pbad = muliu(pbad, q);
93 : }
94 : }
95 3521 : return pbad;
96 : }
97 :
98 : GEN
99 504 : direuler_bad(void *E, GEN (*eval)(void *,GEN,long), GEN a,GEN b,GEN c, GEN Sbad)
100 : {
101 : ulong au, bu, X, sqrtX, n, p;
102 504 : pari_sp av0 = avma;
103 : GEN gp, v, V;
104 : forprime_t T;
105 504 : au = direulertou(a, gceil);
106 497 : bu = direulertou(b, gfloor);
107 490 : X = c ? direulertou(c, gfloor): bu;
108 483 : if (X == 0) return cgetg(1,t_VEC);
109 476 : if (bu > X) bu = X;
110 476 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
111 462 : v = vecsmall_ei(X, 1);
112 462 : V = vec_ei(X, 1);
113 462 : n = 1;
114 462 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
115 427 : p = 1; gp = cgetipos(3); sqrtX = usqrt(X);
116 2758 : while (p <= sqrtX && (p = u_forprime_next(&T)))
117 2352 : if (!Sbad || umodiu(Sbad, p))
118 : {
119 2247 : long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
120 : GEN s;
121 2247 : gp[2] = p; s = eval(E, gp, d);
122 2226 : n = dirmuleuler_small(V, v, n, p, s, d);
123 : }
124 55391 : while ((p = u_forprime_next(&T))) /* sqrt(X) < p <= X */
125 54985 : if (!Sbad || umodiu(Sbad, p))
126 : {
127 : GEN s;
128 54978 : gp[2] = p; s = eval(E, gp, 2); /* s either t_POL or t_SER of val 0 */
129 54978 : if (lg(s) > 3 && !gequal0(gel(s,3)))
130 23072 : dirmuleuler_large(V, p, gel(s,3));
131 : }
132 406 : return gerepilecopy(av0,V);
133 : }
134 :
135 : /* return a t_SER or a truncated t_POL to precision n */
136 : GEN
137 61852 : direuler_factor(GEN s, long n)
138 : {
139 61852 : long t = typ(s);
140 61852 : if (is_scalar_t(t))
141 : {
142 31612 : if (!gequal1(s)) err_direuler(s);
143 31598 : return scalarpol_shallow(s,0);
144 : }
145 30240 : switch(t)
146 : {
147 5467 : case t_POL: break; /* no need to RgXn_red */
148 24458 : case t_RFRAC:
149 : {
150 24458 : GEN p = gel(s,1), q = gel(s,2);
151 24458 : q = RgXn_red_shallow(q,n);
152 24458 : s = RgXn_inv(q, n);
153 24458 : if (typ(p) == t_POL && varn(p) == varn(q))
154 : {
155 28 : p = RgXn_red_shallow(p, n);
156 28 : s = RgXn_mul(s, p, n);
157 : }
158 : else
159 24430 : if (!gequal1(p)) s = RgX_Rg_mul(s, p);
160 24458 : if (!signe(s) || !gequal1(gel(s,2))) err_direuler(s);
161 24444 : break;
162 : }
163 308 : case t_SER:
164 308 : if (!signe(s) || valp(s) || !gequal1(gel(s,2))) err_direuler(s);
165 308 : break;
166 7 : default: pari_err_TYPE("direuler", s);
167 : }
168 30219 : return s;
169 : }
170 :
171 : struct eval_bad
172 : {
173 : void *E;
174 : GEN (*eval)(void *, GEN);
175 : };
176 : static GEN
177 420 : eval_bad(void *E, GEN p, long n)
178 : {
179 420 : struct eval_bad *d = (struct eval_bad*) E;
180 420 : return direuler_factor(d->eval(d->E, p), n);
181 : }
182 : GEN
183 133 : direuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, GEN c)
184 : {
185 : struct eval_bad d;
186 133 : d.E= E; d.eval = eval;
187 133 : return direuler_bad((void*)&d, eval_bad, a, b, c, NULL);
188 : }
189 :
190 : static GEN
191 26614 : primelist(forprime_t *T, GEN Sbad, long n, long *running)
192 : {
193 26614 : GEN P = cgetg(n+1, t_VECSMALL);
194 : long i, j;
195 242816 : for (i = 1, j = 1; i <= n; i++)
196 : {
197 219863 : ulong p = u_forprime_next(T);
198 219863 : if (!p) { *running = 0; break; }
199 216202 : if (Sbad && umodiu(Sbad, p)==0) continue;
200 211701 : uel(P,j++) = p;
201 : }
202 26614 : setlg(P, j);
203 26614 : return P;
204 : }
205 :
206 : GEN
207 3668 : pardireuler(GEN worker, GEN a, GEN b, GEN c, GEN Sbad)
208 : {
209 : ulong au, bu, X, sqrtX, n, snX, nX;
210 3668 : pari_sp av0 = avma;
211 : GEN v, V;
212 : forprime_t T;
213 : struct pari_mt pt;
214 3668 : long running = 1, pending = 0;
215 3668 : au = direulertou(a, gceil);
216 3668 : bu = direulertou(b, gfloor);
217 3668 : X = c ? direulertou(c, gfloor): bu;
218 3668 : if (X == 0) return cgetg(1,t_VEC);
219 3668 : if (bu > X) bu = X;
220 3668 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
221 3661 : v = vecsmall_ei(X, 1);
222 3661 : V = vec_ei(X, 1);
223 3661 : n = 1;
224 3661 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
225 3661 : sqrtX = usqrt(X); snX = uprimepi(sqrtX); nX = uprimepi(X);
226 3661 : if (snX)
227 : {
228 3647 : GEN P = primelist(&T, Sbad, snX, &running);
229 3647 : GEN R = gel(closure_callgenvec(worker, mkvec2(P, utoi(X))), 2);
230 3647 : long i, l = lg(P);
231 17276 : for (i = 1; i < l; i++)
232 : {
233 13629 : GEN s = gel(R,i);
234 13629 : n = dirmuleuler_small(V, v, n, uel(P,i), s, lg(s));
235 : }
236 14 : } else snX = 1;
237 3661 : mt_queue_start_lim(&pt, worker, (nX+snX-1)/snX);
238 29615 : while (running || pending)
239 : {
240 : GEN done;
241 25954 : GEN P = running? primelist(&T, Sbad, snX, &running): NULL;
242 25954 : mt_queue_submit(&pt, 0, P ? mkvec2(P, utoi(X)): NULL);
243 25954 : done = mt_queue_get(&pt, NULL, &pending);
244 25954 : if (done)
245 : {
246 22967 : GEN P = gel(done,1), R = gel(done,2);
247 22967 : long j, l = lg(P);
248 221039 : for (j=1; j<l; j++)
249 : {
250 198072 : GEN F = gel(R,j);
251 198072 : if (degpol(F) && !gequal0(gel(F,3)))
252 126938 : dirmuleuler_large(V, uel(P,j), gel(F,3));
253 : }
254 : }
255 : }
256 3661 : mt_queue_end(&pt);
257 3661 : return gerepilecopy(av0,V);
258 : }
259 :
260 : /********************************************************************/
261 : /** **/
262 : /** DIRPOWERS and DIRPOWERSSUM **/
263 : /** **/
264 : /********************************************************************/
265 :
266 : /* [1^B,...,N^B] */
267 : GEN
268 259 : vecpowuu(long N, ulong B)
269 : {
270 : GEN v;
271 : long p, i;
272 : forprime_t T;
273 :
274 259 : if (B <= 2)
275 : {
276 77 : if (!B) return const_vec(N,gen_1);
277 70 : v = cgetg(N+1, t_VEC); if (N == 0) return v;
278 70 : gel(v,1) = gen_1;
279 70 : if (B == 1)
280 3815 : for (i = 2; i <= N; i++) gel(v,i) = utoipos(i);
281 : else
282 644 : for (i = 2; i <= N; i++) gel(v,i) = sqru(i);
283 70 : return v;
284 : }
285 182 : v = const_vec(N, NULL);
286 182 : u_forprime_init(&T, 3, N);
287 3934 : while ((p = u_forprime_next(&T)))
288 : {
289 : long m, pk, oldpk;
290 3752 : gel(v,p) = powuu(p, B);
291 8365 : for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
292 : {
293 4613 : if (pk != p) gel(v,pk) = mulii(gel(v,oldpk), gel(v,p));
294 20797 : for (m = N/pk; m > 1; m--)
295 16184 : if (gel(v,m) && m%p) gel(v, m*pk) = mulii(gel(v,m), gel(v,pk));
296 : }
297 : }
298 182 : gel(v,1) = gen_1;
299 7602 : for (i = 2; i <= N; i+=2)
300 : {
301 7420 : long vi = vals(i);
302 7420 : gel(v,i) = shifti(gel(v,i >> vi), B * vi);
303 : }
304 182 : return v;
305 : }
306 : /* [1^B,...,N^B] */
307 : GEN
308 2443 : vecpowug(long N, GEN B, long prec)
309 : {
310 2443 : GEN v, pow = NULL;
311 2443 : long p, precp = 2, eB, prec0;
312 : forprime_t T;
313 2443 : if (N == 1) return mkvec(gen_1);
314 2443 : eB = gexpo(B);
315 2443 : prec0 = eB < 5? prec: prec + nbits2extraprec(eB);
316 2443 : u_forprime_init(&T, 2, N);
317 2443 : v = const_vec(N, NULL);
318 2443 : gel(v,1) = gen_1;
319 529634 : while ((p = u_forprime_next(&T)))
320 : {
321 : long m, pk, oldpk;
322 527191 : if (!pow)
323 2443 : pow = gpow(utor(p,prec0), B, prec);
324 : else
325 : {
326 524748 : GEN t = gpow(divru(utor(p,prec0), precp), B, prec);
327 524748 : pow = gmul(pow, t); /* (p / precp)^B * precp^B */
328 : }
329 527191 : precp = p;
330 527191 : gel(v,p) = pow; /* p^B */
331 527191 : if (prec0 != prec) gel(v,p) = gprec_wtrunc(gel(v,p), prec);
332 1093715 : for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
333 : {
334 566524 : if (pk != p) gel(v,pk) = gmul(gel(v,oldpk), gel(v,p));
335 13292209 : for (m = N/pk; m > 1; m--)
336 12725685 : if (gel(v,m) && m%p) gel(v, m*pk) = gmul(gel(v,m), gel(v,pk));
337 : }
338 : }
339 2443 : return v;
340 : }
341 :
342 : GEN
343 266 : dirpowers(long n, GEN x, long prec)
344 : {
345 266 : pari_sp av = avma;
346 : GEN v;
347 266 : if (n <= 0) return cgetg(1, t_VEC);
348 252 : if (typ(x) == t_INT && lgefint(x) <= 3 && signe(x) >= 0)
349 7 : {
350 35 : ulong B = itou(x);
351 35 : v = vecpowuu(n, B);
352 35 : if (B <= 2) return v;
353 : }
354 217 : else v = vecpowug(n, x, prec);
355 224 : return gerepilecopy(av, v);
356 : }
357 :
358 : /* P = prime divisors of (squarefree) n */
359 : static GEN
360 296824 : smallfact(ulong n, GEN P, ulong sq, GEN V)
361 : {
362 296824 : long i, l = lg(P);
363 : ulong p;
364 : GEN c;
365 296824 : if (l == 1) return gen_1;
366 288480 : p = uel(P,l - 1); if (p > sq) return NULL;
367 49861 : c = gel(V, p);
368 66885 : for (i = l-2; i > 0; i--)
369 : {
370 60088 : n /= p; if (n <= sq) return gmul(c, gel(V,n));
371 17024 : p = uel(P, i); c = gmul(c, gel(V,p));
372 : }
373 6797 : return c;
374 : }
375 : /* sum_{n <= N} n^s. */
376 : GEN
377 8372 : dirpowerssum(ulong N, GEN s, long prec)
378 : {
379 8372 : const ulong step = 2048;
380 8372 : pari_sp av = avma, av2;
381 : GEN P, V, W, F, c2, c3, c6, c12, c123, c1234, tmp, S;
382 : forprime_t T;
383 : ulong x1, n, sq, p, precp;
384 : long prec0;
385 :
386 8372 : if (!N) return gen_0;
387 8372 : if (N < 9UL) return gerepileupto(av, RgV_sum(dirpowers(N, s, prec)));
388 8344 : sq = usqrt(N);
389 8344 : V = cgetg(sq+1, t_VEC);
390 8344 : W = cgetg(sq+1, t_VEC);
391 8344 : F = cgetg(sq+1, t_VEC);
392 8344 : prec0 = prec + EXTRAPRECWORD;
393 8344 : s = gprec_w(s, prec0);
394 8344 : gel(V,1) = gel(W,1) = gel(F,1) = gen_1;
395 46037 : for (n = 2; n <= sq; n++)
396 : {
397 37693 : GEN t = divru(utor(n, prec0), n-1);
398 37693 : gel(V,n) = gmul(gel(V,n-1), gpow(t, s, prec0));
399 37693 : gel(W,n) = gadd(gel(W,n-1), gel(V,n));
400 37693 : gel(F,n) = gadd(gel(F,n-1), gsqr(gel(V,n)));
401 : }
402 8344 : c2 = gel(V,2); c3 = gel(V,3); c6 = gmul(c2, c3);
403 8344 : c12 = gaddgs(c2, 1);
404 8344 : c123 = gadd(c12, c3);
405 8344 : c1234 = gadd(gel(F,2), gadd(c2,c3));
406 8344 : precp = 0; tmp = NULL; S = gen_0;
407 8344 : u_forprime_init(&T, sq + 1, N);
408 8344 : av2 = avma;
409 139865 : while ((p = u_forprime_next(&T)))
410 : {
411 131521 : GEN t = utor(p, prec0);
412 131521 : if (precp) t = divru(t, precp);
413 131521 : t = gpow(t, s, prec0);
414 131521 : tmp = precp? gmul(tmp, t): t; /* p^s */
415 131521 : S = gadd(S, gmul(gel(W, N / p), tmp));
416 131521 : precp = p;
417 131521 : if ((p & 0x1ff) == 1) S = gerepileupto(av2, S);
418 : }
419 8344 : P = mkvecsmall2(2, 3);
420 8344 : av2 = avma;
421 8344 : for(x1 = 1;; x1 += step)
422 329 : { /* beware overflow, fuse last two bins (avoid a tiny remainder) */
423 8673 : ulong j, lv, x2 = (N >= 2*step && N - 2*step >= x1)? x1-1 + step: N;
424 8673 : GEN v = vecfactorsquarefreeu_coprime(x1, x2, P);
425 8673 : lv = lg(v);
426 971963 : for (j = 1; j < lv; j++) if (gel(v,j))
427 : {
428 296824 : ulong d = x1-1 + j; /* squarefree, coprime to 6 */
429 : ulong q;
430 296824 : tmp = smallfact(d, gel(v,j), sq, V);
431 296824 : if (!tmp) continue;
432 58205 : q = N / d;
433 58205 : switch(q = N / d)
434 : {
435 15925 : case 1: break;
436 6748 : case 2: tmp = gmul(tmp, c12); break;
437 3864 : case 3: tmp = gmul(tmp, c123); break;
438 5168 : case 4:
439 5168 : case 5: tmp = gmul(tmp, c1234); break;
440 26500 : default:
441 : {
442 26500 : GEN a = gel(F, usqrt(q)), b = gel(F, usqrt(q / 2));
443 26500 : GEN c = gel(F, usqrt(q / 3)), d = gel(F, usqrt(q / 6));
444 26500 : tmp = gmul(tmp, gadd(gadd(a, gmul(c2, b)),
445 : gadd(gmul(c3, c), gmul(c6, d))));
446 : }
447 : }
448 58205 : S = gadd(S, tmp);
449 : }
450 8673 : if (x2 == N) break;
451 329 : S = gerepileupto(av2, S);
452 : }
453 8344 : return gerepilecopy(av, S);
454 : }
455 : GEN
456 35 : dirpowerssum0(GEN N, GEN s, long prec)
457 : {
458 35 : if (typ(N) != t_INT) pari_err_TYPE("dirpowerssum", N);
459 28 : if (signe(N) <= 0) return gen_0;
460 14 : return dirpowerssum(itou(N), s, prec);
461 : }
|
