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; 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 : /********************************************************************/
16 : /** **/
17 : /** Dirichlet series through Euler product **/
18 : /** **/
19 : /********************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : static void
24 28 : err_direuler(GEN x)
25 28 : { pari_err_DOMAIN("direuler","constant term","!=", gen_1,x); }
26 :
27 : /* s = t_POL (tolerate t_SER of valuation 0) of constant term = 1
28 : * d = minimal such that p^d > X
29 : * V indexed by 1..X will contain the a_n
30 : * v[1..n] contains the indices nj such that V[nj] != 0 */
31 : static long
32 28700 : dirmuleuler_small(GEN V, GEN v, long n, ulong p, GEN s, long d)
33 : {
34 28700 : long i, j, m = n, D = minss(d+2, lg(s));
35 28700 : ulong q = 1, X = lg(V)-1;
36 :
37 94724 : for (i = 3, q = p; i < D; i++, q *= p) /* q*p does not overflow */
38 : {
39 66024 : GEN aq = gel(s,i);
40 66024 : if (gequal0(aq)) continue;
41 : /* j = 1 */
42 53753 : gel(V,q) = aq;
43 53753 : v[++n] = q;
44 3268013 : for (j = 2; j <= m; j++)
45 : {
46 3214260 : ulong nj = umuluu_le(uel(v,j), q, X);
47 3214260 : if (!nj) continue;
48 192017 : gel(V,nj) = gmul(aq, gel(V,v[j]));
49 192017 : v[++n] = nj;
50 : }
51 : }
52 28700 : return n;
53 : }
54 :
55 : /* ap != 0 for efficiency, p > sqrt(X) */
56 : static void
57 308798 : dirmuleuler_large(GEN V, ulong p, GEN ap)
58 : {
59 308798 : long j, jp, X = lg(V)-1;
60 308798 : gel(V,p) = ap;
61 1506547 : for (j = 2, jp = 2*p; jp <= X; j++, jp += p) gel(V,jp) = gmul(ap, gel(V,j));
62 308798 : }
63 :
64 : static ulong
65 10269 : direulertou(GEN a, GEN fl(GEN))
66 : {
67 10269 : if (typ(a) != t_INT)
68 : {
69 49 : a = fl(a);
70 28 : if (typ(a) != t_INT) pari_err_TYPE("direuler", a);
71 : }
72 10248 : return signe(a)<=0 ? 0: itou(a);
73 : }
74 :
75 : static GEN
76 3724 : direuler_Sbad(GEN V, GEN v, GEN Sbad, ulong *n)
77 : {
78 3724 : long i, l = lg(Sbad);
79 3724 : ulong X = lg(V)-1;
80 3724 : GEN pbad = gen_1;
81 9646 : for (i = 1; i < l; i++)
82 : {
83 5957 : GEN ai = gel(Sbad,i);
84 : ulong q;
85 5957 : if (typ(ai) != t_VEC || lg(ai) != 3)
86 14 : pari_err_TYPE("direuler [bad primes]",ai);
87 5943 : q = gtou(gel(ai,1));
88 5936 : if (q <= X)
89 : {
90 4809 : long d = ulogint(X, q) + 1;
91 4809 : GEN s = direuler_factor(gel(ai,2), d);
92 4795 : *n = dirmuleuler_small(V, v, *n, q, s, d);
93 4795 : pbad = muliu(pbad, q);
94 : }
95 : }
96 3689 : return pbad;
97 : }
98 :
99 : GEN
100 672 : direuler_bad(void *E, GEN (*eval)(void *,GEN,long), GEN a,GEN b,GEN c, GEN Sbad)
101 : {
102 : ulong au, bu, X, sqrtX, n, p;
103 672 : pari_sp av0 = avma;
104 : GEN gp, v, V;
105 : forprime_t T;
106 672 : au = direulertou(a, gceil);
107 665 : bu = direulertou(b, gfloor);
108 658 : X = c ? direulertou(c, gfloor): bu;
109 651 : if (X == 0) return cgetg(1,t_VEC);
110 644 : if (bu > X) bu = X;
111 644 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
112 630 : v = vecsmall_ei(X, 1);
113 630 : V = vec_ei(X, 1);
114 630 : n = 1;
115 630 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
116 595 : p = 1; gp = cgetipos(3); sqrtX = usqrt(X);
117 8316 : while (p <= sqrtX && (p = u_forprime_next(&T)))
118 7742 : if (!Sbad || umodiu(Sbad, p))
119 : {
120 7637 : long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
121 : GEN s;
122 7637 : gp[2] = p; s = eval(E, gp, d);
123 7616 : n = dirmuleuler_small(V, v, n, p, s, d);
124 : }
125 740201 : while ((p = u_forprime_next(&T))) /* sqrt(X) < p <= X */
126 739627 : if (!Sbad || umodiu(Sbad, p))
127 : {
128 : GEN s;
129 739620 : gp[2] = p; s = eval(E, gp, 2); /* s either t_POL or t_SER of val 0 */
130 739620 : if (lg(s) > 3 && !gequal0(gel(s,3)))
131 139153 : dirmuleuler_large(V, p, gel(s,3));
132 : }
133 574 : return gc_GEN(av0,V);
134 : }
135 :
136 : /* return a t_SER or a truncated t_POL to precision n */
137 : GEN
138 752066 : direuler_factor(GEN s, long n)
139 : {
140 752066 : long t = typ(s);
141 752066 : if (is_scalar_t(t))
142 : {
143 33194 : if (!gequal1(s)) err_direuler(s);
144 33180 : return scalarpol_shallow(s,0);
145 : }
146 718872 : switch(t)
147 : {
148 5712 : case t_POL: break; /* no need to RgXn_red */
149 712845 : case t_RFRAC:
150 : {
151 712845 : GEN p = gel(s,1), q = gel(s,2);
152 712845 : q = RgXn_red_shallow(q,n);
153 712845 : s = RgXn_inv(q, n);
154 712845 : if (typ(p) == t_POL && varn(p) == varn(q))
155 : {
156 28 : p = RgXn_red_shallow(p, n);
157 28 : s = RgXn_mul(s, p, n);
158 : }
159 : else
160 712817 : if (!gequal1(p)) s = RgX_Rg_mul(s, p);
161 712845 : if (!signe(s) || !gequal1(gel(s,2))) err_direuler(s);
162 712831 : break;
163 : }
164 308 : case t_SER:
165 308 : if (!signe(s) || valser(s) || !gequal1(gel(s,2))) err_direuler(s);
166 308 : break;
167 7 : default: pari_err_TYPE("direuler", s);
168 : }
169 718851 : return s;
170 : }
171 :
172 : struct eval_bad
173 : {
174 : void *E;
175 : GEN (*eval)(void *, GEN);
176 : };
177 : static GEN
178 688303 : eval_bad(void *E, GEN p, long n)
179 : {
180 688303 : struct eval_bad *d = (struct eval_bad*) E;
181 688303 : return direuler_factor(d->eval(d->E, p), n);
182 : }
183 : GEN
184 301 : direuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, GEN c)
185 : {
186 : struct eval_bad d;
187 301 : d.E= E; d.eval = eval;
188 301 : return direuler_bad((void*)&d, eval_bad, a, b, c, NULL);
189 : }
190 :
191 : static GEN
192 31157 : primelist(forprime_t *T, GEN Sbad, long n, long *running)
193 : {
194 31157 : GEN P = cgetg(n+1, t_VECSMALL);
195 : long i, j;
196 302631 : for (i = 1, j = 1; i <= n; i++)
197 : {
198 275548 : ulong p = u_forprime_next(T);
199 275548 : if (!p) { *running = 0; break; }
200 271474 : if (Sbad && umodiu(Sbad, p)==0) continue;
201 266791 : uel(P,j++) = p;
202 : }
203 31157 : setlg(P, j);
204 31157 : return P;
205 : }
206 :
207 : GEN
208 4081 : pardireuler(GEN worker, GEN a, GEN b, GEN c, GEN Sbad)
209 : {
210 : ulong au, bu, X, sqrtX, n, snX, nX;
211 4081 : pari_sp av0 = avma;
212 : GEN v, V;
213 : forprime_t T;
214 : struct pari_mt pt;
215 4081 : long running = 1, pending = 0;
216 4081 : au = direulertou(a, gceil);
217 4081 : bu = direulertou(b, gfloor);
218 4081 : X = c ? direulertou(c, gfloor): bu;
219 4081 : if (X == 0) return cgetg(1,t_VEC);
220 4081 : if (bu > X) bu = X;
221 4081 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
222 4074 : v = vecsmall_ei(X, 1);
223 4074 : V = vec_ei(X, 1);
224 4074 : n = 1;
225 4074 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
226 4074 : sqrtX = usqrt(X); snX = uprimepi(sqrtX); nX = uprimepi(X);
227 4074 : if (snX)
228 : {
229 4060 : GEN P = primelist(&T, Sbad, snX, &running);
230 4060 : GEN R = gel(closure_callgenvec(worker, mkvec2(P, utoi(X))), 2);
231 4060 : long i, l = lg(P);
232 20349 : for (i = 1; i < l; i++)
233 : {
234 16289 : GEN s = gel(R,i);
235 16289 : n = dirmuleuler_small(V, v, n, uel(P,i), s, lg(s));
236 : }
237 14 : } else snX = 1;
238 4074 : mt_queue_start_lim(&pt, worker, (nX+snX-1)/snX);
239 34212 : while (running || pending)
240 : {
241 : GEN done;
242 30138 : GEN P = running? primelist(&T, Sbad, snX, &running): NULL;
243 30138 : mt_queue_submit(&pt, 0, P ? mkvec2(P, utoi(X)): NULL);
244 30138 : done = mt_queue_get(&pt, NULL, &pending);
245 30138 : if (done)
246 : {
247 27097 : GEN P = gel(done,1), R = gel(done,2);
248 27097 : long j, l = lg(P);
249 277599 : for (j=1; j<l; j++)
250 : {
251 250502 : GEN F = gel(R,j);
252 250502 : if (degpol(F) && !gequal0(gel(F,3)))
253 169645 : dirmuleuler_large(V, uel(P,j), gel(F,3));
254 : }
255 : }
256 : }
257 4074 : mt_queue_end(&pt);
258 4074 : return gc_GEN(av0,V);
259 : }
260 :
261 : /********************************************************************/
262 : /** **/
263 : /** DIRPOWERS and DIRPOWERSSUM **/
264 : /** **/
265 : /********************************************************************/
266 :
267 : /* [1^B,...,N^B] */
268 : GEN
269 686 : vecpowuu(long N, ulong B)
270 : {
271 : GEN v;
272 : long p, i;
273 : forprime_t T;
274 :
275 686 : if (B <= 8000)
276 : {
277 686 : if (!B) return const_vec(N,gen_1);
278 679 : v = cgetg(N+1, t_VEC); if (N == 0) return v;
279 679 : gel(v,1) = gen_1;
280 679 : if (B == 1)
281 92736 : for (i = 2; i <= N; i++) gel(v,i) = utoipos(i);
282 469 : else if (B == 2)
283 : {
284 : ulong o, s;
285 273 : if (N & HIGHMASK)
286 0 : for (i = 2, o = 3; i <= N; i++, o += 2)
287 0 : gel(v,i) = addiu(gel(v,i-1), o);
288 : else
289 31073 : for (i = 2, s = 1, o = 3; i <= N; i++, s += o, o += 2)
290 30800 : gel(v,i) = utoipos(s + o);
291 : }
292 196 : else if (B == 3)
293 840 : for (i = 2; i <= N; i++) gel(v,i) = powuu(i, B);
294 : else
295 : {
296 182 : long k, Bk, e = expu(N);
297 7553 : for (i = 3; i <= N; i += 2) gel(v,i) = powuu(i, B);
298 1239 : for (k = 1; k <= e; k++)
299 : {
300 1057 : N >>= 1; Bk = B * k;
301 8498 : for (i = 1; i <= N; i += 2) gel(v, i << k) = shifti(gel(v, i), Bk);
302 : }
303 : }
304 679 : return v;
305 : }
306 0 : v = const_vec(N, NULL);
307 0 : u_forprime_init(&T, 3, N);
308 0 : while ((p = u_forprime_next(&T)))
309 : {
310 : long m, pk, oldpk;
311 0 : gel(v,p) = powuu(p, B);
312 0 : for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
313 : {
314 0 : if (pk != p) gel(v,pk) = mulii(gel(v,oldpk), gel(v,p));
315 0 : for (m = N/pk; m > 1; m--)
316 0 : if (gel(v,m) && m%p) gel(v, m*pk) = mulii(gel(v,m), gel(v,pk));
317 : }
318 : }
319 0 : gel(v,1) = gen_1;
320 0 : for (i = 2; i <= N; i+=2)
321 : {
322 0 : long vi = vals(i);
323 0 : gel(v,i) = shifti(gel(v,i >> vi), B * vi);
324 : }
325 0 : return v;
326 : }
327 :
328 : /* does x^s require log(x) ? */
329 : static long
330 22783 : get_needlog(GEN s)
331 : {
332 22783 : switch(typ(s))
333 : {
334 10157 : case t_REAL: return 2; /* yes but not powcx */
335 9492 : case t_COMPLEX: return 1; /* yes using powcx */
336 3134 : default: return 0; /* no */
337 : }
338 : }
339 : /* [1^B,...,N^B] */
340 : GEN
341 21922 : vecpowug(long N, GEN B, long prec)
342 : {
343 21922 : GEN v, logp = NULL;
344 21922 : long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
345 21922 : long p, precp = 2, prec0, prec1, needlog;
346 : forprime_t T;
347 21922 : if (N == 1) return mkvec(gen_1);
348 21901 : if (typ(B) == t_INT && lgefint(B) <= 3 && signe(B) >= 0)
349 168 : return vecpowuu(N, itou(B));
350 21733 : needlog = get_needlog(B);
351 21733 : prec1 = prec0 = prec;
352 21733 : if (needlog == 1) prec1 = powcx_prec(log2((double)N), B, prec);
353 21733 : u_forprime_init(&T, 2, N);
354 21733 : v = const_vec(N, NULL);
355 21733 : gel(v,1) = gen_1;
356 1685820 : while ((p = u_forprime_next(&T)))
357 : {
358 : long m, pk, oldpk;
359 : GEN u;
360 1664087 : gp[2] = p;
361 1664087 : if (needlog)
362 : {
363 205983 : if (!logp)
364 37240 : logp = logr_abs(utor(p, prec1));
365 : else
366 : { /* Assuming p and precp are odd,
367 : * log p = log(precp) + 2 atanh((p - precp) / (p + precp)) */
368 168743 : ulong a = p >> 1, b = precp >> 1; /* p = 2a + 1, precp = 2b + 1 */
369 168743 : GEN z = atanhuu(a - b, a + b + 1, prec1); /* avoid overflow */
370 168743 : shiftr_inplace(z, 1); logp = addrr(logp, z);
371 : }
372 92829 : u = needlog == 1? powcx(gp, logp, B, prec0)
373 205983 : : mpexp(gmul(B, logp));
374 205983 : if (p == 2) logp = NULL; /* reset: precp must be odd */
375 : }
376 : else
377 1458104 : u = gpow(gp, B, prec0);
378 1664087 : precp = p;
379 1664087 : gel(v,p) = u; /* p^B */
380 1664087 : if (prec0 != prec) gel(v,p) = gprec_wtrunc(gel(v,p), prec);
381 3492662 : for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
382 : {
383 1828575 : if (pk != p) gel(v,pk) = gmul(gel(v,oldpk), gel(v,p));
384 47133627 : for (m = N/pk; m > 1; m--)
385 45305052 : if (gel(v,m) && m%p) gel(v, m*pk) = gmul(gel(v,m), gel(v,pk));
386 : }
387 : }
388 21733 : return v;
389 : }
390 :
391 : GEN
392 665 : dirpowers(long n, GEN x, long prec)
393 : {
394 : pari_sp av;
395 : GEN v;
396 665 : if (n <= 0) return cgetg(1, t_VEC);
397 651 : av = avma; v = vecpowug(n, x, prec);
398 651 : if (typ(x) == t_INT && lgefint(x) <= 3 && signe(x) >= 0 && cmpiu(x, 2) <= 0)
399 133 : return v;
400 518 : return gc_GEN(av, v);
401 : }
402 :
403 : /* f is a totally multiplicative function of modulus 0 or 1
404 : * (essentially a Dirichlet character). Compute simultaneously
405 : * sum_{0 < n <= N} f(n)n^s and sum_{0 < n <= N} f(n)n^{-1-conj(s)}
406 : * Warning: s is conjugated, but not f. Main application for Riemann-Siegel,
407 : * where we need R(chi,s) and conj(R(chi,1-conj(s))). */
408 :
409 : static GEN
410 3234 : vecmulsqlv(GEN Q, GEN V)
411 : {
412 : long l, i;
413 : GEN W;
414 3234 : if (typ(V) != t_VEC) return RgV_Rg_mul(Q, V);
415 1239 : l = lg(Q); W = cgetg(l, t_VEC);
416 62244 : for (i = 1; i < l; i++) gel(W, i) = vecmul(gel(Q, i), V);
417 1239 : return W;
418 : }
419 :
420 : /* P = prime divisors of (squarefree) n, V[i] = i^s for i <= sq.
421 : * Return NULL if n is not sq-smooth, else f(n)n^s */
422 : static GEN
423 18216536 : smallfact(ulong n, GEN P, ulong sq, GEN V)
424 : {
425 : long i, l;
426 : ulong p, m, o;
427 : GEN c;
428 18216536 : if (n <= sq) return gel(V,n);
429 18167557 : l = lg(P); m = p = uel(P, l-1); if (p > sq) return NULL;
430 3847551 : for (i = l-2; i > 1; i--, m = o) { p = uel(P,i); o = m*p; if (o > sq) break; }
431 3814244 : c = gel(V,m); n /= m; /* m <= sq, o = m * p > sq */
432 3814244 : if (n > sq) { c = vecmul(c, gel(V,p)); n /= p; }
433 3806016 : return vecmul(c, gel(V,n));
434 : }
435 :
436 : static GEN
437 973 : _Qtor(GEN x, long prec)
438 973 : { return typ(x) == t_FRAC? fractor(x, prec): x; }
439 : static GEN
440 2457 : Qtor(GEN x, long prec)
441 : {
442 2457 : long tx = typ(x);
443 3430 : if (tx == t_VEC || tx == t_COL) pari_APPLY_same(_Qtor(gel(x, i), prec));
444 1589 : return tx == t_FRAC? fractor(x, prec): x;
445 : }
446 :
447 : static GEN
448 238 : vecf(long N, void *E, GEN (*f)(void *, ulong, long), long prec)
449 : {
450 : GEN v;
451 : long n;
452 238 : if (!f) return NULL;
453 140 : v = cgetg(N + 1, t_VEC);
454 33481 : for (n = 1; n <= N; n++) gel(v,n) = f(E, n, prec);
455 140 : return v;
456 : }
457 :
458 : /* Here #V = #F > 0 is small. Analogous to dotproduct, with following
459 : * semantic differences: uses V[1] = 1; V has scalar values but F may have
460 : * vector values */
461 : static GEN
462 301 : naivedirpowerssum(GEN V, GEN F, long prec)
463 : {
464 : GEN S;
465 301 : if (!F) S = RgV_sum(V);
466 : else
467 : {
468 189 : long n, l = lg(V);
469 189 : S = gel(F,1); /* V[1] = 1 */
470 34447 : for (n = 2; n < l; n++) S = gadd(S, gmul(gel(V, n), gel(F, n)));
471 : }
472 301 : return Qtor(S, prec);
473 : }
474 :
475 : /* set *p0 and *p1 to 0 and 1 in the algebra where f takes its values */
476 : static int
477 1092 : mk01(void *E, GEN (*f)(void*,ulong,long), long prec, GEN *p0, GEN *p1)
478 : {
479 1092 : *p0 = gen_0; *p1 = gen_1;
480 1092 : if (!f) return 1;
481 1057 : *p1 = f(E, 1, prec);
482 1057 : if (is_vec_t(typ(*p1)))
483 : {
484 406 : long l = lg(*p1);
485 406 : if (l == 1) { *p0 = *p1 = NULL; return 0; }
486 406 : *p0 = const_vec(l-1, gen_0);
487 : }
488 1057 : return 1;
489 : }
490 : static GEN
491 0 : mktrivial(long both)
492 : {
493 0 : if (!both) return cgetg(1, t_VEC);
494 0 : return mkvec2(cgetg(1,t_VEC), cgetg(1,t_VEC));
495 : }
496 :
497 : static GEN
498 280 : smalldirpowerssum(long N, GEN s, void *E, GEN (*f)(void *, ulong, long),
499 : long both, long prec)
500 : {
501 : GEN F, V, S, SB;
502 280 : if (!N)
503 : {
504 : GEN zerf, onef;
505 42 : if (!mk01(E, f, prec, &zerf, &onef)) return mktrivial(both);
506 42 : return zerf;
507 : }
508 238 : V = vecpowug(N, s, prec);
509 238 : F = vecf(N, E, f, prec);
510 238 : S = naivedirpowerssum(V, F, prec);
511 238 : if (!both) return S;
512 84 : if ((both==2 || !f) && gequalm1(gmul2n(real_i(s), 1)))
513 21 : SB = S;
514 : else
515 : {
516 63 : GEN FB = (both == 1 && F)? conj_i(F): F;
517 : long n;
518 1876 : for (n = 2; n <= N; n++) gel(V, n) = ginv(gmulsg(n, gel(V, n)));
519 63 : SB = naivedirpowerssum(V, FB, prec);
520 : }
521 84 : return mkvec2(S, SB);
522 : }
523 :
524 : static void
525 68509 : v2unpack(GEN v, GEN *pV, GEN *pVB)
526 : {
527 68509 : if (typ(v) == t_COL) { *pV = gel(v,1); *pVB = gel(v,2); }
528 68376 : else { *pV = v; *pVB = NULL; }
529 68509 : }
530 : static GEN
531 69539 : v2pack(GEN V, GEN VB) { return VB? mkcol2(V,VB): V; }
532 :
533 : static GEN
534 1050 : dirpowsuminit(GEN s, GEN onef, GEN zerf, void *E, GEN (*f)(void *, ulong, long), GEN data, long both)
535 : {
536 1050 : long needlog = data[1], prec0 = data[2], prec1 = data[3];
537 1050 : ulong a, b, c, e, q, n, sq = usqrt(data[4]);
538 1050 : GEN V = cgetg(sq+1, t_VEC), W = cgetg(sq+1, t_VEC), VB = NULL, WB = NULL;
539 1050 : GEN Q = cgetg(sq+1, t_VEC), Z = cgetg(sq+1, t_VEC), QB = NULL, ZB = NULL;
540 1050 : GEN logp, c2, Q2, Q3, Q6, c2B = NULL, Q2B = NULL, Q3B = NULL, Q6B = NULL;
541 1050 : long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
542 :
543 1050 : if (both == 1 || (both == 2 && !gequal(real_i(s), gneg(ghalf))))
544 28 : { VB = cgetg(sq+1, t_VEC); WB = cgetg(sq+1, t_VEC); QB = cgetg(sq+1, t_VEC);}
545 1050 : gel(V, 1) = gel(W, 1) = gel(Q, 1) = onef;
546 1050 : if (VB) { gel(VB, 1) = gel(WB, 1) = gel(QB, 1) = onef; }
547 1050 : c2 = gpow(gen_2, s, prec0); if (VB) c2B = ginv(gmul2n(conj_i(c2), 1));
548 1050 : if (f)
549 : {
550 1029 : GEN tmp2 = f(E, 2, prec0);
551 1029 : c2 = gmul(c2, tmp2); if (VB) c2B = gmul(c2B, tmp2);
552 : }
553 1050 : gel(V,2) = c2; /* f(2) 2^s */
554 1050 : gel(W,2) = Qtor(gadd(c2, onef), prec0);
555 1050 : gel(Q,2) = Qtor(gadd(vecsqr(c2), onef), prec0);
556 1050 : if (VB)
557 : {
558 28 : gel(VB, 2) = c2B; gel(WB, 2) = Qtor(gadd(c2B, onef), prec0);
559 28 : gel(QB, 2) = Qtor(gadd(vecsqr(c2B), onef), prec0);
560 : }
561 1050 : logp = NULL;
562 158795 : for (n = 3; n <= sq; n++)
563 : {
564 157745 : GEN u = NULL, uB = NULL, ks = f ? f(E, n, prec0) : gen_1;
565 157745 : if (gequal0(ks)) ks = NULL;
566 157745 : if (odd(n))
567 : {
568 79352 : gp[2] = n;
569 79352 : if (needlog)
570 : {
571 77805 : if (!logp)
572 1029 : logp = logr_abs(utor(n, prec1));
573 : else
574 : { /* log n = log(n-2) + 2 atanh(1 / (n - 1)) */
575 76776 : GEN z = atanhuu(1, n - 1, prec1);
576 76776 : shiftr_inplace(z, 1); logp = addrr(logp, z);
577 : }
578 77805 : if (ks)
579 76230 : u = needlog == 1? powcx(gp, logp, s, prec0) : mpexp(gmul(s, logp));
580 : }
581 1547 : else if (ks) u = gpow(gp, s, prec0);
582 79352 : if (ks)
583 : {
584 77763 : if (VB) uB = gmul(ginv(gmulsg(n, conj_i(u))), ks);
585 77763 : u = gmul(u, ks); /* f(n) n^s */
586 : }
587 : }
588 : else
589 : {
590 78393 : u = vecmul(c2, gel(V, n >> 1));
591 78393 : if (VB) uB = vecmul(c2B, gel(VB, n >> 1));
592 : }
593 157745 : if (ks)
594 : {
595 154693 : gel(V,n) = u; /* f(n) n^s */
596 154693 : gel(W,n) = gadd(gel(W, n-1), gel(V,n)); /*= sum_{i<=n} f(i)i^s */
597 154693 : gel(Q,n) = gadd(gel(Q, n-1), vecsqr(gel(V,n)));/*= sum_{i<=n} f(i^2)i^2s*/
598 154693 : if (VB)
599 : {
600 483 : gel(VB,n) = uB;
601 483 : gel(WB,n) = gadd(gel(WB,n-1), gel(VB,n));
602 483 : gel(QB,n) = gadd(gel(QB,n-1), vecsqr(gel(VB,n)));
603 : }
604 : }
605 : else
606 : {
607 3052 : gel(V,n) = zerf; gel(W,n) = gel(W, n-1); gel(Q,n) = gel(Q, n-1);
608 3052 : if (VB)
609 21 : { gel(VB,n) = zerf; gel(WB,n) = gel(WB, n-1); gel(QB,n) = gel(QB, n-1); }
610 : }
611 : }
612 1050 : Q2 = vecmulsqlv(Q, gel(V,2));
613 1050 : Q3 = vecmulsqlv(Q, gel(V,3));
614 1050 : Q6 = vecmulsqlv(Q, gel(V,6));
615 1050 : if (VB)
616 : {
617 28 : Q2B = vecmulsqlv(QB, gel(VB,2));
618 28 : Q3B = vecmulsqlv(QB, gel(VB,3));
619 28 : Q6B = vecmulsqlv(QB, gel(VB,6));
620 : }
621 : /* a,b,c,e = sqrt(q), sqrt(q/2), sqrt(q/3), sqrt(q/6)
622 : * Z[q] = Q[a] + 2^s Q[b] + 3^s Q[c] + 6^s Q[e], with Q[0] = 0 */
623 1050 : gel(Z, 1) = onef;
624 1050 : gel(Z, 2) = gel(W, 2);
625 1050 : gel(Z, 3) = gel(W, 3);
626 1050 : gel(Z, 4) = gel(Z, 5) = gel(W, 4);
627 1050 : gel(Z, 6) = gel(Z, 7) = gadd(gel(W, 4), gel(V, 6));
628 1050 : if (VB)
629 : {
630 28 : ZB = cgetg(sq+1, t_VEC);
631 28 : gel(ZB, 1) = onef;
632 28 : gel(ZB, 2) = gel(WB, 2);
633 28 : gel(ZB, 3) = gel(WB, 3);
634 28 : gel(ZB, 4) = gel(ZB, 5) = gel(WB, 4);
635 28 : gel(ZB, 6) = gel(ZB, 7) = gadd(gel(WB, 4), gel(VB, 6));
636 : }
637 1050 : a = 2; b = c = e = 1;
638 153545 : for (q = 8; q <= sq; q++)
639 : { /* Gray code: at most one of a,b,c,d differs (by 1) from previous value */
640 152495 : GEN z = gel(Z, q - 1), zB = NULL;
641 : ulong na, nb, nc, ne, na2, nb2, nc2, ne2;
642 152495 : if (VB) zB = gel(ZB, q - 1);
643 152495 : if ((na = usqrt(q)) != a)
644 9191 : { a = na; na2 = na * na; z = gadd(z, gel(V, na2));
645 9191 : if (VB) zB = gadd(zB, gel(VB, na2)); }
646 143304 : else if ((nb = usqrt(q / 2)) != b)
647 6216 : { b = nb; nb2 = 2 * nb * nb; z = gadd(z, gel(V, nb2));
648 6216 : if (VB) zB = gadd(zB, gel(VB, nb2)); }
649 137088 : else if ((nc = usqrt(q / 3)) != c)
650 5418 : { c = nc; nc2 = 3 * nc * nc; z = gadd(z, gel(V, nc2));
651 5418 : if (VB) zB = gadd(zB, gel(VB, nc2)); }
652 131670 : else if ((ne = usqrt(q / 6)) != e)
653 3108 : { e = ne; ne2 = 6 * ne * ne; z = gadd(z, gel(V, ne2));
654 3108 : if (VB) zB = gadd(zB, gel(VB, ne2)); }
655 152495 : gel(Z, q) = z; if (VB) gel(ZB, q) = zB;
656 : }
657 1078 : return v2pack(mkvecn(7, V, W, Q, Q2, Q3, Q6, Z),
658 28 : VB? mkvecn(7, VB, WB, QB, Q2B, Q3B, Q6B, ZB): NULL);
659 : }
660 :
661 : static GEN
662 38301 : sumprimeloop(forprime_t *pT, GEN s, long N, GEN data, GEN S, GEN W, GEN WB,
663 : void *E, GEN (*f)(void *, ulong, long))
664 : {
665 38301 : pari_sp av = avma;
666 38301 : long needlog = data[1], prec0 = data[2], prec1 = data[3];
667 38301 : long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
668 38301 : ulong p, precp = 0;
669 38301 : GEN logp = NULL, SB = WB? S: NULL;
670 :
671 5281977 : while ((p = u_forprime_next(pT)))
672 : {
673 5243396 : GEN u = NULL, ks;
674 5243396 : if (!f) ks = gen_1;
675 : else
676 : {
677 5174509 : ks = f(E, p, prec1);
678 5174271 : if (gequal0(ks)) ks = NULL;
679 5174560 : if (ks && gequal1(ks)) ks = gen_1;
680 : }
681 5243599 : gp[2] = p;
682 5243599 : if (needlog)
683 : {
684 5167383 : if (!logp)
685 30408 : logp = logr_abs(utor(p, prec1));
686 : else
687 : { /* Assuming p and precp are odd,
688 : * log p = log(precp) + 2 atanh((p - precp) / (p + precp)) */
689 5136975 : ulong a = p >> 1, b = precp >> 1; /* p = 2a + 1, precp = 2b + 1 */
690 5136975 : GEN z = atanhuu(a - b, a + b + 1, prec1); /* avoid overflow */
691 5136698 : shiftr_inplace(z, 1); logp = addrr(logp, z);
692 : }
693 10334344 : if (ks) u = needlog == 1? powcx(gp, logp, s, prec0)
694 5167198 : : mpexp(gmul(s, logp));
695 : }
696 76216 : else if (ks) u = gpow(gp, s, prec0);
697 5245453 : if (ks)
698 : {
699 5245478 : S = gadd(S, vecmul(gel(W, N / p), ks == gen_1? u: gmul(ks, u)));
700 5243701 : if (WB)
701 : {
702 2457 : GEN w = gel(WB, N / p);
703 2457 : if (ks != gen_1) w = vecmul(ks, w);
704 2457 : SB = gadd(SB, gdiv(w, gmulsg(p, conj_i(u))));
705 : }
706 : }
707 5243676 : precp = p;
708 5243676 : if ((p & 0x1ff) == 1)
709 : {
710 19565 : if (!logp) (void)gc_all(av, SB? 2: 1, &S, &SB);
711 19285 : else (void)gc_all(av, SB? 3: 2, &S, &logp, &SB);
712 : }
713 : }
714 38182 : return v2pack(S, SB);
715 : }
716 :
717 : static GEN
718 48475 : add4(GEN a, GEN b, GEN c, GEN d) { return gadd(gadd(a,b), gadd(c,d)); }
719 :
720 : static const long step = 2048;
721 : static int
722 29556 : mksqfloop(long N, long x1, GEN R, GEN RB, GEN *pS, GEN *pSB)
723 : {
724 29556 : GEN V = gel(R,1), Q = gel(R,3), Q2 = gel(R,4);
725 29556 : GEN Q3 = gel(R,5), Q6 = gel(R,6), Z = gel(R,7);
726 29556 : GEN v, VB = NULL, QB = NULL, Q2B = NULL, Q3B = NULL, Q6B = NULL, ZB = NULL;
727 29556 : long x2, j, lv, sq = lg(V)-1;
728 :
729 29556 : if (RB)
730 28 : { VB = gel(RB,1); QB = gel(RB,3); Q2B = gel(RB,4);
731 28 : Q3B = gel(RB,5), Q6B = gel(RB,6); ZB = gel(RB,7); }
732 : /* beware overflow, fuse last two bins (avoid a tiny remainder) */
733 29556 : x2 = (N >= 2*step && N - 2*step >= x1)? x1-1 + step: N;
734 29556 : v = vecfactorsquarefreeu_coprime(x1, x2, mkvecsmall2(2, 3));
735 30557 : lv = lg(v);
736 59975025 : for (j = 1; j < lv; j++)
737 59945465 : if (gel(v,j))
738 : {
739 18217302 : ulong d = x1 - 1 + j; /* squarefree, coprime to 6 */
740 18217302 : GEN t = smallfact(d, gel(v,j), sq, V), u;
741 18195518 : GEN tB = NULL, uB = NULL; /* = f(d) d^s */
742 : long a, b, c, e, q;
743 18195518 : if (!t || gequal0(t)) continue;
744 3800273 : if (VB) tB = vecinv(gmulsg(d, conj_i(t)));
745 : /* warning: gives 1/conj(f(d)) d^(-1-conj(s)), equal to
746 : f(d) d^(-1-conj(s)) only if |f(d)|=1. */
747 : /* S += f(d) * d^s * Z[q] */
748 3867379 : q = N / d;
749 3867379 : if (q == 1)
750 : {
751 1474430 : *pS = gadd(*pS, t); if (VB) *pSB = gadd(*pSB, tB);
752 1475662 : continue;
753 : }
754 2392949 : if (q <= sq) { u = gel(Z, q); if (VB) uB = gel(ZB, q); }
755 : else
756 : {
757 48235 : a = usqrt(q); b = usqrt(q / 2); c = usqrt(q / 3); e = usqrt(q / 6);
758 48293 : u = add4(gel(Q,a), gel(Q2,b), gel(Q3,c), gel(Q6,e));
759 48293 : if (VB) uB = add4(gel(QB,a), gel(Q2B,b), gel(Q3B,c), gel(Q6B,e));
760 : }
761 2393007 : *pS = gadd(*pS, vecmul(t, u)); if (VB) *pSB = gadd(*pSB, vecmul(tB, uB));
762 : }
763 29560 : return x2 == N;
764 : }
765 :
766 : static GEN
767 1050 : mkdata(long N, GEN s, long prec)
768 : {
769 1050 : long needlog, prec0, prec1, m = mt_nbthreads(), STEP = maxss(N / (m * m), 1);
770 1050 : prec1 = prec0 = prec + EXTRAPRECWORD;
771 1050 : needlog = get_needlog(s);
772 1050 : if (needlog == 1) prec1 = powcx_prec(log2((double)N), s, prec);
773 1050 : return mkvecsmalln(5, needlog, prec0, prec1, N, STEP);
774 : }
775 :
776 : static GEN
777 13748 : gp_callUp(void *E, ulong x, long prec)
778 : {
779 13748 : long n[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
780 13748 : n[2] = x; return gp_callprec(E, n, prec);
781 : }
782 :
783 : /* both is
784 : * 0: sum_{n<=N}f(n)n^s
785 : * 1: sum for (f,s) and (conj(f),-1-s)
786 : * 2: sum for (f,s) and (f,-1-s), assuming |f(n)| in {0,1} */
787 : static GEN
788 203 : dirpowerssumfun_i(ulong N, GEN s, void *E, GEN (*f)(void *, ulong, long),
789 : long both, long prec)
790 : {
791 : forprime_t T;
792 : pari_sp av;
793 : GEN onef, zerf, R, RB, W, WB, S, SB, data;
794 : ulong x1;
795 :
796 203 : if ((f && N < 49) || (!f && N < 1000))
797 140 : return smalldirpowerssum(N, s, E, f, both, prec);
798 63 : if (!mk01(E, f, prec, &zerf, &onef)) return mktrivial(both);
799 63 : data = mkdata(N, s, prec); s = gprec_w(s, prec + EXTRAPRECWORD);
800 63 : v2unpack(dirpowsuminit(s, onef, zerf, E, f, data, both), &R, &RB);
801 63 : W = gel(R,2); WB = RB? gel(RB,2): NULL;
802 63 : av = avma; u_forprime_init(&T, lg(W), N);
803 63 : v2unpack(sumprimeloop(&T, s, N, data, zerf, W, WB, E, f), &S, &SB);
804 413 : for(x1 = 1;; x1 += step)
805 : {
806 413 : if (mksqfloop(N, x1, R, RB, &S, &SB))
807 63 : return both? mkvec2(S, conj_i(SB? SB: S)): S;
808 350 : (void)gc_all(av, SB? 2: 1, &S, &SB);
809 : }
810 : }
811 : GEN
812 203 : dirpowerssumfun(ulong N, GEN s, void *E, GEN (*f)(void *, ulong, long),
813 : long both, long prec)
814 : {
815 203 : pari_sp av = avma;
816 203 : return gc_GEN(av, dirpowerssumfun_i(N, s, E, f, both, prec));
817 : }
818 :
819 : GEN
820 133 : dirpowerssum(ulong N, GEN s, long both, long prec)
821 133 : { return dirpowerssumfun(N, s, NULL, NULL, both, prec); }
822 :
823 : GEN
824 154 : dirpowerssum0(GEN N, GEN s, GEN f, long both, long prec)
825 : {
826 154 : if (typ(N) != t_INT) pari_err_TYPE("dirpowerssum", N);
827 147 : if (signe(N) <= 0) N = gen_0;
828 147 : if (!f) return dirpowerssum(itou(N), s, both, prec);
829 70 : if (typ(f) != t_CLOSURE) pari_err_TYPE("dirpowerssum", f);
830 70 : return dirpowerssumfun(itou(N), s, (void*)f, gp_callUp, both, prec);
831 : }
832 :
833 : GEN
834 29143 : parsqf_worker(GEN gk, GEN vR, GEN gN)
835 : {
836 29143 : pari_sp av = avma;
837 : GEN R, RB, onef, S, SB;
838 29143 : long k = itou(gk), N = itou(gN), x1 = 1 + step * k;
839 29143 : v2unpack(vR, &R, &RB); onef = gmael(R,1,1);
840 29143 : S = SB = is_vec_t(typ(onef)) ? zerovec(lg(onef) - 1): gen_0;
841 29143 : (void)mksqfloop(N, x1, R, RB, &S, &SB);
842 29147 : return gc_GEN(av, v2pack(S, RB? SB: NULL));
843 : }
844 :
845 : static GEN
846 5351172 : mycallvec(void *f, ulong n, long prec)
847 : {
848 5351172 : GEN F = (GEN)f;
849 5351172 : if (!f) return gen_1;
850 390663 : if (typ(F) == t_CLOSURE) return gp_callUp(f, n, prec);
851 390663 : return gel(F, (n-1) % (lg(F)-1) + 1);
852 : }
853 :
854 : GEN
855 38254 : parsumprimefun_worker(GEN gk, GEN s, GEN zerf, GEN data, GEN vW, GEN f)
856 : {
857 38254 : pari_sp av = avma;
858 : forprime_t T;
859 : GEN W, WB, S;
860 38254 : long k = itou(gk), sq, N = data[4], STEP = data[5];
861 :
862 38253 : v2unpack(vW, &W, &WB); sq = lg(W)-1;
863 38251 : if (isintzero(f)) f = NULL;
864 38248 : u_forprime_init(&T, k * STEP + sq + 1, minss(N, (k + 1) * STEP + sq));
865 38227 : S = sumprimeloop(&T, s, N, data, zerf, W, WB, (void*)f, mycallvec);
866 38287 : return gc_GEN(av, S);
867 : }
868 :
869 : static GEN
870 987 : vR_get_vW(GEN vR)
871 : {
872 : GEN R, RB, W, WB;
873 987 : v2unpack(vR, &R, &RB); W = gel(R,2); WB = RB? gel(RB,2): NULL;
874 987 : return v2pack(W, WB);
875 : }
876 :
877 : static GEN
878 987 : halfconj(long both, GEN V)
879 987 : { return both ? mkvec2(gel(V, 1), gconj(gel(V, 2))) : V; }
880 :
881 : static GEN
882 1127 : pardirpowerssumfun_i(GEN f, ulong N, GEN s, long both, long prec)
883 : {
884 : GEN worker, worker2, data, vR, onef, zerf;
885 :
886 1127 : if ((f && N < 49) || (!f && N < 10000UL))
887 140 : return smalldirpowerssum(N, s, (void*)f, mycallvec, both, prec);
888 987 : if (!mk01((void*)f, mycallvec, prec, &zerf, &onef)) return mktrivial(both);
889 987 : data = mkdata(N, s, prec); s = gprec_w(s, prec + EXTRAPRECWORD);
890 987 : vR = dirpowsuminit(s, onef, zerf, (void*)f, mycallvec, data, both);
891 987 : worker = snm_closure(is_entry("_parsumprimefun_worker"),
892 : mkvecn(5, s, zerf, data, vR_get_vW(vR), f? f: gen_0));
893 987 : worker2 = snm_closure(is_entry("_parsqf_worker"), mkvec2(vR, utoi(N)));
894 1974 : return halfconj(both,
895 987 : gadd(parsum(gen_0, utoipos((N-1) / data[5]), worker),
896 987 : parsum(gen_0, utoipos(maxss((N-1) / step - 1, 0)), worker2)));
897 : }
898 :
899 : GEN
900 1127 : pardirpowerssumfun(GEN f, ulong N, GEN s, long both, long prec)
901 : {
902 1127 : pari_sp av = avma;
903 1127 : return gc_GEN(av, pardirpowerssumfun_i(f, N, s, both, prec));
904 : }
905 : GEN
906 0 : pardirpowerssum0(GEN N, GEN s, GEN f, long both, long prec)
907 : {
908 0 : if (typ(N) != t_INT) pari_err_TYPE("pardirpowerssum", N);
909 0 : return pardirpowerssumfun(f, itou(N), s, both, prec);
910 : }
911 : GEN
912 0 : pardirpowerssum(ulong N, GEN s, long prec)
913 0 : { return pardirpowerssumfun(NULL, N, s, 0, prec); }
|
