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 24 : err_direuler(GEN x)
25 24 : { 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 24618 : dirmuleuler_small(GEN V, GEN v, long n, ulong p, GEN s, long d)
33 : {
34 24618 : long i, j, m = n, D = minss(d+2, lg(s));
35 24618 : ulong q = 1, X = lg(V)-1;
36 :
37 81227 : for (i = 3, q = p; i < D; i++, q *= p) /* q*p does not overflow */
38 : {
39 56609 : GEN aq = gel(s,i);
40 56609 : if (gequal0(aq)) continue;
41 : /* j = 1 */
42 46079 : gel(V,q) = aq;
43 46079 : v[++n] = q;
44 2801152 : for (j = 2; j <= m; j++)
45 : {
46 2755073 : ulong nj = umuluu_le(uel(v,j), q, X);
47 2755073 : if (!nj) continue;
48 164579 : gel(V,nj) = gmul(aq, gel(V,v[j]));
49 164579 : v[++n] = nj;
50 : }
51 : }
52 24618 : return n;
53 : }
54 :
55 : /* ap != 0 for efficiency, p > sqrt(X) */
56 : static void
57 264729 : dirmuleuler_large(GEN V, ulong p, GEN ap)
58 : {
59 264729 : long j, jp, X = lg(V)-1;
60 264729 : gel(V,p) = ap;
61 1291437 : for (j = 2, jp = 2*p; jp <= X; j++, jp += p) gel(V,jp) = gmul(ap, gel(V,j));
62 264729 : }
63 :
64 : static ulong
65 8810 : direulertou(GEN a, GEN fl(GEN))
66 : {
67 8810 : if (typ(a) != t_INT)
68 : {
69 42 : a = fl(a);
70 24 : if (typ(a) != t_INT) pari_err_TYPE("direuler", a);
71 : }
72 8792 : return signe(a)<=0 ? 0: itou(a);
73 : }
74 :
75 : static GEN
76 3198 : direuler_Sbad(GEN V, GEN v, GEN Sbad, ulong *n)
77 : {
78 3198 : long i, l = lg(Sbad);
79 3198 : ulong X = lg(V)-1;
80 3198 : GEN pbad = gen_1;
81 8292 : for (i = 1; i < l; i++)
82 : {
83 5124 : GEN ai = gel(Sbad,i), gq;
84 : ulong q;
85 5124 : if (typ(ai) != t_VEC || lg(ai) != 3)
86 12 : pari_err_TYPE("direuler [bad primes]",ai);
87 5112 : gq = gel(ai,1);
88 5112 : if (typ(gq) != t_INT || signe(gq) <= 0)
89 6 : pari_err_TYPE("direuler [prime expected]",gq);
90 5106 : if (lgefint(gq) == 3 && (q = uel(gq,2)) <= X)
91 : {
92 4134 : long d = ulogint(X, q) + 1;
93 4134 : GEN s = direuler_factor(gel(ai,2), d);
94 4122 : *n = dirmuleuler_small(V, v, *n, q, s, d);
95 4122 : pbad = muliu(pbad, q);
96 : }
97 : }
98 3168 : return pbad;
99 : }
100 :
101 : GEN
102 574 : direuler_bad(void *E, GEN (*eval)(void *,GEN,long), GEN a,GEN b,GEN c, GEN Sbad)
103 : {
104 : ulong au, bu, X, sqrtX, n, p;
105 574 : pari_sp av0 = avma;
106 : GEN gp, v, V;
107 : forprime_t T;
108 574 : au = direulertou(a, gceil);
109 568 : bu = direulertou(b, gfloor);
110 562 : X = c ? direulertou(c, gfloor): bu;
111 556 : if (X == 0) return cgetg(1,t_VEC);
112 550 : if (bu > X) bu = X;
113 550 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
114 538 : v = vecsmall_ei(X, 1);
115 538 : V = vec_ei(X, 1);
116 538 : n = 1;
117 538 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
118 508 : p = 1; gp = cgetipos(3); sqrtX = usqrt(X);
119 7120 : while (p <= sqrtX && (p = u_forprime_next(&T)))
120 6630 : if (!Sbad || umodiu(Sbad, p))
121 : {
122 6540 : long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
123 : GEN s;
124 6540 : gp[2] = p; s = eval(E, gp, d);
125 6522 : n = dirmuleuler_small(V, v, n, p, s, d);
126 : }
127 634447 : while ((p = u_forprime_next(&T))) /* sqrt(X) < p <= X */
128 633957 : if (!Sbad || umodiu(Sbad, p))
129 : {
130 : GEN s;
131 633951 : gp[2] = p; s = eval(E, gp, 2); /* s either t_POL or t_SER of val 0 */
132 633951 : if (lg(s) > 3 && !gequal0(gel(s,3)))
133 119265 : dirmuleuler_large(V, p, gel(s,3));
134 : }
135 490 : return gc_GEN(av0,V);
136 : }
137 :
138 : /* return a t_SER or a truncated t_POL to precision n */
139 : GEN
140 644625 : direuler_factor(GEN s, long n)
141 : {
142 644625 : long t = typ(s);
143 644625 : if (is_scalar_t(t))
144 : {
145 28452 : if (!gequal1(s)) err_direuler(s);
146 28440 : return scalarpol_shallow(s,0);
147 : }
148 616173 : switch(t)
149 : {
150 4896 : case t_POL: break; /* no need to RgXn_red */
151 611007 : case t_RFRAC:
152 : {
153 611007 : GEN p = gel(s,1), q = gel(s,2);
154 611007 : q = RgXn_red_shallow(q,n);
155 611007 : s = RgXn_inv(q, n);
156 611007 : if (typ(p) == t_POL && varn(p) == varn(q))
157 : {
158 24 : p = RgXn_red_shallow(p, n);
159 24 : s = RgXn_mul(s, p, n);
160 : }
161 : else
162 610983 : if (!gequal1(p)) s = RgX_Rg_mul(s, p);
163 611007 : if (!signe(s) || !gequal1(gel(s,2))) err_direuler(s);
164 610995 : break;
165 : }
166 264 : case t_SER:
167 264 : if (!signe(s) || valser(s) || !gequal1(gel(s,2))) err_direuler(s);
168 264 : break;
169 6 : default: pari_err_TYPE("direuler", s);
170 : }
171 616155 : return s;
172 : }
173 :
174 : struct eval_bad
175 : {
176 : void *E;
177 : GEN (*eval)(void *, GEN);
178 : };
179 : static GEN
180 589959 : eval_bad(void *E, GEN p, long n)
181 : {
182 589959 : struct eval_bad *d = (struct eval_bad*) E;
183 589959 : return direuler_factor(d->eval(d->E, p), n);
184 : }
185 : GEN
186 256 : direuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, GEN c)
187 : {
188 : struct eval_bad d;
189 256 : d.E= E; d.eval = eval;
190 256 : return direuler_bad((void*)&d, eval_bad, a, b, c, NULL);
191 : }
192 :
193 : static GEN
194 26748 : primelist(forprime_t *T, GEN Sbad, long n, long *running)
195 : {
196 26748 : GEN P = cgetg(n+1, t_VECSMALL);
197 : long i, j;
198 259590 : for (i = 1, j = 1; i <= n; i++)
199 : {
200 236340 : ulong p = u_forprime_next(T);
201 236340 : if (!p) { *running = 0; break; }
202 232842 : if (Sbad && umodiu(Sbad, p)==0) continue;
203 228816 : uel(P,j++) = p;
204 : }
205 26748 : setlg(P, j);
206 26748 : return P;
207 : }
208 :
209 : GEN
210 3504 : pardireuler(GEN worker, GEN a, GEN b, GEN c, GEN Sbad)
211 : {
212 : ulong au, bu, X, sqrtX, n, snX, nX;
213 3504 : pari_sp av0 = avma;
214 : GEN v, V;
215 : forprime_t T;
216 : struct pari_mt pt;
217 3504 : long running = 1, pending = 0;
218 3504 : au = direulertou(a, gceil);
219 3504 : bu = direulertou(b, gfloor);
220 3504 : X = c ? direulertou(c, gfloor): bu;
221 3504 : if (X == 0) return cgetg(1,t_VEC);
222 3504 : if (bu > X) bu = X;
223 3504 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
224 3498 : v = vecsmall_ei(X, 1);
225 3498 : V = vec_ei(X, 1);
226 3498 : n = 1;
227 3498 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
228 3498 : sqrtX = usqrt(X); snX = uprimepi(sqrtX); nX = uprimepi(X);
229 3498 : if (snX)
230 : {
231 3486 : GEN P = primelist(&T, Sbad, snX, &running);
232 3486 : GEN R = gel(closure_callgenvec(worker, mkvec2(P, utoi(X))), 2);
233 3486 : long i, l = lg(P);
234 17460 : for (i = 1; i < l; i++)
235 : {
236 13974 : GEN s = gel(R,i);
237 13974 : n = dirmuleuler_small(V, v, n, uel(P,i), s, lg(s));
238 : }
239 12 : } else snX = 1;
240 3498 : mt_queue_start_lim(&pt, worker, (nX+snX-1)/snX);
241 26760 : while (running || pending)
242 : {
243 : GEN done;
244 23262 : GEN P = running? primelist(&T, Sbad, snX, &running): NULL;
245 23262 : mt_queue_submit(&pt, 0, P ? mkvec2(P, utoi(X)): NULL);
246 23262 : done = mt_queue_get(&pt, NULL, &pending);
247 23262 : if (done)
248 : {
249 23262 : GEN P = gel(done,1), R = gel(done,2);
250 23262 : long j, l = lg(P);
251 238104 : for (j=1; j<l; j++)
252 : {
253 214842 : GEN F = gel(R,j);
254 214842 : if (degpol(F) && !gequal0(gel(F,3)))
255 145464 : dirmuleuler_large(V, uel(P,j), gel(F,3));
256 : }
257 : }
258 : }
259 3498 : mt_queue_end(&pt);
260 3498 : return gc_GEN(av0,V);
261 : }
262 :
263 : /********************************************************************/
264 : /** **/
265 : /** DIRPOWERS and DIRPOWERSSUM **/
266 : /** **/
267 : /********************************************************************/
268 :
269 : /* [1^B,...,N^B] */
270 : GEN
271 579 : vecpowuu(long N, ulong B)
272 : {
273 : GEN v;
274 : long p, i;
275 : forprime_t T;
276 :
277 579 : if (B <= 8000)
278 : {
279 579 : if (!B) return const_vec(N,gen_1);
280 573 : v = cgetg(N+1, t_VEC); if (N == 0) return v;
281 573 : gel(v,1) = gen_1;
282 573 : if (B == 1)
283 79488 : for (i = 2; i <= N; i++) gel(v,i) = utoipos(i);
284 393 : else if (B == 2)
285 : {
286 : ulong o, s;
287 232 : if (N & HIGHMASK)
288 0 : for (i = 2, o = 3; i <= N; i++, o += 2)
289 0 : gel(v,i) = addiu(gel(v,i-1), o);
290 : else
291 26546 : for (i = 2, s = 1, o = 3; i <= N; i++, s += o, o += 2)
292 26314 : gel(v,i) = utoipos(s + o);
293 : }
294 161 : else if (B == 3)
295 720 : for (i = 2; i <= N; i++) gel(v,i) = powuu(i, B);
296 : else
297 : {
298 149 : long k, Bk, e = expu(N);
299 6135 : for (i = 3; i <= N; i += 2) gel(v,i) = powuu(i, B);
300 1012 : for (k = 1; k <= e; k++)
301 : {
302 863 : N >>= 1; Bk = B * k;
303 6909 : for (i = 1; i <= N; i += 2) gel(v, i << k) = shifti(gel(v, i), Bk);
304 : }
305 : }
306 573 : return v;
307 : }
308 0 : v = const_vec(N, NULL);
309 0 : u_forprime_init(&T, 3, N);
310 0 : while ((p = u_forprime_next(&T)))
311 : {
312 : long m, pk, oldpk;
313 0 : gel(v,p) = powuu(p, B);
314 0 : for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
315 : {
316 0 : if (pk != p) gel(v,pk) = mulii(gel(v,oldpk), gel(v,p));
317 0 : for (m = N/pk; m > 1; m--)
318 0 : if (gel(v,m) && m%p) gel(v, m*pk) = mulii(gel(v,m), gel(v,pk));
319 : }
320 : }
321 0 : gel(v,1) = gen_1;
322 0 : for (i = 2; i <= N; i+=2)
323 : {
324 0 : long vi = vals(i);
325 0 : gel(v,i) = shifti(gel(v,i >> vi), B * vi);
326 : }
327 0 : return v;
328 : }
329 :
330 : /* does x^s require log(x) ? */
331 : static long
332 17435 : get_needlog(GEN s)
333 : {
334 17435 : switch(typ(s))
335 : {
336 7259 : case t_REAL: return 2; /* yes but not powcx */
337 7471 : case t_COMPLEX: return 1; /* yes using powcx */
338 2705 : default: return 0; /* no */
339 : }
340 : }
341 : /* [1^B,...,N^B] */
342 : GEN
343 16696 : vecpowug(long N, GEN B, long prec)
344 : {
345 16696 : GEN v, logp = NULL;
346 16696 : long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
347 16696 : long p, precp = 2, prec0, prec1, needlog;
348 : forprime_t T;
349 16696 : if (N == 1) return mkvec(gen_1);
350 16679 : if (typ(B) == t_INT && lgefint(B) <= 3 && signe(B) >= 0)
351 144 : return vecpowuu(N, itou(B));
352 16535 : needlog = get_needlog(B);
353 16535 : prec1 = prec0 = prec;
354 16535 : if (needlog == 1) prec1 = powcx_prec(log2((double)N), B, prec);
355 16535 : u_forprime_init(&T, 2, N);
356 16535 : v = const_vec(N, NULL);
357 16535 : gel(v,1) = gen_1;
358 1422172 : while ((p = u_forprime_next(&T)))
359 : {
360 : long m, pk, oldpk;
361 : GEN u;
362 1405637 : gp[2] = p;
363 1405637 : if (needlog)
364 : {
365 154082 : if (!logp)
366 27696 : logp = logr_abs(utor(p, prec1));
367 : else
368 : { /* Assuming p and precp are odd,
369 : * log p = log(precp) + 2 atanh((p - precp) / (p + precp)) */
370 126386 : ulong a = p >> 1, b = precp >> 1; /* p = 2a + 1, precp = 2b + 1 */
371 126386 : GEN z = atanhuu(a - b, a + b + 1, prec1); /* avoid overflow */
372 126386 : shiftr_inplace(z, 1); logp = addrr(logp, z);
373 : }
374 73280 : u = needlog == 1? powcx(gp, logp, B, prec0)
375 154082 : : mpexp(gmul(B, logp));
376 154082 : if (p == 2) logp = NULL; /* reset: precp must be odd */
377 : }
378 : else
379 1251555 : u = gpow(gp, B, prec0);
380 1405637 : precp = p;
381 1405637 : gel(v,p) = u; /* p^B */
382 1405637 : if (prec0 != prec) gel(v,p) = gprec_wtrunc(gel(v,p), prec);
383 2941495 : for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
384 : {
385 1535858 : if (pk != p) gel(v,pk) = gmul(gel(v,oldpk), gel(v,p));
386 40279446 : for (m = N/pk; m > 1; m--)
387 38743588 : if (gel(v,m) && m%p) gel(v, m*pk) = gmul(gel(v,m), gel(v,pk));
388 : }
389 : }
390 16535 : return v;
391 : }
392 :
393 : GEN
394 570 : dirpowers(long n, GEN x, long prec)
395 : {
396 : pari_sp av;
397 : GEN v;
398 570 : if (n <= 0) return cgetg(1, t_VEC);
399 558 : av = avma; v = vecpowug(n, x, prec);
400 558 : if (typ(x) == t_INT && lgefint(x) <= 3 && signe(x) >= 0 && cmpiu(x, 2) <= 0)
401 114 : return v;
402 444 : return gc_GEN(av, v);
403 : }
404 :
405 : /* f is a totally multiplicative function of modulus 0 or 1
406 : * (essentially a Dirichlet character). Compute simultaneously
407 : * sum_{0 < n <= N} f(n)n^s and sum_{0 < n <= N} f(n)n^{-1-conj(s)}
408 : * Warning: s is conjugated, but not f. Main application for Riemann-Siegel,
409 : * where we need R(chi,s) and conj(R(chi,1-conj(s))). */
410 :
411 : static GEN
412 2772 : vecmulsqlv(GEN Q, GEN V)
413 : {
414 : long l, i;
415 : GEN W;
416 2772 : if (typ(V) != t_VEC) return RgV_Rg_mul(Q, V);
417 1062 : l = lg(Q); W = cgetg(l, t_VEC);
418 53352 : for (i = 1; i < l; i++) gel(W, i) = vecmul(gel(Q, i), V);
419 1062 : return W;
420 : }
421 :
422 : /* P = prime divisors of (squarefree) n, V[i] = i^s for i <= sq.
423 : * Return NULL if n is not sq-smooth, else f(n)n^s */
424 : static GEN
425 15804066 : smallfact(ulong n, GEN P, ulong sq, GEN V)
426 : {
427 : long i, l;
428 : ulong p, m, o;
429 : GEN c;
430 15804066 : if (n <= sq) return gel(V,n);
431 15762084 : l = lg(P); m = p = uel(P, l-1); if (p > sq) return NULL;
432 3370530 : for (i = l-2; i > 1; i--, m = o) { p = uel(P,i); o = m*p; if (o > sq) break; }
433 3341970 : c = gel(V,m); n /= m; /* m <= sq, o = m * p > sq */
434 3341970 : if (n > sq) { c = vecmul(c, gel(V,p)); n /= p; }
435 3341970 : return vecmul(c, gel(V,n));
436 : }
437 :
438 : static GEN
439 834 : _Qtor(GEN x, long prec)
440 834 : { return typ(x) == t_FRAC? fractor(x, prec): x; }
441 : static GEN
442 2099 : Qtor(GEN x, long prec)
443 : {
444 2099 : long tx = typ(x);
445 2933 : if (tx == t_VEC || tx == t_COL) pari_APPLY_same(_Qtor(gel(x, i), prec));
446 1355 : return tx == t_FRAC? fractor(x, prec): x;
447 : }
448 :
449 : static GEN
450 197 : vecf(long N, void *E, GEN (*f)(void *, ulong, long), long prec)
451 : {
452 : GEN v;
453 : long n;
454 197 : if (!f) return NULL;
455 120 : v = cgetg(N + 1, t_VEC);
456 28698 : for (n = 1; n <= N; n++) gel(v,n) = f(E, n, prec);
457 120 : return v;
458 : }
459 :
460 : /* Here #V = #F > 0 is small. Analogous to dotproduct, with following
461 : * semantic differences: uses V[1] = 1; V has scalar values but F may have
462 : * vector values */
463 : static GEN
464 251 : naivedirpowerssum(GEN V, GEN F, long prec)
465 : {
466 : GEN S;
467 251 : if (!F) S = RgV_sum(V);
468 : else
469 : {
470 162 : long n, l = lg(V);
471 162 : S = gel(F,1); /* V[1] = 1 */
472 29526 : for (n = 2; n < l; n++) S = gadd(S, gmul(gel(V, n), gel(F, n)));
473 : }
474 251 : return Qtor(S, prec);
475 : }
476 :
477 : /* set *p0 and *p1 to 0 and 1 in the algebra where f takes its values */
478 : static int
479 936 : mk01(void *E, GEN (*f)(void*,ulong,long), long prec, GEN *p0, GEN *p1)
480 : {
481 936 : *p0 = gen_0; *p1 = gen_1;
482 936 : if (!f) return 1;
483 906 : *p1 = f(E, 1, prec);
484 906 : if (is_vec_t(typ(*p1)))
485 : {
486 348 : long l = lg(*p1);
487 348 : if (l == 1) { *p0 = *p1 = NULL; return 0; }
488 348 : *p0 = const_vec(l-1, gen_0);
489 : }
490 906 : return 1;
491 : }
492 : static GEN
493 0 : mktrivial(long both)
494 : {
495 0 : if (!both) return cgetg(1, t_VEC);
496 0 : return mkvec2(cgetg(1,t_VEC), cgetg(1,t_VEC));
497 : }
498 :
499 : static GEN
500 233 : smalldirpowerssum(long N, GEN s, void *E, GEN (*f)(void *, ulong, long),
501 : long both, long prec)
502 : {
503 : GEN F, V, S, SB;
504 233 : if (!N)
505 : {
506 : GEN zerf, onef;
507 36 : if (!mk01(E, f, prec, &zerf, &onef)) return mktrivial(both);
508 36 : return zerf;
509 : }
510 197 : V = vecpowug(N, s, prec);
511 197 : F = vecf(N, E, f, prec);
512 197 : S = naivedirpowerssum(V, F, prec);
513 197 : if (!both) return S;
514 72 : if ((both==2 || !f) && gequalm1(gmul2n(real_i(s), 1)))
515 18 : SB = S;
516 : else
517 : {
518 54 : GEN FB = (both == 1 && F)? conj_i(F): F;
519 : long n;
520 1608 : for (n = 2; n <= N; n++) gel(V, n) = ginv(gmulsg(n, gel(V, n)));
521 54 : SB = naivedirpowerssum(V, FB, prec);
522 : }
523 72 : return mkvec2(S, SB);
524 : }
525 :
526 : static void
527 26784 : v2unpack(GEN v, GEN *pV, GEN *pVB)
528 : {
529 26784 : if (typ(v) == t_COL) { *pV = gel(v,1); *pVB = gel(v,2); }
530 26730 : else { *pV = v; *pVB = NULL; }
531 26784 : }
532 : static GEN
533 27630 : v2pack(GEN V, GEN VB) { return VB? mkcol2(V,VB): V; }
534 :
535 : static GEN
536 900 : dirpowsuminit(GEN s, GEN onef, GEN zerf, void *E, GEN (*f)(void *, ulong, long), GEN data, long both)
537 : {
538 900 : long needlog = data[1], prec0 = data[2], prec1 = data[3];
539 900 : ulong a, b, c, e, q, n, sq = usqrt(data[4]);
540 900 : GEN V = cgetg(sq+1, t_VEC), W = cgetg(sq+1, t_VEC), VB = NULL, WB = NULL;
541 900 : GEN Q = cgetg(sq+1, t_VEC), Z = cgetg(sq+1, t_VEC), QB = NULL, ZB = NULL;
542 900 : GEN logp, c2, Q2, Q3, Q6, c2B = NULL, Q2B = NULL, Q3B = NULL, Q6B = NULL;
543 900 : long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
544 :
545 900 : if (both == 1 || (both == 2 && !gequal(real_i(s), gneg(ghalf))))
546 24 : { VB = cgetg(sq+1, t_VEC); WB = cgetg(sq+1, t_VEC); QB = cgetg(sq+1, t_VEC);}
547 900 : gel(V, 1) = gel(W, 1) = gel(Q, 1) = onef;
548 900 : if (VB) { gel(VB, 1) = gel(WB, 1) = gel(QB, 1) = onef; }
549 900 : c2 = gpow(gen_2, s, prec0); if (VB) c2B = ginv(gmul2n(conj_i(c2), 1));
550 900 : if (f)
551 : {
552 882 : GEN tmp2 = f(E, 2, prec0);
553 882 : c2 = gmul(c2, tmp2); if (VB) c2B = gmul(c2B, tmp2);
554 : }
555 900 : gel(V,2) = c2; /* f(2) 2^s */
556 900 : gel(W,2) = Qtor(gadd(c2, onef), prec0);
557 900 : gel(Q,2) = Qtor(gadd(vecsqr(c2), onef), prec0);
558 900 : if (VB)
559 : {
560 24 : gel(VB, 2) = c2B; gel(WB, 2) = Qtor(gadd(c2B, onef), prec0);
561 24 : gel(QB, 2) = Qtor(gadd(vecsqr(c2B), onef), prec0);
562 : }
563 900 : logp = NULL;
564 136110 : for (n = 3; n <= sq; n++)
565 : {
566 135210 : GEN u = NULL, uB = NULL, ks = f ? f(E, n, prec0) : gen_1;
567 135210 : if (gequal0(ks)) ks = NULL;
568 135210 : if (odd(n))
569 : {
570 68016 : gp[2] = n;
571 68016 : if (needlog)
572 : {
573 66690 : if (!logp)
574 882 : logp = logr_abs(utor(n, prec1));
575 : else
576 : { /* log n = log(n-2) + 2 atanh(1 / (n - 1)) */
577 65808 : GEN z = atanhuu(1, n - 1, prec1);
578 65808 : shiftr_inplace(z, 1); logp = addrr(logp, z);
579 : }
580 66690 : if (ks)
581 65340 : u = needlog == 1? powcx(gp, logp, s, prec0) : mpexp(gmul(s, logp));
582 : }
583 1326 : else if (ks) u = gpow(gp, s, prec0);
584 68016 : if (ks)
585 : {
586 66654 : if (VB) uB = gmul(ginv(gmulsg(n, conj_i(u))), ks);
587 66654 : u = gmul(u, ks); /* f(n) n^s */
588 : }
589 : }
590 : else
591 : {
592 67194 : u = vecmul(c2, gel(V, n >> 1));
593 67194 : if (VB) uB = vecmul(c2B, gel(VB, n >> 1));
594 : }
595 135210 : if (ks)
596 : {
597 132594 : gel(V,n) = u; /* f(n) n^s */
598 132594 : gel(W,n) = gadd(gel(W, n-1), gel(V,n)); /*= sum_{i<=n} f(i)i^s */
599 132594 : gel(Q,n) = gadd(gel(Q, n-1), vecsqr(gel(V,n)));/*= sum_{i<=n} f(i^2)i^2s*/
600 132594 : if (VB)
601 : {
602 414 : gel(VB,n) = uB;
603 414 : gel(WB,n) = gadd(gel(WB,n-1), gel(VB,n));
604 414 : gel(QB,n) = gadd(gel(QB,n-1), vecsqr(gel(VB,n)));
605 : }
606 : }
607 : else
608 : {
609 2616 : gel(V,n) = zerf; gel(W,n) = gel(W, n-1); gel(Q,n) = gel(Q, n-1);
610 2616 : if (VB)
611 18 : { gel(VB,n) = zerf; gel(WB,n) = gel(WB, n-1); gel(QB,n) = gel(QB, n-1); }
612 : }
613 : }
614 900 : Q2 = vecmulsqlv(Q, gel(V,2));
615 900 : Q3 = vecmulsqlv(Q, gel(V,3));
616 900 : Q6 = vecmulsqlv(Q, gel(V,6));
617 900 : if (VB)
618 : {
619 24 : Q2B = vecmulsqlv(QB, gel(VB,2));
620 24 : Q3B = vecmulsqlv(QB, gel(VB,3));
621 24 : Q6B = vecmulsqlv(QB, gel(VB,6));
622 : }
623 : /* a,b,c,e = sqrt(q), sqrt(q/2), sqrt(q/3), sqrt(q/6)
624 : * Z[q] = Q[a] + 2^s Q[b] + 3^s Q[c] + 6^s Q[e], with Q[0] = 0 */
625 900 : gel(Z, 1) = onef;
626 900 : gel(Z, 2) = gel(W, 2);
627 900 : gel(Z, 3) = gel(W, 3);
628 900 : gel(Z, 4) = gel(Z, 5) = gel(W, 4);
629 900 : gel(Z, 6) = gel(Z, 7) = gadd(gel(W, 4), gel(V, 6));
630 900 : if (VB)
631 : {
632 24 : ZB = cgetg(sq+1, t_VEC);
633 24 : gel(ZB, 1) = onef;
634 24 : gel(ZB, 2) = gel(WB, 2);
635 24 : gel(ZB, 3) = gel(WB, 3);
636 24 : gel(ZB, 4) = gel(ZB, 5) = gel(WB, 4);
637 24 : gel(ZB, 6) = gel(ZB, 7) = gadd(gel(WB, 4), gel(VB, 6));
638 : }
639 900 : a = 2; b = c = e = 1;
640 131610 : for (q = 8; q <= sq; q++)
641 : { /* Gray code: at most one of a,b,c,d differs (by 1) from previous value */
642 130710 : GEN z = gel(Z, q - 1), zB = NULL;
643 : ulong na, nb, nc, ne, na2, nb2, nc2, ne2;
644 130710 : if (VB) zB = gel(ZB, q - 1);
645 130710 : if ((na = usqrt(q)) != a)
646 7878 : { a = na; na2 = na * na; z = gadd(z, gel(V, na2));
647 7878 : if (VB) zB = gadd(zB, gel(VB, na2)); }
648 122832 : else if ((nb = usqrt(q / 2)) != b)
649 5328 : { b = nb; nb2 = 2 * nb * nb; z = gadd(z, gel(V, nb2));
650 5328 : if (VB) zB = gadd(zB, gel(VB, nb2)); }
651 117504 : else if ((nc = usqrt(q / 3)) != c)
652 4644 : { c = nc; nc2 = 3 * nc * nc; z = gadd(z, gel(V, nc2));
653 4644 : if (VB) zB = gadd(zB, gel(VB, nc2)); }
654 112860 : else if ((ne = usqrt(q / 6)) != e)
655 2664 : { e = ne; ne2 = 6 * ne * ne; z = gadd(z, gel(V, ne2));
656 2664 : if (VB) zB = gadd(zB, gel(VB, ne2)); }
657 130710 : gel(Z, q) = z; if (VB) gel(ZB, q) = zB;
658 : }
659 924 : return v2pack(mkvecn(7, V, W, Q, Q2, Q3, Q6, Z),
660 24 : VB? mkvecn(7, VB, WB, QB, Q2B, Q3B, Q6B, ZB): NULL);
661 : }
662 :
663 : static GEN
664 900 : sumprimeloop(forprime_t *pT, GEN s, long N, GEN data, GEN S, GEN W, GEN WB,
665 : void *E, GEN (*f)(void *, ulong, long))
666 : {
667 900 : pari_sp av = avma;
668 900 : long needlog = data[1], prec0 = data[2], prec1 = data[3];
669 900 : long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
670 900 : ulong p, precp = 0;
671 900 : GEN logp = NULL, SB = WB? S: NULL;
672 :
673 4501146 : while ((p = u_forprime_next(pT)))
674 : {
675 4500246 : GEN u = NULL, ks;
676 4500246 : if (!f) ks = gen_1;
677 : else
678 : {
679 4441200 : ks = f(E, p, prec1);
680 4441200 : if (gequal0(ks)) ks = NULL;
681 4441200 : if (ks && gequal1(ks)) ks = gen_1;
682 : }
683 4500246 : gp[2] = p;
684 4500246 : if (needlog)
685 : {
686 4434918 : if (!logp)
687 882 : logp = logr_abs(utor(p, prec1));
688 : else
689 : { /* Assuming p and precp are odd,
690 : * log p = log(precp) + 2 atanh((p - precp) / (p + precp)) */
691 4434036 : ulong a = p >> 1, b = precp >> 1; /* p = 2a + 1, precp = 2b + 1 */
692 4434036 : GEN z = atanhuu(a - b, a + b + 1, prec1); /* avoid overflow */
693 4434036 : shiftr_inplace(z, 1); logp = addrr(logp, z);
694 : }
695 8869830 : if (ks) u = needlog == 1? powcx(gp, logp, s, prec0)
696 4434912 : : mpexp(gmul(s, logp));
697 : }
698 65328 : else if (ks) u = gpow(gp, s, prec0);
699 4500246 : if (ks)
700 : {
701 4500240 : S = gadd(S, vecmul(gel(W, N / p), ks == gen_1? u: gmul(ks, u)));
702 4500240 : if (WB)
703 : {
704 2106 : GEN w = gel(WB, N / p);
705 2106 : if (ks != gen_1) w = vecmul(ks, w);
706 2106 : SB = gadd(SB, gdiv(w, gmulsg(p, conj_i(u))));
707 : }
708 : }
709 4500246 : precp = p;
710 4500246 : if ((p & 0x1ff) == 1)
711 : {
712 16770 : if (!logp) (void)gc_all(av, SB? 2: 1, &S, &SB);
713 16530 : else (void)gc_all(av, SB? 3: 2, &S, &logp, &SB);
714 : }
715 : }
716 900 : return v2pack(S, SB);
717 : }
718 :
719 : static GEN
720 41550 : add4(GEN a, GEN b, GEN c, GEN d) { return gadd(gadd(a,b), gadd(c,d)); }
721 :
722 : static const long step = 2048;
723 : static int
724 25338 : mksqfloop(long N, long x1, GEN R, GEN RB, GEN *pS, GEN *pSB)
725 : {
726 25338 : GEN V = gel(R,1), Q = gel(R,3), Q2 = gel(R,4);
727 25338 : GEN Q3 = gel(R,5), Q6 = gel(R,6), Z = gel(R,7);
728 25338 : GEN v, VB = NULL, QB = NULL, Q2B = NULL, Q3B = NULL, Q6B = NULL, ZB = NULL;
729 25338 : long x2, j, lv, sq = lg(V)-1;
730 :
731 25338 : if (RB)
732 24 : { VB = gel(RB,1); QB = gel(RB,3); Q2B = gel(RB,4);
733 24 : Q3B = gel(RB,5), Q6B = gel(RB,6); ZB = gel(RB,7); }
734 : /* beware overflow, fuse last two bins (avoid a tiny remainder) */
735 25338 : x2 = (N >= 2*step && N - 2*step >= x1)? x1-1 + step: N;
736 25338 : v = vecfactorsquarefreeu_coprime(x1, x2, mkvecsmall2(2, 3));
737 25338 : lv = lg(v);
738 52026342 : for (j = 1; j < lv; j++)
739 52001004 : if (gel(v,j))
740 : {
741 15804066 : ulong d = x1 - 1 + j; /* squarefree, coprime to 6 */
742 15804066 : GEN t = smallfact(d, gel(v,j), sq, V), u;
743 15804066 : GEN tB = NULL, uB = NULL; /* = f(d) d^s */
744 : long a, b, c, e, q;
745 15804066 : if (!t || gequal0(t)) continue;
746 3347334 : if (VB) tB = vecinv(gmulsg(d, conj_i(t)));
747 : /* warning: gives 1/conj(f(d)) d^(-1-conj(s)), equal to
748 : f(d) d^(-1-conj(s)) only if |f(d)|=1. */
749 : /* S += f(d) * d^s * Z[q] */
750 3347334 : q = N / d;
751 3347334 : if (q == 1)
752 : {
753 1275018 : *pS = gadd(*pS, t); if (VB) *pSB = gadd(*pSB, tB);
754 1275018 : continue;
755 : }
756 2072316 : if (q <= sq) { u = gel(Z, q); if (VB) uB = gel(ZB, q); }
757 : else
758 : {
759 41394 : a = usqrt(q); b = usqrt(q / 2); c = usqrt(q / 3); e = usqrt(q / 6);
760 41394 : u = add4(gel(Q,a), gel(Q2,b), gel(Q3,c), gel(Q6,e));
761 41394 : if (VB) uB = add4(gel(QB,a), gel(Q2B,b), gel(Q3B,c), gel(Q6B,e));
762 : }
763 2072316 : *pS = gadd(*pS, vecmul(t, u)); if (VB) *pSB = gadd(*pSB, vecmul(tB, uB));
764 : }
765 25338 : return x2 == N;
766 : }
767 :
768 : static GEN
769 900 : mkdata(long N, GEN s, long prec)
770 : {
771 900 : long needlog, prec0, prec1, m = mt_nbthreads(), STEP = maxss(N / (m * m), 1);
772 900 : prec1 = prec0 = prec + EXTRAPRECWORD;
773 900 : needlog = get_needlog(s);
774 900 : if (needlog == 1) prec1 = powcx_prec(log2((double)N), s, prec);
775 900 : return mkvecsmalln(5, needlog, prec0, prec1, N, STEP);
776 : }
777 :
778 : static GEN
779 11784 : gp_callUp(void *E, ulong x, long prec)
780 : {
781 11784 : long n[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
782 11784 : n[2] = x; return gp_callprec(E, n, prec);
783 : }
784 :
785 : /* both is
786 : * 0: sum_{n<=N}f(n)n^s
787 : * 1: sum for (f,s) and (conj(f),-1-s)
788 : * 2: sum for (f,s) and (f,-1-s), assuming |f(n)| in {0,1} */
789 : static GEN
790 167 : dirpowerssumfun_i(ulong N, GEN s, void *E, GEN (*f)(void *, ulong, long),
791 : long both, long prec)
792 : {
793 : forprime_t T;
794 : pari_sp av;
795 : GEN onef, zerf, R, RB, W, WB, S, SB, data;
796 : ulong x1;
797 :
798 167 : if ((f && N < 49) || (!f && N < 1000))
799 113 : return smalldirpowerssum(N, s, E, f, both, prec);
800 54 : if (!mk01(E, f, prec, &zerf, &onef)) return mktrivial(both);
801 54 : data = mkdata(N, s, prec); s = gprec_w(s, prec + EXTRAPRECWORD);
802 54 : v2unpack(dirpowsuminit(s, onef, zerf, E, f, data, both), &R, &RB);
803 54 : W = gel(R,2); WB = RB? gel(RB,2): NULL;
804 54 : av = avma; u_forprime_init(&T, lg(W), N);
805 54 : v2unpack(sumprimeloop(&T, s, N, data, zerf, W, WB, E, f), &S, &SB);
806 354 : for(x1 = 1;; x1 += step)
807 : {
808 354 : if (mksqfloop(N, x1, R, RB, &S, &SB))
809 54 : return both? mkvec2(S, conj_i(SB? SB: S)): S;
810 300 : (void)gc_all(av, SB? 2: 1, &S, &SB);
811 : }
812 : }
813 : GEN
814 167 : dirpowerssumfun(ulong N, GEN s, void *E, GEN (*f)(void *, ulong, long),
815 : long both, long prec)
816 : {
817 167 : pari_sp av = avma;
818 167 : return gc_GEN(av, dirpowerssumfun_i(N, s, E, f, both, prec));
819 : }
820 :
821 : GEN
822 107 : dirpowerssum(ulong N, GEN s, long both, long prec)
823 107 : { return dirpowerssumfun(N, s, NULL, NULL, both, prec); }
824 :
825 : GEN
826 132 : dirpowerssum0(GEN N, GEN s, GEN f, long both, long prec)
827 : {
828 132 : if (typ(N) != t_INT) pari_err_TYPE("dirpowerssum", N);
829 126 : if (signe(N) <= 0) N = gen_0;
830 126 : if (!f) return dirpowerssum(itou(N), s, both, prec);
831 60 : if (typ(f) != t_CLOSURE) pari_err_TYPE("dirpowerssum", f);
832 60 : return dirpowerssumfun(itou(N), s, (void*)f, gp_callUp, both, prec);
833 : }
834 :
835 : GEN
836 24984 : parsqf_worker(GEN gk, GEN vR, GEN gN)
837 : {
838 24984 : pari_sp av = avma;
839 : GEN R, RB, onef, S, SB;
840 24984 : long k = itou(gk), N = itou(gN), x1 = 1 + step * k;
841 24984 : v2unpack(vR, &R, &RB); onef = gmael(R,1,1);
842 24984 : S = SB = is_vec_t(typ(onef)) ? zerovec(lg(onef) - 1): gen_0;
843 24984 : (void)mksqfloop(N, x1, R, RB, &S, &SB);
844 24984 : return gc_GEN(av, v2pack(S, RB? SB: NULL));
845 : }
846 :
847 : static GEN
848 4592760 : mycallvec(void *f, ulong n, long prec)
849 : {
850 4592760 : GEN F = (GEN)f;
851 4592760 : if (!f) return gen_1;
852 335280 : if (typ(F) == t_CLOSURE) return gp_callUp(f, n, prec);
853 335280 : return gel(F, (n-1) % (lg(F)-1) + 1);
854 : }
855 :
856 : GEN
857 846 : parsumprimefun_worker(GEN gk, GEN s, GEN zerf, GEN data, GEN vW, GEN f)
858 : {
859 846 : pari_sp av = avma;
860 : forprime_t T;
861 : GEN W, WB, S;
862 846 : long k = itou(gk), sq, N = data[4], STEP = data[5];
863 :
864 846 : v2unpack(vW, &W, &WB); sq = lg(W)-1;
865 846 : if (isintzero(f)) f = NULL;
866 846 : u_forprime_init(&T, k * STEP + sq + 1, minss(N, (k + 1) * STEP + sq));
867 846 : S = sumprimeloop(&T, s, N, data, zerf, W, WB, (void*)f, mycallvec);
868 846 : return gc_GEN(av, S);
869 : }
870 :
871 : static GEN
872 846 : vR_get_vW(GEN vR)
873 : {
874 : GEN R, RB, W, WB;
875 846 : v2unpack(vR, &R, &RB); W = gel(R,2); WB = RB? gel(RB,2): NULL;
876 846 : return v2pack(W, WB);
877 : }
878 :
879 : static GEN
880 846 : halfconj(long both, GEN V)
881 846 : { return both ? mkvec2(gel(V, 1), gconj(gel(V, 2))) : V; }
882 :
883 : static GEN
884 966 : pardirpowerssumfun_i(GEN f, ulong N, GEN s, long both, long prec)
885 : {
886 : GEN worker, worker2, data, vR, onef, zerf;
887 :
888 966 : if ((f && N < 49) || (!f && N < 10000UL))
889 120 : return smalldirpowerssum(N, s, (void*)f, mycallvec, both, prec);
890 846 : if (!mk01((void*)f, mycallvec, prec, &zerf, &onef)) return mktrivial(both);
891 846 : data = mkdata(N, s, prec); s = gprec_w(s, prec + EXTRAPRECWORD);
892 846 : vR = dirpowsuminit(s, onef, zerf, (void*)f, mycallvec, data, both);
893 846 : worker = snm_closure(is_entry("_parsumprimefun_worker"),
894 : mkvecn(5, s, zerf, data, vR_get_vW(vR), f? f: gen_0));
895 846 : worker2 = snm_closure(is_entry("_parsqf_worker"), mkvec2(vR, utoi(N)));
896 1692 : return halfconj(both,
897 846 : gadd(parsum(gen_0, utoipos((N-1) / data[5]), worker),
898 846 : parsum(gen_0, utoipos(maxss((N-1) / step - 1, 0)), worker2)));
899 : }
900 :
901 : GEN
902 966 : pardirpowerssumfun(GEN f, ulong N, GEN s, long both, long prec)
903 : {
904 966 : pari_sp av = avma;
905 966 : return gc_GEN(av, pardirpowerssumfun_i(f, N, s, both, prec));
906 : }
907 : GEN
908 0 : pardirpowerssum0(GEN N, GEN s, GEN f, long both, long prec)
909 : {
910 0 : if (typ(N) != t_INT) pari_err_TYPE("pardirpowerssum", N);
911 0 : return pardirpowerssumfun(f, itou(N), s, both, prec);
912 : }
913 : GEN
914 0 : pardirpowerssum(ulong N, GEN s, long prec)
915 0 : { return pardirpowerssumfun(NULL, N, s, 0, prec); }
|
