Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - qfisom.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25685-2ac6d95b9c) Lines: 870 890 97.8 %
Date: 2020-08-06 06:06:08 Functions: 55 55 100.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13