Line data Source code
1 : /* Copyright (C) 2000, 2012 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 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : static long
19 15102 : conginlist(GEN L, GEN g, void *E, long (*in)(void *, GEN ))
20 : {
21 15102 : pari_sp av = avma;
22 15102 : long i, l = lg(L);
23 15102 : GEN gi = ginv(g);
24 738558 : for (i = 1; i < l; i++)
25 736530 : if (in(E, gmul(gel(L,i), gi))) break;
26 15102 : return gc_long(av, i);
27 : }
28 :
29 : static GEN
30 19176 : normalise(GEN M)
31 : {
32 19176 : long sd = signe(gcoeff(M,2,2));
33 19176 : if (sd < 0 || (!sd && signe(gcoeff(M,1,2)) < 0)) M = ZM_neg(M);
34 19176 : return M;
35 : }
36 :
37 : static void
38 3150 : filln(GEN V, long n, long a, long c)
39 : {
40 : long i, j;
41 15054 : for (j = a + 1, i = 1; i < n; i++)
42 : { /* j != a (mod n) */
43 11904 : gel(V,i) = mkvecsmall2(c, j);
44 11904 : if (++j > n) j = 1;
45 : }
46 3150 : }
47 : /* set v[k+1..k+n-1] or (k == l) append to v; 0 <= a < n */
48 : static GEN
49 3150 : vec_insertn(GEN v, long n, long k, long a, long c)
50 : {
51 3150 : long i, j, l = lg(v), L = l + n-1;
52 3150 : GEN V = cgetg(L, t_VEC);
53 3150 : if (k == l)
54 : {
55 0 : for (i = 1; i < l; i++) gel(V,i) = gel(v,i);
56 0 : filln(V + i-1, n, a, c);
57 : }
58 : else
59 : {
60 162786 : for (i = 1; i <= k; i++) gel(V,i) = gel(v,i);
61 3150 : filln(V + i-1, n, a, c);
62 3150 : i += n - 1;
63 86430 : for (j = k + 1; j < l; j++) gel(V,i++) = gel(v,j);
64 : }
65 3150 : return V;
66 : }
67 : /* append the [c,L[i]], i=1..#L to v */
68 : static GEN
69 6300 : vec_appendL(GEN v, GEN L, long c)
70 : {
71 6300 : long i, j, lv, l = lg(L);
72 : GEN w;
73 6300 : if (l == 1) return v;
74 6216 : lv = lg(v); w = cgetg(lv + l -1, typ(v));
75 248712 : for (i = 1; i < lv; i++) gel(w,i) = gel(v,i);
76 21270 : for (j = 1; j < l; i++, j++) gel(w,i) = mkvecsmall2(c, L[j]);
77 6216 : return w;
78 : }
79 : #define newcoset(g, k, a) \
80 : { \
81 : long _c = lg(C); \
82 : C = vec_append(C, g); \
83 : M = vec_append(M, zero_zv(n)); \
84 : L3= vec_appendL(L3, list3, _c); \
85 : L = vec_appendL(L, list, _c); \
86 : B = vec_insertn(B, n, k, a % n, _c); \
87 : }
88 :
89 : static long
90 2022 : _isin2(GEN L, long m, long a)
91 : {
92 2022 : pari_sp av = avma;
93 2022 : long k = RgV_isin(L, mkvecsmall2(m,a));
94 2022 : return gc_long(av, k? k: lg(L));
95 : }
96 : static void
97 24258 : get2(GEN x, long *a, long *b) { *a = x[1]; *b = x[2]; }
98 :
99 : static GEN
100 576 : denval(GEN g)
101 : {
102 576 : GEN a = gcoeff(g,1,1), c = gcoeff(g,2,1);
103 576 : return signe(c)? denom_i(gdiv(a,c)): gen_0;
104 : }
105 : /* M * S, S = [0,1;-1,0] */
106 : static GEN
107 384 : mulS(GEN g)
108 : {
109 384 : GEN a = gcoeff(g,1,1), b = gcoeff(g,1,2);
110 384 : GEN c = gcoeff(g,2,1), d = gcoeff(g,2,2);
111 384 : retmkmat22(negi(b), a, negi(d), c);
112 : }
113 : /* remove extra scales and reduce ast to involution */
114 : static GEN
115 90 : rectify(GEN V, GEN ast, GEN gam)
116 : {
117 90 : long n = lg(V)-1, n1, i, def, m, dec;
118 : GEN V1, a1, g1, d, inj;
119 : pari_sp av;
120 :
121 9288 : for(i = 1, def = 0; i <= n; i++)
122 9198 : if (ast[ast[i]] != i) def++;
123 90 : def /= 3;
124 :
125 90 : if (!def) return mkvec3(V, ast, gam);
126 6 : n1 = n + def;
127 6 : g1 = cgetg(n1+1, t_VEC);
128 6 : V1 = cgetg(n1+1, t_VEC);
129 6 : a1 = cgetg(n1+1, t_VECSMALL);
130 6 : d = cgetg(def+1, t_VECSMALL);
131 6 : av = avma;
132 7890 : for (i = m = 1; i <= n; i++)
133 : {
134 7884 : long i2 = ast[i], i3 = ast[i2];
135 7884 : if (i2 > i && i3 > i)
136 : {
137 192 : GEN d1 = denval(ZM_mul(gel(gam,i), gel(V,ast[i])));
138 192 : GEN d2 = denval(ZM_mul(gel(gam,i2), gel(V,ast[i2])));
139 192 : GEN d3 = denval(ZM_mul(gel(gam,i3), gel(V,ast[i3])));
140 192 : if (cmpii(d1,d2) <= 0)
141 90 : d[m++] = cmpii(d1,d3) <= 0? i: i3;
142 : else
143 102 : d[m++] = cmpii(d2,d3) <= 0? i2: i3;
144 : }
145 : }
146 6 : set_avma(av); inj = zero_zv(n);
147 198 : for (i = 1; i <= def; i++) inj[d[i]] = 1;
148 7890 : for (i = 1, dec = 0; i <= n; i++) { dec += inj[i]; inj[i] = i + dec; }
149 7890 : for (i = 1; i <= n; i++)
150 7884 : if (ast[ast[i]] == i)
151 : {
152 7308 : gel(g1, inj[i]) = gel(gam,i);
153 7308 : gel(V1, inj[i]) = gel(V,i);
154 7308 : a1[inj[i]] = inj[ast[i]];
155 : }
156 198 : for (i = 1; i <= def; i++)
157 : {
158 192 : long a = d[i], b = ast[a], c = ast[b];
159 : GEN igc;
160 :
161 192 : gel(V1, inj[b]) = gel(V, b);
162 192 : gel(g1, inj[b]) = normalise(SL2_inv_shallow(gel(gam,a)));
163 192 : a1[inj[b]] = inj[a]-1;
164 :
165 192 : gel(V1, inj[c]) = gel(V, c);
166 192 : gel(g1, inj[c]) = gel(gam, c);
167 192 : a1[inj[c]] = inj[a];
168 :
169 192 : gel(V1, inj[a]-1) = normalise(ZM_mul(gel(gam,a), mulS(gel(V,b))));
170 192 : gel(g1, inj[a]-1) = gel(gam, a);
171 192 : a1[inj[a]-1] = inj[b];
172 :
173 192 : igc = SL2_inv_shallow(gel(gam,c));
174 192 : gel(V1, inj[a]) = normalise(ZM_mul(igc, mulS(gel(V,c))));
175 192 : gel(g1, inj[a]) = normalise(igc);
176 192 : a1[inj[a]] = inj[c];
177 : }
178 6 : return mkvec3(V1, a1, g1);
179 : }
180 : static GEN
181 15054 : vecpop(GEN v)
182 : {
183 15054 : long l = lg(v);
184 15054 : *v++ = evaltyp(t_VEC)|_evallg(1); /* stackdummy */
185 15054 : *v = evaltyp(t_VEC)|_evallg(l-1);
186 15054 : return v;
187 : }
188 :
189 : GEN
190 96 : msfarey(GEN F, void *E, long (*in)(void *, GEN), GEN *pCM)
191 : {
192 96 : pari_sp av = avma, av2, av3;
193 96 : GEN V = gel(F,1), ast = gel(F,2), gam = gel(F,3), V2, ast2, gam2;
194 : GEN C, M, L3, L, B, g, list3, list, perm, v2;
195 96 : long n = lg(gam)-1, i, k, m, a, l, c, c3;
196 :
197 96 : list = cgetg(n+1, t_VECSMALL);
198 96 : list3 = cgetg(n+1, t_VECSMALL);
199 666 : for (i = c = c3 = 1; i <= n; i++)
200 : {
201 : long t;
202 570 : if (ast[i] == i)
203 174 : t = !isintzero(gtrace(gel(gam,i)));
204 : else
205 396 : t = ast[ast[i]] != i;
206 570 : if (t) list3[c3++] = i; else list[c++] = i;
207 : }
208 96 : setlg(list, c); setlg(list3, c3);
209 96 : if (typ(ast) == t_VEC) ast = ZV_to_zv(ast);
210 96 : av2 = avma;
211 96 : C = M = L = L3 = cgetg(1, t_VEC);
212 96 : B = mkvec(mkvecsmall2(1,1));
213 96 : newcoset(matid(2),1,1);
214 11436 : while(lg(L)-1 + lg(L3)-1)
215 : {
216 15054 : while(lg(L3)-1)
217 : {
218 3714 : get2(gel(L3,1), &m,&a); L3 = vecpop(L3);
219 3714 : av3 = avma;
220 3714 : g = ZM_mul(gel(C,m), gel(gam,a));
221 3714 : k = conginlist(C, g, E, in);
222 3714 : gel(M,m)[a] = k;
223 3714 : if (k < lg(C)) set_avma(av3);
224 : else
225 : {
226 1032 : k = _isin2(B, m, a);
227 1032 : newcoset(g, k, ast[a]);
228 1032 : newcoset(ZM_mul(g,gel(gam,ast[a])), k+n-1, ast[ast[a]]);
229 1032 : B = vecsplice(B, k);
230 : }
231 : }
232 11340 : get2(gel(L,1), &m,&a); L = vecpop(L);
233 11340 : if (gc_needed(av,2))
234 : {
235 0 : if (DEBUGMEM>1) pari_warn(warnmem,"msfarey, #L = %ld", lg(L)-1);
236 0 : (void)gc_all(av2, 4, &C, &M, &L, &B); L3 = cgetg(1, t_VEC);
237 : }
238 11340 : av3 = avma;
239 11340 : g = ZM_mul(gel(C,m), gel(gam,a));
240 11340 : k = conginlist(C, g, E, in);
241 11340 : gel(M,m)[a] = k; /* class of C[m]*gam[a] */
242 11340 : if (k < lg(C)) set_avma(av3);
243 : else
244 : {
245 990 : k = _isin2(B, m, a);
246 990 : newcoset(g,k,ast[a]);
247 990 : B = vecsplice(B,k);
248 : }
249 : }
250 96 : vecvecsmall_sort_inplace(B, &perm);
251 96 : l = lg(B);
252 96 : V2 = cgetg(l, t_VEC);
253 96 : gam2 = cgetg(l, t_VEC);
254 96 : ast2 = cgetg(l, t_VECSMALL);
255 96 : v2 = cgetg(3, t_VECSMALL);
256 9294 : for (i = 1; i < l; i++)
257 : {
258 9204 : long r, j = perm[i];
259 : GEN ig;
260 9204 : get2(gel(B,i), &m,&a);
261 9204 : r = gel(M,m)[a]; ig = SL2_inv_shallow(gel(C,r));
262 9204 : gel(V2, j) = normalise(ZM_mul(gel(C,m), gel(V,a)));
263 9204 : gel(gam2, j) = normalise(ZM_mul(ZM_mul(gel(C,m), gel(gam,a)), ig));
264 9204 : v2[1] = r; v2[2] = ast[a]; k = vecvecsmall_search(B,v2);
265 9204 : if (k < 0)
266 6 : pari_err(e_MISC, "msfarey: H is not a subgroup of PSL_2(Z)");
267 9198 : ast2[j] = perm[k];
268 : }
269 90 : F = rectify(V2, ast2, gam2);
270 90 : if (pCM) *pCM = mkvec2(C,M);
271 90 : return gc_all(av, pCM? 2: 1, &F, pCM);
272 : }
273 :
274 : GEN
275 12 : mscosets(GEN G, void *E, long (*in)(void *, GEN))
276 : {
277 12 : pari_sp av = avma;
278 : GEN g, L, M;
279 12 : long n = lg(G)-1, i, m, k;
280 12 : g = gel(G,1);
281 12 : L = mkvec(typ(g) == t_VECSMALL? identity_perm(lg(g)-1): gdiv(g,g));
282 12 : M = mkvec(zero_zv(n));
283 30 : for (m = 1; m < lg(L); m++)
284 66 : for (i = 1; i <= n; i++)
285 : {
286 48 : g = gmul(gel(L,m), gel(G,i));
287 48 : mael(M, m, i) = k = conginlist(L, g, E, in);
288 48 : if (k > lg(L)-1) { L = vec_append(L,g); M = vec_append(M, zero_zv(n)); }
289 48 : if (gc_needed(av,2))
290 : {
291 0 : if (DEBUGMEM>1) pari_warn(warnmem,"mscosets, #L = %ld", lg(L)-1);
292 0 : (void)gc_all(av, 2, &M, &L);
293 : }
294 : }
295 12 : return gc_GEN(av, mkvec2(L, M));
296 : }
297 :
298 : int
299 2046 : checkfarey_i(GEN F)
300 : {
301 : GEN V, ast, gam;
302 2046 : if (typ(F) != t_VEC || lg(F) < 4) return 0;
303 2046 : V = gel(F,1);
304 2046 : ast = gel(F,2);
305 2046 : gam = gel(F,3);
306 2046 : if (typ(V) != t_VEC
307 2046 : || (typ(ast) != t_VECSMALL && (typ(ast) != t_VEC || !RgV_is_ZV(ast)))
308 96 : || typ(gam) != t_VEC
309 2046 : || lg(V) != lg(ast) || lg(ast) != lg(gam)) return 0;
310 96 : return 1;
311 : }
312 : static int
313 120 : check_inH(GEN inH)
314 : {
315 120 : return (typ(inH) == t_CLOSURE && closure_arity(inH) == 1
316 240 : && !closure_is_variadic(inH));
317 : }
318 : GEN
319 96 : msfarey0(GEN F, GEN code, GEN *pCM)
320 : {
321 96 : if (!checkfarey_i(F)) pari_err_TYPE("msfarey", F);
322 96 : if (!check_inH(code)) pari_err_TYPE("msfarey", code);
323 96 : return msfarey(F, (void*)code, gp_callbool, pCM);
324 : }
325 : GEN
326 24 : mscosets0(GEN V, GEN code)
327 : {
328 24 : if (typ(V) != t_VEC) pari_err_TYPE("mscosets", V);
329 24 : if (!check_inH(code)) pari_err_TYPE("mscosets", code);
330 18 : if (lg(V) == 1) pari_err_TYPE("mscosets [trivial group]", V);
331 12 : return mscosets(V, (void*)code, gp_callbool);
332 : }
|
