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 39166 : Hspec(GEN N)
26 : {
27 39166 : long v2 = Z_lvalrem(N, 2, &N), v2f = v2 >> 1;
28 : GEN t;
29 39166 : if (odd(v2)) { v2f--; N = shifti(N,3); }
30 34973 : 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 39166 : t = addui(3, muliu(subiu(int2n(v2f+1), 3), 2 - kroiu(N,2)));
33 39166 : return mulii(t, hclassno6(N));
34 : }
35 :
36 : static GEN
37 80614 : tauprime_i(ulong t, GEN p2_7, GEN p_9, GEN p, ulong tin)
38 : {
39 80614 : 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 80614 : h = ((t&1UL) == tin)? hclassno6(D): Hspec(D);
42 80610 : a = mulii(powiu(t2,3), addii(p2_7, mulii(t2, subii(shifti(t2,2), p_9))));
43 80610 : return mulii(a, h);
44 : }
45 : GEN
46 93 : ramanujantau_worker(GEN T, GEN p2_7, GEN p_9, GEN p)
47 : {
48 93 : ulong tin = mod4(p) == 3? 1: 0;
49 93 : long i, l = lg(T);
50 93 : GEN s = gen_0;
51 3093 : for (i = 1; i < l; i++) s = addii(s, tauprime_i(T[i], p2_7, p_9, p, tin));
52 93 : return s;
53 : }
54 :
55 : /* B <- {a + k * m : k = 0, ..., (b-a)/m)} */
56 : static void
57 453 : arithprogset(GEN B, ulong a, ulong b, ulong m)
58 : {
59 : long k;
60 27687 : for (k = 1; a <= b; a += m, k++) B[k] = a;
61 453 : setlg(B, k);
62 453 : }
63 : /* sum_{1 <= t <= N} f(t), f has integer values */
64 : static GEN
65 9 : parsum_u(ulong N, GEN worker)
66 : {
67 9 : long a, r, pending = 0, m = usqrt(N);
68 : struct pari_mt pt;
69 9 : GEN v, s = gen_0;
70 : pari_sp av;
71 :
72 9 : mt_queue_start_lim(&pt, worker, m);
73 9 : v = mkvec(cgetg(m + 2, t_VECSMALL)); av = avma;
74 597 : for (r = 1, a = 1; r <= m || pending; r++)
75 : {
76 : long workid;
77 : GEN done;
78 588 : if (r <= m) { arithprogset(gel(v,1), a, N, m); a++; }
79 588 : mt_queue_submit(&pt, 0, r <= m? v: NULL);
80 588 : done = mt_queue_get(&pt, &workid, &pending);
81 588 : if (done) s = gc_INT(av, addii(s, done));
82 : }
83 9 : mt_queue_end(&pt); return s;
84 : }
85 :
86 : static int
87 11448 : tau_parallel(GEN n) { return mt_nbthreads() > 1 && expi(n) > 18; }
88 :
89 : /* Ramanujan tau function for p prime */
90 : static GEN
91 14909 : tauprime(GEN p)
92 : {
93 14909 : pari_sp av = avma;
94 : GEN s, p2, p2_7, p_9, T;
95 : ulong lim;
96 :
97 14909 : if (absequaliu(p, 2)) return utoineg(24);
98 : /* p > 2 */
99 11402 : p2 = sqri(p); p2_7 = mului(7, p2); p_9 = mului(9, p);
100 11402 : lim = itou(sqrtint(p));
101 11402 : if (tau_parallel(p))
102 : {
103 3 : GEN worker = snm_closure(is_entry("_ramanujantau_worker"),
104 : mkvec3(p2_7, p_9, p));
105 3 : s = parsum_u(lim, worker);
106 : }
107 : else
108 : {
109 11399 : pari_sp av2 = avma;
110 11399 : ulong tin = mod4(p) == 3? 1: 0, t;
111 11399 : s = gen_0;
112 89009 : for (t = 1; t <= lim; t++)
113 : {
114 77614 : s = addii(s, tauprime_i(t, p2_7, p_9, p, tin));
115 77610 : if (!(t & 255)) s = gc_INT(av2, s);
116 : }
117 : }
118 : /* 28p^3 - 28p^2 - 90p - 35 */
119 11398 : T = subii(shifti(mulii(p2_7, subiu(p,1)), 2), addiu(mului(90,p), 35));
120 11398 : s = shifti(diviuexact(s, 3), 6);
121 11398 : return gc_INT(av, subii(mulii(mulii(p2, p), T), addui(1, s)));
122 : }
123 :
124 : static GEN
125 105133 : taugen_n_i(ulong t, GEN G, GEN n4)
126 : {
127 105133 : GEN t2 = sqru(t);
128 105133 : return mulii(mfrhopol_eval(G, t2), hclassno6(subii(n4, t2)));
129 : }
130 : GEN
131 360 : taugen_n_worker(GEN T, GEN G, GEN n4)
132 : {
133 360 : long i, l = lg(T);
134 360 : GEN s = gen_0;
135 24594 : for (i = 1; i < l; i++) s = addii(s, taugen_n_i(T[i], G, n4));
136 360 : return s;
137 : }
138 :
139 : static GEN
140 46 : taugen_n(GEN n, GEN G)
141 : {
142 46 : GEN S, r, n4 = shifti(n, 2);
143 46 : ulong t, lim = itou(sqrtremi(n4, &r));
144 :
145 46 : if (r == gen_0) lim--;
146 46 : G = ZX_unscale(G, n);
147 46 : if (tau_parallel(n))
148 : {
149 6 : GEN worker = snm_closure(is_entry("_taugen_n_worker"), mkvec2(G, n4));
150 6 : S = parsum_u(lim, worker);
151 : }
152 : else
153 : {
154 40 : pari_sp av2 = avma;
155 40 : S = gen_0;
156 80939 : for (t = 1; t <= lim; t++)
157 : {
158 80899 : S = addii(S, taugen_n_i(t, G, n4));
159 80899 : if (!(t & 255)) S = gc_INT(av2, S);
160 : }
161 : }
162 46 : S = addii(shifti(S,1), mulii(leading_coeff(G), hclassno6(n4)));
163 46 : return gdivgu(S, 12);
164 : }
165 :
166 : /* ell != 12 */
167 : static GEN
168 20 : newtrace(GEN fan, GEN n, long ell)
169 : {
170 20 : pari_sp av = avma;
171 20 : GEN D = divisors(fan), G = mfrhopol(ell-2), T = taugen_n(n, G);
172 20 : long i, l = lg(D);
173 :
174 60 : for (i = 1; i < l; i++)
175 : {
176 60 : GEN d = gel(D, i), q;
177 60 : long c = cmpii(sqri(d), n);
178 60 : if (c > 0) break;
179 40 : q = powiu(d, ell-1);
180 40 : 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 20 : 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 7088 : ramanujantau(GEN n, long ell)
206 : {
207 7088 : pari_sp av = avma;
208 7088 : GEN T, P, E, G, F = check_arith_all(n, "ramanujantau");
209 : long j, lP;
210 :
211 7081 : if (ell < 12 || ell == 14 || odd(ell)) return gen_0;
212 7081 : if (!F)
213 : {
214 7046 : if (signe(n) <= 0) return gen_0;
215 7039 : F = Z_factor(n); P = gel(F,1);
216 : }
217 : else
218 : {
219 35 : P = gel(F,1);
220 35 : if (lg(P) == 1 || signe(gel(P,1)) <= 0) return gen_0;
221 14 : n = typ(n) == t_VEC? gel(n,1): NULL;
222 : }
223 7053 : if (ell > 26 || ell == 24) return newtrace(F, n? n: factorback(F), ell);
224 : /* dimension 1: tau is multiplicative */
225 7033 : E = gel(F,2); lP = lg(P); T = gen_1;
226 7033 : G = ell == 12? NULL: mfrhopol(ell - 2);
227 21964 : for (j = 1; j < lP; j++)
228 : {
229 14935 : GEN p = gel(P,j), q = powiu(p, ell-1), t, t1, t0 = gen_1;
230 14935 : long k, e = itou(gel(E,j));
231 26 : t1 = t = G? subsi(-1, taugen_n(p, G))
232 14935 : : tauprime(p);
233 20188 : for (k = 1; k < e; k++)
234 : {
235 5257 : GEN t2 = subii(mulii(t, t1), mulii(q, t0));
236 5257 : t0 = t1; t1 = t2;
237 : }
238 14931 : T = mulii(T, t1);
239 : }
240 : #ifdef DEBUG
241 : checkellcong(T, n, ell);
242 : #endif
243 7029 : return gc_INT(av, T);
244 : }
|
