Line data Source code
1 : /* Copyright (C) 2018 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 : /** L-functions: values at integers of L-functions **/
17 : /** of primitive quadratic characters **/
18 : /********************************************************************/
19 : #include "pari.h"
20 : #include "paripriv.h"
21 :
22 : static GEN
23 812 : RCpol(long k, long t, GEN c)
24 : {
25 812 : GEN P = cgetg(k+3, t_POL);
26 : long l;
27 :
28 812 : gel(P,k+2) = c;
29 2891 : for(l = 0; l < k; l++)
30 : {
31 2079 : c = diviiexact(mulii(c, muluu(2*k-1 - 2*l, k-l)), mulss(l+1, 2*l-t));
32 2079 : gel(P,k-l+1) = c;
33 : }
34 812 : P[1] = evalsigne(1) | evalvarn(0); return P;
35 : }
36 : static GEN
37 329 : vecRCpol(long r, long d)
38 : {
39 329 : long k, K = d - 1, t = 2*r - 3;
40 329 : GEN v = cgetg(d + 1, t_VEC), c = int2n(2*K);
41 812 : for (k = 0; k <= K; k++)
42 : { /* c = 2^(2K) binomial(n/2,k), an integer */
43 812 : gel(v,k+1) = RCpol(k, t, c);
44 812 : if (k == K) break;
45 483 : c = diviuexact(muliu(c, t - 2*k), 2*k + 2);
46 : }
47 329 : return v;
48 : }
49 : static GEN
50 95739 : euler_sumdiv(GEN q, long v)
51 : {
52 95739 : GEN u = addui(1, q);
53 123536 : for (; v > 1; v--) u = addui(1, mulii(q, u));
54 95739 : return u;
55 : }
56 :
57 : /* [p^{k-1},p^{k-3},...,p^{k-2(d-1)-1}] * (s/p), s = 1 or -1 */
58 : static GEN
59 6587 : vpowp(long k, long d, long p, long s)
60 : {
61 6587 : GEN v = cgetg(d + 1, t_VEC), p2 = sqru(p);
62 : long j;
63 6587 : gel(v, d) = powuu(p, k - 2*d + 1);
64 6587 : if (s == -1 && (p & 3L) == 3) togglesign_safe(&gel(v,d));
65 95739 : for (j = d-1; j >= 1; j--) gel(v, j) = mulii(p2, gel(v, j+1));
66 6587 : return v;
67 : }
68 : static GEN
69 259 : usumdivk_0_all(long k, long d)
70 : {
71 259 : GEN v = cgetg(d + 1, t_COL);
72 : long j;
73 259 : constbern(k >> 1);
74 917 : for (j = 1; j <= d; j++)
75 : {
76 658 : long n = k + 2 - 2*j;
77 658 : gel(v,j) = gdivgs(bernfrac(n), - (n << 1));
78 : }
79 259 : return v;
80 : }
81 : static GEN
82 3115 : usumdivk_fact_all(GEN fa, long k, long d)
83 : {
84 : GEN res, P, E, pow;
85 : long i, j, l;
86 3115 : res = cgetg(d + 1, t_COL);
87 3115 : P = gel(fa, 1); l = lg(P);
88 3115 : E = gel(fa, 2); pow = cgetg(l, t_VEC);
89 8400 : for (i = 1; i < l; i++) gel(pow, i) = vpowp(k, d, P[i], 1);
90 45500 : for (j = 1; j <= d; j++)
91 : {
92 42385 : GEN v = cgetg(l, t_VEC);
93 135254 : for (i = 1; i < l; i++) gel(v,i) = euler_sumdiv(gmael(pow,i,j), E[i]);
94 42385 : gel(res, j) = ZV_prod(v);
95 : }
96 3115 : return res;
97 : }
98 :
99 : /* Hadamard product */
100 : static GEN
101 3801 : RgV_mul(GEN a, GEN b)
102 : {
103 3801 : long j, l = lg(a);
104 3801 : GEN v = cgetg(l, t_COL);
105 57771 : for (j = 1; j < l; j++) gel(v,j) = gmul(gel(a,j), gel(b,j));
106 3801 : return v;
107 : }
108 : static GEN
109 1512 : RgV_multwist(GEN a, GEN P, long k, long dim, long d, long v2, long N4)
110 : {
111 1512 : GEN v = cgetg(dim+1, t_COL);
112 : long j;
113 4879 : for (j = 1; j <= d; j++)
114 : {
115 : GEN z;
116 3367 : gel(v,j) = z = gmul(gel(a,j), gel(P,j));
117 3367 : if (j + d <= dim)
118 : {
119 2086 : if (N4 == 3) z = negi(z);
120 2086 : if (v2) z = shifti(z, (k - 2*j + 1)*v2);
121 2086 : gel(v, j + d) = z;
122 : }
123 : }
124 1512 : return v;
125 : }
126 :
127 : /* r = k - 2*j, 0<=j<d, factor s=an+b, 0<=s<lim. Check if n starts at 0 or 1
128 : * P(D,(an+b)^2), (D-s^2)/N = (D-b^2)/N - 2abn/N - a^2n^2/N and guarantee
129 : * N | D-b^2, N | 2ab, and N | a^2 (except N=8, D odd):
130 : * N=4: a=2, b=0,1\equiv D: D = 0,1 mod 4.
131 : * N=8: a=4, b=2 if D/4 odd, 0 if D/4 even: D = 0 mod 4 or 1 mod 8
132 : * N=12: a=6, b=3 if D odd, 0 if D even: D = 0,1 mod 4
133 : * N=-12: a=6, b=5,1 if D odd, 4,2 if D even: D = 0,1 mod 4
134 : * N=16: a=8, b=7,1 if D = 1 mod 16, 5,3 if D = 9 mod 16: D = 1 mod 8 */
135 : /* Cost: O( sqrt(D)/a d^3 log(D) ) */
136 : static GEN
137 1386 : sigsum(long k, long d, long a, long b, long D, long N, GEN vs, GEN vP)
138 : {
139 : pari_sp av;
140 1386 : GEN S, keep0 = NULL, vPD = RgXV_rescale(vP, stoi(D));
141 1386 : long D2, n, c1, c2, s, lim = usqrt(labs(D));
142 :
143 1386 : D2 = (D - b*b)/N; c1 = (2*a*b)/N; c2 = (a*a)/N;
144 1386 : av = avma; S = zerocol(d);
145 5187 : for (s = b, n = 0; s <= lim; s += a, n++)
146 : {
147 3801 : long Ds = c2 ? D2 - n*(c2*n + c1) : D2 - ((n*(n+1)) >> 1);
148 3801 : GEN v, P = gsubst(vPD, 0, utoi(s*s));
149 3801 : if (vs)
150 1260 : v = gel(vs, Ds+1);
151 : else
152 2541 : v = Ds? usumdivk_fact_all(factoru(Ds), k, d)
153 2541 : : usumdivk_0_all(k,d);
154 3801 : v = RgV_mul(v, P);
155 3801 : if (!s) keep0 = gclone(v); else S = gadd(S, v);
156 3801 : if (gc_needed(av, 1)) S = gerepileupto(av, S);
157 : }
158 1386 : S = gmul2n(S, 1);
159 1386 : if (keep0) { S = gadd(S, keep0); gunclone(keep0); }
160 1386 : return S;
161 : }
162 :
163 : static GEN
164 119 : sigsum4(long k, long d, long D, GEN vs, GEN vP)
165 119 : { return sigsum(k, d, 2, odd(D), D, 4, vs, vP); }
166 :
167 : /* D != 5 (mod 8) */
168 : static GEN
169 210 : sigsum8(long k, long d, long D, GEN vs, GEN vP)
170 : {
171 210 : if (D&1L) return gmul2n(sigsum(k, d, 2, 1, D, 8, vs, vP), -1);
172 210 : return sigsum(k, d, 4, 2*odd(D >> 2), D, 8, vs, vP);
173 : }
174 :
175 : /* D = 0 (mod 3) */
176 : static GEN
177 273 : sigsum12(long k, long d, long D, GEN vs, GEN vP)
178 273 : { return sigsum(k, d, 6, 3*odd(D), D, 12, vs, vP); }
179 :
180 : /* D = 1 (mod 3) */
181 : static GEN
182 35 : sigsumm12(long k, long d, long D, GEN vs, GEN vP)
183 : {
184 35 : long fl = odd(D);
185 35 : GEN res = sigsum(k, d, 6, 4 + fl, D, 12, vs, vP);
186 35 : res = gadd(res, sigsum(k, d, 6, 2 - fl, D, 12, vs, vP));
187 35 : return gmul2n(res, -1);
188 : }
189 :
190 : /* D = 1 (mod 8) */
191 : static GEN
192 357 : sigsum16(long k, long d, long D, GEN vs, GEN vP)
193 : {
194 357 : long fl = (D&15L) == 1;
195 357 : GEN res = sigsum(k, d, 8, 5 + 2*fl, D, 16, vs, vP);
196 357 : return gadd(res, sigsum(k, d, 8, 3 - 2*fl, D, 16, vs, vP));
197 : }
198 :
199 : /* N = 4 (as above), 8 (factor (1+(D/2))), 12 (factor (1+(D/3))),
200 : 16 (only D=1 mod 8). */
201 : static GEN
202 259 : Dpos(long d, long N, long B)
203 : {
204 259 : GEN vD = cgetg(maxss(B, d) + 1, t_VECSMALL);
205 : long D, step, c;
206 259 : switch(N)
207 : {
208 49 : case 4: D = 5; step = 1; break;
209 63 : case 8: D = 8; step = 4; break;
210 56 : case 12: D = 12; step = 3; break;
211 77 : case 16: D = 17; step = 8; break;
212 14 : default: D = 13; step = 3; break; /* -12 */
213 : }
214 1547 : for (c = 1; c <= d || D <= B; D += step)
215 1288 : if (sisfundamental(D)) vD[c++] = D;
216 259 : setlg(vD, c); return vD;
217 : }
218 :
219 : typedef GEN (*SIGMA_F)(long,long,long,GEN,GEN);
220 : static SIGMA_F
221 259 : get_S_even(long N)
222 : {
223 259 : switch(N) {
224 49 : case 4: return sigsum4;
225 63 : case 8: return sigsum8;
226 56 : case 12:return sigsum12;
227 77 : case 16:return sigsum16;
228 14 : default:return sigsumm12; /* -12 */
229 : }
230 : }
231 :
232 : static GEN
233 357 : mfDcoefs(GEN F, GEN vD, long d)
234 : {
235 357 : long l = lg(vD), i;
236 357 : GEN v = mfcoefs(F, vD[l-1], d), w = cgetg(l, t_COL);
237 357 : if (d == 4)
238 252 : for (i = 1; i < l; i++) gel(w, i) = gel(v, (vD[i]>>2)+1);
239 : else
240 854 : for (i = 1; i < l; i++) gel(w, i) = gel(v, vD[i]+1);
241 357 : return w;
242 : }
243 :
244 : static GEN
245 329 : myinverseimage(GEN M, GEN R, GEN *pden)
246 : {
247 329 : GEN c = Q_remove_denom(QM_gauss_i(M, R, 1), pden);/* M*res / den = R */
248 329 : if (!c) pari_err_BUG("theta brackets");
249 329 : return c;
250 : }
251 :
252 : static GEN Lfeq(long D, long k);
253 : static GEN
254 329 : Hcol(GEN k, long r, GEN vD, long d, long N2)
255 : {
256 329 : long i, l = lg(vD);
257 : GEN v;
258 329 : if (r < 5)
259 : {
260 203 : v = mfDcoefs(mfEH(k),vD,d);
261 763 : for (i = 1; i < l; i++)
262 560 : if (N2 != 1 && vD[i] % N2) gel(v,i) = gmul2n(gel(v,i), 1);
263 203 : return v;
264 : }
265 126 : v = cgetg(l, t_COL);
266 896 : for (i = 1; i < l; i++)
267 : {
268 770 : pari_sp av = avma;
269 770 : GEN c = Lfeq(odd(r)? -vD[i]: vD[i], r); /* fundamental */
270 770 : if (N2 != 1 && vD[i] % N2) c = gmul2n(c, 1);
271 770 : gel(v, i) = gerepileupto(av, c);
272 : }
273 126 : return v;
274 : }
275 :
276 : /***********************************************************/
277 : /* Modular form method using Half-Integral Weight forms */
278 : /* Case D > 0 */
279 : /***********************************************************/
280 : static long
281 259 : dimeven(long r, long N)
282 : {
283 259 : switch(N)
284 : {
285 49 : case 4: return r / 6 + 1;
286 70 : case 12: return r / 3 + 1;
287 140 : default: return r / 4 + 1;
288 : }
289 : }
290 : static long
291 259 : muleven(long N) { return (N == 4)? 1: 2; }
292 :
293 : /* L(\chi_D, 1-r) for D > 0 and r > 0 even. */
294 : static GEN
295 259 : modulareven(long D, long r, long N0)
296 : {
297 259 : long B, d, i, l, N = labs(N0);
298 259 : GEN V, vs, R, M, C, den, L, vP, vD, k = uutoQ(2*r+1, 2);
299 259 : SIGMA_F S = get_S_even(N0);
300 :
301 259 : d = dimeven(r, N);
302 259 : B = muleven(N) * mfsturmNgk(N, k);
303 259 : vD = Dpos(d, N0, B);
304 259 : vP = vecRCpol(r, d);
305 259 : l = lg(vD); B = vD[l-1] / N; V = vecfactoru_i(1, B);
306 259 : vs = cgetg(B+2, t_VEC); gel(vs,1) = usumdivk_0_all(r, d);
307 833 : for (i = 1; i <= B; i++) gel(vs, i+1) = usumdivk_fact_all(gel(V,i), r, d);
308 259 : M = cgetg(l, t_MAT);
309 994 : for (i = 1; i < l; i++)
310 : {
311 735 : pari_sp av = avma;
312 735 : gel(M,i) = gerepileupto(av, S(r, d, vD[i], vs, vP));
313 : }
314 259 : M = shallowtrans(M);
315 259 : if (r == 2*d)
316 : { /* r = 2 or (r = 4 and N = 4) */
317 154 : GEN v = mfDcoefs(mfderiv(mfTheta(NULL), d+1), vD, 1);
318 154 : gel(M, d) = gadd(gel(M, d), gdivgu(v, N*(2*d - 1)));
319 : }
320 259 : R = Hcol(k, r, vD, 1, (N == 8 || N0 == 12)? N >> 2: 1);
321 : /* Cost is O(d^2) * bitsize(result) ~ O(d^3.8) [heuristic] */
322 259 : C = myinverseimage(M, R, &den);
323 :
324 : /* Cost: O( sqrt(D)/c d^3 log(D) ), c from findNeven */
325 259 : L = RgV_dotproduct(C, S(r, lg(C)-1, D, NULL, vP));
326 259 : return den? gdiv(L, den): L;
327 : }
328 :
329 : /***********************************************************/
330 : /* Modular form method using Half-Integral Weight forms */
331 : /* Case D < 0 */
332 : /***********************************************************/
333 :
334 : static long
335 70 : dimodd(long r, long kro, long N)
336 : {
337 70 : switch(N)
338 : {
339 0 : case 1: switch (kro)
340 : {
341 0 : case -1:return (r + 3) >> 2;
342 0 : case 0: return (r + 2)/3;
343 0 : case 1: return (r + 1) >> 2;
344 : }
345 7 : case 3: return kro? (r + 1) >> 1: ((r << 1) + 2)/3;
346 28 : case 5: switch (kro)
347 : {
348 0 : case -1:return (3*r + 2) >> 2;
349 28 : case 0: return r;
350 0 : case 1: return (3*r - 1) >> 2;
351 : }
352 7 : case 6: return kro == 1 ? (r + 1) >> 1 : r;
353 28 : default: return r;
354 : }
355 : }
356 :
357 : static GEN
358 70 : Dneg(long n, long kro, long d, long N)
359 : {
360 70 : GEN vD = cgetg(maxss(n, d) + 1, t_VECSMALL);
361 70 : long D, c, step, N2 = odd(N)? N: N>> 1;
362 70 : switch(kro)
363 : {
364 21 : case -1: D = -3; step = 8; break;
365 14 : case 1: D = -7; step = 8; break;
366 35 : default: D = -8; step = 4; break;
367 : }
368 1694 : for (c = 1; D >= -n || c <= d; D -= step)
369 1624 : if (kross(-D, N2) != -1 && sisfundamental(D)) vD[c++] = -D;
370 70 : setlg(vD, c); return vD;
371 : }
372 :
373 : static GEN
374 35 : div4(GEN V)
375 : {
376 35 : long l = lg(V), i;
377 35 : GEN W = cgetg(l, t_VECSMALL);
378 329 : for (i = 1; i < l; i++) W[i] = V[i] >> 2;
379 35 : return W;
380 : }
381 :
382 : static GEN
383 1498 : usumdivktwist_fact_all(GEN fa, long k, long d)
384 : {
385 1498 : GEN V, P, E, pow, res = cgetg(d + 1, t_VEC);
386 : long i, j, l;
387 :
388 1498 : P = gel(fa, 1); l = lg(P);
389 1498 : E = gel(fa, 2);
390 1498 : if (l > 1 && P[1] == 2) { l--; P++; E++; } /* odd part */
391 1498 : pow = cgetg(l, t_VEC);
392 2800 : for (i = 1; i < l; i++) gel(pow, i) = vpowp(k, d, P[i], -1);
393 1498 : V = cgetg(l, t_VEC);
394 4795 : for (j = 1; j <= d; j++)
395 : {
396 6167 : for (i = 1; i < l; i++) gel(V,i) = euler_sumdiv(gmael(pow,i,j), E[i]);
397 3297 : gel(res, j) = ZV_prod(V);
398 : }
399 1498 : return res;
400 : }
401 :
402 : static long
403 70 : mulodd(long N, long kro)
404 : {
405 70 : if (N == 1 || N == 2) return 1;
406 56 : if (kro != 1) return kro? 5: 7;
407 0 : if (N == 3) return 4;
408 0 : if (N == 5) return 5;
409 0 : return 2;
410 : }
411 :
412 : /* Cost: O( sqrt(D)/a d^3 log(D) ) */
413 : static GEN
414 1015 : sigsumtwist(long k, long dim, long a, long b, long Da, long N, GEN vs, GEN vP)
415 : {
416 1015 : GEN vPD, S = zerocol(dim), keep0 = NULL;
417 1015 : long D2, n, c1, c2, s, lim = usqrt(Da), d;
418 : pari_sp av;
419 :
420 1015 : if (N > 2 && kross(Da, N == 6 ? 3 : N) == -1) return S;
421 1015 : d = (dim + 1) >> 1;
422 1015 : vPD = RgXV_rescale(vP, stoi(Da));
423 1015 : D2 = (Da - b*b)/N; c1 = (2*a*b)/N; c2 = (a*a)/N;
424 1015 : av = avma;
425 2527 : for (s = b, n = 0; s <= lim; s += a, n++)
426 : {
427 1512 : long v2, D4, Ds2, Ds = D2 - n*(c2*n + c1); /* (Da - s^2) / N */
428 : GEN v, P;
429 1512 : if (!Ds) continue;
430 1512 : v2 = vals(Ds); Ds2 = Ds >> v2; D4 = Ds2 & 3L; /* (Ds/2^oo) mod 4 */
431 1512 : if (vs)
432 1323 : v = gel(vs, Ds+1);
433 : else
434 189 : v = usumdivktwist_fact_all(factoru(Ds2), k, d);
435 1512 : P = gsubst(vPD, 0, utoi(s*s));
436 1512 : v = RgV_multwist(v, P, k, dim, d, v2, D4);
437 1512 : if (!s) keep0 = gclone(v); else S = gadd(S, v);
438 1512 : if (gc_needed(av, 1)) S = gerepileupto(av, S);
439 : }
440 1015 : S = gmul2n(S, 1);
441 1015 : if (keep0) { S = gadd(S, keep0); gunclone(keep0); }
442 1015 : return gmul2n(S, -2*(d-1));
443 : }
444 :
445 : /* Da = |D|; [sum sigma_r^(1)(Da-s^2), sum sigma_r^(2)(Da-s^2)], N = 1 */
446 : static GEN
447 0 : sigsumtwist11(long k, long dim, long Da, long N, GEN vs, GEN vP)
448 0 : { return sigsumtwist(k, dim, 1, 0, Da, N, vs, vP); }
449 :
450 : /* Da = |D| or |D|/4 */
451 : /* [sum sigma_r^(1)((Da-s^2)/N), sum sigma_r^(2)((Da-s^2)/N)] */
452 : /* Case N|Da; N not necessarily prime. */
453 : static GEN
454 161 : sigsumtwist12p0(long k, long dim, long Da, long N, GEN vs, GEN vP)
455 161 : { return sigsumtwist(k, dim, N, 0, Da, N, vs, vP); }
456 :
457 : /* [sum sigma_r^(1)((Da-s^2)/p), sum sigma_r^(2)((Da-s^2)/p)] */
458 : /* Case p\nmid Da */
459 : /* p = 3: s = +-1 mod 3;
460 : * p = 5: s = +-1 mod 5 if Da = 1 mod 5, s = +-2 mod 5 if Da = 2 mod 5;
461 : * p = 7: s=+-1, +-2, +-3 if Da=1,4,2 mod 7;
462 : * p = 6: s=+-1, +-2, +-3 if Da=1,4,3 mod 6 */
463 : static GEN
464 504 : sigsumtwist12pt(long k, long dim, long Da, long N, GEN vs, GEN vP)
465 : {
466 504 : long t = Da%N, e = 0;
467 : GEN res;
468 504 : if (t == 1) e = 1;
469 210 : else if (t == 4) e = 2;
470 63 : else if (t == 2 || t == 3) e = 3;
471 504 : res = sigsumtwist(k, dim, N, N-e, Da, N, vs, vP);
472 504 : if (N-e != e) res = gadd(res, sigsumtwist(k, dim, N, e, Da, N, vs, vP));
473 504 : return res;
474 : }
475 :
476 : static GEN
477 63 : sigsumtwist12_6(long r, long dim, long Da, long N, GEN vs, GEN vP)
478 : {
479 63 : if (Da%12 == 6) return sigsumtwist12p0(r, dim, Da, N, vs, vP);
480 42 : return sigsumtwist12pt(r, dim, Da, N, vs, vP);
481 : }
482 : static GEN
483 602 : sigsumtwist12_N(long r, long dim, long Da, long N, GEN vs, GEN vP)
484 : {
485 602 : if (Da%N == 0) return sigsumtwist12p0(r, dim, Da, N, vs, vP);
486 462 : return sigsumtwist12pt(r, dim, Da, N, vs, vP);
487 : }
488 :
489 : typedef GEN (*SIGMA_Fodd)(long,long,long,long,GEN,GEN);
490 : static SIGMA_Fodd
491 70 : get_S_odd(long N)
492 : {
493 70 : if (N == 1) return sigsumtwist11;
494 70 : if (N == 6) return sigsumtwist12_6;
495 63 : return sigsumtwist12_N;
496 : }
497 :
498 : /* L(\chi_D, 1-r) for D < 0 and r > 0 odd. */
499 : static GEN
500 70 : modularodd(long D, long r, long N0)
501 : {
502 70 : long B, d, i, l, dim, kro = kross(D, 2), Da = labs(D), N = labs(N0);
503 70 : GEN V, vs, R, M, C, den, L, vP, vD, vD4, k = uutoQ(2*r+1, 2);
504 70 : SIGMA_Fodd S = get_S_odd(N);
505 :
506 70 : dim = dimodd(r, kro, N); d = (dim + 1) >> 1;
507 70 : vP = vecRCpol(r, d);
508 70 : B = mulodd(N, kro) * mfsturmNgk(4*N, k);
509 70 : vD = Dneg(B, kro, dim + 5, N);
510 70 : vD4 = kro ? vD : div4(vD);
511 70 : l = lg(vD); B = vD4[l-1] / N; V = vecfactoru_i(1, B);
512 70 : vs = cgetg(B+2, t_VEC); gel(vs,1) = NULL; /* unused */
513 1379 : for (i = 1; i <= B; i++) gel(vs,i+1) = usumdivktwist_fact_all(gel(V,i), r, d);
514 70 : M = cgetg(l, t_MAT);
515 665 : for (i = 1; i < l; i++)
516 : {
517 595 : pari_sp av = avma;
518 595 : gel(M,i) = gerepileupto(av, S(r, dim, vD4[i], N, vs, vP));
519 : }
520 70 : M = shallowtrans(M);
521 70 : R = Hcol(k, r, vD, kro? 1: 4, odd(N)? N: N >>1);
522 : /* Cost O(d^2) * bitsize(result) ~ O(d^3.7) [heuristic] */
523 70 : C = myinverseimage(M, R, &den);
524 :
525 70 : if (!kro) Da >>= 2;
526 : /* Cost: O( sqrt(D)/c d^3 log(D) ), c from findNodd */
527 70 : L = RgV_dotproduct(C, S(r, lg(C)-1, Da, N, NULL, vP));
528 70 : if (N0 < 0 && (N0 != -6 || Da%3)) den = den? shifti(den,1): gen_2;
529 70 : return den? gdiv(L, den): L;
530 : }
531 :
532 : /********************************************************/
533 : /* Using the Full Functional Equation */
534 : /********************************************************/
535 : /* prod_p (1 - (D/p)p^(-k))
536 : * Cost O( D/log(D) (k log(kD))^mu ), mu = multiplication exponent */
537 : static GEN
538 1848 : Linv(long D, long k, ulong den)
539 : {
540 : pari_sp av;
541 1848 : long s, bit, lim, Da = labs(D), prec;
542 1848 : double km = k - 1, B = (k-0.5) * log(km*Da/17.079) + 12; /* 17.079 ~ 2Pi e */
543 : forprime_t iter;
544 : ulong p;
545 : GEN P, Q;
546 1848 : if (den) B += log((double)den);
547 1848 : bit = maxss((long)(B * k)/(M_LN2 * km), 32) + 32;
548 1848 : prec = nbits2prec(bit);
549 1848 : lim = (long)exp( (B-log(km)) / km ); /* ~ D / (2Pi e) */
550 1848 : u_forprime_init(&iter, 3, lim); av = avma;
551 1848 : s = kross(D, 2);
552 1848 : if (!s) P = real_1(prec);
553 : else
554 : {
555 1113 : Q = real2n(-k, nbits2prec(bit - k));
556 1113 : P = (s == 1)? subir(gen_1, Q): addir(gen_1, Q);
557 : }
558 105742 : while ((p = u_forprime_next(&iter)))
559 : {
560 : long bitnew;
561 : GEN Q;
562 103894 : s = kross(D, p); if (!s) continue;
563 101724 : bitnew = (long)(bit - k * log2(p));
564 101724 : Q = divrr(P, rpowuu(p, k, nbits2prec(maxss(64, bitnew))));
565 101724 : P = s == 1? subrr(P, Q): addrr(P, Q);
566 101724 : if (gc_needed(av,1)) P = gerepileuptoleaf(av, P);
567 : }
568 1848 : return P;
569 : }
570 :
571 : static GEN
572 1848 : myround(GEN z, ulong d)
573 : {
574 : long e;
575 1848 : if (d) z = mulru(z, d);
576 1848 : z = grndtoi(z, &e); if (e >= -4) pari_err_BUG("lfunquad");
577 1848 : return d? Qdiviu(z, d): z;
578 : }
579 :
580 : /* D != 1, k > 2; L(\chi_D, 1-k) using func. eq. */
581 : static GEN
582 1848 : Lfeq(long D, long k)
583 : {
584 : GEN z, res;
585 1848 : long Da, prec, den = 0;
586 :
587 1848 : if ((D > 0 && odd(k)) || (D < 0 && !odd(k))) return gen_0;
588 1848 : Da = labs(D);
589 1848 : if (Da & 3)
590 : {
591 1113 : long d = (Da - 1) >> 1, kd = k / d;
592 1113 : if (odd(kd) && !(k % d) && uisprime(Da)) den = kd * Da;
593 : }
594 735 : else if (Da == 4) den = 2;
595 1848 : z = Linv(D, k, den); prec = realprec(z);
596 1848 : z = mulrr(z, powrs(divru(Pi2n(1, prec), Da), k));
597 1848 : if (Da != 4) { z = mulrr(z, sqrtr_abs(utor(Da,prec))); shiftr_inplace(z,-1); }
598 1848 : res = divrr(mpfactr(k-1, prec), z);
599 1848 : if (odd(k/2)) togglesign(res);
600 1848 : return myround(res, den);
601 : }
602 :
603 : /* heuristic */
604 : static long
605 1407 : usefeq(long D, long k, double c)
606 : {
607 1407 : if (k == 2) return 0;
608 1281 : if (D < 0) { k = 2*k; D = -D; }
609 1281 : return sqrt(D*c) <= k;
610 : }
611 :
612 : static long
613 644 : findNeven(long D, double *c)
614 : {
615 644 : long r = D%3;
616 644 : if (!r) { *c = 3; return 12; }
617 553 : if ((D&7L) == 1) { *c = 2; return 16; }
618 476 : if (!odd(D)) { *c = 2; return 8; }
619 280 : if (r == 1) { *c = 1.5; return -12; }
620 189 : *c = 1; return 4;
621 : }
622 : static long
623 763 : findNodd(long D, long k, double *c)
624 : {
625 763 : long Dmod8 = D&7L, r;
626 763 : if (log(k) > 0.7 * log((double)-D)) { *c = 1; return odd(D)? 2: 1; }
627 343 : if (D%7 == 0 && Dmod8 == 5) { *c = 3.5; return 7; }
628 343 : if (D%6 == 0) { *c = 3; return 6; }
629 315 : if (D%5 == 0) { *c = 2.5; return 5; }
630 294 : if (D%3 == 0) { *c = 1.5; return 3; }
631 245 : if (Dmod8 == 5)
632 : {
633 63 : r = smodss(D, 7);
634 63 : if (r!=1 && r!=2 && r!=4) { *c = 7./6; return -7; }
635 : }
636 182 : if (smodss(D, 3) != 1 && !odd(D)) { *c = 1.5; return -6; }
637 182 : r = smodss(D, 5); if (r != 2 && r != 3) { *c = 5./4; return -5; }
638 70 : *c = 1; return 2;
639 : }
640 :
641 : /* k <= 0 */
642 : static GEN
643 3241 : lfunquadneg_i(long D, long k)
644 : {
645 : double c;
646 : long N;
647 :
648 3241 : if (D == 1) return k == 0 ? gneg(ghalf) : gdivgs(bernfrac(1-k), k-1);
649 2625 : if (!sisfundamental(D)) pari_err_TYPE("lfunquad [D not fundamental]",stoi(D));
650 2625 : if (k == 0) return D < 0? hclassno(stoi(-D)): gen_0;
651 2576 : if ((D > 0 && !odd(k)) || (D < 0 && odd(k))) return gen_0;
652 1498 : if (D == -4) return gmul2n(eulerfrac(-k), -1);
653 1407 : k = 1 - k;
654 1407 : N = D < 0? findNodd(D, k, &c): findNeven(D, &c);
655 1407 : if (usefeq(D, k, c)) return Lfeq(D, k);
656 329 : return D < 0? modularodd(D,k,N): modulareven(D,k,N);
657 : }
658 : /* need k <= 0 and D fundamental */
659 : GEN
660 3241 : lfunquadneg(long D, long k)
661 3241 : { pari_sp av = avma; return gerepileupto(av, lfunquadneg_i(D, k)); }
|
