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