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 23247 : dirmuleuler_small(GEN V, GEN v, long n, ulong p, GEN s, long d)
33 : {
34 23247 : long i, j, m = n, D = minss(d+2, lg(s));
35 23247 : ulong q = 1, X = lg(V)-1;
36 :
37 79982 : for (i = 3, q = p; i < D; i++, q *= p) /* q*p does not overflow */
38 : {
39 56735 : GEN aq = gel(s,i);
40 56735 : if (gequal0(aq)) continue;
41 : /* j = 1 */
42 48559 : gel(V,q) = aq;
43 48559 : v[++n] = q;
44 2283547 : for (j = 2; j <= m; j++)
45 : {
46 2234988 : ulong nj = umuluu_le(uel(v,j), q, X);
47 2234988 : if (!nj) continue;
48 167566 : gel(V,nj) = gmul(aq, gel(V,v[j]));
49 167566 : v[++n] = nj;
50 : }
51 : }
52 23247 : return n;
53 : }
54 :
55 : /* ap != 0 for efficiency, p > sqrt(X) */
56 : static void
57 193249 : dirmuleuler_large(GEN V, ulong p, GEN ap)
58 : {
59 193249 : long j, jp, X = lg(V)-1;
60 193249 : gel(V,p) = ap;
61 789404 : for (j = 2, jp = 2*p; jp <= X; j++, jp += p) gel(V,jp) = gmul(ap, gel(V,j));
62 193249 : }
63 :
64 : static ulong
65 9877 : direulertou(GEN a, GEN fl(GEN))
66 : {
67 9877 : 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 9856 : return signe(a)<=0 ? 0: itou(a);
73 : }
74 :
75 : static GEN
76 3696 : direuler_Sbad(GEN V, GEN v, GEN Sbad, ulong *n)
77 : {
78 3696 : long i, l = lg(Sbad);
79 3696 : ulong X = lg(V)-1;
80 3696 : GEN pbad = gen_1;
81 9562 : for (i = 1; i < l; i++)
82 : {
83 5901 : GEN ai = gel(Sbad,i);
84 : ulong q;
85 5901 : if (typ(ai) != t_VEC || lg(ai) != 3)
86 14 : pari_err_TYPE("direuler [bad primes]",ai);
87 5887 : q = gtou(gel(ai,1));
88 5880 : if (q <= X)
89 : {
90 4774 : long d = ulogint(X, q) + 1;
91 4774 : GEN s = direuler_factor(gel(ai,2), d);
92 4760 : *n = dirmuleuler_small(V, v, *n, q, s, d);
93 4760 : pbad = muliu(pbad, q);
94 : }
95 : }
96 3661 : return pbad;
97 : }
98 :
99 : GEN
100 504 : 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 504 : pari_sp av0 = avma;
104 : GEN gp, v, V;
105 : forprime_t T;
106 504 : au = direulertou(a, gceil);
107 497 : bu = direulertou(b, gfloor);
108 490 : X = c ? direulertou(c, gfloor): bu;
109 483 : if (X == 0) return cgetg(1,t_VEC);
110 476 : if (bu > X) bu = X;
111 476 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
112 462 : v = vecsmall_ei(X, 1);
113 462 : V = vec_ei(X, 1);
114 462 : n = 1;
115 462 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
116 427 : p = 1; gp = cgetipos(3); sqrtX = usqrt(X);
117 2765 : while (p <= sqrtX && (p = u_forprime_next(&T)))
118 2359 : if (!Sbad || umodiu(Sbad, p))
119 : {
120 2254 : long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
121 : GEN s;
122 2254 : gp[2] = p; s = eval(E, gp, d);
123 2233 : n = dirmuleuler_small(V, v, n, p, s, d);
124 : }
125 57533 : while ((p = u_forprime_next(&T))) /* sqrt(X) < p <= X */
126 57127 : if (!Sbad || umodiu(Sbad, p))
127 : {
128 : GEN s;
129 57120 : gp[2] = p; s = eval(E, gp, 2); /* s either t_POL or t_SER of val 0 */
130 57120 : if (lg(s) > 3 && !gequal0(gel(s,3)))
131 23611 : dirmuleuler_large(V, p, gel(s,3));
132 : }
133 406 : return gerepilecopy(av0,V);
134 : }
135 :
136 : /* return a t_SER or a truncated t_POL to precision n */
137 : GEN
138 64148 : direuler_factor(GEN s, long n)
139 : {
140 64148 : long t = typ(s);
141 64148 : if (is_scalar_t(t))
142 : {
143 33194 : if (!gequal1(s)) err_direuler(s);
144 33180 : return scalarpol_shallow(s,0);
145 : }
146 30954 : switch(t)
147 : {
148 5705 : case t_POL: break; /* no need to RgXn_red */
149 24934 : case t_RFRAC:
150 : {
151 24934 : GEN p = gel(s,1), q = gel(s,2);
152 24934 : q = RgXn_red_shallow(q,n);
153 24934 : s = RgXn_inv(q, n);
154 24934 : 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 24906 : if (!gequal1(p)) s = RgX_Rg_mul(s, p);
161 24934 : if (!signe(s) || !gequal1(gel(s,2))) err_direuler(s);
162 24920 : break;
163 : }
164 308 : case t_SER:
165 308 : if (!signe(s) || valp(s) || !gequal1(gel(s,2))) err_direuler(s);
166 308 : break;
167 7 : default: pari_err_TYPE("direuler", s);
168 : }
169 30933 : return s;
170 : }
171 :
172 : struct eval_bad
173 : {
174 : void *E;
175 : GEN (*eval)(void *, GEN);
176 : };
177 : static GEN
178 420 : eval_bad(void *E, GEN p, long n)
179 : {
180 420 : struct eval_bad *d = (struct eval_bad*) E;
181 420 : return direuler_factor(d->eval(d->E, p), n);
182 : }
183 : GEN
184 133 : direuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, GEN c)
185 : {
186 : struct eval_bad d;
187 133 : d.E= E; d.eval = eval;
188 133 : return direuler_bad((void*)&d, eval_bad, a, b, c, NULL);
189 : }
190 :
191 : static GEN
192 31066 : primelist(forprime_t *T, GEN Sbad, long n, long *running)
193 : {
194 31066 : GEN P = cgetg(n+1, t_VECSMALL);
195 : long i, j;
196 302386 : for (i = 1, j = 1; i <= n; i++)
197 : {
198 275366 : ulong p = u_forprime_next(T);
199 275366 : if (!p) { *running = 0; break; }
200 271320 : if (Sbad && umodiu(Sbad, p)==0) continue;
201 266672 : uel(P,j++) = p;
202 : }
203 31066 : setlg(P, j);
204 31066 : return P;
205 : }
206 :
207 : GEN
208 4053 : pardireuler(GEN worker, GEN a, GEN b, GEN c, GEN Sbad)
209 : {
210 : ulong au, bu, X, sqrtX, n, snX, nX;
211 4053 : pari_sp av0 = avma;
212 : GEN v, V;
213 : forprime_t T;
214 : struct pari_mt pt;
215 4053 : long running = 1, pending = 0;
216 4053 : au = direulertou(a, gceil);
217 4053 : bu = direulertou(b, gfloor);
218 4053 : X = c ? direulertou(c, gfloor): bu;
219 4053 : if (X == 0) return cgetg(1,t_VEC);
220 4053 : if (bu > X) bu = X;
221 4053 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
222 4046 : v = vecsmall_ei(X, 1);
223 4046 : V = vec_ei(X, 1);
224 4046 : n = 1;
225 4046 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
226 4046 : sqrtX = usqrt(X); snX = uprimepi(sqrtX); nX = uprimepi(X);
227 4046 : if (snX)
228 : {
229 4032 : GEN P = primelist(&T, Sbad, snX, &running);
230 4032 : GEN R = gel(closure_callgenvec(worker, mkvec2(P, utoi(X))), 2);
231 4032 : long i, l = lg(P);
232 20286 : for (i = 1; i < l; i++)
233 : {
234 16254 : GEN s = gel(R,i);
235 16254 : n = dirmuleuler_small(V, v, n, uel(P,i), s, lg(s));
236 : }
237 14 : } else snX = 1;
238 4046 : mt_queue_start_lim(&pt, worker, (nX+snX-1)/snX);
239 34614 : while (running || pending)
240 : {
241 : GEN done;
242 30568 : GEN P = running? primelist(&T, Sbad, snX, &running): NULL;
243 30568 : mt_queue_submit(&pt, 0, P ? mkvec2(P, utoi(X)): NULL);
244 30568 : done = mt_queue_get(&pt, NULL, &pending);
245 30568 : if (done)
246 : {
247 27034 : GEN P = gel(done,1), R = gel(done,2);
248 27034 : long j, l = lg(P);
249 277452 : for (j=1; j<l; j++)
250 : {
251 250418 : GEN F = gel(R,j);
252 250418 : if (degpol(F) && !gequal0(gel(F,3)))
253 169638 : dirmuleuler_large(V, uel(P,j), gel(F,3));
254 : }
255 : }
256 : }
257 4046 : mt_queue_end(&pt);
258 4046 : return gerepilecopy(av0,V);
259 : }
260 :
261 : /********************************************************************/
262 : /** **/
263 : /** DIRPOWERS and DIRPOWERSSUM **/
264 : /** **/
265 : /********************************************************************/
266 :
267 : /* [1^B,...,N^B] */
268 : GEN
269 644 : vecpowuu(long N, ulong B)
270 : {
271 : GEN v;
272 : long p, i;
273 : forprime_t T;
274 :
275 644 : if (B <= 8000)
276 : {
277 644 : if (!B) return const_vec(N,gen_1);
278 637 : v = cgetg(N+1, t_VEC); if (N == 0) return v;
279 637 : gel(v,1) = gen_1;
280 637 : if (B == 1)
281 92106 : for (i = 2; i <= N; i++) gel(v,i) = utoipos(i);
282 462 : else if (B == 2)
283 : {
284 : ulong o, s;
285 266 : 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 32145 : for (i = 2, s = 1, o = 3; i <= N; i++, s += o, o += 2)
290 31879 : 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 637 : 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 n^s require log(x) ? */
329 : static long
330 11625 : get_needlog(GEN s)
331 : {
332 11625 : switch(typ(s))
333 : {
334 294 : case t_REAL: return 2; /* yes but not powcx */
335 8316 : case t_COMPLEX: return 1; /* yes using powcx */
336 3015 : default: return 0; /* no */
337 : }
338 : }
339 : /* [1^B,...,N^B] */
340 : GEN
341 11751 : vecpowug(long N, GEN B, long prec)
342 : {
343 11751 : GEN v, logp = NULL;
344 11751 : long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
345 11751 : long p, precp = 2, prec0, prec1, needlog;
346 : forprime_t T;
347 11751 : if (N == 1) return mkvec(gen_1);
348 11744 : if (typ(B) == t_INT && lgefint(B) <= 3 && signe(B) >= 0)
349 133 : return vecpowuu(N, itou(B));
350 11611 : needlog = get_needlog(B);
351 11611 : prec1 = prec0 = prec;
352 11611 : if (needlog == 1) prec1 = powcx_prec(log2((double)N), B, prec);
353 11611 : u_forprime_init(&T, 2, N);
354 11611 : v = const_vec(N, NULL);
355 11611 : gel(v,1) = gen_1;
356 758716 : while ((p = u_forprime_next(&T)))
357 : {
358 : long m, pk, oldpk;
359 : GEN u;
360 747105 : gp[2] = p;
361 747105 : if (needlog)
362 : {
363 90613 : if (!logp)
364 17220 : 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 73393 : ulong a = p >> 1, b = precp >> 1; /* p = 2a + 1, precp = 2b + 1 */
369 73393 : GEN z = atanhuu(a - b, a + b + 1, prec1); /* avoid overflow */
370 73393 : shiftr_inplace(z, 1); logp = addrr(logp, z);
371 : }
372 87089 : u = needlog == 1? powcx(gp, logp, B, prec0)
373 90613 : : mpexp(gmul(B, logp));
374 90613 : if (p == 2) logp = NULL; /* reset: precp must be odd */
375 : }
376 : else
377 656492 : u = gpow(gp, B, prec0);
378 747105 : precp = p;
379 747105 : gel(v,p) = u; /* p^B */
380 747105 : if (prec0 != prec) gel(v,p) = gprec_wtrunc(gel(v,p), prec);
381 1598707 : for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
382 : {
383 851602 : if (pk != p) gel(v,pk) = gmul(gel(v,oldpk), gel(v,p));
384 16711815 : for (m = N/pk; m > 1; m--)
385 15860213 : if (gel(v,m) && m%p) gel(v, m*pk) = gmul(gel(v,m), gel(v,pk));
386 : }
387 : }
388 11611 : return v;
389 : }
390 :
391 : GEN
392 602 : dirpowers(long n, GEN x, long prec)
393 : {
394 : pari_sp av;
395 : GEN v;
396 602 : if (n <= 0) return cgetg(1, t_VEC);
397 588 : av = avma; v = vecpowug(n, x, prec);
398 588 : if (typ(x) == t_INT && lgefint(x) <= 3 && signe(x) >= 0 && cmpiu(x, 2) <= 0)
399 105 : return v;
400 483 : return gerepilecopy(av, v);
401 : }
402 :
403 : /* P = prime divisors of (squarefree) n, V[i] = i^s for i <= sq.
404 : * Return NULL if n is not sq-smooth, else f(n)n^s */
405 : static GEN
406 234087 : smallfact(ulong n, GEN P, ulong sq, GEN V)
407 : {
408 : long i, l;
409 : ulong p, m, o;
410 : GEN c;
411 234087 : if (n <= sq) return gel(V,n);
412 233198 : l = lg(P); m = p = uel(P, l-1); if (p > sq) return NULL;
413 46739 : for (i = l-2; i > 1; i--, m = o) { p = uel(P,i); o = m*p; if (o > sq) break; }
414 46431 : c = gel(V,m); n /= m; /* m <= sq, o = m * p > sq */
415 46431 : if (n > sq) { c = gmul(c, gel(V,p)); n /= p; }
416 46431 : return gmul(c, gel(V,n));
417 : }
418 : static GEN
419 105 : Qtor(GEN x, long prec)
420 105 : { return typ(x) == t_FRAC? fractor(x, prec): x; }
421 : /* sum_{n <= N} f(n)n^s. */
422 : GEN
423 91 : dirpowerssumfun(ulong N, GEN s, void *E, GEN (*f)(void *, ulong, long),
424 : long prec)
425 : {
426 91 : const ulong step = 2048;
427 91 : pari_sp av = avma, av2;
428 : GEN P, V, W, Q, c2, Q2, Q3, Q6, S, Z, logp;
429 : forprime_t T;
430 91 : long gp[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
431 : ulong a, b, c, e, q, x1, n, sq, p, precp;
432 : long prec0, prec1, needlog;
433 :
434 91 : if (!N) return gen_0;
435 91 : if (f)
436 : {
437 14 : if (N < 49)
438 : {
439 7 : V = vecpowug(N, s, prec);
440 7 : S = gen_1;
441 28 : for (n = 2; n <= N; n++)
442 21 : S = gadd(S, gmul(gel(V,n), f(E, n, prec)));
443 7 : return gerepileupto(av, Qtor(S, prec));
444 : }
445 : }
446 77 : else if (N < 1000UL)
447 70 : return gerepileupto(av, Qtor(RgV_sum(vecpowug(N, s, prec)), prec));
448 14 : sq = usqrt(N);
449 14 : V = cgetg(sq+1, t_VEC);
450 14 : W = cgetg(sq+1, t_VEC);
451 14 : Q = cgetg(sq+1, t_VEC);
452 14 : prec1 = prec0 = prec + EXTRAPRECWORD;
453 14 : s = gprec_w(s, prec0);
454 14 : needlog = get_needlog(s);
455 14 : if (needlog == 1) prec1 = powcx_prec(log2((double)N), s, prec);
456 14 : gel(V,1) = gel(W,1) = gel(Q,1) = gen_1;
457 14 : c2 = gpow(gen_2, s, prec0);
458 14 : if (f) c2 = gmul(c2, f(E, 2, prec));
459 14 : gel(V,2) = c2; /* f(2) 2^s */
460 14 : gel(W,2) = Qtor(gaddgs(c2, 1), prec0);
461 14 : gel(Q,2) = Qtor(gaddgs(gsqr(c2), 1), prec0);
462 14 : logp = NULL;
463 2898 : for (n = 3; n <= sq; n++)
464 : {
465 : GEN u;
466 2884 : if (odd(n))
467 : {
468 1442 : gp[2] = n;
469 1442 : if (needlog)
470 : {
471 0 : if (!logp)
472 0 : logp = logr_abs(utor(n, prec1));
473 : else
474 : { /* log n = log(n-2) + 2 atanh(1 / (n - 1)) */
475 0 : GEN z = atanhuu(1, n - 1, prec1);
476 0 : shiftr_inplace(z, 1); logp = addrr(logp, z);
477 : }
478 0 : u = needlog == 1? powcx(gp, logp, s, prec0)
479 0 : : mpexp(gmul(s, logp));
480 :
481 : }
482 : else
483 1442 : u = gpow(gp, s, prec0);
484 1442 : if (f) u = gmul(u, f(E, n, prec0)); /* f(n) n^s */
485 : }
486 : else
487 1442 : u = gmul(c2, gel(V, n>> 1));
488 2884 : gel(V,n) = u; /* f(n) n^s */
489 2884 : gel(W,n) = gadd(gel(W,n-1), gel(V,n)); /* = sum_{i<=n} f(i)i^s */
490 2884 : gel(Q,n) = gadd(gel(Q,n-1), gsqr(gel(V,n))); /* = sum_{i<=n} f(i^2)i^2s */
491 : }
492 14 : Q2 = RgV_Rg_mul(Q, gel(V,2));
493 14 : Q3 = RgV_Rg_mul(Q, gel(V,3));
494 14 : Q6 = RgV_Rg_mul(Q, gel(V,6));
495 14 : precp = 0; logp = NULL; S = gen_0;
496 14 : u_forprime_init(&T, sq + 1, N);
497 14 : av2 = avma;
498 75131 : while ((p = u_forprime_next(&T)))
499 : {
500 : GEN u;
501 75117 : gp[2] = p;
502 75117 : if (needlog)
503 : {
504 0 : if (!logp)
505 0 : logp = logr_abs(utor(p, prec1));
506 : else
507 : { /* log p = log(precp) + 2 atanh((p - precp) / (p + precp)) */
508 0 : ulong a = p >> 1, b = precp >> 1; /* p = 2a + 1, precp = 2b + 1 */
509 0 : GEN z = atanhuu(a - b, a + b + 1, prec1); /* avoid overflow */
510 0 : shiftr_inplace(z, 1); logp = addrr(logp, z);
511 : }
512 0 : u = needlog == 1? powcx(gp, logp, s, prec0)
513 0 : : mpexp(gmul(s, logp));
514 : }
515 : else
516 75117 : u = gpow(gp, s, prec0);
517 75117 : if (f) u = gmul(u, f(E, p, prec0)); /* f(p) p^s */
518 75117 : S = gadd(S, gmul(gel(W, N / p), u));
519 75117 : precp = p;
520 75117 : if ((p & 0x1ff) == 1)
521 : {
522 280 : if (!logp) S = gerepileupto(av2,S);
523 0 : else gerepileall(av2, 2, &S, &logp);
524 : }
525 : }
526 14 : P = mkvecsmall2(2, 3);
527 14 : Z = cgetg(sq+1, t_VEC);
528 : /* a,b,c,e = sqrt(q), sqrt(q/2), sqrt(q/3), sqrt(q/6)
529 : * Z[q] = Q[a] + 2^s Q[b] + 3^s Q[c] + 6^s Q[e], with Q[0] = 0 */
530 14 : gel(Z, 1) = gen_1;
531 14 : gel(Z, 2) = gel(W, 2);
532 14 : gel(Z, 3) = gel(W, 3);
533 14 : gel(Z, 4) = gel(Z, 5) = gel(W,4);
534 14 : gel(Z, 6) = gel(Z, 7) = gadd(gel(W,4), gel(V,6));
535 14 : a = 2; b = c = e = 1;
536 2828 : for (q = 8; q <= sq; q++)
537 : { /* Gray code: at most one of a,b,c,d differs (by 1) from previous value */
538 2814 : GEN z = gel(Z, q - 1);
539 : ulong na, nb, nc, ne;
540 2814 : if ((na = usqrt(q)) != a)
541 161 : { a = na; z = gadd(z, gel(V, na * na)); }
542 2653 : else if ((nb = usqrt(q / 2)) != b)
543 119 : { b = nb; z = gadd(z, gel(V, 2 * nb * nb)); }
544 2534 : else if ((nc = usqrt(q / 3)) != c)
545 91 : { c = nc; z = gadd(z, gel(V, 3 * nc * nc)); }
546 2443 : else if ((ne = usqrt(q / 6)) != e)
547 63 : { e = ne; z = gadd(z, gel(V, 6 * ne * ne)); }
548 2814 : gel(Z,q) = z;
549 : }
550 14 : av2 = avma;
551 14 : for(x1 = 1;; x1 += step)
552 350 : { /* beware overflow, fuse last two bins (avoid a tiny remainder) */
553 364 : ulong j, lv, x2 = (N >= 2*step && N - 2*step >= x1)? x1-1 + step: N;
554 364 : GEN v = vecfactorsquarefreeu_coprime(x1, x2, P);
555 364 : lv = lg(v);
556 770364 : for (j = 1; j < lv; j++) if (gel(v,j))
557 : {
558 234087 : ulong d = x1-1 + j; /* squarefree, coprime to 6 */
559 234087 : GEN t = smallfact(d, gel(v,j), sq, V), u; /* = d^s */
560 234087 : if (!t) continue;
561 : /* S += d^s * Z[q] */
562 47320 : q = N / d;
563 47320 : if (q == 1) { S = gadd(S, t); continue; }
564 30184 : if (q <= sq) u = gel(Z, q);
565 : else
566 : {
567 889 : a = usqrt(q); b = usqrt(q / 2); c = usqrt(q / 3); e = usqrt(q / 6);
568 889 : u = gadd(gadd(gel(Q,a), gel(Q2,b)), gadd(gel(Q3,c), gel(Q6,e)));
569 : }
570 30184 : S = gadd(S, gmul(t, u));
571 : }
572 364 : if (x2 == N) break;
573 350 : S = gerepileupto(av2, S);
574 : }
575 14 : return gerepileupto(av, S);
576 : }
577 : GEN
578 77 : dirpowerssum(ulong N, GEN s, long prec)
579 77 : { return dirpowerssumfun(N, s, NULL, NULL, prec); }
580 : static GEN
581 8799 : gp_callUp(void *E, ulong x, long prec)
582 : {
583 8799 : long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
584 8799 : court[2] = x; return gp_callprec(E, court, prec);
585 : }
586 : GEN
587 56 : dirpowerssum0(GEN N, GEN s, GEN f, long prec)
588 : {
589 56 : if (typ(N) != t_INT) pari_err_TYPE("dirpowerssum", N);
590 49 : if (signe(N) <= 0) return gen_0;
591 35 : if (!f) return dirpowerssum(itou(N), s, prec);
592 14 : if (typ(f) != t_CLOSURE) pari_err_TYPE("dirpowerssum", f);
593 14 : return dirpowerssumfun(itou(N), s, (void*)f, gp_callUp, prec);
594 : }
|
