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 21 : err_direuler(GEN x)
24 21 : { 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 19670 : dirmuleuler_small(GEN V, GEN v, long n, ulong p, GEN s, long d)
32 : {
33 19670 : long i, j, m = n, D = minss(d+2, lg(s));
34 19670 : ulong q = 1, X = lg(V)-1;
35 :
36 67088 : for (i = 3, q = p; i < D; i++, q *= p) /* q*p does not overflow */
37 : {
38 47418 : GEN aq = gel(s,i);
39 47418 : if (gequal0(aq)) continue;
40 : /* j = 1 */
41 40201 : gel(V,q) = aq;
42 40201 : v[++n] = q;
43 1663221 : for (j = 2; j <= m; j++)
44 : {
45 1623020 : ulong nj = umuluu_le(uel(v,j), q, X);
46 1623020 : if (!nj) continue;
47 125622 : gel(V,nj) = gmul(aq, gel(V,v[j]));
48 125622 : v[++n] = nj;
49 : }
50 : }
51 19670 : return n;
52 : }
53 :
54 : /* ap != 0 for efficiency, p > sqrt(X) */
55 : static void
56 141449 : dirmuleuler_large(GEN V, ulong p, GEN ap)
57 : {
58 141449 : long j, jp, X = lg(V)-1;
59 141449 : gel(V,p) = ap;
60 141449 : for (j = 2, jp = 2*p; jp <= X; j++, jp += p) gel(V,jp) = gmul(ap, gel(V,j));
61 141449 : }
62 :
63 : static ulong
64 8771 : direulertou(GEN a, GEN fl(GEN))
65 : {
66 8771 : 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 8750 : return signe(a)<=0 ? 0: itou(a);
72 : }
73 :
74 : static GEN
75 3458 : direuler_Sbad(GEN V, GEN v, GEN Sbad, ulong *n)
76 : {
77 3458 : long i, l = lg(Sbad);
78 3458 : ulong X = lg(V)-1;
79 3458 : GEN pbad = gen_1;
80 9107 : for (i = 1; i < l; i++)
81 : {
82 5656 : GEN ai = gel(Sbad,i);
83 : ulong q;
84 5656 : if (typ(ai) != t_VEC || lg(ai) != 3)
85 7 : pari_err_TYPE("direuler [bad primes]",ai);
86 5649 : q = gtou(gel(ai,1));
87 5649 : if (q <= X)
88 : {
89 4543 : long d = ulogint(X, q) + 1;
90 4543 : GEN s = direuler_factor(gel(ai,2), d);
91 4543 : *n = dirmuleuler_small(V, v, *n, q, s, d);
92 4543 : pbad = muliu(pbad, q);
93 : }
94 : }
95 3451 : return pbad;
96 : }
97 :
98 : GEN
99 455 : 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 455 : pari_sp av0 = avma;
103 : GEN gp, v, V;
104 : forprime_t T;
105 455 : au = direulertou(a, gceil);
106 448 : bu = direulertou(b, gfloor);
107 441 : X = c ? direulertou(c, gfloor): bu;
108 434 : if (X == 0) return cgetg(1,t_VEC);
109 427 : if (bu > X) bu = X;
110 427 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
111 413 : v = vecsmall_ei(X, 1);
112 413 : V = vec_ei(X, 1);
113 413 : n = 1;
114 413 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
115 406 : p = 1; gp = cgetipos(3); sqrtX = usqrt(X);
116 3101 : while (p <= sqrtX && (p = u_forprime_next(&T)))
117 2310 : if (!Sbad || umodiu(Sbad, p))
118 : {
119 2205 : long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
120 : GEN s;
121 2205 : gp[2] = p; s = eval(E, gp, d);
122 2184 : n = dirmuleuler_small(V, v, n, p, s, d);
123 : }
124 55734 : while ((p = u_forprime_next(&T))) /* sqrt(X) < p <= X */
125 54964 : if (!Sbad || umodiu(Sbad, p))
126 : {
127 : GEN s;
128 54957 : gp[2] = p; s = eval(E, gp, 2); /* s either t_POL or t_SER of val 0 */
129 54957 : if (lg(s) > 3 && !gequal0(gel(s,3)))
130 23051 : dirmuleuler_large(V, p, gel(s,3));
131 : }
132 385 : return gerepilecopy(av0,V);
133 : }
134 :
135 : /* return a t_SER or a truncated t_POL to precision n */
136 : GEN
137 61705 : direuler_factor(GEN s, long n)
138 : {
139 61705 : long t = typ(s);
140 61705 : if (is_scalar_t(t))
141 : {
142 31605 : if (!gequal1(s)) err_direuler(s);
143 31598 : return scalarpol_shallow(s,0);
144 : }
145 30100 : switch(t)
146 : {
147 5467 : case t_POL: break; /* no need to RgXn_red */
148 : case t_RFRAC:
149 : {
150 24325 : GEN p = gel(s,1), q = gel(s,2);
151 24325 : q = RgXn_red_shallow(q,n);
152 24325 : s = RgXn_inv(q, n);
153 24325 : 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 24297 : if (!gequal1(p)) s = RgX_Rg_mul(s, p);
160 24325 : if (!signe(s) || !gequal1(gel(s,2))) err_direuler(s);
161 24311 : break;
162 : }
163 : case t_SER:
164 308 : if (!signe(s) || valp(s) || !gequal1(gel(s,2))) err_direuler(s);
165 308 : break;
166 0 : default: pari_err_TYPE("direuler", s);
167 : }
168 30086 : return s;
169 : }
170 :
171 : struct eval_bad
172 : {
173 : void *E;
174 : GEN (*eval)(void *, GEN);
175 : };
176 : static GEN
177 399 : eval_bad(void *E, GEN p, long n)
178 : {
179 399 : struct eval_bad *d = (struct eval_bad*) E;
180 399 : return direuler_factor(d->eval(d->E, p), n);
181 : }
182 : GEN
183 126 : direuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, GEN c)
184 : {
185 : struct eval_bad d;
186 126 : d.E= E; d.eval = eval;
187 126 : return direuler_bad((void*)&d, eval_bad, a, b, c, NULL);
188 : }
189 :
190 : static GEN
191 25480 : primelist(forprime_t *T, GEN Sbad, long n, long *running)
192 : {
193 25480 : GEN P = cgetg(n+1, t_VECSMALL);
194 : long i, j;
195 229488 : for (i = 1, j = 1; i <= n; i++)
196 : {
197 207585 : ulong p = u_forprime_next(T);
198 207585 : if (!p) { *running = 0; break; }
199 204008 : if (Sbad && umodiu(Sbad, p)==0) continue;
200 199577 : uel(P,j++) = p;
201 : }
202 25480 : setlg(P, j);
203 25480 : return P;
204 : }
205 :
206 : GEN
207 3584 : pardireuler(GEN worker, GEN a, GEN b, GEN c, GEN Sbad)
208 : {
209 : ulong au, bu, X, sqrtX, n, snX, nX;
210 3584 : pari_sp av0 = avma;
211 : GEN v, V;
212 : forprime_t T;
213 : struct pari_mt pt;
214 3584 : long running = 1, pending = 0;
215 3584 : au = direulertou(a, gceil);
216 3584 : bu = direulertou(b, gfloor);
217 3584 : X = c ? direulertou(c, gfloor): bu;
218 3584 : if (X == 0) return cgetg(1,t_VEC);
219 3584 : if (bu > X) bu = X;
220 3584 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
221 3577 : v = vecsmall_ei(X, 1);
222 3577 : V = vec_ei(X, 1);
223 3577 : n = 1;
224 3577 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, &n);
225 3577 : sqrtX = usqrt(X); snX = uprimepi(sqrtX); nX = uprimepi(X);
226 3577 : if (snX)
227 : {
228 3563 : GEN P = primelist(&T, Sbad, snX, &running);
229 3563 : GEN R = gel(closure_callgenvec(worker, mkvec2(P, utoi(X))), 2);
230 3563 : long i, l = lg(P);
231 16506 : for (i = 1; i < l; i++)
232 : {
233 12943 : GEN s = gel(R,i);
234 12943 : n = dirmuleuler_small(V, v, n, uel(P,i), s, lg(s));
235 : }
236 14 : } else snX = 1;
237 3577 : mt_queue_start_lim(&pt, worker, (nX+snX-1)/snX);
238 31912 : while (running || pending)
239 : {
240 : GEN done;
241 24758 : GEN P = running? primelist(&T, Sbad, snX, &running): NULL;
242 24758 : mt_queue_submit(&pt, 0, P ? mkvec2(P, utoi(X)): NULL);
243 24758 : done = mt_queue_get(&pt, NULL, &pending);
244 24758 : if (done)
245 : {
246 21917 : GEN P = gel(done,1), R = gel(done,2);
247 21917 : long j, l = lg(P);
248 208551 : for (j=1; j<l; j++)
249 : {
250 186634 : GEN F = gel(R,j);
251 186634 : if (degpol(F) && !gequal0(gel(F,3)))
252 118398 : dirmuleuler_large(V, uel(P,j), gel(F,3));
253 : }
254 : }
255 : }
256 3577 : mt_queue_end(&pt);
257 3577 : return gerepilecopy(av0,V);
258 : }
|
