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 17619 : conginlist(GEN L, GEN g, void *E, long (*in)(void *, GEN ))
20 : {
21 17619 : pari_sp av = avma;
22 17619 : long i, l = lg(L);
23 17619 : GEN gi = ginv(g);
24 861651 : for (i = 1; i < l; i++)
25 859285 : if (in(E, gmul(gel(L,i), gi))) break;
26 17619 : return gc_long(av, i);
27 : }
28 :
29 : static GEN
30 22372 : normalise(GEN M)
31 : {
32 22372 : long sd = signe(gcoeff(M,2,2));
33 22372 : if (sd < 0 || (!sd && signe(gcoeff(M,1,2)) < 0)) M = ZM_neg(M);
34 22372 : return M;
35 : }
36 :
37 : static void
38 3675 : filln(GEN V, long n, long a, long c)
39 : {
40 : long i, j;
41 17563 : for (j = a + 1, i = 1; i < n; i++)
42 : { /* j != a (mod n) */
43 13888 : gel(V,i) = mkvecsmall2(c, j);
44 13888 : if (++j > n) j = 1;
45 : }
46 3675 : }
47 : /* set v[k+1..k+n-1] or (k == l) append to v; 0 <= a < n */
48 : static GEN
49 3675 : vec_insertn(GEN v, long n, long k, long a, long c)
50 : {
51 3675 : long i, j, l = lg(v), L = l + n-1;
52 3675 : GEN V = cgetg(L, t_VEC);
53 3675 : 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 189917 : for (i = 1; i <= k; i++) gel(V,i) = gel(v,i);
61 3675 : filln(V + i-1, n, a, c);
62 3675 : i += n - 1;
63 100835 : for (j = k + 1; j < l; j++) gel(V,i++) = gel(v,j);
64 : }
65 3675 : return V;
66 : }
67 : /* append the [c,L[i]], i=1..#L to v */
68 : static GEN
69 7350 : vec_appendL(GEN v, GEN L, long c)
70 : {
71 7350 : long i, j, lv, l = lg(L);
72 : GEN w;
73 7350 : if (l == 1) return v;
74 7252 : lv = lg(v); w = cgetg(lv + l -1, typ(v));
75 290164 : for (i = 1; i < lv; i++) gel(w,i) = gel(v,i);
76 24815 : for (j = 1; j < l; i++, j++) gel(w,i) = mkvecsmall2(c, L[j]);
77 7252 : 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 2359 : _isin2(GEN L, long m, long a)
91 : {
92 2359 : pari_sp av = avma;
93 2359 : long k = RgV_isin(L, mkvecsmall2(m,a));
94 2359 : return gc_long(av, k? k: lg(L));
95 : }
96 : static void
97 28301 : get2(GEN x, long *a, long *b) { *a = x[1]; *b = x[2]; }
98 :
99 : static GEN
100 672 : denval(GEN g)
101 : {
102 672 : GEN a = gcoeff(g,1,1), c = gcoeff(g,2,1);
103 672 : return signe(c)? denom_i(gdiv(a,c)): gen_0;
104 : }
105 : /* M * S, S = [0,1;-1,0] */
106 : static GEN
107 448 : mulS(GEN g)
108 : {
109 448 : GEN a = gcoeff(g,1,1), b = gcoeff(g,1,2);
110 448 : GEN c = gcoeff(g,2,1), d = gcoeff(g,2,2);
111 448 : retmkmat22(negi(b), a, negi(d), c);
112 : }
113 : /* remove extra scales and reduce ast to involution */
114 : static GEN
115 105 : rectify(GEN V, GEN ast, GEN gam)
116 : {
117 105 : long n = lg(V)-1, n1, i, def, m, dec;
118 : GEN V1, a1, g1, d, inj;
119 : pari_sp av;
120 :
121 10836 : for(i = 1, def = 0; i <= n; i++)
122 10731 : if (ast[ast[i]] != i) def++;
123 105 : def /= 3;
124 :
125 105 : if (!def) return mkvec3(V, ast, gam);
126 7 : n1 = n + def;
127 7 : g1 = cgetg(n1+1, t_VEC);
128 7 : V1 = cgetg(n1+1, t_VEC);
129 7 : a1 = cgetg(n1+1, t_VECSMALL);
130 7 : d = cgetg(def+1, t_VECSMALL);
131 7 : av = avma;
132 9205 : for (i = m = 1; i <= n; i++)
133 : {
134 9198 : long i2 = ast[i], i3 = ast[i2];
135 9198 : if (i2 > i && i3 > i)
136 : {
137 224 : GEN d1 = denval(ZM_mul(gel(gam,i), gel(V,ast[i])));
138 224 : GEN d2 = denval(ZM_mul(gel(gam,i2), gel(V,ast[i2])));
139 224 : GEN d3 = denval(ZM_mul(gel(gam,i3), gel(V,ast[i3])));
140 224 : if (cmpii(d1,d2) <= 0)
141 105 : d[m++] = cmpii(d1,d3) <= 0? i: i3;
142 : else
143 119 : d[m++] = cmpii(d2,d3) <= 0? i2: i3;
144 : }
145 : }
146 7 : set_avma(av); inj = zero_zv(n);
147 231 : for (i = 1; i <= def; i++) inj[d[i]] = 1;
148 9205 : for (i = 1, dec = 0; i <= n; i++) { dec += inj[i]; inj[i] = i + dec; }
149 9205 : for (i = 1; i <= n; i++)
150 9198 : if (ast[ast[i]] == i)
151 : {
152 8526 : gel(g1, inj[i]) = gel(gam,i);
153 8526 : gel(V1, inj[i]) = gel(V,i);
154 8526 : a1[inj[i]] = inj[ast[i]];
155 : }
156 231 : for (i = 1; i <= def; i++)
157 : {
158 224 : long a = d[i], b = ast[a], c = ast[b];
159 : GEN igc;
160 :
161 224 : gel(V1, inj[b]) = gel(V, b);
162 224 : gel(g1, inj[b]) = normalise(SL2_inv_shallow(gel(gam,a)));
163 224 : a1[inj[b]] = inj[a]-1;
164 :
165 224 : gel(V1, inj[c]) = gel(V, c);
166 224 : gel(g1, inj[c]) = gel(gam, c);
167 224 : a1[inj[c]] = inj[a];
168 :
169 224 : gel(V1, inj[a]-1) = normalise(ZM_mul(gel(gam,a), mulS(gel(V,b))));
170 224 : gel(g1, inj[a]-1) = gel(gam, a);
171 224 : a1[inj[a]-1] = inj[b];
172 :
173 224 : igc = SL2_inv_shallow(gel(gam,c));
174 224 : gel(V1, inj[a]) = normalise(ZM_mul(igc, mulS(gel(V,c))));
175 224 : gel(g1, inj[a]) = normalise(igc);
176 224 : a1[inj[a]] = inj[c];
177 : }
178 7 : return mkvec3(V1, a1, g1);
179 : }
180 : static GEN
181 17563 : vecpop(GEN v)
182 : {
183 17563 : long l = lg(v);
184 17563 : *v++ = evaltyp(t_VEC)|_evallg(1); /* stackdummy */
185 17563 : *v = evaltyp(t_VEC)|_evallg(l-1);
186 17563 : return v;
187 : }
188 :
189 : GEN
190 112 : msfarey(GEN F, void *E, long (*in)(void *, GEN), GEN *pCM)
191 : {
192 112 : pari_sp av = avma, av2, av3;
193 112 : 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 112 : long n = lg(gam)-1, i, k, m, a, l, c, c3;
196 :
197 112 : list = cgetg(n+1, t_VECSMALL);
198 112 : list3 = cgetg(n+1, t_VECSMALL);
199 777 : for (i = c = c3 = 1; i <= n; i++)
200 : {
201 : long t;
202 665 : if (ast[i] == i)
203 203 : t = !isintzero(gtrace(gel(gam,i)));
204 : else
205 462 : t = ast[ast[i]] != i;
206 665 : if (t) list3[c3++] = i; else list[c++] = i;
207 : }
208 112 : setlg(list, c); setlg(list3, c3);
209 112 : if (typ(ast) == t_VEC) ast = ZV_to_zv(ast);
210 112 : av2 = avma;
211 112 : C = M = L = L3 = cgetg(1, t_VEC);
212 112 : B = mkvec(mkvecsmall2(1,1));
213 112 : newcoset(matid(2),1,1);
214 13342 : while(lg(L)-1 + lg(L3)-1)
215 : {
216 17563 : while(lg(L3)-1)
217 : {
218 4333 : get2(gel(L3,1), &m,&a); L3 = vecpop(L3);
219 4333 : av3 = avma;
220 4333 : g = ZM_mul(gel(C,m), gel(gam,a));
221 4333 : k = conginlist(C, g, E, in);
222 4333 : gel(M,m)[a] = k;
223 4333 : if (k < lg(C)) set_avma(av3);
224 : else
225 : {
226 1204 : k = _isin2(B, m, a);
227 1204 : newcoset(g, k, ast[a]);
228 1204 : newcoset(ZM_mul(g,gel(gam,ast[a])), k+n-1, ast[ast[a]]);
229 1204 : B = vecsplice(B, k);
230 : }
231 : }
232 13230 : get2(gel(L,1), &m,&a); L = vecpop(L);
233 13230 : 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 13230 : av3 = avma;
239 13230 : g = ZM_mul(gel(C,m), gel(gam,a));
240 13230 : k = conginlist(C, g, E, in);
241 13230 : gel(M,m)[a] = k; /* class of C[m]*gam[a] */
242 13230 : if (k < lg(C)) set_avma(av3);
243 : else
244 : {
245 1155 : k = _isin2(B, m, a);
246 1155 : newcoset(g,k,ast[a]);
247 1155 : B = vecsplice(B,k);
248 : }
249 : }
250 112 : vecvecsmall_sort_inplace(B, &perm);
251 112 : l = lg(B);
252 112 : V2 = cgetg(l, t_VEC);
253 112 : gam2 = cgetg(l, t_VEC);
254 112 : ast2 = cgetg(l, t_VECSMALL);
255 112 : v2 = cgetg(3, t_VECSMALL);
256 10843 : for (i = 1; i < l; i++)
257 : {
258 10738 : long r, j = perm[i];
259 : GEN ig;
260 10738 : get2(gel(B,i), &m,&a);
261 10738 : r = gel(M,m)[a]; ig = SL2_inv_shallow(gel(C,r));
262 10738 : gel(V2, j) = normalise(ZM_mul(gel(C,m), gel(V,a)));
263 10738 : gel(gam2, j) = normalise(ZM_mul(ZM_mul(gel(C,m), gel(gam,a)), ig));
264 10738 : v2[1] = r; v2[2] = ast[a]; k = vecvecsmall_search(B,v2);
265 10738 : if (k < 0)
266 7 : pari_err(e_MISC, "msfarey: H is not a subgroup of PSL_2(Z)");
267 10731 : ast2[j] = perm[k];
268 : }
269 105 : F = rectify(V2, ast2, gam2);
270 105 : if (pCM) *pCM = mkvec2(C,M);
271 105 : return gc_all(av, pCM? 2: 1, &F, pCM);
272 : }
273 :
274 : GEN
275 14 : mscosets(GEN G, void *E, long (*in)(void *, GEN))
276 : {
277 14 : pari_sp av = avma;
278 : GEN g, L, M;
279 14 : long n = lg(G)-1, i, m, k;
280 14 : g = gel(G,1);
281 14 : L = mkvec(typ(g) == t_VECSMALL? identity_perm(lg(g)-1): gdiv(g,g));
282 14 : M = mkvec(zero_zv(n));
283 35 : for (m = 1; m < lg(L); m++)
284 77 : for (i = 1; i <= n; i++)
285 : {
286 56 : g = gmul(gel(L,m), gel(G,i));
287 56 : mael(M, m, i) = k = conginlist(L, g, E, in);
288 56 : if (k > lg(L)-1) { L = vec_append(L,g); M = vec_append(M, zero_zv(n)); }
289 56 : 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 14 : return gc_GEN(av, mkvec2(L, M));
296 : }
297 :
298 : int
299 2387 : checkfarey_i(GEN F)
300 : {
301 : GEN V, ast, gam;
302 2387 : if (typ(F) != t_VEC || lg(F) < 4) return 0;
303 2387 : V = gel(F,1);
304 2387 : ast = gel(F,2);
305 2387 : gam = gel(F,3);
306 2387 : if (typ(V) != t_VEC
307 2387 : || (typ(ast) != t_VECSMALL && (typ(ast) != t_VEC || !RgV_is_ZV(ast)))
308 112 : || typ(gam) != t_VEC
309 2387 : || lg(V) != lg(ast) || lg(ast) != lg(gam)) return 0;
310 112 : return 1;
311 : }
312 : static int
313 140 : check_inH(GEN inH)
314 : {
315 140 : return (typ(inH) == t_CLOSURE && closure_arity(inH) == 1
316 280 : && !closure_is_variadic(inH));
317 : }
318 : GEN
319 112 : msfarey0(GEN F, GEN code, GEN *pCM)
320 : {
321 112 : if (!checkfarey_i(F)) pari_err_TYPE("msfarey", F);
322 112 : if (!check_inH(code)) pari_err_TYPE("msfarey", code);
323 112 : return msfarey(F, (void*)code, gp_callbool, pCM);
324 : }
325 : GEN
326 28 : mscosets0(GEN V, GEN code)
327 : {
328 28 : if (typ(V) != t_VEC) pari_err_TYPE("mscosets", V);
329 28 : if (!check_inH(code)) pari_err_TYPE("mscosets", code);
330 21 : if (lg(V) == 1) pari_err_TYPE("mscosets [trivial group]", V);
331 14 : return mscosets(V, (void*)code, gp_callbool);
332 : }
|
