Line data Source code
1 : /* Copyright (C) 2000 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 : #include "pari.h"
15 : #include "paripriv.h"
16 :
17 : /******************************************************************/
18 : /* */
19 : /* RAMANUJAN's TAU FUNCTION */
20 : /* */
21 : /******************************************************************/
22 : /* 4|N > 0, not fundamental at 2; 6 * Hurwitz class number in level 2,
23 : * equal to 6*(H(N)+2H(N/4)), H=qfbhclassno */
24 : static GEN
25 26250 : Hspec(GEN N)
26 : {
27 26250 : long v2 = Z_lvalrem(N, 2, &N), v2f = v2 >> 1;
28 : GEN t;
29 26250 : if (odd(v2)) { v2f--; N = shifti(N,3); }
30 23255 : else if (mod4(N)!=3) { v2f--; N = shifti(N,2); }
31 : /* N fundamental at 2, v2f = v2(f) s.t. N = f^2 D, D fundamental */
32 26250 : t = addui(3, muliu(subiu(int2n(v2f+1), 3), 2 - kroiu(N,2)));
33 26250 : return mulii(t, hclassno6(N));
34 : }
35 :
36 : static GEN
37 54130 : tauprime_i(ulong t, GEN p2_7, GEN p_9, GEN p, ulong tin)
38 : {
39 54130 : GEN h, a, t2 = sqru(t), D = shifti(subii(p, t2), 2); /* 4(p-t^2) */
40 : /* t mod 2 != tin <=> D not fundamental at 2 */
41 54130 : h = ((t&1UL) == tin)? hclassno6(D): Hspec(D);
42 54130 : a = mulii(powiu(t2,3), addii(p2_7, mulii(t2, subii(shifti(t2,2), p_9))));
43 54130 : return mulii(a, h);
44 : }
45 : GEN
46 0 : ramanujantau_worker(GEN T, GEN p2_7, GEN p_9, GEN p)
47 : {
48 0 : ulong tin = mod4(p) == 3? 1: 0;
49 0 : long i, l = lg(T);
50 0 : GEN s = gen_0;
51 0 : for (i = 1; i < l; i++) s = addii(s, tauprime_i(T[i], p2_7, p_9, p, tin));
52 0 : return s;
53 : }
54 :
55 : /* B <- {a + k * m : k = 0, ..., (b-a)/m)} */
56 : static void
57 0 : arithprogset(GEN B, ulong a, ulong b, ulong m)
58 : {
59 : long k;
60 0 : for (k = 1; a <= b; a += m, k++) B[k] = a;
61 0 : setlg(B, k);
62 0 : }
63 : /* sum_{1 <= t <= N} f(t), f has integer values */
64 : static GEN
65 0 : parsum_u(ulong N, GEN worker)
66 : {
67 0 : long a, r, pending = 0, m = usqrt(N);
68 : struct pari_mt pt;
69 0 : GEN v, s = gen_0;
70 : pari_sp av;
71 :
72 0 : mt_queue_start_lim(&pt, worker, m);
73 0 : v = mkvec(cgetg(m + 2, t_VECSMALL)); av = avma;
74 0 : for (r = 1, a = 1; r <= m || pending; r++)
75 : {
76 : long workid;
77 : GEN done;
78 0 : if (r <= m) { arithprogset(gel(v,1), a, N, m); a++; }
79 0 : mt_queue_submit(&pt, 0, r <= m? v: NULL);
80 0 : done = mt_queue_get(&pt, &workid, &pending);
81 0 : if (done) s = gc_INT(av, addii(s, done));
82 : }
83 0 : mt_queue_end(&pt); return s;
84 : }
85 :
86 : static int
87 8160 : tau_parallel(GEN n) { return mt_nbthreads() > 1 && expi(n) > 18; }
88 :
89 : /* Ramanujan tau function for p prime */
90 : static GEN
91 10645 : tauprime(GEN p)
92 : {
93 10645 : pari_sp av = avma;
94 : GEN s, p2, p2_7, p_9, T;
95 : ulong lim;
96 :
97 10645 : if (absequaliu(p, 2)) return utoineg(24);
98 : /* p > 2 */
99 8140 : p2 = sqri(p); p2_7 = mului(7, p2); p_9 = mului(9, p);
100 8140 : lim = itou(sqrtint(p));
101 8140 : if (tau_parallel(p))
102 : {
103 0 : GEN worker = snm_closure(is_entry("_ramanujantau_worker"),
104 : mkvec3(p2_7, p_9, p));
105 0 : s = parsum_u(lim, worker);
106 : }
107 : else
108 : {
109 8140 : pari_sp av2 = avma;
110 8140 : ulong tin = mod4(p) == 3? 1: 0, t;
111 8140 : s = gen_0;
112 62270 : for (t = 1; t <= lim; t++)
113 : {
114 54130 : s = addii(s, tauprime_i(t, p2_7, p_9, p, tin));
115 54130 : if (!(t & 255)) s = gc_INT(av2, s);
116 : }
117 : }
118 : /* 28p^3 - 28p^2 - 90p - 35 */
119 8140 : T = subii(shifti(mulii(p2_7, subiu(p,1)), 2), addiu(mului(90,p), 35));
120 8140 : s = shifti(diviuexact(s, 3), 6);
121 8140 : return gc_INT(av, subii(mulii(mulii(p2, p), T), addui(1, s)));
122 : }
123 :
124 : static GEN
125 40445 : taugen_n_i(ulong t, GEN G, GEN n4)
126 : {
127 40445 : GEN t2 = sqru(t);
128 40445 : return mulii(mfrhopol_eval(G, t2), hclassno6(subii(n4, t2)));
129 : }
130 : GEN
131 0 : taugen_n_worker(GEN T, GEN G, GEN n4)
132 : {
133 0 : long i, l = lg(T);
134 0 : GEN s = gen_0;
135 0 : for (i = 1; i < l; i++) s = addii(s, taugen_n_i(T[i], G, n4));
136 0 : return s;
137 : }
138 :
139 : static GEN
140 20 : taugen_n(GEN n, GEN G)
141 : {
142 20 : GEN S, r, n4 = shifti(n, 2);
143 20 : ulong t, lim = itou(sqrtremi(n4, &r));
144 :
145 20 : if (r == gen_0) lim--;
146 20 : G = ZX_unscale(G, n);
147 20 : if (tau_parallel(n))
148 : {
149 0 : GEN worker = snm_closure(is_entry("_taugen_n_worker"), mkvec2(G, n4));
150 0 : S = parsum_u(lim, worker);
151 : }
152 : else
153 : {
154 20 : pari_sp av2 = avma;
155 20 : S = gen_0;
156 40465 : for (t = 1; t <= lim; t++)
157 : {
158 40445 : S = addii(S, taugen_n_i(t, G, n4));
159 40445 : if (!(t & 255)) S = gc_INT(av2, S);
160 : }
161 : }
162 20 : S = addii(shifti(S,1), mulii(leading_coeff(G), hclassno6(n4)));
163 20 : return gdivgu(S, 12);
164 : }
165 :
166 : /* ell != 12 */
167 : static GEN
168 10 : newtrace(GEN fan, GEN n, long ell)
169 : {
170 10 : pari_sp av = avma;
171 10 : GEN D = divisors(fan), G = mfrhopol(ell-2), T = taugen_n(n, G);
172 10 : long i, l = lg(D);
173 :
174 30 : for (i = 1; i < l; i++)
175 : {
176 30 : GEN d = gel(D, i), q;
177 30 : long c = cmpii(sqri(d), n);
178 30 : if (c > 0) break;
179 20 : q = powiu(d, ell-1);
180 20 : if (c < 0) T = gadd(T, q);
181 : else /* d^2 = n */
182 : {
183 0 : T = gadd(T, gmul2n(q, -1));
184 0 : T = gsub(T, gdivgu(mulii(diviiexact(q,d), mfrhopol_eval(G, utoipos(4))), 12));
185 0 : break;
186 : }
187 : }
188 10 : return gc_INT(av, negi(T));
189 : }
190 :
191 : #ifdef DEBUG
192 : static void
193 : checkellcong(GEN T, GEN n, long ell)
194 : {
195 : long V[] = { 0, 691, 0, 3617, 43867, 283*617, 131*593, 0, 657931 };
196 : if (typ(T) != t_INT
197 : || umodiu(subii(T, sumdivk(n, ell-1)), V[ell / 2 - 5]))
198 : pari_err_BUG("ramanujantau");
199 : }
200 : #endif
201 :
202 : /* Ramanujan tau function for weights ell = 12, 16, 18, 20, 22, 26,
203 : * return 0 for <= 0 */
204 : GEN
205 5050 : ramanujantau(GEN n, long ell)
206 : {
207 5050 : pari_sp av = avma;
208 5050 : GEN T, P, E, G, F = check_arith_all(n, "ramanujantau");
209 : long j, lP;
210 :
211 5045 : if (ell < 12 || ell == 14 || odd(ell)) return gen_0;
212 5045 : if (!F)
213 : {
214 5020 : if (signe(n) <= 0) return gen_0;
215 5015 : F = Z_factor(n); P = gel(F,1);
216 : }
217 : else
218 : {
219 25 : P = gel(F,1);
220 25 : if (lg(P) == 1 || signe(gel(P,1)) <= 0) return gen_0;
221 10 : n = typ(n) == t_VEC? gel(n,1): NULL;
222 : }
223 5025 : if (ell > 26 || ell == 24) return newtrace(F, n? n: factorback(F), ell);
224 : /* dimension 1: tau is multiplicative */
225 5015 : E = gel(F,2); lP = lg(P); T = gen_1;
226 5015 : G = ell == 12? NULL: mfrhopol(ell - 2);
227 15670 : for (j = 1; j < lP; j++)
228 : {
229 10655 : GEN p = gel(P,j), q = powiu(p, ell-1), t, t1, t0 = gen_1;
230 10655 : long k, e = itou(gel(E,j));
231 10 : t1 = t = G? subsi(-1, taugen_n(p, G))
232 10655 : : tauprime(p);
233 14410 : for (k = 1; k < e; k++)
234 : {
235 3755 : GEN t2 = subii(mulii(t, t1), mulii(q, t0));
236 3755 : t0 = t1; t1 = t2;
237 : }
238 10655 : T = mulii(T, t1);
239 : }
240 : #ifdef DEBUG
241 : checkellcong(T, n, ell);
242 : #endif
243 5015 : return gc_INT(av, T);
244 : }
|
