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 23114 : dirmuleuler_small(GEN V, GEN v, long n, ulong p, GEN s, long d)
32 : {
33 23114 : long i, j, m = n, D = minss(d+2, lg(s));
34 23114 : ulong q = 1, X = lg(V)-1;
35 :
36 73738 : for (i = 3, q = p; i < D; i++, q *= p) /* q*p does not overflow */
37 : {
38 50624 : GEN aq = gel(s,i);
39 50624 : if (gequal0(aq)) continue;
40 : /* j = 1 */
41 43078 : gel(V,q) = aq;
42 43078 : v[++n] = q;
43 1789991 : for (j = 2; j <= m; j++)
44 : {
45 1746913 : ulong nj = umuluu_le(uel(v,j), q, X);
46 1746913 : if (!nj) continue;
47 135905 : gel(V,nj) = gmul(aq, gel(V,v[j]));
48 135905 : v[++n] = nj;
49 : }
50 : }
51 23114 : return n;
52 : }
53 :
54 : /* ap != 0 for efficiency, p > sqrt(X) */
55 : static void
56 144627 : dirmuleuler_large(GEN V, ulong p, GEN ap)
57 : {
58 144627 : long j, jp, X = lg(V)-1;
59 144627 : gel(V,p) = ap;
60 144627 : for (j = 2, jp = 2*p; jp <= X; j++, jp += p) gel(V,jp) = gmul(ap, gel(V,j));
61 144627 : }
62 :
63 : static ulong
64 8988 : direulertou(GEN a, GEN fl(GEN))
65 : {
66 8988 : 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 8967 : return signe(a)<=0 ? 0: itou(a);
72 : }
73 :
74 : GEN
75 4102 : direuler_bad(void *E, GEN (*eval)(void *,GEN,long), GEN a,GEN b,GEN c, GEN Sbad)
76 : {
77 : ulong au, bu, X, sqrtX, n, p;
78 4102 : pari_sp av0 = avma;
79 : GEN gp, v, V;
80 : forprime_t T;
81 : long i;
82 4102 : au = direulertou(a, gceil);
83 4095 : bu = direulertou(b, gfloor);
84 4088 : X = c ? direulertou(c, gfloor): bu;
85 4081 : if (X == 0) return cgetg(1,t_VEC);
86 4074 : if (bu > X) bu = X;
87 4074 : if (!u_forprime_init(&T, au, bu)) { avma = av0; return mkvec(gen_1); }
88 4053 : v = vecsmall_ei(X, 1);
89 4053 : V = vec_ei(X, 1);
90 4053 : n = 1;
91 4053 : if (Sbad)
92 : {
93 3423 : long l = lg(Sbad);
94 3423 : GEN pbad = gen_1;
95 9030 : for (i = 1; i < l; i++)
96 : {
97 5607 : GEN ai = gel(Sbad,i);
98 : ulong q;
99 5607 : if (typ(ai) != t_VEC || lg(ai) != 3)
100 0 : pari_err_TYPE("direuler [bad primes]",ai);
101 5607 : q = gtou(gel(ai,1));
102 5607 : if (q <= X)
103 : {
104 4501 : long d = ulogint(X, q) + 1;
105 4501 : GEN s = direuler_factor(gel(ai,2), d);
106 4501 : n = dirmuleuler_small(V, v, n, q, s, d);
107 4501 : pbad = muliu(pbad, q);
108 : }
109 : }
110 3423 : Sbad = pbad;
111 : }
112 4053 : p = 1; gp = cgetipos(3); sqrtX = usqrt(X);
113 30051 : while (p <= sqrtX && (p = u_forprime_next(&T)))
114 21966 : if (!Sbad || umodiu(Sbad, p))
115 : {
116 18634 : long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
117 : GEN s;
118 18634 : gp[2] = p; s = eval(E, gp, d);
119 18613 : n = dirmuleuler_small(V, v, n, p, s, d);
120 : }
121 271425 : while ((p = u_forprime_next(&T))) /* sqrt(X) < p <= X */
122 263361 : if (!Sbad || umodiu(Sbad, p))
123 : {
124 : GEN s;
125 262192 : gp[2] = p; s = eval(E, gp, 2);
126 262192 : if (degpol(s) && !gequal0(gel(s,3))) dirmuleuler_large(V, p, gel(s,3));
127 : }
128 4032 : return gerepilecopy(av0,V);
129 : }
130 :
131 : /* return a t_SER or a truncated t_POL to precision n */
132 : GEN
133 86744 : direuler_factor(GEN s, long n)
134 : {
135 86744 : long t = typ(s);
136 86744 : if (is_scalar_t(t))
137 : {
138 49791 : if (!gequal1(s)) err_direuler(s);
139 49784 : return scalarpol_shallow(s,0);
140 : }
141 36953 : switch(t)
142 : {
143 7203 : case t_POL: break; /* no need to RgXn_red */
144 : case t_RFRAC:
145 : {
146 29750 : GEN p = gel(s,1), q = gel(s,2);
147 29750 : q = RgXn_red_shallow(q,n);
148 29750 : s = RgXn_inv(q, n);
149 29750 : if (typ(p) == t_POL && varn(p) == varn(q))
150 : {
151 28 : p = RgXn_red_shallow(p, n);
152 28 : s = RgXn_mul(s, p, n);
153 : }
154 : else
155 29722 : if (!gequal1(p)) s = RgX_Rg_mul(s, p);
156 29750 : if (!signe(s) || !gequal1(gel(s,2))) err_direuler(s);
157 29736 : break;
158 : }
159 : case t_SER:
160 0 : if (!signe(s) || valp(s) || !gequal1(gel(s,2))) err_direuler(s);
161 0 : break;
162 0 : default: pari_err_TYPE("direuler", s);
163 : }
164 36939 : return s;
165 : }
166 :
167 : struct eval_bad
168 : {
169 : void *E;
170 : GEN (*eval)(void *, GEN);
171 : };
172 : static GEN
173 399 : eval_bad(void *E, GEN p, long n)
174 : {
175 399 : struct eval_bad *d = (struct eval_bad*) E;
176 399 : return direuler_factor(d->eval(d->E, p), n);
177 : }
178 : GEN
179 126 : direuler(void *E, GEN (*eval)(void *, GEN), GEN a, GEN b, GEN c)
180 : {
181 : struct eval_bad d;
182 126 : d.E= E; d.eval = eval;
183 126 : return direuler_bad((void*)&d, eval_bad, a, b, c, NULL);
184 : }
|
