1 /* Copyright (C) 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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
14
15 #include "pari.h"
16 #include "paripriv.h"
17
18 /*To be moved elsewhere */
19
20 static long
Z_trunc(GEN z)21 Z_trunc(GEN z)
22 {
23 long s = signe(z);
24 return s==0 ? 0: (long)(s*mod2BIL(z));
25 }
26
27 static GEN
ZV_trunc_to_zv(GEN z)28 ZV_trunc_to_zv(GEN z)
29 {
30 long i, l = lg(z);
31 GEN x = cgetg(l, t_VECSMALL);
32 for (i=1; i<l; i++) x[i] = Z_trunc(gel(z,i));
33 return x;
34 }
35
36 /* returns scalar product of vectors x and y with respect to Gram-matrix F */
37 static long
scp(GEN x,GEN F,GEN y)38 scp(GEN x, GEN F, GEN y)
39 {
40 long i, j, n = lg(F)-1;
41 ulong sum = 0;
42 for (i = 1; i <= n; ++i)
43 {
44 ulong xi = uel(x,i);
45 if (xi)
46 {
47 GEN Fi = gel(F, i);
48 for (j = 1; j <= n; ++j) sum += xi * uel(Fi,j) * uel(y,j);
49 }
50 }
51 return sum;
52 }
53
54 static GEN
ZM_trunc_to_zm(GEN z)55 ZM_trunc_to_zm(GEN z)
56 {
57 long i, l = lg(z);
58 GEN x = cgetg(l,t_MAT);
59 for (i=1; i<l; i++) gel(x,i) = ZV_trunc_to_zv(gel(z,i));
60 return x;
61 }
62
63 static GEN
zm_divmod(GEN A,GEN B,ulong p)64 zm_divmod(GEN A, GEN B, ulong p)
65 {
66 pari_sp av = avma;
67 GEN Ap = zm_to_Flm(A,p), Bp = zm_to_Flm(B,p);
68 GEN C = Flm_center(Flm_mul(Flm_inv(Ap, p), Bp, p), p, p>>1);
69 return gerepileupto(av, C);
70 }
71
72 static int
zv_canon(GEN V)73 zv_canon(GEN V)
74 {
75 long l = lg(V), j, k;
76 for (j = 1; j < l && V[j] == 0; ++j);
77 if (j < l && V[j] < 0)
78 {
79 for (k = j; k < l; ++k) V[k] = -V[k];
80 return -1;
81 }
82 return 1;
83 }
84 static GEN
ZM_to_zm_canon(GEN V)85 ZM_to_zm_canon(GEN V)
86 {
87 GEN W = ZM_to_zm(V);
88 long i, l = lg(W);
89 for (i=1; i<l; i++) (void)zv_canon(gel(W,i));
90 return W;
91 }
92
93 static GEN
zm_apply_zm(GEN M,GEN U)94 zm_apply_zm(GEN M, GEN U)
95 { return zm_mul(zm_transpose(U),zm_mul(M, U)); }
96
97 static GEN
zmV_apply_zm(GEN v,GEN U)98 zmV_apply_zm(GEN v, GEN U)
99 {
100 long i, l;
101 GEN V = cgetg_copy(v, &l);
102 for (i=1; i<l; i++) gel(V,i) = zm_apply_zm(gel(v,i), U);
103 return V;
104 }
105
106 /********************************************************************/
107 /** **/
108 /** QFAUTO/QFISOM ported from B. Souvignier ISOM program **/
109 /** **/
110 /********************************************************************/
111
112 /* This is a port by Bill Allombert of the program ISOM by Bernt Souvignier
113 which implements
114 W. PLESKEN, B. SOUVIGNIER, Computing Isometries of Lattices,
115 Journal of Symbolic Computation, Volume 24, Issues 3-4, September 1997,
116 Pages 327-334, ISSN 0747-7171, 10.1006/jsco.1996.0130.
117 (http://www.sciencedirect.com/science/article/pii/S0747717196901303)
118
119 We thank Professor Souvignier for giving us permission to port his code.
120 */
121
122 struct group
123 {
124 GEN ord, ng, nsg, g;
125 };
126
127 struct fingerprint
128 {
129 GEN diag, per, e;
130 };
131
132 struct qfauto
133 {
134 long dim;
135 GEN F, U, V, W, v;
136 ulong p;
137 };
138
139 struct qfcand
140 {
141 long cdep;
142 GEN comb, bacher_pol;
143 };
144
145 static long
possible(GEN F,GEN Ftr,GEN V,GEN W,GEN per,long I,long J)146 possible(GEN F, GEN Ftr, GEN V, GEN W, GEN per, long I, long J)
147 {
148 long i, j, k, count = 0, n = lg(W)-1, f = lg(F)-1;
149
150 for (j = 1; j <= n; j++)
151 {
152 GEN Wj = gel(W,j), Vj = gel(V,j);
153 i = I+1;
154 /* check vector length */
155 for (k = 1; k <= f && i > I && Wj[k] == mael3(F,k,J,J); k++)
156 for (i = 1; i <= I; i++) /* check scalar products with basis vectors */
157 if (zv_dotproduct(Vj,gmael(Ftr,k,per[i])) != mael(gel(F,k),J,per[i]))
158 break;
159 if (k == f+1 && i > I) ++count;
160 /* same for the negative vector */
161 i = I+1;
162 for (k = 1; k <= f && i > I && Wj[k] == mael3(F,k,J,J); k++)
163 for (i = 1; i <= I ; i++)
164 if (zv_dotproduct(Vj,gmael(Ftr,k,per[i])) != -mael(gel(F,k),J,per[i]))
165 break;
166 if (k == f+1 && i > I) ++count;
167 }
168 return count;
169 }
170
171 static void
fingerprint(struct fingerprint * fp,struct qfauto * qf)172 fingerprint(struct fingerprint *fp, struct qfauto *qf)
173 {
174 pari_sp av;
175 GEN V=qf->V, W=qf->W, F=qf->F, Mf, Ftr;
176 long i, j, k, min, dim = qf->dim, n = lg(V)-1, f = lg(F)-1;
177 fp->per = identity_perm(dim);
178 fp->e = cgetg(dim+1, t_VECSMALL);
179 fp->diag = cgetg(dim+1, t_VECSMALL);
180 av = avma;
181 Ftr = cgetg(f+1,t_VEC);
182 for (i = 1; i <= f; i++) gel(Ftr,i) = zm_transpose(gel(F,i));
183 Mf = zero_Flm_copy(dim, dim);
184 /* the first row of the fingerprint has as entry nr. i the number of
185 vectors, which have the same length as the i-th basis-vector with
186 respect to every invariant form */
187 for (j = 1; j <= n; j++)
188 {
189 GEN Wj = gel(W,j);
190 for (i = 1; i <= dim; i++)
191 {
192 for (k = 1; k <= f && Wj[k] == mael3(F,k,i,i); ++k);
193 if (k == f+1) mael(Mf,1,i) += 2;
194 }
195 }
196 for (i = 1; i <= dim-1; ++i)
197 { /* a minimal entry != 0 in the i-th row is chosen */
198 min = i;
199 for (j = i+1; j <= dim; j++)
200 if (mael(Mf,i,fp->per[j]) < mael(Mf,i,fp->per[min])) min = j;
201 lswap(fp->per[i],fp->per[min]);
202 /* the column below the minimal entry is set to 0 */
203 for (j = i+1; j <= dim; j++) mael(Mf,j,fp->per[i]) = 0;
204 /* compute the row i+1 of the fingerprint */
205 for (j = i+1; j <= dim; j++)
206 mael(Mf,i+1,fp->per[j]) = possible(F, Ftr, V, W, fp->per, i, fp->per[j]);
207 }
208 for (i = 1; i <= dim; i++)
209 {
210 fp->diag[i] = mael(Mf,i,fp->per[i]); /* only diag(f) is needed later */
211 fp->e[i] = vecvecsmall_search(V,vecsmall_ei(dim,fp->per[i]),0);
212 if (!fp->e[i]) pari_err_BUG("qfisom, standard basis vector not found");
213 }
214 set_avma(av);
215 }
216
217 /* The Bacher polynomial for v[I] with scalar product S is defined as follows:
218 * let L be the vectors with same length as v[I] and scalar product S with v[I];
219 * for each vector w in L let nw be the number of pairs (y,z) of vectors in L,
220 * such that all scalar products between w,y and z are S, then the Bacher
221 * polynomial is the sum over the w in list of the monomials X^nw */
222 static GEN
bacher(long I,long S,struct qfauto * qf)223 bacher(long I, long S, struct qfauto *qf)
224 {
225 pari_sp av = avma;
226 GEN V=qf->V, W=qf->W, Fv=gel(qf->v,1), list, listxy, counts, vI, coef;
227 long i, j, k, nlist, nxy, sum, mind, maxd, n = lg(V)-1;
228
229 I = labs(I); /* Bacher polynomials of v[I] and -v[I] are equal */
230 vI = gel(V,I);
231 list = zero_Flv(2*n); /* vectors that have scalar product S with v[I] */
232 nlist = 0;
233 for (i = 1; i <= n; ++i)
234 if (mael(W,i,1) == mael(W,I,1))
235 {
236 long s = zv_dotproduct(vI, gel(Fv,i));
237 if (s == S) list[++nlist] = i;
238 if (-s == S) list[++nlist] = -i;
239 }
240 /* there are nlist vectors that have scalar product S with v[I] */
241 sum = nlist;
242 if (nlist==0) retmkvec2(mkvecsmall3(0,0,0),cgetg(1,t_VEC));
243 counts = cgetg(nlist+1, t_VECSMALL);
244 listxy = cgetg(nlist+1, t_VECSMALL);
245 for (i = 1; i <= nlist; ++i)
246 {
247 long S1;
248 /* listxy is the list of the nxy vectors from list that have scalar
249 product S with v[list[i]] */
250 for (j = 1; j <= nlist; ++j) listxy[j] = 0;
251 nxy = 0;
252 S1 = list[i] > 0 ? S : -S;
253 for (j = 1; j <= nlist; ++j)
254 {
255 long S2 = list[j] > 0 ? S1 : -S1;
256 /* note: for i > 0 is v[-i] = -v[i] */
257 if (zv_dotproduct(gel(V,labs(list[i])), gel(Fv,labs(list[j]))) == S2)
258 listxy[++nxy] = list[j];
259 }
260 /* counts[i] is the number of pairs for the vector v[list[i]] */
261 counts[i] = 0;
262 for (j = 1; j <= nxy; ++j)
263 {
264 long S1 = listxy[j] > 0 ? S : -S;
265 GEN Vj = gel(V, labs(listxy[j]));
266 for (k = j+1; k <= nxy; ++k)
267 {
268 long S2 = listxy[k] > 0 ? S1 : -S1, lk = labs(listxy[k]);
269 if (zv_dotproduct(Vj, gel(Fv,lk)) == S2) counts[i]++;
270 }
271 }
272 }
273 /* maxd = max degree of the Bacher-polynomial, mind = min degree */
274 mind = maxd = counts[1];
275 for (i = 2; i <= nlist; i++)
276 {
277 if (counts[i] > maxd) maxd = counts[i];
278 else if (counts[i] < mind) mind = counts[i];
279 }
280 coef = zero_Flv(maxd - mind + 1);
281 for (i = 1; i <= nlist; i++) coef[1+counts[i] - mind] += 1;
282 if (DEBUGLEVEL)
283 err_printf("QFIsom: mind=%ld maxd=%ld sum=%ld\n",mind,maxd,sum);
284 /* Bacher polynomial = sum_{mind <= i <= maxd} coef[i - mind] * X^i */
285 return gerepilecopy(av, mkvec2(mkvecsmall3(sum, mind, maxd), coef));
286 }
287
288 static GEN
init_bacher(long bachdep,struct fingerprint * fp,struct qfauto * qf)289 init_bacher(long bachdep, struct fingerprint *fp, struct qfauto *qf)
290 {
291 GEN z = cgetg(bachdep+1,t_VEC);
292 long i;
293 for (i=1;i<=bachdep;i++)
294 {
295 long bachscp = mael(qf->W,fp->e[i],1) / 2;
296 gel(z,i) = bacher(fp->e[i], bachscp, qf);
297 }
298 return z;
299 }
300
301 /* checks, whether the vector v[I] has Bacher polynomial pol */
302 static long
bachcomp(GEN pol,long I,GEN V,GEN W,GEN Fv)303 bachcomp(GEN pol, long I, GEN V, GEN W, GEN Fv)
304 {
305 pari_sp av = avma;
306 GEN co, list, listxy, vI, coef = gel(pol,2);
307 long i, j, k, nlist, nxy, count;
308 const long n = lg(V)-1, S = mael(W,I,1) / 2;
309 long sum = mael(pol,1,1), mind = mael(pol,1,2), maxd = mael(pol,1,3);
310 vI = gel(V,I);
311 list = zero_Flv(sum);
312 nlist = 0; /* nlist should be equal to pol.sum */
313 for (i = 1; i <= n && nlist <= sum; i++)
314 if (mael(W,i,1) == mael(W,I,1))
315 {
316 long s = zv_dotproduct(vI, gel(Fv,i));
317 if (s == S)
318 {
319 if (nlist < sum) list[nlist+1] = i;
320 nlist++;
321 }
322 if (-s == S)
323 {
324 if (nlist < sum) list[nlist+1] = -i;
325 nlist++;
326 }
327 }
328 /* the number of vectors with scalar product S is already different */
329 if (nlist != sum) return gc_long(av,0);
330 if (nlist == 0) return gc_long(av,1);
331 /* listxy is the list of the nxy vectors from list that have scalar product S
332 with v[list[i]] */
333 listxy = cgetg(nlist+1,t_VECSMALL);
334 co = zero_Flv(maxd - mind + 1);
335 for (i = 1; i <= nlist; i++)
336 {
337 long S1 = list[i] > 0 ? S : -S;
338 GEN Vi = gel(V, labs(list[i]));
339 for (j = 1; j <= nlist; j++) listxy[j] = 0;
340 nxy = 0;
341 for (j = 1; j <= nlist; j++)
342 {
343 long S2 = list[j] > 0 ? S1 : -S1;
344 if (zv_dotproduct(Vi, gel(Fv,labs(list[j]))) == S2)
345 listxy[++nxy] = list[j];
346 }
347 count = 0; /* number of pairs */
348 for (j = 1; j <= nxy && count <= maxd; j++)
349 {
350 long S1 = listxy[j] > 0 ? S : -S;
351 GEN Vj = gel(V, labs(listxy[j]));
352 for (k = j+1; k <= nxy && count <= maxd; k++)
353 {
354 long S2 = listxy[k] > 0 ? S1 : -S1, lk = labs(listxy[k]);
355 if (zv_dotproduct(Vj, gel(Fv,lk)) == S2) count++;
356 }
357 }
358 /* Bacher polynomials can not be equal */
359 if (count < mind || count > maxd ||
360 co[count-mind+1] >= coef[count-mind+1]) return gc_long(av, 0);
361 co[count-mind+1]++;
362 }
363 return gc_long(av, 1); /* Bacher-polynomials are equal */
364 }
365
366 static GEN
checkvecs(GEN V,GEN F,GEN norm)367 checkvecs(GEN V, GEN F, GEN norm)
368 {
369 long i, j, k, n = lg(V)-1, f = lg(F)-1;
370 GEN W = cgetg(n+1, t_MAT);
371 j = 0;
372 for (i = 1; i <= n; i++)
373 {
374 GEN normvec = cgetg(f+1, t_VECSMALL), Vi = gel(V,i);
375 for (k = 1; k <= f; ++k) normvec[k] = scp(Vi, gel(F,k), Vi);
376 if (!vecvecsmall_search(norm,normvec,0)) ++j;
377 else
378 {
379 gel(V,i-j) = Vi;
380 gel(W,i-j) = normvec;
381 }
382 }
383 setlg(V, n+1-j);
384 setlg(W, n+1-j); return W;
385 }
386
387 static long
operate(long nr,GEN A,GEN V)388 operate(long nr, GEN A, GEN V)
389 {
390 pari_sp av = avma;
391 long im,eps;
392 GEN w = zm_zc_mul(A,gel(V,labs(nr)));
393 eps = zv_canon(w);
394 if (nr < 0) eps = -eps; /* -w */
395 im = vecvecsmall_search(V,w,0);
396 if (!im) pari_err_BUG("qfauto, image of vector not found");
397 return gc_long(av, eps * im);
398 }
399
400 static GEN
orbit(GEN pt,long ipt,long npt,GEN H,GEN V)401 orbit(GEN pt, long ipt, long npt, GEN H, GEN V)
402 {
403 pari_sp av = avma;
404 long i, cnd, im, n = lg(V)-1, nH = lg(H)-1, no = npt;
405 GEN flag = zero_Flv(2*n+1)+n+1; /* need negative indices */
406 GEN orb = cgetg(2*n+1,t_VECSMALL);
407 for (i = 1; i <= npt; ++i)
408 {
409 orb[i] = pt[ipt+i];
410 flag[pt[ipt+i]] = 1;
411 }
412 for (cnd=1; cnd <= no; ++cnd)
413 for (i = 1; i <= nH; ++i)
414 {
415 im = operate(orb[cnd], gel(H,i), V);
416 /* image is a new point in the orbit */
417 if (flag[im] == 0) { orb[++no] = im; flag[im] = 1; }
418 }
419 setlg(orb,no+1); return gerepileuptoleaf(av, orb);
420 }
421
422 /* return the length of the orbit of pt under the first nG matrices in G */
423 static long
orbitlen(long pt,long orblen,GEN G,long nG,GEN V)424 orbitlen(long pt, long orblen, GEN G, long nG, GEN V)
425 {
426 pari_sp av = avma;
427 long i, len, cnd, n = lg(V)-1;
428 GEN orb, flag;
429 /* if flag[i + n+1] = 1, -n <= i <= n, then i is already in the orbit */
430 flag = zero_F2v(2*n + 1);
431 orb = zero_Flv(orblen); orb[1] = pt;
432 F2v_set(flag,pt+n+1);
433 len = 1;
434 for (cnd = 1; cnd <= len && len < orblen; cnd++)
435 for (i = 1; i <= nG && len < orblen; i++)
436 {
437 long im = operate(orb[cnd], gel(G,i), V);
438 /* image is a new point in the orbit */
439 if (!F2v_coeff(flag,im+n+1)) { orb[++len] = im; F2v_set(flag,im+n+1); }
440 }
441 return gc_long(av, len);
442 }
443
444 /* delete the elements in orb2 from orb1, an entry 0 marks the end of the
445 * list, returns the length of orb1 */
446 static long
orbdelete(GEN orb1,GEN orb2)447 orbdelete(GEN orb1, GEN orb2)
448 {
449 long i, j, len, l1 = lg(orb1)-1, l2 = lg(orb2)-1;
450 for (i = 1; i <= l1 && orb1[i] != 0; ++i);
451 len = i - 1;
452 for (i = 1; i <= l2 && orb2[i] != 0; ++i)
453 {
454 long o2i = orb2[i];
455 for (j = 1; j <= len && orb1[j] != o2i; ++j);
456 /* orb1[j] = orb2[i], hence delete orb1[j] from orb1 */
457 if (j <= len) { orb1[j] = orb1[len]; orb1[len--] = 0; }
458 }
459 return len;
460 }
461
462 static long
orbsubtract(GEN Cs,GEN pt,long ipt,long npt,GEN H,GEN V,long * len)463 orbsubtract(GEN Cs, GEN pt, long ipt, long npt, GEN H, GEN V, long *len)
464 {
465 pari_sp av = avma;
466 GEN orb = orbit(pt, ipt, npt, H, V);
467 if (len) *len = lg(orb)-1;
468 return gc_long(av, orbdelete(Cs, orb));
469 }
470
471 /* Generates the matrix X whose per[i]-th row is the vector V[x[i]] */
472 static GEN
matgen(GEN x,GEN per,GEN V)473 matgen(GEN x, GEN per, GEN V)
474 {
475 long i, j, n = lg(x)-1;
476 GEN X = cgetg(n+1,t_MAT);
477 for (i = 1; i <= n; i++)
478 {
479 GEN Xp = cgetg(n+1,t_VECSMALL);
480 long xi = x[i];
481 for (j = 1; j <= n; j++) Xp[j] = xi > 0? mael(V,xi,j): -mael(V,-xi,j);
482 gel(X, per[i]) = Xp;
483 }
484 return X;
485 }
486 /* x1 corresponds to an element X1 mapping some vector e on p1, x2 to an
487 * element X2 mapping e on p2 and G is a generator mapping p1 on p2, then
488 * S = X1*G*X2^-1 stabilizes e */
489 static GEN
stabil(GEN x1,GEN x2,GEN per,GEN G,GEN V,ulong p)490 stabil(GEN x1, GEN x2, GEN per, GEN G, GEN V, ulong p)
491 {
492 pari_sp av = avma;
493 long i, n = lg(x1)-1;
494 GEN XG, X2, x = cgetg(n+1,t_VECSMALL);
495 for (i = 1; i <= n; i++) x[i] = operate(x1[i], G, V);
496 /* XG is the composite mapping of the matrix corresponding to x1 and G */
497 XG = matgen(x, per, V);
498 X2 = matgen(x2, per, V);
499 return gerepileupto(av, zm_divmod(X2,XG,p));
500 }
501
502 /* computes the orbit of fp.e[I] under the generators in G->g[I]...G->g[n-1]
503 * and elements stabilizing fp.e[I], has some heuristic break conditions, the
504 * generators in G->g[i] stabilize fp.e[0]...fp.e[i-1] but not fp.e[i],
505 * G->ng[i] is the number of generators in G->g[i], the first G->nsg[i] of
506 * which are elements which are obtained as stabilizer elements in
507 * <G->g[0],...,G->g[i-1]>, G->ord[i] is the orbit length of fp.e[i] under
508 * <G->g[i],...,G->g[n-1]> */
509 static void
stab(long I,struct group * G,struct fingerprint * fp,GEN V,ulong p)510 stab(long I, struct group *G, struct fingerprint *fp, GEN V, ulong p)
511 {
512 GEN orb, w, flag, H, Hj, S;
513 long len, cnd, i, j, k, l, im, nH, fail;
514 long Maxfail, Rest, dim = lg(fp->diag)-1, n = lg(V)-1;
515 /* Heuristic break conditions for the computation of stabilizer elements:
516 * it is too expensive to calculate all the stabilizer generators, which are
517 * obtained from the orbit, since this is highly redundant. On the other hand
518 * every new generator which enlarges the group is cheaper than one obtained
519 * from the backtrack, after Maxfail subsequent stabilizer elements, that do
520 * not enlarge the group, Rest more elements are calculated even if they
521 * leave the group unchanged, since it turned out that this is often useful
522 * in the following steps. Increasing the parameters may decrease the number
523 * of generators for the group while increasing the running time. */
524 for (Rest = 0, i = I; i <= dim; i++)
525 if (fp->diag[i] > 1 && G->ord[i] < fp->diag[i]) ++Rest;
526 for (Maxfail = Rest, i = 1; i <= dim; i++)
527 if (fp->diag[i] > 1) ++Maxfail;
528 for (nH = 0, i = I; i <= dim; i++)
529 nH += G->ng[i];
530 /* generators of the group in which the stabilizer is computed */
531 H = cgetg(nH+1,t_MAT);
532 Hj = cgetg(nH+2,t_MAT);
533 for (i = I, k = 0; i <= dim; i++)
534 for (j = 1; j <= G->ng[i]; j++) gel(H,++k) = gmael(G->g,i,j);
535 /* in w[V.n+i] an element is stored that maps fp.e[I] on v[i] */
536 w = cgetg(2*n+2,t_VEC);
537 orb = zero_Flv(2*n); /* the orbit of fp.e[I] */
538 flag = zero_F2v(2*n+1); /* if flag[i + V.n], then i is already in orbit */
539 orb[1] = fp->e[I];
540 F2v_set(flag,orb[1]+n+1);
541 gel(w,orb[1]+n+1) = cgetg(dim+1,t_VECSMALL);
542 for (i = 1; i <= dim; i++) mael(w,orb[1]+n+1,i) = fp->e[i];
543 cnd = len = 1;
544 fail = 0; /* number of successive failures */
545 while (cnd <= len && fail < Maxfail+Rest)
546 {
547 for (i = 1; i <= nH && fail < Maxfail+Rest; ++i)
548 {
549 if (fail >= Maxfail)
550 /* already Maxfail successive failures: apply a random generator to a
551 * random point of the orbit to get Rest more stabilizer elements */
552 {
553 cnd = 1+(long)random_Fl(len);
554 i = 1+(long)random_Fl(nH);
555 }
556 im = operate(orb[cnd], gel(H,i), V);
557 if (F2v_coeff(flag,im+n+1) == 0)
558 { /* found new element, appended to the orbit; an element mapping
559 * fp.e[I] to im is stored in w[im+V.n] */
560 GEN wim;
561 orb[++len] = im;
562 F2v_set(flag,im+n+1);
563 wim = cgetg(dim+1,t_VECSMALL);
564 gel(w,im+n+1) = wim;
565 for (j = 1; j <= dim; ++j)
566 wim[j] = operate(mael(w,orb[cnd]+n+1,j), gel(H,i), V);
567 }
568 else /* image was already in the orbit */
569 { /* j = first index where images of old and new element mapping e[I] on
570 * im differ */
571 for (j = I; j <= dim; j++)
572 if (operate(mael(w,orb[cnd]+n+1,j), gel(H,i), V) != mael(w,im+n+1,j))
573 break;
574 if (j <= dim && (G->ord[j] < fp->diag[j] || fail >= Maxfail))
575 {
576 GEN wo = gel(w,orb[cnd]+n+1);
577 long tmplen, nHj = 1;
578 /* new stabilizer element S = w[orb[cnd]+V.n] * H[i] * (w[im+V.n])^-1 */
579 S = stabil(wo, gel(w,im+n+1), fp->per, gel(H,i), V, p);
580 gel(Hj,1) = S;
581 for (k = j; k <= dim; k++)
582 for (l = 1; l <= G->ng[k]; l++) gel(Hj, ++nHj) = gmael(G->g,k,l);
583 tmplen = orbitlen(fp->e[j], fp->diag[j], Hj, nHj, V);
584 if (tmplen > G->ord[j] || fail >= Maxfail)
585 /* the new stabilizer element S either enlarges the orbit of e[j]
586 or it is one of the additional elements after MAXFAIL failures */
587 {
588 GEN Ggj;
589 G->ord[j] = tmplen;
590 G->ng[j]++;
591 G->nsg[j]++;
592 /* allocate memory for new generator */
593 gel(G->g,j) = vec_lengthen(gel(G->g,j),G->ng[j]);
594 Ggj = gel(G->g,j);
595 /* new generator is inserted as stabilizer element nsg[j]-1 */
596 for (k = G->ng[j]; k > G->nsg[j]; k--) gel(Ggj,k) = gel(Ggj,k-1);
597 gel(Ggj,G->nsg[j]) = S;
598 nH++;
599 H = vec_lengthen(H, nH);
600 Hj = vec_lengthen(Hj, nH+1);
601 gel(H,nH) = gel(Ggj,G->nsg[j]); /* append new generator to H */
602 if (fail < Maxfail)
603 fail = 0; /* number of failures is reset to 0 */
604 else
605 ++fail;
606 }
607 else /* new stabilizer element S does not enlarge the orbit of e[j] */
608 ++fail;
609 }
610 else if ((j < dim && fail < Maxfail) || (j == dim && fail >= Maxfail))
611 ++fail;
612 /* if S is the identity and fail < Maxfail, nothing is done */
613 }
614 }
615 if (fail < Maxfail) ++cnd;
616 }
617 }
618
619 /* tests, whether x[1],...,x[I-1] is a partial automorphism, using scalar
620 * product combinations and Bacher-polynomials depending on the chosen options,
621 * puts the candidates for x[I] (i.e. the vectors vec such that the scalar
622 * products of x[1],...,x[I-1],vec are correct) on CI, returns their number
623 * (should be fp.diag[I]) */
624 static long
qfisom_candidates_novec(GEN CI,long I,GEN x,struct qfauto * qf,struct qfauto * qff,struct fingerprint * fp)625 qfisom_candidates_novec(GEN CI, long I, GEN x, struct qfauto *qf,
626 struct qfauto *qff, struct fingerprint *fp)
627 {
628 pari_sp av = avma;
629 long i, j, k, okp, okm, nr, fail;
630 GEN vec, F =qf->F, V=qff->V, W=qff->W, v=qff->v;
631 long n = lg(V)-1, f = lg(F)-1;
632 vec = cgetg(I,t_VECSMALL);
633 for (i = 1; i <= fp->diag[I]; i++) CI[i] = 0; /* list for the candidates */
634 nr = fail = 0;
635 for (j = 1; j <= n && fail == 0; j++)
636 {
637 GEN Vj = gel(V,j), Wj = gel(W, j);
638 okp = okm = 0;
639 for (i = 1; i <= f; i++)
640 {
641 GEN FAiI = gmael(F,i,fp->per[I]), FFvi = gel(v,i);
642 /* vec is the vector of scalar products of V.v[j] with the first I base
643 vectors x[0]...x[I-1] */
644 for (k = 1; k < I; k++)
645 {
646 long xk = x[k];
647 if (xk > 0)
648 vec[k] = zv_dotproduct(Vj, gel(FFvi,xk));
649 else
650 vec[k] = -zv_dotproduct(Vj, gel(FFvi,-xk));
651 }
652 for (k = 1; k < I && vec[k] == FAiI[fp->per[k]]; k++);
653 if (k == I && Wj[i] == FAiI[fp->per[I]]) ++okp;
654 /* V.v[j] is a candidate for x[I] with respect to the form F.A[i] */
655 for (k = 1; k < I && vec[k] == -FAiI[fp->per[k]]; k++);
656 if (k == I && Wj[i] == FAiI[fp->per[I]]) ++okm;
657 /* -V.v[j] is a candidate for x[I] with respect to the form F.A[i] */
658 }
659 if (okp == f) /* V.v[j] is a candidate for x[I] */
660 {
661 if (nr < fp->diag[I])
662 CI[++nr] = j;
663 else
664 fail = 1; /* too many candidates */
665 }
666 if (okm == f) /* -V.v[j] is a candidate for x[I] */
667 {
668 if (nr < fp->diag[I])
669 CI[++nr] = -j;
670 else
671 fail = 1; /* too many candidates */
672 }
673 }
674 return gc_long(av, fail? 0: nr);
675 }
676
677 static long
qfisom_candidates(GEN CI,long I,GEN x,struct qfauto * qf,struct qfauto * qff,struct fingerprint * fp,struct qfcand * qfcand)678 qfisom_candidates(GEN CI, long I, GEN x, struct qfauto *qf,
679 struct qfauto *qff, struct fingerprint *fp, struct qfcand *qfcand)
680 {
681 pari_sp av = avma;
682 GEN vec, xvec, xbase, Fxbase, scpvec, com, list, trans, ccoef, cF;
683 GEN F =qf->F, V=qff->V, W=qff->W, v=qff->v, FF= qff->F;
684 long i, j, k, okp, okm, nr, nc, vj, rank, num;
685 long dim = qf->dim, n = lg(V)-1, f = lg(F)-1;
686 long DEP = qfcand->cdep, len = f * DEP;
687 if (I >= 2 && I <= lg(qfcand->bacher_pol))
688 {
689 long t = labs(x[I-1]);
690 GEN bpolI = gel(qfcand->bacher_pol,I-1);
691 if (bachcomp(bpolI, t, V, W, gel(v,1)) == 0) return 0;
692 }
693 if (I==1 || DEP ==0) return qfisom_candidates_novec(CI,I,x,qf,qff,fp);
694 vec = cgetg(I,t_VECSMALL);
695 scpvec = cgetg(len+1,t_VECSMALL);
696 com = gel(qfcand->comb,I-1);
697 list=gel(com,1); trans = gel(com,2); ccoef = gel(com,3); cF=gel(com,4);
698 rank = lg(trans)-1;
699 nc = lg(list)-2;
700 /* xvec is the list of vector sums which are computed with respect to the
701 partial basis in x */
702 xvec = zero_Flm_copy(dim,nc+1);
703 /* xbase should be a basis for the lattice generated by the vectors in xvec,
704 it is obtained via the transformation matrix comb[I-1].trans */
705 xbase = cgetg(rank+1,t_MAT);
706 for (i = 1; i <= rank; ++i)
707 gel(xbase,i) = cgetg(dim+1,t_VECSMALL);
708 Fxbase = cgetg(rank+1,t_MAT); /* product of a form F with the base xbase */
709 for (i = 1; i <= rank; ++i) gel(Fxbase,i) = cgetg(dim+1,t_VECSMALL);
710 for (i = 1; i <= fp->diag[I]; ++i) CI[i] = 0; /* list for the candidates */
711 nr = 0;
712 for (j = 1; j <= n; ++j)
713 {
714 long sign;
715 GEN Vj = gel(V,j), Wj = gel(W, j);
716 okp = okm = 0;
717 for (k = 1; k <= len; ++k) scpvec[k] = 0;
718 for (i = 1; i <= f; ++i)
719 {
720 GEN FAiI = gmael(F,i,fp->per[I]), FFvi = gel(v,i);
721 /* vec is the vector of scalar products of V.v[j] with the first I base
722 vectors x[0]...x[I-1] */
723 for (k = 1; k < I; ++k)
724 {
725 long xk = x[k];
726 if (xk > 0)
727 vec[k] = zv_dotproduct(Vj, gel(FFvi,xk));
728 else
729 vec[k] = -zv_dotproduct(Vj, gel(FFvi,-xk));
730 }
731 for (k = 1; k < I && vec[k] == FAiI[fp->per[k]]; ++k);
732 if (k == I && Wj[i] == FAiI[fp->per[I]]) ++okp;
733 /* V.v[j] is a candidate for x[I] with respect to the form F.A[i] */
734 for (k = 1; k < I && vec[k] == -FAiI[fp->per[k]]; ++k);
735 if (k == I && Wj[i] == FAiI[fp->per[I]]) ++okm;
736 /* -V.v[j] is a candidate for x[I] with respect to the form F.A[i] */
737 for (k = I-1; k >= 1 && k > I-1-DEP; --k)
738 scpvec[(i-1)*DEP+I-k] = vec[k];
739 }
740 /* check, whether the scalar product combination scpvec is contained in the
741 list comb[I-1].list */
742 if (!zv_equal0(scpvec))
743 {
744 sign = zv_canon(scpvec);
745 num = vecvecsmall_search(list,scpvec,0);
746 if (!num) return gc_long(av,0); /* x[0..I-1] not a partial automorphism */
747 else
748 {
749 GEN xnum = gel(xvec,num);
750 for (k = 1; k <= dim; ++k) xnum[k] += sign * Vj[k];
751 }
752 }
753 if (okp == f) /* V.v[j] is a candidate for x[I] */
754 {
755 if (nr >= fp->diag[I]) return gc_long(av,0); /* too many candidates */
756 CI[++nr] = j;
757 }
758 if (okm == f) /* -V.v[j] is a candidate for x[I] */
759 {
760 if (nr >= fp->diag[I]) return gc_long(av,0); /* too many candidates */
761 CI[++nr] = -j;
762 }
763 }
764 if (nr == fp->diag[I])
765 { /* compute the basis of the lattice generated by the vectors in xvec via
766 the transformation matrix comb[I-1].trans */
767 for (i = 1; i <= rank; ++i)
768 {
769 GEN comtri = gel(trans,i);
770 for (j = 1; j <= dim; ++j)
771 {
772 long xbij = 0;
773 for (k = 1; k <= nc+1; ++k) xbij += comtri[k] * mael(xvec,k,j);
774 mael(xbase,i,j) = xbij;
775 }
776 }
777 /* check, whether the base xbase has the right scalar products */
778 for (i = 1; i <= f; ++i)
779 {
780 for (j = 1; j <= rank; ++j)
781 for (k = 1; k <= dim; ++k)
782 mael(Fxbase,j,k) = zv_dotproduct(gmael(FF,i,k), gel(xbase,j));
783 for (j = 1; j <= rank; ++j)
784 for (k = 1; k <= j; ++k) /* a scalar product is wrong ? */
785 if (zv_dotproduct(gel(xbase,j), gel(Fxbase,k)) != mael3(cF,i,j,k))
786 return gc_long(av, 0);
787 }
788 for (i = 1; i <= nc+1; ++i)
789 {
790 GEN comcoi = gel(ccoef,i);
791 for (j = 1; j <= dim; ++j)
792 {
793 vj = 0;
794 for (k = 1; k <= rank; ++k)
795 vj += comcoi[k] * mael(xbase,k,j);
796 if (vj != mael(xvec,i,j)) return gc_long(av,0); /* an entry is wrong */
797 }
798 }
799 }
800 return gc_long(av, nr);
801 }
802
803 static long
aut(long step,GEN x,GEN C,struct group * G,struct qfauto * qf,struct fingerprint * fp,struct qfcand * cand)804 aut(long step, GEN x, GEN C, struct group *G, struct qfauto *qf,
805 struct fingerprint *fp, struct qfcand *cand)
806 {
807 long dim = qf->dim;
808 GEN orb;
809 /* found new automorphism */
810 if (step == dim && mael(C,dim,1)) { x[dim] = mael(C,dim,1); return 1; }
811 orb = cgetg(2,t_VECSMALL);
812 while (mael(C,step,1))
813 {
814 long nbc;
815 /* choose the image of the base-vector nr. step */
816 x[step] = mael(C,step,1);
817 /* check, whether x[0..step] is a partial automorphism and compute
818 the candidates for x[step+1] */
819 nbc = qfisom_candidates(gel(C,step+1), step+1, x, qf, qf, fp, cand);
820 if (nbc == fp->diag[step+1]
821 && aut(step+1, x, C, G, qf, fp, cand)) return 1;
822 orb[1] = x[step];
823 /* delete chosen vector from list of candidates */
824 (void)orbdelete(gel(C,step), orb);
825 }
826 return 0;
827 }
828
829 /* search new automorphisms until on all levels representatives for all orbits
830 * have been tested */
831 static void
autom(struct group * G,struct qfauto * qf,struct fingerprint * fp,struct qfcand * cand)832 autom(struct group *G, struct qfauto *qf, struct fingerprint *fp,
833 struct qfcand *cand)
834 {
835 long i, j, step, im, nC, found, tries, nbad;
836 GEN x, bad, H, V = qf->V;
837 long dim = qf->dim, n = lg(V)-1, STAB = 1;
838 GEN C = cgetg(dim+1,t_VEC);
839 /* C[i] is the list of candidates for the image of the i-th base-vector */
840 for (i = 1; i <= dim; ++i)
841 gel(C,i) = cgetg(fp->diag[i]+1, t_VECSMALL);
842 /* x is the new base i.e. x[i] is the index in V.v of the i-th base-vector */
843 x = cgetg(dim+1, t_VECSMALL);
844 for (step = STAB; step <= dim; ++step)
845 {
846 long nH = 0;
847 if (DEBUGLEVEL) err_printf("QFAuto: Step %ld/%ld\n",step,dim);
848 for (nH = 0, i = step; i <= dim; ++i)
849 nH += G->ng[i];
850 H = cgetg(nH+1,t_VEC);
851 for (nH = 0, i = step; i <= dim; ++i)
852 for (j = 1; j <= G->ng[i]; ++j)
853 gel(H,++nH) = gmael(G->g,i,j);
854 bad = zero_Flv(2*n);
855 nbad = 0;
856 /* the first step base-vectors are fixed */
857 for (i = 1; i < step; ++i)
858 x[i] = fp->e[i];
859 /* compute the candidates for x[step] */
860 if (fp->diag[step] > 1)
861 /* if fp.diag[step] > 1 compute the candidates for x[step] */
862 nC = qfisom_candidates(gel(C,step), step, x, qf, qf, fp, cand);
863 else
864 /* if fp.diag[step] == 1, fp.e[step] is the only candidate */
865 {
866 mael(C,step,1) = fp->e[step];
867 nC = 1;
868 }
869 /* delete the orbit of the step-th base-vector from the candidates */
870 nC = orbsubtract(gel(C,step), fp->e, step-1, 1, H, V, &(G->ord[step]));
871 while (nC > 0 && (im = mael(C,step,1)) != 0)
872 {
873 found = 0;
874 /* tries vector V.v[im] as image of the step-th base-vector */
875 x[step] = im;
876 if (step < dim)
877 {
878 long nbc;
879 /* check, whether x[0]...x[step] is a partial basis and compute the
880 candidates for x[step+1] */
881 nbc = qfisom_candidates(gel(C,step+1), step+1, x, qf, qf, fp, cand);
882 if (nbc == fp->diag[step+1])
883 /* go into the recursion */
884 found = aut(step+1, x, C, G, qf, fp, cand);
885 else
886 found = 0;
887 }
888 else
889 found = 1;
890 if (!found) /* x[0..step] can not be continued to an automorphism */
891 { /* delete the orbit of im from the candidates for x[step] */
892 nC = orbsubtract(gel(C,step),mkvecsmall(im), 0, 1, H, V, NULL);
893 bad[++nbad] = im;
894 }
895 else
896 { /* new generator has been found */
897 GEN Gstep;
898 ++G->ng[step];
899 /* append new generator to G->g[step] */
900 Gstep = vec_lengthen(gel(G->g,step),G->ng[step]);
901 gel(Gstep,G->ng[step]) = matgen(x, fp->per, V);
902 gel(G->g,step) = Gstep;
903 ++nH;
904 H = cgetg(nH+1, t_VEC);
905 for (nH = 0, i = step; i <= dim; i++)
906 for (j = 1; j <= G->ng[i]; j++) gel(H,++nH) = gmael(G->g,i,j);
907 nC = orbsubtract(gel(C,step), fp->e, step-1, 1, H, V, &(G->ord[step]));
908 nC = orbsubtract(gel(C,step), bad, 0, nbad, H, V, NULL);
909 }
910 }
911 /* test, whether on step STAB some generators may be omitted */
912 if (step == STAB) for (tries = G->nsg[step]; tries <= G->ng[step]; tries++)
913 {
914 nH = 0;
915 for (j = 1; j < tries; j++) gel(H,++nH) = gmael(G->g,step,j);
916 for (j = tries+1; j < G->ng[step]; j++) gel(H,++nH) = gmael(G->g,step,j);
917 for (i = step+1; i <= dim; i++)
918 for (j = 1; j <= G->ng[i]; j++) gel(H,++nH) = gmael(G->g,i,j);
919 if (orbitlen(fp->e[step], G->ord[step], H, nH, V) == G->ord[step])
920 { /* generator g[step][tries] can be omitted */
921 G->ng[step]--;
922 for (i = tries; i < G->ng[step]; ++i)
923 gmael(G->g,step,i) = gmael(G->g,step,i+1);
924 tries--;
925 }
926 }
927 /* calculate stabilizer elements fixing basis-vectors fp.e[0..step] */
928 if (step < dim && G->ord[step] > 1) stab(step, G, fp, V, qf->p);
929 }
930 }
931
932 #define MAXENTRY (1L<<((BITS_IN_LONG-2)>>1))
933 #define MAXNORM (1L<<(BITS_IN_LONG-2))
934
935 static long
zm_maxdiag(GEN A)936 zm_maxdiag(GEN A)
937 {
938 long dim = lg(A)-1, max = coeff(A,1,1), i;
939 for (i = 2; i <= dim; ++i)
940 if (coeff(A,i,i) > max) max = coeff(A,i,i);
941 return max;
942 }
943
944 static GEN
init_qfauto(GEN F,GEN U,long max,struct qfauto * qf,GEN norm,GEN minvec)945 init_qfauto(GEN F, GEN U, long max, struct qfauto *qf, GEN norm, GEN minvec)
946 {
947 GEN V = minvec ? ZM_to_zm_canon(minvec)
948 : gel(minim_zm(zm_to_ZM(gel(F,1)), stoi(max), NULL), 3);
949 GEN W, v;
950 long i, j, k, n = lg(V)-1, f = lg(F)-1, dim = lg(gel(F,1))-1;
951 for (i = 1; i <= n; i++)
952 {
953 GEN Vi = gel(V,i);
954 for (k = 1; k <= dim; k++)
955 {
956 long l = labs(Vi[k]);
957 if (l > max) max = l;
958 }
959 }
960 if (max > MAXENTRY) pari_err_OVERFLOW("qfisom [lattice too large]");
961 qf->p = unextprime(2*max+1);
962 V = vecvecsmall_sort_uniq(V);
963 if (!norm)
964 {
965 norm = cgetg(dim+1,t_VEC);
966 for (i = 1; i <= dim; i++)
967 {
968 GEN Ni = cgetg(f+1,t_VECSMALL);
969 for (k = 1; k <= f; k++) Ni[k] = mael3(F,k,i,i);
970 gel(norm,i) = Ni;
971 }
972 norm = vecvecsmall_sort_uniq(norm);
973 }
974 W = checkvecs(V, F, norm);
975 v = cgetg(f+1,t_VEC);
976 /* the product of the maximal entry in the short vectors with the maximal
977 entry in v[i] should not exceed MAXNORM to avoid overflow */
978 max = MAXNORM / max;
979 for (i = 1; i <= f; ++i)
980 {
981 GEN Fi = gel(F,i), vi = cgetg(n+1,t_MAT);
982 gel(v,i) = vi;
983 for (j = 1; j <= n; ++j)
984 {
985 GEN Vj = gel(V,j), vij = cgetg(dim+1, t_VECSMALL);
986 gel(vi,j) = vij;
987 for (k = 1; k <= dim; ++k)
988 {
989 vij[k] = zv_dotproduct(gel(Fi,k), Vj);
990 if (labs(vij[k]) > max) pari_err_OVERFLOW("qfisom [lattice too large]");
991 }
992 }
993 }
994 qf->dim = dim; qf->F = F; qf->V = V; qf->W = W; qf->v = v; qf->U = U;
995 return norm;
996 }
997
998 static void
init_qfgroup(struct group * G,struct fingerprint * fp,struct qfauto * qf)999 init_qfgroup(struct group *G, struct fingerprint *fp, struct qfauto *qf)
1000 {
1001 GEN H, M, V = qf->V;
1002 long nH, i, j, k, dim = qf->dim;
1003 G->ng = zero_Flv(dim+1);
1004 G->nsg = zero_Flv(dim+1);
1005 G->ord = cgetg(dim+1,t_VECSMALL);
1006 G->g = cgetg(dim+1,t_VEC);
1007 for (i = 1; i <= dim; ++i) gel(G->g,i) = mkvec(gen_0);
1008 M = matid_Flm(dim);
1009 gmael(G->g,1,1) = M;
1010 G->ng[1] = 1;
1011 /* -Id is always an automorphism */
1012 for (i = 1; i <= dim; i++) mael(M,i,i) = -1;
1013 nH = 0;
1014 for (i = 1; i <= dim; i++) nH += G->ng[i];
1015 H = cgetg(nH+1,t_MAT);
1016 /* calculate the orbit lengths under the automorphisms known so far */
1017 for (i = 1; i <= dim; ++i)
1018 {
1019 if (G->ng[i] > 0)
1020 {
1021 nH = 0;
1022 for (j = i; j <= dim; j++)
1023 for (k = 1; k <= G->ng[j]; k++) gel(H,++nH) = gmael(G->g,j,k);
1024 G->ord[i] = orbitlen(fp->e[i], fp->diag[i], H, nH, V);
1025 }
1026 else
1027 G->ord[i] = 1;
1028 }
1029 }
1030
1031 /* calculates the scalar products of the vector w with the base vectors
1032 * v[b[I]] down to v[b[I-dep+1]] with respect to all invariant forms and puts
1033 * them on scpvec */
1034 static GEN
scpvector(GEN w,GEN b,long I,long dep,GEN v)1035 scpvector(GEN w, GEN b, long I, long dep, GEN v)
1036 {
1037 long i, j, n = lg(v)-1;
1038 GEN scpvec = zero_Flv(dep*n);
1039 for (i = I; i >= 1 && i > I-dep; i--)
1040 {
1041 long bi = b[i];
1042 if (bi > 0)
1043 for (j = 1; j <= n; j++)
1044 scpvec[1+(j-1)*dep + I-i] = zv_dotproduct(w, gmael(v,j,bi));
1045 else
1046 for (j = 1; j <= n; j++)
1047 scpvec[1+(j-1)*dep + I-i] = -zv_dotproduct(w, gmael(v,j,-bi));
1048 }
1049 return scpvec;
1050 }
1051
1052 /* computes the list of scalar product combinations of the vectors
1053 * in V.v with the basis-vectors in b */
1054 static GEN
scpvecs(GEN * pt_vec,long I,GEN b,long dep,struct qfauto * qf)1055 scpvecs(GEN *pt_vec, long I, GEN b, long dep, struct qfauto *qf)
1056 {
1057 GEN list, vec, V = qf->V, F = qf->F, v = qf->v;
1058 long n = lg(V)-1, dim = lg(gel(F,1))-1, len = (lg(F)-1)*dep;
1059 long j;
1060 /* the first vector in the list is the 0-vector and is not counted */
1061 list = mkvec(zero_Flv(len));
1062 vec = mkvec(zero_Flv(dim));
1063 for (j = 1; j <= n; ++j)
1064 {
1065 GEN Vvj = gel(V,j), scpvec = scpvector(Vvj, b, I, dep, v);
1066 long i, nr, sign;
1067 if (zv_equal0(scpvec)) continue;
1068 sign = zv_canon(scpvec);
1069 nr = vecvecsmall_search(list, scpvec, 0);
1070 if (nr > 0) /* scpvec already in list */
1071 {
1072 GEN vecnr = gel(vec,nr);
1073 for (i = 1; i <= dim; i++) vecnr[i] += sign * Vvj[i];
1074 }
1075 else /* scpvec is a new scalar product combination */
1076 {
1077 nr = vecvecsmall_search(list, scpvec, 1);
1078 list = vec_insert(list,nr,scpvec);
1079 vec = vec_insert(vec,nr,sign < 0 ? zv_neg(Vvj) : zv_copy(Vvj));
1080 }
1081 }
1082 settyp(list,t_MAT);
1083 settyp(vec,t_MAT);
1084 *pt_vec = vec; return list;
1085 }
1086
1087 /* com->F[i] is the Gram-matrix of the basis b with respect to F.A[i] */
1088 static GEN
scpforms(GEN b,struct qfauto * qf)1089 scpforms(GEN b, struct qfauto *qf)
1090 {
1091 GEN F = qf->F;
1092 long i, j, k, n = lg(F)-1, dim = lg(gel(F,1))-1, nb = lg(b)-1;
1093 GEN gram = cgetg(n+1, t_VEC), Fbi = cgetg(nb+1, t_MAT);
1094 /* list of products of F.A[i] with the vectors in b */
1095 for (j = 1; j <= nb; j++) gel(Fbi, j) = cgetg(dim+1, t_VECSMALL);
1096 for (i = 1; i <= n; i++)
1097 {
1098 GEN FAi = gel(F,i);
1099 gel(gram, i) = cgetg(nb+1, t_MAT);
1100 for (j = 1; j <= nb; j++)
1101 for (k = 1; k <= dim; k++)
1102 mael(Fbi,j,k) = zv_dotproduct(gel(FAi,k), gel(b,j));
1103 for (j = 1; j <= nb; j++)
1104 {
1105 GEN Gij = cgetg(nb+1, t_VECSMALL);
1106 for (k = 1; k <= nb; k++) Gij[k] = zv_dotproduct(gel(b,j), gel(Fbi,k));
1107 gmael(gram,i,j) = Gij;
1108 }
1109 }
1110 return gram;
1111 }
1112
1113 static GEN
gen_comb(long cdep,GEN A,GEN e,struct qfauto * qf,long lim)1114 gen_comb(long cdep, GEN A, GEN e, struct qfauto *qf, long lim)
1115 {
1116 long i, dim = lg(A)-1;
1117 GEN comb = cgetg(dim+1,t_VEC);
1118 for (i = 1; i <= dim; ++i)
1119 {
1120 pari_sp av = avma;
1121 GEN trans, ccoef, cF, B, BI, sumveclist, sumvecbase;
1122 GEN list = scpvecs(&sumveclist, i, e, cdep, qf);
1123 GEN M = zm_to_ZM(sumveclist);
1124 GEN T = lllgramint(qf_apply_ZM(A,M));
1125 if (lim && lg(T)-1>=lim) return NULL;
1126 B = ZM_mul(M,T);
1127 BI = RgM_inv(B);
1128 sumvecbase = ZM_trunc_to_zm(B);
1129 trans = ZM_trunc_to_zm(T);
1130 ccoef = ZM_trunc_to_zm(RgM_mul(BI,M));
1131 cF = scpforms(sumvecbase, qf);
1132 gel(comb,i) = gerepilecopy(av, mkvec4(list, trans, ccoef, cF));
1133 }
1134 return comb;
1135 }
1136
1137 static void
init_comb(struct qfcand * cand,GEN A,GEN e,struct qfauto * qf)1138 init_comb(struct qfcand *cand, GEN A, GEN e, struct qfauto *qf)
1139 {
1140 long dim = lg(A)-1;
1141 GEN Am = zm_to_ZM(A);
1142 for (cand->cdep = 1; ; cand->cdep++)
1143 {
1144 cand->comb = gen_comb(cand->cdep, Am, e, qf, (dim+1)>>1);
1145 if (!cand->comb) break;
1146 }
1147 cand->cdep= maxss(1, cand->cdep-1);
1148 cand->comb = gen_comb(cand->cdep, Am, e, qf, 0);
1149 }
1150
1151 static void
init_flags(struct qfcand * cand,GEN A,struct fingerprint * fp,struct qfauto * qf,GEN flags)1152 init_flags(struct qfcand *cand, GEN A, struct fingerprint *fp,
1153 struct qfauto *qf, GEN flags)
1154 {
1155 if (!flags)
1156 {
1157 init_comb(cand, A, fp->e, qf);
1158 cand->bacher_pol = init_bacher(0, fp, qf);
1159 }
1160 else
1161 {
1162 long cdep, bach;
1163 if (typ(flags)!=t_VEC || lg(flags)!=3) pari_err_TYPE("qfisominit",flags);
1164 cdep = gtos(gel(flags,1));
1165 bach = minss(gtos(gel(flags,2)), lg(fp->e)-1);
1166 if (cdep<0 || bach<0) pari_err_FLAG("qfisom");
1167 cand->cdep = cdep;
1168 cand->comb = cdep? gen_comb(cdep, zm_to_ZM(A), fp->e, qf, 0): NULL;
1169 cand->bacher_pol = init_bacher(bach, fp, qf);
1170 }
1171 }
1172
1173 static GEN
gen_group(struct group * G,GEN U)1174 gen_group(struct group *G, GEN U)
1175 {
1176 GEN V, o = gen_1;
1177 long i, j, n, dim = lg(G->ord)-1;
1178 for (i = 1; i <= dim; i++) o = muliu(o, G->ord[i]);
1179 for (i = n = 1; i <= dim; i++) n += G->ng[i] - G->nsg[i];
1180 V = cgetg(n, t_VEC);
1181 for (i = n = 1; i <= dim; ++i)
1182 for (j=G->nsg[i]+1; j<=G->ng[i]; j++)
1183 gel(V,n++) = U ? zm_mul(gel(U,1), zm_mul(gmael(G->g,i,j), gel(U,2)))
1184 : gmael(G->g,i,j);
1185 return mkvec2(o, V);
1186 }
1187
1188 static long
is_qfisom(GEN F)1189 is_qfisom(GEN F)
1190 {
1191 return (lg(F)==6 && typ(F)==t_VEC && typ(gel(F,1))==t_VEC
1192 && typ(gel(F,3))==t_VEC && typ(gel(F,4))==t_VEC);
1193 }
1194
1195 static GEN
unpack_qfisominit(GEN F,GEN * norm,struct qfauto * qf,struct fingerprint * fp,struct qfcand * cand)1196 unpack_qfisominit(GEN F, GEN *norm, struct qfauto *qf,
1197 struct fingerprint *fp, struct qfcand *cand)
1198 {
1199 GEN QF = gel(F,3);
1200 qf->F = gel(QF,1);
1201 qf->V = gel(QF,2);
1202 qf->W = gel(QF,3);
1203 qf->v = gel(QF,4);
1204 qf->p = itou(gel(QF,5));
1205 qf->U = lg(QF)>6 ? gel(QF,6):NULL;
1206 QF = gel(F,4);
1207 fp->diag = gel(QF,1);
1208 fp->per = gel(QF,2);
1209 fp->e = gel(QF,3);
1210 QF = gel(F,5);
1211 cand->cdep =itos(gel(QF,1));
1212 cand->comb = gel(QF,2);
1213 cand->bacher_pol = gel(QF,3);
1214 *norm = gel(F,2);
1215 qf->dim = lg(gmael(F,1,1))-1;
1216 return qf->F;
1217 }
1218
1219 static GEN
qfisom_bestmat(GEN A,long * pt_max)1220 qfisom_bestmat(GEN A, long *pt_max)
1221 {
1222 long max = zm_maxdiag(A), max2;
1223 GEN A2, A1 = zm_to_ZM(A), U = lllgramint(A1);
1224 if (lg(U) != lg(A1))
1225 pari_err_DOMAIN("qfisom","form","is not",
1226 strtoGENstr("positive definite"), A1);
1227 A2 = ZM_to_zm(qf_apply_ZM(A1, U));
1228 max2 = zm_maxdiag(A2);
1229 if (max2 >= max) { *pt_max = max; return NULL; }
1230 *pt_max = max2; return mkvec2(ZM_to_zm(U), ZM_to_zm(ZM_inv(U,NULL)));
1231 }
1232
1233 static GEN
init_qfisom(GEN F,struct fingerprint * fp,struct qfcand * cand,struct qfauto * qf,GEN flags,long * max,GEN minvec)1234 init_qfisom(GEN F, struct fingerprint *fp, struct qfcand *cand,
1235 struct qfauto *qf, GEN flags, long *max, GEN minvec)
1236 {
1237 GEN U, A, norm;
1238 if (is_qfisom(F))
1239 {
1240 F = unpack_qfisominit(F, &norm, qf, fp, cand);
1241 A = gel(F,1);
1242 *max = zm_maxdiag(A);
1243 if (flags)
1244 init_flags(cand, A, fp, qf, flags);
1245 }
1246 else
1247 {
1248 if (lg(F)<2) pari_err_TYPE("qfisom",F);
1249 A = gel(F,1);
1250 if (lg(A)<2) pari_err_TYPE("qfisom",A);
1251 if (!minvec)
1252 {
1253 U = qfisom_bestmat(A, max);
1254 if (DEBUGLEVEL) err_printf("QFIsom: max=%ld\n",*max);
1255 if (U) F = zmV_apply_zm(F, gel(U,1));
1256 } else
1257 {
1258 *max = zm_maxdiag(A); U = NULL;
1259 if (typ(minvec)==t_VEC && lg(minvec)==4 && typ(gel(minvec,2))==t_INT)
1260 {
1261 long n = itos(gel(minvec,2));
1262 if (n != *max)
1263 pari_err_DOMAIN("qfisominit","m[2]","!=",stoi(*max),stoi(n));
1264 minvec = gel(minvec, 3);
1265 }
1266 if (typ(minvec)!=t_MAT || lg(gel(minvec,1))!=lg(A))
1267 pari_err_TYPE("qfisominit",minvec);
1268 }
1269 norm = init_qfauto(F, U, *max, qf, NULL, minvec);
1270 fingerprint(fp, qf);
1271 if (DEBUGLEVEL) err_printf("QFIsom: fp=%Ps\n",fp->diag);
1272 init_flags(cand, A, fp, qf, flags);
1273 }
1274 return norm;
1275 }
1276
1277 GEN
qfauto(GEN F,GEN flags)1278 qfauto(GEN F, GEN flags)
1279 {
1280 pari_sp av = avma;
1281 struct fingerprint fp;
1282 struct group G;
1283 struct qfcand cand;
1284 struct qfauto qf;
1285 long max;
1286 (void)init_qfisom(F, &fp, &cand, &qf, flags, &max, NULL);
1287 init_qfgroup(&G, &fp, &qf);
1288 autom(&G, &qf, &fp, &cand);
1289 return gerepilecopy(av, gen_group(&G, qf.U));
1290 }
1291
1292 static GEN
qf_to_zmV(GEN F)1293 qf_to_zmV(GEN F)
1294 {
1295 switch(typ(F))
1296 {
1297 case t_MAT: return RgM_is_ZM(F)? mkvec(ZM_to_zm(F)): NULL;
1298 case t_VEC: return RgV_is_ZMV(F)? ZMV_to_zmV(F): NULL;
1299 }
1300 return NULL;
1301 }
1302
1303 GEN
qfauto0(GEN x,GEN flags)1304 qfauto0(GEN x, GEN flags)
1305 {
1306 pari_sp av = avma;
1307 GEN F, G;
1308 if (is_qfisom(x)) F = x;
1309 else
1310 {
1311 F = qf_to_zmV(x);
1312 if (!F) pari_err_TYPE("qfauto",x);
1313 }
1314 G = qfauto(F, flags);
1315 return gerepilecopy(av, mkvec2(gel(G,1), zmV_to_ZMV(gel(G,2))));
1316 }
1317 /* computes the orbit of V.v[pt] under the generators G[0],...,G[nG-1] and
1318 * elements stabilizing V.v[pt], which are stored in H, returns the number of
1319 * generators in H */
1320 static GEN
isostab(long pt,GEN G,GEN V,long Maxfail,ulong p)1321 isostab(long pt, GEN G, GEN V, long Maxfail, ulong p)
1322 {
1323 pari_sp av = avma;
1324 long i, im, nH, fail, len, cnd, orblen, rpt;
1325 long dim = lg(gel(V,1))-1, n = lg(V)-1, nG = lg(G)-1;
1326 GEN w, flag, orb, H;
1327 H = cgetg(2,t_VEC); /* generators of the stabilizer of V.v[pt] */
1328 nH = 0;
1329 w = cgetg(2*n+2,t_MAT); /* w[i+V.n] is a matrix that maps V.v[pt] on V.v[i] */
1330 orb = zero_Flv(2*n);
1331 orblen = 1; /* length of the orbit of a random vector in V.v */
1332 flag = zero_Flv(2*n+1); /* if flag[i+V.n] = 1, then i is already in orbit */
1333 orb[1] = pt;
1334 flag[orb[1]+n+1] = 1;
1335 /* w[pt+V.n] is the Identity */
1336 gel(w,orb[1]+n+1) = matid_Flm(dim);
1337 cnd = len = 1;
1338 fail = 0; /* number of successive failures */
1339 while (cnd <= len && fail < Maxfail)
1340 {
1341 for (i = 1; i <= nG && fail < Maxfail; i++)
1342 {
1343 im = operate(orb[cnd], gel(G,i), V);
1344 if (flag[im+n+1] == 0)
1345 { /* new element is found, append to the orbit and store an element
1346 * mapping V.v[pt] to im in w[im+V.n] */
1347 orb[++len] = im;
1348 flag[im+n+1] = 1;
1349 gel(w,im+n+1) = zm_mul(gel(G,i), gel(w,orb[cnd]+n+1));
1350 }
1351 else
1352 { /* image was already in orbit */
1353 GEN B = zm_mul(gel(G,i), gel(w,orb[cnd]+n+1));
1354 /* check whether the old and the new element mapping pt on im differ */
1355 if (!zvV_equal(B, gel(w,im+n+1)))
1356 {
1357 long tmplen;
1358 gel(H,nH+1) = zm_divmod(gel(w,im+n+1),B,p);
1359 rpt = 1+(long)random_Fl(n);
1360 tmplen = orbitlen(rpt, 2*n, H, nH+1, V);
1361 while (tmplen < orblen)
1362 { /* orbit of this vector is shorter than a previous one:
1363 * choose new random vector */
1364 rpt = 1+(long)random_Fl(n);
1365 tmplen = orbitlen(rpt, 2*n, H, nH+1, V);
1366 }
1367 if (tmplen > orblen)
1368 { /* new stabilizer element H[nH] enlarges the group generated by H */
1369 orblen = tmplen;
1370 H = vec_lengthen(H, (++nH)+1); /* allocate for new generator */
1371 fail = 0;
1372 }
1373 else /* new stabilizer element does not enlarge random vector orbit */
1374 ++fail;
1375 }
1376 /* if H[nH] is the identity, do nothing */
1377 }
1378 }
1379 ++cnd;
1380 }
1381 setlg(H, nH+1); return gerepilecopy(av, H);
1382 }
1383
1384 /* the heart of the program: the recursion */
1385 static long
iso(long step,GEN x,GEN C,struct qfauto * qf,struct qfauto * qff,struct fingerprint * fp,GEN G,struct qfcand * cand)1386 iso(long step, GEN x, GEN C, struct qfauto *qf, struct qfauto *qff,
1387 struct fingerprint *fp, GEN G, struct qfcand *cand)
1388 {
1389 long dim = qf->dim;
1390 /* found isomorphism */
1391 if (step == dim && mael(C,dim,1)) { x[dim] = mael(C,dim,1); return 1; }
1392 while (mael(C,step,1))
1393 {
1394 long nbc;
1395 /* choose the image of the base-vector nr. step */
1396 x[step] = mael(C,step,1);
1397 /* check whether x[0]...x[step] is a partial automorphism and compute the
1398 candidates for x[step+1] */
1399 nbc = qfisom_candidates(gel(C,step+1), step+1, x, qf, qff, fp, cand);
1400 if (nbc == fp->diag[step+1])
1401 { /* go deeper into the recursion */
1402 long i, Maxfail = 0;
1403 GEN H;
1404 /* heuristic value of Maxfail for break condition in isostab */
1405 for (i = 1; i <= step; ++i)
1406 if (fp->diag[i] > 1) Maxfail++;
1407 for (; i <= dim; ++i)
1408 if (fp->diag[i] > 1) Maxfail += 2;
1409 /* compute the stabilizer H of x[step] in G */
1410 H = isostab(x[step], G, qff->V, Maxfail,qff->p);
1411 if (iso(step+1, x, C, qf, qff, fp, H, cand)) return 1;
1412 }
1413 /* delete the orbit of the chosen vector from the list of candidates */
1414 orbsubtract(gel(C,step), x, step-1, 1, G, qff->V, NULL);
1415 }
1416 return 0;
1417 }
1418
1419 /* search for an isometry */
1420 static GEN
isometry(struct qfauto * qf,struct qfauto * qff,struct fingerprint * fp,GEN G,struct qfcand * cand)1421 isometry(struct qfauto *qf, struct qfauto *qff, struct fingerprint *fp, GEN G,
1422 struct qfcand *cand)
1423
1424 {
1425 long i, dim = qf->dim;
1426 GEN x, C = cgetg(dim+1,t_VEC);
1427 /* C[i] is the list of candidates for the image of the i-th base-vector */
1428 for (i = 1; i <= dim; ++i) gel(C,i) = cgetg(fp->diag[i]+1, t_VECSMALL);
1429 x = cgetg(dim+1, t_VECSMALL);
1430 /* compute the candidates for x[1] */
1431 qfisom_candidates(gel(C,1), 1, x, qf, qff, fp, cand);
1432 return iso(1, x, C, qf, qff, fp, G, cand)? matgen(x, fp->per, qff->V): NULL;
1433 }
1434
1435 GEN
qfisominit(GEN F,GEN flags,GEN minvec)1436 qfisominit(GEN F, GEN flags, GEN minvec)
1437 {
1438 pari_sp av = avma;
1439 struct fingerprint fp;
1440 struct qfauto qf;
1441 struct qfcand cand;
1442 long max;
1443 GEN norm = init_qfisom(F, &fp, &cand, &qf, flags, &max, minvec);
1444 return gerepilecopy(av, mkvec5(F, norm,
1445 mkvecn(qf.U?6:5, qf.F, qf.V, qf.W, qf.v, utoi(qf.p), qf.U),
1446 mkvec3(fp.diag, fp.per, fp.e),
1447 mkvec3(stoi(cand.cdep),cand.comb?cand.comb:cgetg(1,t_VEC),
1448 cand.bacher_pol)));
1449 }
1450
1451 GEN
qfisominit0(GEN x,GEN flags,GEN minvec)1452 qfisominit0(GEN x, GEN flags, GEN minvec)
1453 {
1454 pari_sp av = avma;
1455 GEN F = qf_to_zmV(x);
1456 if (!F) pari_err_TYPE("qfisom",x);
1457 return gerepileupto(av, qfisominit(F, flags, minvec));
1458 }
1459
1460 GEN
qfisom(GEN F,GEN FF,GEN flags,GEN G)1461 qfisom(GEN F, GEN FF, GEN flags, GEN G)
1462 {
1463 pari_sp av = avma;
1464 struct fingerprint fp;
1465 GEN res, detf, detff;
1466 struct qfauto qf, qff;
1467 struct qfcand cand;
1468 long max;
1469 GEN norm = init_qfisom(F, &fp, &cand, &qf, flags, &max, NULL);
1470 init_qfauto(FF, NULL, max, &qff, norm, NULL);
1471 detf = ZM_det(zm_to_ZM(gel(qf.F,1)));
1472 detff = ZM_det(zm_to_ZM(gel(qff.F,1)));
1473 if (lg(qf.W)!=lg(qff.W) || !equalii(detf, detff)
1474 || !zvV_equal(vecvecsmall_sort_shallow(qf.W),
1475 vecvecsmall_sort_shallow(qff.W))) return gc_const(av,gen_0);
1476 if (!G) G = mkvec(scalar_Flm(-1, qff.dim));
1477 res = isometry(&qf, &qff, &fp, G, &cand);
1478 if (!res) return gc_const(av, gen_0);
1479 return gerepilecopy(av, zm_to_ZM(qf.U? zm_mul(res,gel(qf.U, 2)): res));
1480 }
1481
1482 static GEN
check_qfauto(GEN G)1483 check_qfauto(GEN G)
1484 {
1485 if (typ(G)==t_VEC && lg(G)==3 && typ(gel(G,1))==t_INT) G = gel(G,2);
1486 return qf_to_zmV(G);
1487 }
1488
1489 GEN
qfisom0(GEN x,GEN y,GEN flags,GEN G)1490 qfisom0(GEN x, GEN y, GEN flags, GEN G)
1491 {
1492 pari_sp av = avma;
1493 GEN F, FF;
1494 if (is_qfisom(x)) F = x;
1495 else
1496 {
1497 F = qf_to_zmV(x);
1498 if (!F) pari_err_TYPE("qfisom",x);
1499 }
1500 FF = qf_to_zmV(y);
1501 if (!FF) pari_err_TYPE("qfisom",y);
1502 if (G) G = check_qfauto(G);
1503 return gerepileupto(av, qfisom(F, FF, flags, G));
1504 }
1505
1506 static GEN
ZM_to_GAP(GEN M)1507 ZM_to_GAP(GEN M)
1508 {
1509 pari_sp av = avma;
1510 long i, j, c, rows = nbrows(M), cols = lg(M)-1;
1511 GEN comma = strtoGENstr(", "), bra = strtoGENstr("["), ket = strtoGENstr("]");
1512 GEN s = cgetg(2*rows*cols+2*rows+2,t_VEC);
1513 gel(s,1) = bra; c=2;
1514 for (i = 1; i <= rows; i++)
1515 {
1516 if (i > 1) gel(s,c++) = comma;
1517 gel(s,c++) = bra;
1518 for (j = 1; j <= cols; j++)
1519 {
1520 if (j > 1) gel(s,c++) = comma;
1521 gel(s,c++) = GENtoGENstr(gcoeff(M,i,j));
1522 }
1523 gel(s,c++) = ket;
1524 }
1525 gel(s,c++) = ket;
1526 return gerepilecopy(av, shallowconcat1(s));
1527 }
1528
1529 GEN
qfautoexport(GEN G,long flag)1530 qfautoexport(GEN G, long flag)
1531 {
1532 pari_sp av = avma;
1533 long i, lgen, c = 2;
1534 GEN gen, str, comma = strtoGENstr(", ");
1535 if (typ(G)!=t_VEC || lg(G)!=3) pari_err_TYPE("qfautoexport", G);
1536 if (flag!=0 && flag!=1) pari_err_FLAG("qfautoexport");
1537 gen = gel(G,2); lgen = lg(gen)-1;
1538 str = cgetg(2+2*lgen,t_VEC);
1539 /* in GAP or MAGMA the matrix group is called BG */
1540 if (flag == 0)
1541 gel(str,1) = strtoGENstr("Group(");
1542 else
1543 {
1544 long dim = lg(gmael(gen,1,1))-1;
1545 gel(str,1) = gsprintf("MatrixGroup<%d, Integers() |",dim);
1546 }
1547 for (i = 1; i <= lgen; i++)
1548 {
1549 if (i!=1) gel(str,c++) = comma;
1550 gel(str,c++) = ZM_to_GAP(gel(gen,i));
1551 }
1552 gel(str,c++) = strtoGENstr(flag ? ">":")");
1553 return gerepilecopy(av, shallowconcat1(str));
1554 }
1555
1556 GEN
qforbits(GEN G,GEN V)1557 qforbits(GEN G, GEN V)
1558 {
1559 pari_sp av = avma;
1560 GEN gen, w, W, p, v, orb, o;
1561 long i, j, n, ng, nborbits = 0;
1562 gen = check_qfauto(G);
1563 if (!gen) pari_err_TYPE("qforbits", G);
1564 if (typ(V)==t_VEC && lg(V)==4
1565 && typ(gel(V,1))==t_INT && typ(gel(V,2))==t_INT) V = gel(V,3);
1566 if (typ(V)!=t_MAT || !RgM_is_ZM(V)) pari_err_TYPE("qforbits", V);
1567 n = lg(V)-1; ng = lg(gen)-1;
1568 W = ZM_to_zm_canon(V);
1569 p = vecvecsmall_indexsort(W);
1570 v = vecpermute(W, p);
1571 w = zero_zv(n);
1572 orb = cgetg(n+1, t_VEC);
1573 o = cgetg(n+1, t_VECSMALL);
1574 if (lg(v) != lg(V)) return gen_0;
1575 for (i = 1; i <= n; i++)
1576 {
1577 long cnd, no = 1;
1578 GEN T;
1579 if (w[i]) continue;
1580 w[i] = ++nborbits;
1581 o[1] = i;
1582 for (cnd=1; cnd <= no; ++cnd)
1583 for (j=1; j <= ng; j++)
1584 {
1585 GEN Vij = zm_zc_mul(gel(gen, j), gel(v, o[cnd]));
1586 long k;
1587 (void) zv_canon(Vij);
1588 k = vecvecsmall_search(v, Vij, 0);
1589 if (k == 0) return gc_const(av, gen_0);
1590 if (w[k] == 0) { o[++no] = k; w[k] = nborbits; }
1591 }
1592 gel(orb, nborbits) = T = cgetg(no+1, t_VEC);
1593 for (j=1; j<=no; j++) gel(T,j) = gel(V,p[o[j]]);
1594 }
1595 setlg(orb, nborbits+1); return gerepilecopy(av, orb);
1596 }
1597