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