Code coverage tests

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

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

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

Generated by: LCOV version 1.9