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