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 20230 : dirmuleuler_small(GEN V, GEN v, long n, ulong p, GEN s, long d)
32 : {
33 20230 : long i, j, m = n, D = minss(d+2, lg(s));
34 20230 : ulong q = 1, X = lg(V)-1;
35 :
36 68698 : for (i = 3, q = p; i < D; i++, q *= p) /* q*p does not overflow */
37 : {
38 48468 : GEN aq = gel(s,i);
39 48468 : if (gequal0(aq)) continue;
40 : /* j = 1 */
41 40922 : gel(V,q) = aq;
42 40922 : v[++n] = q;
43 1682919 : for (j = 2; j <= m; j++)
44 : {
45 1641997 : ulong nj = umuluu_le(uel(v,j), q, X);
46 1641997 : if (!nj) continue;
47 125657 : gel(V,nj) = gmul(aq, gel(V,v[j]));
48 125657 : v[++n] = nj;
49 : }
50 : }
51 20230 : return n;
52 : }
53 :
54 : /* ap != 0 for efficiency, p > sqrt(X) */
55 : static void
56 146909 : dirmuleuler_large(GEN V, ulong p, GEN ap)
57 : {
58 146909 : long j, jp, X = lg(V)-1;
59 146909 : gel(V,p) = ap;
60 146909 : for (j = 2, jp = 2*p; jp <= X; j++, jp += p) gel(V,jp) = gmul(ap, gel(V,j));
61 146909 : }
62 :
63 : static ulong
64 9030 : direulertou(GEN a, GEN fl(GEN))
65 : {
66 9030 : 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 9009 : return signe(a)<=0 ? 0: itou(a);
72 : }
73 :
74 : static GEN
75 3444 : direuler_Sbad(GEN V, GEN v, GEN Sbad, ulong X, ulong *n)
76 : {
77 3444 : long i, l = lg(Sbad);
78 3444 : GEN pbad = gen_1;
79 9072 : for (i = 1; i < l; i++)
80 : {
81 5628 : GEN ai = gel(Sbad,i);
82 : ulong q;
83 5628 : if (typ(ai) != t_VEC || lg(ai) != 3)
84 0 : pari_err_TYPE("direuler [bad primes]",ai);
85 5628 : q = gtou(gel(ai,1));
86 5628 : if (q <= X)
87 : {
88 4522 : long d = ulogint(X, q) + 1;
89 4522 : GEN s = direuler_factor(gel(ai,2), d);
90 4522 : *n = dirmuleuler_small(V, v, *n, q, s, d);
91 4522 : pbad = muliu(pbad, q);
92 : }
93 : }
94 3444 : return pbad;
95 : }
96 :
97 : GEN
98 560 : direuler_bad(void *E, GEN (*eval)(void *,GEN,long), GEN a,GEN b,GEN c, GEN Sbad)
99 : {
100 : ulong au, bu, X, sqrtX, n, p;
101 560 : pari_sp av0 = avma;
102 : GEN gp, v, V;
103 : forprime_t T;
104 560 : au = direulertou(a, gceil);
105 553 : bu = direulertou(b, gfloor);
106 546 : X = c ? direulertou(c, gfloor): bu;
107 539 : if (X == 0) return cgetg(1,t_VEC);
108 532 : if (bu > X) bu = X;
109 532 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
110 518 : v = vecsmall_ei(X, 1);
111 518 : V = vec_ei(X, 1);
112 518 : n = 1;
113 518 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, X, &n);
114 518 : p = 1; gp = cgetipos(3); sqrtX = usqrt(X);
115 3997 : while (p <= sqrtX && (p = u_forprime_next(&T)))
116 2982 : if (!Sbad || umodiu(Sbad, p))
117 : {
118 2856 : long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
119 : GEN s;
120 2856 : gp[2] = p; s = eval(E, gp, d);
121 2835 : n = dirmuleuler_small(V, v, n, p, s, d);
122 : }
123 80388 : while ((p = u_forprime_next(&T))) /* sqrt(X) < p <= X */
124 79394 : if (!Sbad || umodiu(Sbad, p))
125 : {
126 : GEN s;
127 79387 : gp[2] = p; s = eval(E, gp, 2);
128 79387 : if (degpol(s) && !gequal0(gel(s,3))) dirmuleuler_large(V, p, gel(s,3));
129 : }
130 497 : return gerepilecopy(av0,V);
131 : }
132 :
133 : /* return a t_SER or a truncated t_POL to precision n */
134 : GEN
135 86765 : direuler_factor(GEN s, long n)
136 : {
137 86765 : long t = typ(s);
138 86765 : if (is_scalar_t(t))
139 : {
140 49791 : if (!gequal1(s)) err_direuler(s);
141 49784 : return scalarpol_shallow(s,0);
142 : }
143 36974 : switch(t)
144 : {
145 7224 : case t_POL: break; /* no need to RgXn_red */
146 : case t_RFRAC:
147 : {
148 29750 : GEN p = gel(s,1), q = gel(s,2);
149 29750 : q = RgXn_red_shallow(q,n);
150 29750 : s = RgXn_inv(q, n);
151 29750 : if (typ(p) == t_POL && varn(p) == varn(q))
152 : {
153 28 : p = RgXn_red_shallow(p, n);
154 28 : s = RgXn_mul(s, p, n);
155 : }
156 : else
157 29722 : if (!gequal1(p)) s = RgX_Rg_mul(s, p);
158 29750 : if (!signe(s) || !gequal1(gel(s,2))) err_direuler(s);
159 29736 : break;
160 : }
161 : case t_SER:
162 0 : if (!signe(s) || valp(s) || !gequal1(gel(s,2))) err_direuler(s);
163 0 : break;
164 0 : default: pari_err_TYPE("direuler", s);
165 : }
166 36960 : return s;
167 : }
168 :
169 : struct eval_bad
170 : {
171 : void *E;
172 : GEN (*eval)(void *, GEN);
173 : };
174 : static GEN
175 399 : eval_bad(void *E, GEN p, long n)
176 : {
177 399 : struct eval_bad *d = (struct eval_bad*) E;
178 399 : return direuler_factor(d->eval(d->E, p), n);
179 : }
180 : GEN
181 126 : direuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, GEN c)
182 : {
183 : struct eval_bad d;
184 126 : d.E= E; d.eval = eval;
185 126 : return direuler_bad((void*)&d, eval_bad, a, b, c, NULL);
186 : }
187 :
188 : static GEN
189 25326 : primelist(forprime_t *T, GEN Sbad, long n, long *running)
190 : {
191 25326 : GEN P = cgetg(n+1, t_VECSMALL);
192 : long i, j;
193 228403 : for (i = 1, j = 1; i <= n; i++)
194 : {
195 206633 : ulong p = u_forprime_next(T);
196 206633 : if (!p) { *running = 0; break; }
197 203077 : if (Sbad && umodiu(Sbad, p)==0) continue;
198 198688 : uel(P,j++) = p;
199 : }
200 25326 : setlg(P, j);
201 25326 : return P;
202 : }
203 :
204 : GEN
205 3563 : pardireuler(GEN worker, GEN a, GEN b, GEN c, GEN Sbad)
206 : {
207 : ulong au, bu, X, sqrtX, n, snX, nX;
208 3563 : pari_sp av0 = avma;
209 : GEN v, V;
210 : forprime_t T;
211 : struct pari_mt pt;
212 3563 : long running = 1, pending = 0;
213 3563 : au = direulertou(a, gceil);
214 3563 : bu = direulertou(b, gfloor);
215 3563 : X = c ? direulertou(c, gfloor): bu;
216 3563 : if (X == 0) return cgetg(1,t_VEC);
217 3563 : if (bu > X) bu = X;
218 3563 : if (!u_forprime_init(&T, au, bu)) { set_avma(av0); return mkvec(gen_1); }
219 3556 : v = vecsmall_ei(X, 1);
220 3556 : V = vec_ei(X, 1);
221 3556 : n = 1;
222 3556 : if (Sbad) Sbad = direuler_Sbad(V, v, Sbad, X, &n);
223 3556 : sqrtX = usqrt(X); snX = uprimepi(sqrtX); nX = uprimepi(X);
224 3556 : if (snX)
225 : {
226 3542 : GEN P = primelist(&T, Sbad, snX, &running);
227 3542 : GEN R = gel(closure_callgenvec(worker, mkvec2(P, utoi(X))), 2);
228 3542 : long i, l = lg(P);
229 16415 : for (i = 1; i < l; i++)
230 : {
231 12873 : GEN s = gel(R,i);
232 12873 : n = dirmuleuler_small(V, v, n, uel(P,i), s, lg(s));
233 : }
234 14 : } else snX = 1;
235 3556 : mt_queue_start_lim(&pt, worker, (nX+snX-1)/snX);
236 31716 : while (running || pending)
237 : {
238 : GEN done;
239 24604 : GEN P = running? primelist(&T, Sbad, snX, &running): NULL;
240 24604 : mt_queue_submit(&pt, 0, P ? mkvec2(P, utoi(X)): NULL);
241 24604 : done = mt_queue_get(&pt, NULL, &pending);
242 24604 : if (done)
243 : {
244 21784 : GEN P = gel(done,1), R = gel(done,2);
245 21784 : long j, l = lg(P);
246 207599 : for (j=1; j<l; j++)
247 : {
248 185815 : GEN F = gel(R,j);
249 185815 : if (degpol(F) && !gequal0(gel(F,3)))
250 117579 : dirmuleuler_large(V, uel(P,j), gel(F,3));
251 : }
252 : }
253 : }
254 3556 : mt_queue_end(&pt);
255 3556 : return gerepilecopy(av0,V);
256 : }
|
