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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 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.8.0 lcov report (development 16937-4bd9b4e) Lines: 917 959 95.6 %
Date: 2014-10-24 Functions: 49 49 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 611 708 86.3 %

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

Generated by: LCOV version 1.9