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