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.14.0 lcov report (development 27775-aca467eab2) Lines: 871 888 98.1 %
Date: 2022-07-03 07:33:15 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; 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    50840244 :     for (k = j; k < l; ++k) V[k] = -V[k];
      82    10163531 :     return -1;
      83             :   }
      84    12517624 :   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_apply_ZM(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         154 :   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_apply_ZM(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             : }

Generated by: LCOV version 1.13