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.10.0 lcov report (development 20459-9710128) Lines: 985 1027 95.9 %
Date: 2017-03-30 05:32:39 Functions: 54 54 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2012  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /*To be moved elsewhere */
      18             : 
      19             : static long
      20      354174 : Z_trunc(GEN z)
      21             : {
      22      354174 :   long s = signe(z);
      23      354174 :   return s==0 ? 0: (long)(s*mod2BIL(z));
      24             : }
      25             : 
      26             : static GEN
      27       50598 : ZV_trunc_to_zv(GEN z)
      28             : {
      29       50598 :   long i, l = lg(z);
      30       50598 :   GEN x = cgetg(l, t_VECSMALL);
      31       50598 :   for (i=1; i<l; i++) x[i] = Z_trunc(gel(z,i));
      32       50598 :   return x;
      33             : }
      34             : 
      35             : /* returns scalar product of vectors x and y with respect to Gram-matrix F */
      36             : static long
      37      173538 : scp(GEN x, GEN F, GEN y)
      38             : {
      39             :   long i, j;
      40      173538 :   ulong sum = 0;
      41      173538 :   long n = lg(F)-1;
      42     2893644 :   for (i = 1; i <= n; ++i)
      43             :   {
      44     2720106 :     ulong xi = uel(x,i);
      45     2720106 :     if (xi)
      46             :     {
      47     1118196 :       GEN Fi = gel(F, i);
      48    18683406 :       for (j = 1; j <= n; ++j)
      49    17565210 :         sum += xi * uel(Fi,j) * uel(y,j);
      50             :     }
      51             :   }
      52      173538 :   return sum;
      53             : }
      54             : 
      55             : static GEN
      56       10512 : ZM_trunc_to_zm(GEN z)
      57             : {
      58       10512 :   long i, l = lg(z);
      59       10512 :   GEN x = cgetg(l,t_MAT);
      60       10512 :   for (i=1; i<l; i++) gel(x,i) = ZV_trunc_to_zv(gel(z,i));
      61       10512 :   return x;
      62             : }
      63             : 
      64             : static GEN
      65        6768 : zm_divmod(GEN A, GEN B, ulong p)
      66             : {
      67        6768 :   pari_sp av = avma;
      68        6768 :   GEN Ap = zm_to_Flm(A,p), Bp = zm_to_Flm(B,p);
      69        6768 :   GEN C = Flm_center(Flm_mul(Flm_inv(Ap, p), Bp, p), p, p>>1);
      70        6768 :   return gerepileupto(av, C);
      71             : }
      72             : 
      73             : static int
      74    11910882 : zv_canon(GEN V)
      75             : {
      76    11910882 :   long l = lg(V), j, k;
      77    11910882 :   for (j = 1; j < l  &&  V[j] == 0; ++j);
      78    11910882 :   if (j < l  &&  V[j] < 0)
      79             :   {
      80    19401648 :     for (k = j; k < l; ++k)
      81    13778796 :       V[k] = -V[k];
      82     5622852 :     return -1;
      83             :   }
      84     6288030 :   return 1;
      85             : }
      86             : 
      87             : static GEN
      88         174 : ZM_to_zm_canon(GEN V)
      89             : {
      90         174 :   GEN W = ZM_to_zm(V);
      91         174 :   long i, l = lg(W);
      92      174000 :   for(i=1; i<l; i++)
      93      173826 :     (void)zv_canon(gel(W,i));
      94         174 :   return W;
      95             : }
      96             : 
      97             : static GEN
      98          24 : zm_apply_zm(GEN M, GEN U)
      99             : {
     100          24 :   return zm_mul(zm_transpose(U),zm_mul(M, U));
     101             : }
     102             : 
     103             : static GEN
     104          24 : zmV_apply_zm(GEN v, GEN U)
     105             : {
     106             :   long i, l;
     107          24 :   GEN V = cgetg_copy(v, &l);
     108          24 :   for(i=1; i<l; i++) gel(V,i) = zm_apply_zm(gel(v,i), U);
     109          24 :   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        9096 : possible(GEN F, GEN Ftr, GEN V, GEN W, GEN per, long I, long J)
     159             : {
     160             :   long i, j, k, count;
     161        9096 :   long n = lg(W)-1, f = lg(F)-1;
     162        9096 :   count = 0;
     163    14179464 :   for (j = 1; j <= n; ++j)
     164             :   {
     165    14170368 :     GEN Wj = gel(W,j), Vj = gel(V,j);
     166    14170368 :     i = I+1;
     167             :     /* check the length of the vector */
     168    26604096 :     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    18144888 :       for (i = 1; i <= I; ++i)
     171    17627874 :         if (zv_dotproduct(Vj,gmael(Ftr,k,per[i])) != coeff(gel(F,k),J,per[i]))
     172    11916714 :           break;
     173    14170368 :     if (k == f+1  &&  i > I)
     174      517014 :       ++count;
     175             :     /* the same for the negative vector */
     176    14170368 :     i = I+1;
     177    26604096 :     for (k = 1; k <= f  &&  i > I  &&  Wj[k] == mael3(F,k,J,J); ++k)
     178    18538116 :       for (i = 1; i <= I ; ++i)
     179    17924298 :         if (zv_dotproduct(Vj,gmael(Ftr,k,per[i])) != -coeff(gel(F,k),J,per[i]))
     180    11819910 :           break;
     181    14170368 :     if (k == f+1  &&  i > I)
     182      613818 :       ++count;
     183             :   }
     184        9096 :   return count;
     185             : }
     186             : 
     187             : static void
     188         120 : fingerprint(struct fingerprint *fp, struct qfauto *qf)
     189             : {
     190             :   pari_sp av;
     191         120 :   GEN V=qf->V, W=qf->W, F=qf->F;
     192             :   GEN Mf;
     193             :   long i, j, k, min;
     194         120 :   long dim = qf->dim, n = lg(V)-1, f = lg(F)-1;
     195             :   GEN Ftr;
     196         120 :   fp->per = identity_perm(dim);
     197         120 :   fp->e = cgetg(dim+1, t_VECSMALL);
     198         120 :   fp->diag = cgetg(dim+1, t_VECSMALL);
     199         120 :   av = avma;
     200         120 :   Ftr = cgetg(f+1,t_VEC);
     201         240 :   for (i = 1; i <= f; ++i)
     202         120 :     gel(Ftr,i) = zm_transpose(gel(F,i));
     203         120 :   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      124446 :   for (j = 1; j <= n; ++j)
     208             :   {
     209      124326 :     GEN Wj = gel(W,j);
     210     2057700 :     for (i = 1; i <= dim; ++i)
     211             :     {
     212     1933374 :       for (k = 1; k <= f  &&  Wj[k] == mael3(F,k,i,i); ++k);
     213     1933374 :       if (k == f+1)
     214     1266366 :         mael(Mf,1,i) += 2;
     215             :     }
     216             :   }
     217        1422 :   for (i = 1; i <= dim-1; ++i)
     218             :   {
     219             :     /* a minimal entry != 0 in the i-th row is chosen */
     220        1302 :     min = i;
     221       10398 :     for (j = i+1; j <= dim; ++j)
     222             :     {
     223        9096 :       if (mael(Mf,i,fp->per[j]) < mael(Mf,i,fp->per[min]))
     224         594 :         min = j;
     225             :     }
     226        1302 :     lswap(fp->per[i],fp->per[min]);
     227             :     /* the column below the minimal entry is set to 0 */
     228       10398 :     for (j = i+1; j <= dim; ++j)
     229        9096 :       mael(Mf,j,fp->per[i]) = 0;
     230             :     /* compute the row i+1 of the fingerprint */
     231       10398 :     for (j = i+1; j <= dim; ++j)
     232        9096 :       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        1542 :   for (i = 1; i <= dim; ++i)
     236        1422 :     fp->diag[i] = mael(Mf,i,fp->per[i]);
     237        1542 :   for (i = 1; i <= dim; ++i)
     238             :   {
     239        1422 :     fp->e[i] = vecvecsmall_search(V,vecsmall_ei(dim,fp->per[i]),0);
     240        1422 :     if (!fp->e[i])
     241           0 :       pari_err_BUG("qfisom, standard basis vector not found");
     242             :   }
     243         120 :   avma = av;
     244         120 : }
     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          36 : bacher(long I, long S, struct qfauto *qf)
     255             : {
     256          36 :   pari_sp av=avma;
     257          36 :   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          36 :   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          36 :   I = labs(I);
     266          36 :   vI = gel(V,I);
     267             :   /* list of vectors that have scalar product S with v[I] */
     268          36 :   list = zero_Flv(2*n);
     269          36 :   nlist = 0;
     270       53724 :   for (i = 1; i <= n; ++i)
     271       53688 :     if (mael(W,i,1) == mael(W,I,1))
     272             :     {
     273        4920 :       long s = zv_dotproduct(vI, gel(Fv,i));
     274        4920 :       if (s == S)
     275         696 :         list[++nlist] = i;
     276        4920 :       if (-s == S)
     277         288 :         list[++nlist] = -i;
     278             :     }
     279             :   /* there are nlist vectors that have scalar product S with v[I] */
     280          36 :   sum = nlist;
     281          36 :   if (nlist==0) retmkvec2(mkvecsmall3(0,0,0),cgetg(1,t_VEC));
     282          12 :   counts = cgetg(nlist+1, t_VECSMALL);
     283          12 :   listxy = cgetg(nlist+1, t_VECSMALL);
     284         996 :   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       81672 :     for (j = 1; j <= nlist; ++j)
     290       80688 :       listxy[j] = 0;
     291         984 :     nxy = 0;
     292         984 :     S1 = list[i] > 0 ? S : -S;
     293       81672 :     for (j = 1; j <= nlist; ++j)
     294             :     {
     295       80688 :       long S2 = list[j] > 0 ? S1 : -S1;
     296             :       /* note: for i > 0 is v[-i] = -v[i] */
     297       80688 :       if (zv_dotproduct(gel(V,labs(list[i])), gel(Fv,labs(list[j]))) == S2)
     298       25920 :         listxy[++nxy] = list[j];
     299             :     }
     300             :     /* counts[i] is the number of pairs for the vector v[list[i]] */
     301         984 :     counts[i] = 0;
     302       26904 :     for (j = 1; j <= nxy; ++j)
     303             :     {
     304       25920 :       long S1 = listxy[j] > 0 ? S : -S;
     305      362880 :       for (k = j+1; k <= nxy; ++k)
     306             :       {
     307      336960 :         long S2 = listxy[k] > 0 ? S1 : -S1;
     308      336960 :         long lj = labs(listxy[j]), lk = labs(listxy[k]);
     309      336960 :         if (zv_dotproduct(gel(V,lj), gel(Fv,lk)) == S2)
     310      155520 :           counts[i] += 1;
     311             :       }
     312             :     }
     313             :   }
     314             :    /* maxd is the maximal degree of the Bacher-polynomial,
     315             :       mind the minimal degree */
     316          12 :   maxd = counts[1];
     317          12 :   mind = counts[1];
     318         984 :   for (i = 2; i <= nlist; ++i)
     319             :   {
     320         972 :     if (counts[i] > maxd)
     321           0 :       maxd = counts[i];
     322         972 :     else if (counts[i] < mind)
     323          12 :       mind = counts[i];
     324             :   }
     325          12 :   coef = zero_Flv(maxd - mind + 1);
     326         996 :   for (i = 1; i <= nlist; ++i)
     327         984 :     coef[1+counts[i] - mind] += 1;
     328          12 :   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          12 :   return gerepilecopy(av, mkvec2(mkvecsmall3(sum, mind, maxd),coef));
     333             : }
     334             : 
     335             : static GEN
     336         120 : init_bacher(long bachdep, struct fingerprint *fp, struct qfauto *qf)
     337             : {
     338         120 :   GEN z = cgetg(bachdep+1,t_VEC);
     339             :   long i;
     340         156 :   for (i=1;i<=bachdep;i++)
     341             :   {
     342          36 :     long bachscp = mael(qf->W,fp->e[i],1) / 2;
     343          36 :     gel(z,i) = bacher(fp->e[i], bachscp, qf);
     344             :   }
     345         120 :   return z;
     346             : }
     347             : 
     348             : /* checks, whether the vector v[I] has the Bacher-polynomial pol  */
     349             : static long
     350         132 : bachcomp(GEN pol, long I, GEN V, GEN W, GEN Fv)
     351             : {
     352         132 :   pari_sp av = avma;
     353             :   GEN co, list, listxy, vI;
     354             :   long i, j, k;
     355             :   long nlist, nxy, count;
     356         132 :   const long n = lg(V)-1, S = mael(W,I,1) / 2;
     357         132 :   long sum = mael(pol,1,1), mind = mael(pol,1,2), maxd = mael(pol,1,3);
     358         132 :   GEN coef = gel(pol,2);
     359         132 :   vI = gel(V,I);
     360         132 :   list = zero_Flv(sum);
     361             :   /* nlist should be equal to pol.sum */
     362         132 :   nlist = 0;
     363      117588 :   for (i = 1; i <= n  &&  nlist <= sum; ++i)
     364             :   {
     365      117456 :     if (mael(W,i,1) == mael(W,I,1))
     366             :     {
     367       19044 :       long s = zv_dotproduct(vI, gel(Fv,i));
     368       19044 :       if (s == S)
     369             :       {
     370        3096 :         if (nlist < sum)
     371        3084 :           list[nlist+1] = i;
     372        3096 :         nlist++;
     373             :       }
     374       19044 :       if (-s == S)
     375             :       {
     376         876 :         if (nlist < sum)
     377         852 :           list[nlist+1] = -i;
     378         876 :         nlist++;
     379             :       }
     380             :     }
     381             :   }
     382         132 :   if (nlist != sum)
     383             :   {
     384             :     /* the number of vectors with scalar product S is already different */
     385          36 :     avma=av; return 0;
     386             :   }
     387          96 :   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          48 :   listxy = cgetg(nlist+1,t_VECSMALL);
     391          48 :   co = zero_Flv(maxd - mind + 1);
     392        3984 :   for (i = 1; i <= nlist; ++i)
     393             :   {
     394        3936 :     long S1 = list[i] > 0 ? S : -S;
     395      326688 :     for (j = 1; j <= nlist; ++j)
     396      322752 :       listxy[j] = 0;
     397        3936 :     nxy = 0;
     398      326688 :     for (j = 1; j <= nlist; ++j)
     399             :     {
     400      322752 :       long S2 = list[j] > 0 ? S1 : -S1;
     401      322752 :       if (zv_dotproduct(gel(V,labs(list[i])), gel(Fv,labs(list[j]))) == S2)
     402      103680 :         listxy[++nxy] = list[j];
     403             :     }
     404             :     /* count is the number of pairs */
     405        3936 :     count = 0;
     406      107616 :     for (j = 1; j <= nxy  &&  count <= maxd; ++j)
     407             :     {
     408      103680 :       long S1 = listxy[j] > 0 ? S : -S;
     409     1451520 :       for (k = j+1; k <= nxy  &&  count <= maxd; ++k)
     410             :       {
     411     1347840 :         long S2 = listxy[k] > 0 ? S1 : -S1;
     412     1347840 :         long lj = labs(listxy[j]), lk = labs(listxy[k]);
     413     1347840 :         if (zv_dotproduct(gel(V,lj), gel(Fv,lk)) == S2)
     414      622080 :           count++;
     415             :       }
     416             :     }
     417        7872 :     if (count < mind  ||  count > maxd  ||
     418        3936 :         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        3936 :       co[count-mind+1]++;
     428             :   }
     429             :   /* the Bacher-polynomials are equal */
     430          48 :   avma = av;
     431          48 :   return 1;
     432             : }
     433             : 
     434             : static GEN
     435         156 : checkvecs(GEN V, GEN F, GEN norm)
     436             : {
     437             :   long i, j, k;
     438         156 :   long n = lg(V)-1, f = lg(F)-1;
     439         156 :   GEN W = cgetg(n+1, t_MAT);
     440         156 :   j = 0;
     441      173694 :   for (i = 1; i <= n; ++i)
     442             :   {
     443      173538 :     GEN normvec = cgetg(f+1, t_VECSMALL);
     444      173538 :     GEN Vi = gel(V,i);
     445      347076 :     for (k = 1; k <= f; ++k)
     446      173538 :       normvec[k] = scp(Vi, gel(F,k), Vi);
     447      173538 :     if (!vecvecsmall_search(norm,normvec,0))
     448           0 :       ++j;
     449             :     else
     450             :     {
     451      173538 :       gel(V,i-j) = Vi;
     452      173538 :       gel(W,i-j) = normvec;
     453             :     }
     454             :   }
     455         156 :   setlg(V, n+1-j);
     456         156 :   setlg(W, n+1-j);
     457         156 :   return W;
     458             : }
     459             : 
     460             : static long
     461     1412364 : operate(long nr, GEN A, GEN V)
     462             : {
     463     1412364 :   pari_sp av = avma;
     464             :   long im,eps;
     465     1412364 :   GEN w = zm_zc_mul(A,gel(V,labs(nr)));
     466     1412364 :   eps = zv_canon(w);
     467     1412364 :   if (nr < 0) eps = -eps; /* -w */
     468     1412364 :   im = vecvecsmall_search(V,w,0);
     469     1412364 :   if (!im) pari_err_BUG("qfauto, image of vector not found");
     470     1412364 :   avma = av;
     471     1412364 :   return eps*im;
     472             : }
     473             : 
     474             : static GEN
     475        2406 : orbit(GEN pt, long ipt, long npt, GEN H, GEN V)
     476             : {
     477        2406 :   pari_sp av = avma;
     478             :   long i, cnd, im;
     479        2406 :   long n = lg(V)-1, nH = lg(H)-1, no = npt;
     480        2406 :   GEN flag = zero_Flv(2*n+1)+n+1; /*We need negative indices*/
     481        2406 :   GEN orb = cgetg(2*n+1,t_VECSMALL);
     482        4500 :   for (i = 1; i <= npt; ++i)
     483             :   {
     484        2094 :     orb[i] = pt[ipt+i];
     485        2094 :     flag[pt[ipt+i]] = 1;
     486             :   }
     487       48750 :   for (cnd=1; cnd <= no; ++cnd)
     488      282774 :     for (i = 1; i <= nH; ++i)
     489             :     {
     490      236430 :       im = operate(orb[cnd], gel(H,i), V);
     491      236430 :       if (flag[im] == 0)
     492             :         /* the image is a new point in the orbit */
     493             :       {
     494       44250 :         orb[++no] = im;
     495       44250 :         flag[im] = 1;
     496             :       }
     497             :     }
     498        2406 :   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        7212 : orbitlen(long pt, long orblen, GEN G, long nG, GEN V)
     505             : {
     506        7212 :   pari_sp av = avma;
     507             :   long i, len, cnd;
     508        7212 :   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        7212 :   flag = zero_Flv(2*n + 1)+n+1;
     513        7212 :   orb  = zero_Flv(orblen);
     514        7212 :   orb[1] = pt;
     515        7212 :   flag[pt] = 1;
     516        7212 :   len = 1;
     517      154872 :   for(cnd = 1; cnd <= len  &&  len < orblen; ++cnd)
     518      976188 :     for (i = 1; i <= nG  &&  len < orblen; ++i)
     519             :     {
     520      828528 :       long im = operate(orb[cnd], gel(G,i), V);
     521      828528 :       if (flag[im] == 0)
     522             :         /* the image is a new point in the orbit */
     523             :       {
     524      144360 :         orb[++len] = im;
     525      144360 :         flag[im] = 1;
     526             :       }
     527             :     }
     528        7212 :   avma = av;
     529        7212 :   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        6618 : orbdelete(GEN orb1, GEN orb2)
     537             : {
     538             :   long i, j, len;
     539        6618 :   long l1 = lg(orb1)-1;
     540        6618 :   long l2 = lg(orb2)-1;
     541        6618 :   for (i = 1; i <= l1  &&  orb1[i] != 0; ++i);
     542        6618 :   len = i - 1;
     543       57174 :   for (i = 1; i <= l2  &&  orb2[i] != 0; ++i)
     544             :   {
     545       50556 :     long o2i = orb2[i];
     546       50556 :     for (j = 1; j <= len  &&  orb1[j] != o2i; ++j);
     547             :     /* orb1[j] = orb2[i], hence delete orb1[j] from orb1 */
     548       50556 :     if (j <= len)
     549             :     {
     550       40314 :       orb1[j] = orb1[len];
     551       40314 :       orb1[len--] = 0;
     552             :     }
     553             :   }
     554        6618 :   return len;
     555             : }
     556             : 
     557             : static long
     558        2406 : orbsubtract(GEN Cs, GEN pt, long ipt, long npt, GEN H, GEN V, long *len)
     559             : {
     560        2406 :   pari_sp av = avma;
     561             :   long nC;
     562        2406 :   GEN orb = orbit(pt, ipt, npt, H, V);
     563        2406 :   if (len) *len = lg(orb)-1;
     564        2406 :   nC = orbdelete(Cs, orb);
     565        2406 :   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       14004 : matgen(GEN x, GEN per, GEN V)
     573             : {
     574             :   long i, j;
     575       14004 :   long dim = lg(x)-1;
     576       14004 :   GEN X = cgetg(dim+1,t_MAT);
     577      212010 :   for (i = 1; i <= dim; ++i)
     578             :   {
     579      198006 :     long xi = x[i];
     580      198006 :     GEN Xp = cgetg(dim+1,t_VECSMALL);
     581     3077892 :     for (j = 1; j <= dim; ++j)
     582     2879886 :       Xp[j] = xi > 0? mael(V,xi,j): -mael(V,-xi,j);
     583      198006 :     gel(X,per[i]) = Xp;
     584             :   }
     585       14004 :   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        6768 : stabil(GEN x1, GEN x2, GEN per, GEN G, GEN V, ulong p)
     594             : {
     595        6768 :   pari_sp av = avma;
     596             :   long i;
     597             :   GEN x, XG, X2;
     598        6768 :   long dim = lg(x1)-1;
     599        6768 :   x = cgetg(dim+1,t_VECSMALL);
     600      102474 :   for (i = 1; i <= dim; ++i)
     601       95706 :     x[i] = operate(x1[i], G, V);
     602             :   /* XG is the composite mapping of the matrix corresponding to x1 and G */
     603        6768 :   XG = matgen(x, per, V);
     604        6768 :   X2 = matgen(x2, per, V);
     605        6768 :   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         624 : 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         624 :   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        5766 :   for (Rest = 0, i = I; i <= dim; ++i)
     637        5142 :     if (fp->diag[i] > 1  &&  G->ord[i] < fp->diag[i])
     638        2634 :       ++Rest;
     639        9216 :   for (Maxfail = Rest, i = 1; i <= dim; ++i)
     640        8592 :     if (fp->diag[i] > 1)
     641        6864 :       ++Maxfail;
     642        5766 :   for (nH = 0, i = I; i <= dim; ++i)
     643        5142 :     nH += G->ng[i];
     644             :   /* H are the generators of the group in which the stabilizer is computed */
     645         624 :   H = cgetg(nH+1,t_MAT);
     646         624 :   Hj = cgetg(nH+2,t_MAT);
     647         624 :   k = 0;
     648        5766 :   for (i = I; i <= dim; ++i)
     649       11880 :     for (j = 1; j <= G->ng[i]; ++j)
     650        6738 :       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         624 :   w = cgetg(2*n+2,t_VEC);
     653             :   /* orb contains the orbit of fp.e[I] */
     654         624 :   orb = zero_Flv(2*n);
     655             :   /* if flag[i + V.n] = 1, then the point i is already in the orbit */
     656         624 :   flag = zero_Flv(2*n+1);
     657         624 :   orb[1] = fp->e[I];
     658         624 :   flag[orb[1]+n+1] = 1;
     659         624 :   gel(w,orb[1]+n+1) = cgetg(dim+1,t_VECSMALL);
     660        9216 :   for (i = 1; i <= dim; ++i)
     661        8592 :     mael(w,orb[1]+n+1,i) = fp->e[i];
     662         624 :   cnd = 1;
     663         624 :   len = 1;
     664             :   /* fail is the number of successive failures */
     665         624 :   fail = 0;
     666        4488 :   while (cnd <= len && fail < Maxfail+Rest)
     667             :   {
     668       30924 :     for (i = 1; i <= nH  &&  fail < Maxfail+Rest; ++i)
     669             :     {
     670       27684 :       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        3966 :         cnd = 1+(long)random_Fl(len);
     676        3966 :         i = 1+(long)random_Fl(nH);
     677             :       }
     678       27684 :       im = operate(orb[cnd], gel(H,i), V);
     679       27684 :       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       10908 :         orb[++len] = im;
     685       10908 :         flag[im+n+1] = 1;
     686       10908 :         wim = cgetg(dim+1,t_VECSMALL);
     687       10908 :         gel(w,im+n+1) = wim;
     688      146502 :         for (j = 1; j <= dim; ++j)
     689      135594 :           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       92412 :         for (j = I; j <= dim; j++)
     697       88350 :           if (operate(mael(w,orb[cnd]+n+1,j), gel(H,i), V) != mael(w,im+n+1,j))
     698       12714 :             break;
     699       29490 :         if (j <= dim  &&
     700       19776 :             (G->ord[j] < fp->diag[j] || fail >= Maxfail))
     701        6768 :         {
     702        6768 :           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        6768 :           S = stabil(wo, gel(w,im+n+1), fp->per, gel(H,i), V, p);
     705        6768 :           gel(Hj,1) = S;
     706        6768 :           nHj = 1;
     707       71676 :           for (k = j; k <= dim; ++k)
     708      107748 :             for (l = 1; l <= G->ng[k]; ++l)
     709       42840 :               gel(Hj,++nHj) = gmael(G->g,k,l);
     710        6768 :           tmplen = orbitlen(fp->e[j], fp->diag[j], Hj, nHj, V);
     711        6768 :           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        2862 :           {
     715             :             GEN Ggj;
     716        2862 :             G->ord[j] = tmplen;
     717        2862 :             G->ng[j]++;
     718        2862 :             G->nsg[j]++;
     719             :             /* allocate memory for the new generator */
     720        2862 :             gel(G->g,j) = vec_lengthen(gel(G->g,j),G->ng[j]);
     721        2862 :             Ggj = gel(G->g,j);
     722             :             /* the new generator is inserted as stabilizer element
     723             :                nr. nsg[j]-1 */
     724        2862 :             for (k = G->ng[j]; k > G->nsg[j]; --k)
     725           0 :               gel(Ggj,k) = gel(Ggj,k-1);
     726        2862 :             gel(Ggj,G->nsg[j]) = S;
     727        2862 :             nH++;
     728        2862 :             H = vec_lengthen(H, nH);
     729        2862 :             Hj = vec_lengthen(Hj, nH+1);
     730             :             /* the new generator is appended to H */
     731        2862 :             gel(H,nH) = gel(Ggj,G->nsg[j]);
     732             :             /* the number of failures is reset to 0 */
     733        2862 :             if (fail < Maxfail)
     734         714 :               fail = 0;
     735             :             else
     736        2148 :               ++fail;
     737             :           }
     738             :           else
     739             :             /*the new stabilizer element S does not enlarge the orbit of e[j]*/
     740        3906 :             ++fail;
     741             :         }
     742       10008 :         else if ((j < dim  &&  fail < Maxfail)  ||
     743          36 :             (j == dim  &&  fail >= Maxfail))
     744        5910 :           ++fail;
     745             :         /* if S is the identity and fail < Maxfail, nothing is done */
     746             :       }
     747             :     }
     748        3240 :     if (fail < Maxfail)
     749        2646 :       ++cnd;
     750             :   }
     751         624 : }
     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         120 : qfisom_candidates_novec(GEN CI, long I, GEN x, struct qfauto *qf,
     761             :        struct qfauto *qff, struct fingerprint *fp)
     762             : {
     763         120 :   pari_sp av = avma;
     764             :   long i, j, k, okp, okm, nr, fail;
     765             :   GEN vec;
     766         120 :   GEN F =qf->F, V=qff->V, W=qff->W, v=qff->v;
     767         120 :   long n = lg(V)-1, f = lg(F)-1;
     768         120 :   vec = cgetg(I,t_VECSMALL);
     769             :   /* CI is the list for the candidates */
     770       29316 :   for (i = 1; i <= fp->diag[I]; ++i)
     771       29196 :     CI[i] = 0;
     772         120 :   nr = 0;
     773         120 :   fail = 0;
     774      124446 :   for (j = 1; j <= n  &&  fail == 0; ++j)
     775             :   {
     776      124326 :     GEN Vj = gel(V,j), Wj = gel(W, j);
     777      124326 :     okp = 0;
     778      124326 :     okm = 0;
     779      248652 :     for (i = 1; i <= f; ++i)
     780             :     {
     781      124326 :       GEN FAiI = gmael(F,i,fp->per[I]);
     782      124326 :       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      124326 :       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      124326 :       for (k = 1; k < I  &&  vec[k] == FAiI[fp->per[k]]; ++k);
     794      124326 :       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       14598 :         ++okp;
     797      124326 :       for (k = 1; k < I  &&  vec[k] == -FAiI[fp->per[k]]; ++k);
     798      124326 :       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       14598 :         ++okm;
     801             :     }
     802      124326 :     if (okp == f)
     803             :       /* V.v[j] is a candidate for x[I] */
     804             :     {
     805       14598 :       if (nr < fp->diag[I])
     806       14598 :         CI[++nr] = j;
     807             :       else
     808             :         /* there are too many candidates */
     809           0 :         fail = 1;
     810             :     }
     811      124326 :     if (okm == f)
     812             :       /* -V.v[j] is a candidate for x[I] */
     813             :     {
     814       14598 :       if (nr < fp->diag[I])
     815       14598 :         CI[++nr] = -j;
     816             :       else
     817             :         /* there are too many candidates */
     818           0 :         fail = 1;
     819             :     }
     820             :   }
     821         120 :   if (fail == 1)
     822           0 :     nr = 0;
     823         120 :   avma = av;
     824         120 :   return nr;
     825             : }
     826             : 
     827             : static long
     828       10224 : qfisom_candidates(GEN CI, long I, GEN x, struct qfauto *qf,
     829             :       struct qfauto *qff, struct fingerprint *fp, struct qfcand *qfcand)
     830             : {
     831       10224 :   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       10224 :   GEN F =qf->F, V=qff->V, W=qff->W, v=qff->v, FF= qff->F;
     838       10224 :   long dim = qf->dim, n = lg(V)-1, f = lg(F)-1;
     839             :   long nc;
     840       10224 :   long DEP = qfcand->cdep, len = f * DEP;
     841       10224 :   if (I >= 2  &&  I <= lg(qfcand->bacher_pol))
     842             :   {
     843         132 :     long t = labs(x[I-1]);
     844         132 :     GEN bpolI = gel(qfcand->bacher_pol,I-1);
     845         132 :     if (bachcomp(bpolI, t, V, W, gel(v,1)) == 0) return 0;
     846             :   }
     847       10188 :   if (I==1 || DEP ==0)
     848         120 :     return qfisom_candidates_novec(CI,I,x,qf,qff,fp);
     849       10068 :   vec = cgetg(I,t_VECSMALL);
     850       10068 :   scpvec = cgetg(len+1,t_VECSMALL);
     851       10068 :   com = gel(qfcand->comb,I-1);
     852       10068 :   list=gel(com,1); trans = gel(com,2); ccoef = gel(com,3); cF=gel(com,4);
     853       10068 :   rank = lg(trans)-1;
     854       10068 :   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       10068 :   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       10068 :   xbase = cgetg(rank+1,t_MAT);
     861       39426 :   for (i = 1; i <= rank; ++i)
     862       29358 :     gel(xbase,i) = cgetg(dim+1,t_VECSMALL);
     863             :   /* Fxbase is the product of a form F with the base xbase */
     864       10068 :   Fxbase = cgetg(rank+1,t_MAT);
     865       39426 :   for (i = 1; i <= rank; ++i)
     866       29358 :     gel(Fxbase,i) = cgetg(dim+1,t_VECSMALL);
     867             :   /* CI is the list for the candidates */
     868       69198 :   for (i = 1; i <= fp->diag[I]; ++i)
     869       59130 :     CI[i] = 0;
     870       10068 :   nr = 0;
     871       10068 :   fail = 0;
     872    10238544 :   for (j = 1; j <= n  &&  fail == 0; ++j)
     873             :   {
     874             :     long sign;
     875    10228476 :     GEN Vj = gel(V,j), Wj = gel(W, j);
     876    10228476 :     okp = 0;
     877    10228476 :     okm = 0;
     878    26670072 :     for (k = 1; k <= len; ++k)
     879    16441596 :       scpvec[k] = 0;
     880    20456952 :     for (i = 1; i <= f; ++i)
     881             :     {
     882    10228476 :       GEN FAiI = gmael(F,i,fp->per[I]);
     883    10228476 :       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    91286676 :       for (k = 1; k < I; ++k)
     887             :       {
     888    81058200 :         long xk = x[k];
     889    81058200 :         if (xk > 0)
     890    49315188 :           vec[k] = zv_dotproduct(Vj, gel(FFvi,xk));
     891             :         else
     892    31743012 :           vec[k] = -zv_dotproduct(Vj, gel(FFvi,-xk));
     893             :       }
     894    10228476 :       for (k = 1; k < I  &&  vec[k] == FAiI[fp->per[k]]; ++k);
     895    10228476 :       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       25986 :         ++okp;
     898    10228476 :       for (k = 1; k < I  &&  vec[k] == -FAiI[fp->per[k]]; ++k);
     899    10228476 :       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       17376 :         ++okm;
     902    26400072 :       for (k = I-1; k >= 1  &&  k > I-1-DEP; --k)
     903    16171596 :         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    10228476 :     if (!zv_equal0(scpvec))
     908             :     {
     909     7152420 :       sign = zv_canon(scpvec);
     910     7152420 :       num = vecvecsmall_search(list,scpvec,0);
     911     7152420 :       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     7152420 :         GEN xnum = gel(xvec,num);
     920   114574212 :         for (k = 1; k <= dim; ++k)
     921   107421792 :           xnum[k] += sign * Vj[k];
     922             :       }
     923             :     }
     924    10228476 :     if (okp == f)
     925             :       /* V.v[j] is a candidate for x[I] */
     926             :     {
     927       25986 :       if (nr < fp->diag[I])
     928       25956 :         CI[++nr] = j;
     929             :       else
     930             :         /* there are too many candidates */
     931          30 :         fail = 1;
     932             :     }
     933    10228476 :     if (okm == f)
     934             :       /* -V.v[j] is a candidate for x[I] */
     935             :     {
     936       17376 :       if (nr < fp->diag[I])
     937       17334 :         CI[++nr] = -j;
     938             :       else
     939             :         /* there are too many candidates */
     940          42 :         fail = 1;
     941             :     }
     942             :   }
     943       10068 :   if (fail == 1)
     944          72 :     nr = 0;
     945       10068 :   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       26502 :     for (i = 1; i <= rank; ++i)
     950             :     {
     951       19770 :       GEN comtri = gel(trans,i);
     952      302250 :       for (j = 1; j <= dim; ++j)
     953             :       {
     954      282480 :         long xbij = 0;
     955     4132464 :         for (k = 1; k <= nc+1; ++k)
     956     3849984 :           xbij += comtri[k] * mael(xvec,k,j);
     957      282480 :         mael(xbase,i,j) = xbij;
     958             :       }
     959             :     }
     960             :     /* check, whether the base xbase has the right scalar products */
     961       13464 :     for (i = 1; i <= f &&  nr > 0; ++i)
     962             :     {
     963       26502 :       for (j = 1; j <= rank; ++j)
     964      302250 :         for (k = 1; k <= dim; ++k)
     965      282480 :           mael(Fxbase,j,k) = zv_dotproduct(gmael(FF,i,k), gel(xbase,j));
     966       26502 :       for (j = 1; j <= rank  &&  nr > 0; ++j)
     967       67914 :         for (k = 1; k <= j  &&  nr > 0; ++k)
     968             :         {
     969       48144 :           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       76230 :     for (i = 1; i <= nc+1  &&  nr > 0; ++i)
     975             :     {
     976       69498 :       GEN comcoi = gel(ccoef,i);
     977     1057674 :       for (j = 1; j <= dim && nr > 0; ++j)
     978             :       {
     979      988176 :         vj = 0;
     980     4838160 :         for (k = 1; k <= rank; ++k)
     981     3849984 :           vj += comcoi[k] * mael(xbase,k,j);
     982      988176 :         if (vj != mael(xvec,i,j))
     983             :           /* an entry is wrong */
     984           0 :           nr = 0;
     985             :       }
     986             :     }
     987             :   }
     988       10068 :   avma = av;
     989       10068 :   return nr;
     990             : }
     991             : 
     992             : static long
     993        5628 : aut(long step, GEN x, GEN C, struct group *G, struct qfauto *qf,
     994             :                              struct fingerprint *fp, struct qfcand *cand)
     995             : {
     996        5628 :   long dim = qf->dim;
     997        5628 :   GEN orb = cgetg(2,t_VECSMALL);
     998        5628 :   while (mael(C,step,1) != 0)
     999             :   {
    1000        8436 :     if (step < dim)
    1001             :     {
    1002             :       long nbc;
    1003             :       /* choose the image of the base-vector nr. step */
    1004        8040 :       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        8040 :       nbc = qfisom_candidates(gel(C,step+1), step+1, x, qf, qf, fp, cand);
    1008        8040 :       if (nbc == fp->diag[step+1])
    1009             :         /* go deeper into the recursion */
    1010        5124 :         if (aut(step+1, x, C, G, qf, fp, cand))
    1011        3828 :           return 1;
    1012        4212 :       orb[1] = x[step];
    1013             :       /* delete the chosen vector from the list of candidates */
    1014        4212 :       (void)orbdelete(gel(C,step), orb);
    1015             :     }
    1016             :     else
    1017             :     {
    1018             :       /* a new automorphism is found */
    1019         396 :       gel(x,dim) = gmael(C,dim,1);
    1020         396 :       return 1;
    1021             :     }
    1022             :   }
    1023        1404 :   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          84 : 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          84 :   GEN V = qf->V;
    1037          84 :   long dim = qf->dim, n = lg(V)-1;
    1038          84 :   long STAB = 1;
    1039          84 :   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        1062 :   for (i = 1; i <= dim; ++i)
    1042         978 :     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          84 :   x = cgetg(dim+1, t_VECSMALL);
    1045        1062 :   for (step = STAB; step <= dim; ++step)
    1046             :   {
    1047         978 :     if(DEBUGLEVEL) err_printf("QFAuto: Step %ld/%ld\n",step,dim);
    1048         978 :     nH = 0;
    1049        8052 :     for (nH = 0, i = step; i <= dim; ++i)
    1050        7074 :       nH += G->ng[i];
    1051         978 :     H = cgetg(nH+1,t_VEC);
    1052        8052 :     for (nH = 0, i = step; i <= dim; ++i)
    1053       15960 :       for (j = 1; j <= G->ng[i]; ++j)
    1054        8886 :         gel(H,++nH) = gmael(G->g,i,j);
    1055         978 :     bad = zero_Flv(2*n);
    1056         978 :     nbad = 0;
    1057             :     /* the first step base-vectors are fixed */
    1058        7074 :     for (i = 1; i < step; ++i)
    1059        6096 :       x[i] = fp->e[i];
    1060             :     /* compute the candidates for x[step] */
    1061         978 :     if (fp->diag[step] > 1)
    1062             :       /* if fp.diag[step] > 1 compute the candidates for x[step] */
    1063         780 :       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         198 :       mael(C,step,1) = fp->e[step];
    1068         198 :       nC = 1;
    1069             :     }
    1070             :     /* delete the orbit of the step-th base-vector from the candidates */
    1071         978 :     nC = orbsubtract(gel(C,step), fp->e, step-1, 1, H, V, &(G->ord[step]));
    1072        2856 :     while (nC > 0  &&  (im = mael(C,step,1)) != 0)
    1073             :     {
    1074         900 :       found = 0;
    1075             :       /* tries vector V.v[im] as image of the step-th base-vector */
    1076         900 :       x[step] = im;
    1077         900 :       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         864 :         nbc = qfisom_candidates(gel(C,step+1), step+1, x, qf, qf, fp, cand);
    1083         864 :         if (nbc == fp->diag[step+1])
    1084             :           /* go into the recursion */
    1085         504 :           found = aut(step+1, x, C, G, qf, fp, cand);
    1086             :         else
    1087         360 :           found = 0;
    1088             :       }
    1089             :       else
    1090          36 :         found = 1;
    1091         900 :       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         468 :         nC = orbsubtract(gel(C,step),mkvecsmall(im), 0, 1, H, V, NULL);
    1096         468 :         bad[++nbad] = im;
    1097             :       }
    1098             :       else
    1099             :         /* a new generator has been found */
    1100             :       {
    1101             :         GEN Gstep;
    1102         432 :         ++G->ng[step];
    1103             :         /* append the new generator to G->g[step] */
    1104         432 :         Gstep = vec_lengthen(gel(G->g,step),G->ng[step]);
    1105         432 :         gel(Gstep,G->ng[step]) = matgen(x, fp->per, V);
    1106         432 :         gel(G->g,step) = Gstep;
    1107         432 :         ++nH;
    1108         432 :         H = cgetg(nH+1, t_VEC);
    1109        5088 :         for (nH = 0, i = step; i <= dim; ++i)
    1110        8040 :           for (j = 1; j <= G->ng[i]; ++j)
    1111        3384 :             gel(H,++nH) = gmael(G->g,i,j);
    1112         432 :         nC = orbsubtract(gel(C,step), fp->e, step-1, 1, H, V, &(G->ord[step]));
    1113         432 :         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         978 :     if (step == STAB)
    1118         444 :       for (tries = G->nsg[step]; tries <= G->ng[step]; ++tries)
    1119             :       {
    1120         360 :         nH = 0;
    1121         702 :         for (j = 1; j < tries; ++j)
    1122         342 :           gel(H,++nH) = gmael(G->g,step,j);
    1123         702 :         for (j = tries+1; j < G->ng[step]; ++j)
    1124         342 :           gel(H,++nH) = gmael(G->g,step,j);
    1125        4332 :         for (i = step+1; i <= dim; ++i)
    1126        3972 :           for (j = 1; j <= G->ng[i]; ++j)
    1127           0 :             gel(H,++nH) = gmael(G->g,i,j);
    1128         360 :         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         978 :     if (step < dim  &&  G->ord[step] > 1)
    1138             :       /* calculate stabilizer elements fixing the basis-vectors
    1139             :          fp.e[0]...fp.e[step] */
    1140         624 :       stab(step, G, fp, V, qf->p);
    1141             :   }
    1142          84 : }
    1143             : 
    1144             : #define MAXENTRY (1L<<((BITS_IN_LONG-2)>>1))
    1145             : #define MAXNORM (1L<<(BITS_IN_LONG-2))
    1146             : 
    1147             : static long
    1148         288 : zm_maxdiag(GEN A)
    1149             : {
    1150         288 :   long dim = lg(A)-1;
    1151         288 :   long max = coeff(A,1,1);
    1152             :   long i;
    1153        3432 :   for (i = 2; i <= dim; ++i)
    1154        3144 :     if (coeff(A,i,i) > max)
    1155         132 :       max = coeff(A,i,i);
    1156         288 :   return max;
    1157             : }
    1158             : 
    1159             : static GEN
    1160         156 : 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         156 :   GEN M = minvec? minvec: gel(minim(zm_to_ZM(gel(F,1)), stoi(max), NULL), 3);
    1165         156 :   GEN V = ZM_to_zm_canon(M);
    1166         156 :   long n = lg(V)-1, f = lg(F)-1, dim = lg(gel(F,1))-1;
    1167      173694 :   for (i = 1; i <= n; ++i)
    1168             :   {
    1169      173538 :     GEN Vi = gel(V,i);
    1170     2893644 :     for (k = 1; k <= dim; ++k)
    1171             :     {
    1172     2720106 :       long l = labs(Vi[k]);
    1173     2720106 :       if (l > max)
    1174           0 :         max = l;
    1175             :     }
    1176             :   }
    1177         156 :   if (max > MAXENTRY) pari_err_OVERFLOW("qfisom [lattice too large]");
    1178         156 :   qf->p = unextprime(2*max+1);
    1179         156 :   V = vecvecsmall_sort_uniq(V);
    1180         156 :   if (!norm)
    1181             :   {
    1182         120 :     norm = cgetg(dim+1,t_VEC);
    1183        1542 :     for (i = 1; i <= dim; ++i)
    1184             :     {
    1185        1422 :       GEN Ni = cgetg(f+1,t_VECSMALL);
    1186        2844 :       for (k = 1; k <= f; ++k)
    1187        1422 :         Ni[k] = mael3(F,k,i,i);
    1188        1422 :       gel(norm,i) = Ni;
    1189             :     }
    1190         120 :     norm = vecvecsmall_sort_uniq(norm);
    1191             :   }
    1192         156 :   W = checkvecs(V, F, norm);
    1193         156 :   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         156 :   max = MAXNORM / max;
    1197         312 :   for (i = 1; i <= f; ++i)
    1198             :   {
    1199         156 :     GEN Fi = gel(F,i), vi;
    1200         156 :     vi = cgetg(n+1,t_MAT);
    1201         156 :     gel(v,i) = vi;
    1202      173694 :     for (j = 1; j <= n; ++j)
    1203             :     {
    1204      173538 :       GEN Vj = gel(V,j);
    1205      173538 :       GEN vij = cgetg(dim+1, t_VECSMALL);
    1206      173538 :       gel(vi,j) = vij;
    1207     2893644 :       for (k = 1; k <= dim; ++k)
    1208             :       {
    1209     2720106 :         vij[k] = zv_dotproduct(gel(Fi,k), Vj);
    1210     2720106 :         if (labs(vij[k]) > max) pari_err_OVERFLOW("qfisom [lattice too large]");
    1211             :       }
    1212             :     }
    1213             :   }
    1214         156 :   qf->dim = dim; qf->F = F; qf->V = V; qf->W = W; qf->v = v; qf->U = U;
    1215         156 :   return norm;
    1216             : }
    1217             : 
    1218             : static void
    1219          84 : init_qfgroup(struct group *G, struct fingerprint *fp, struct qfauto *qf)
    1220             : {
    1221          84 :   GEN H, M, V = qf->V;
    1222             :   long nH;
    1223             :   long i, j, k;
    1224          84 :   long dim = qf->dim;
    1225          84 :   G->ng  = zero_Flv(dim+1);
    1226          84 :   G->nsg = zero_Flv(dim+1);
    1227          84 :   G->ord = cgetg(dim+1,t_VECSMALL);
    1228          84 :   G->g = cgetg(dim+1,t_VEC);
    1229        1062 :   for (i = 1; i <= dim; ++i)
    1230         978 :     gel(G->g,i) = mkvec(gen_0);
    1231          84 :   M = matid_Flm(dim);
    1232          84 :   gmael(G->g,1,1) = M;
    1233          84 :   G->ng[1] = 1;
    1234             :   /* -Id is always an automorphism */
    1235        1062 :   for (i = 1; i <= dim; ++i)
    1236         978 :     mael(M,i,i) = -1;
    1237          84 :   nH = 0;
    1238        1062 :   for (i = 1; i <= dim; ++i)
    1239         978 :     nH += G->ng[i];
    1240          84 :   H = cgetg(nH+1,t_MAT);
    1241             :   /* calculate the orbit lengths under the automorphisms known so far */
    1242        1062 :   for (i = 1; i <= dim; ++i)
    1243             :   {
    1244         978 :     if (G->ng[i] > 0)
    1245             :     {
    1246          84 :       nH = 0;
    1247        1062 :       for (j = i; j <= dim; ++j)
    1248        1062 :         for (k = 1; k <= G->ng[j]; ++k)
    1249          84 :           gel(H,++nH) = gmael(G->g,j,k);
    1250          84 :       G->ord[i] = orbitlen(fp->e[i], fp->diag[i], H, nH, V);
    1251             :     }
    1252             :     else
    1253         894 :       G->ord[i] = 1;
    1254             :   }
    1255          84 : }
    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     4637484 : scpvector(GEN w, GEN b, long I, long dep, GEN v)
    1262             : {
    1263     4637484 :   long  i, j, n = lg(v)-1;
    1264     4637484 :   GEN scpvec = zero_Flv(dep*n);
    1265    12145008 :   for (i = I; i >= 1  &&  i > I-dep; --i)
    1266             :   {
    1267     7507524 :     long bi = b[i];
    1268     7507524 :     if (bi > 0)
    1269    15015048 :       for (j = 1; j <= n; ++j)
    1270     7507524 :         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     4637484 :   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        3588 : 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        3588 :   GEN V = qf->V, F = qf->F, v = qf->v;
    1287        3588 :   long n = lg(V)-1;
    1288        3588 :   long dim = lg(gel(F,1))-1;
    1289        3588 :   long len = (lg(F)-1)*dep;
    1290             :   /* the first vector in the list is the 0-vector and is not counted */
    1291        3588 :   list = mkmat(zero_Flv(len));
    1292        3588 :   vec  = mkmat(zero_Flv(dim));
    1293     4641072 :   for (j = 1; j <= n; ++j)
    1294             :   {
    1295     4637484 :     GEN Vvj = gel(V,j);
    1296     4637484 :     GEN scpvec = scpvector(Vvj, b, I, dep, v);
    1297     4637484 :     if (zv_equal0(scpvec))
    1298     1468092 :       nr = -1;
    1299             :     else
    1300             :     {
    1301     3169392 :       sign = zv_canon(scpvec);
    1302     3169392 :       nr = vecvecsmall_search(list,scpvec,0);
    1303             :     }
    1304             :     /* scpvec is already in list */
    1305     4637484 :     if (nr > 0)
    1306             :     {
    1307     3136974 :       vecnr = gel(vec,nr);
    1308    52123854 :       for (i = 1; i <= dim; ++i)
    1309    48986880 :         vecnr[i] += sign * Vvj[i];
    1310             :     }
    1311             :     /* scpvec is a new scalar product combination */
    1312     1500510 :     else if (nr==0)
    1313             :     {
    1314       32418 :       nr = vecvecsmall_search(list,scpvec,1);
    1315       32418 :       list=vec_insert(list,nr,scpvec);
    1316       32418 :       vec=vec_insert(vec,nr,sign < 0 ? zv_neg(Vvj) : zv_copy(Vvj));
    1317             :     }
    1318             :   }
    1319        3588 :   settyp(list,t_MAT);
    1320        3588 :   *pt_vec = vec;
    1321        3588 :   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        3504 : scpforms(GEN b, struct qfauto *qf)
    1327             : {
    1328             :   long i, j, k;
    1329        3504 :   GEN F = qf->F;
    1330        3504 :   long n = lg(F)-1, dim = lg(gel(F,1))-1;
    1331        3504 :   long nb = lg(b)-1;
    1332        3504 :   GEN gram = cgetg(n+1, t_VEC);
    1333             :   /* Fbi is the list of products of F.A[i] with the vectors in b */
    1334        3504 :   GEN Fbi = cgetg(nb+1, t_MAT);
    1335       12216 :   for (j = 1; j <= nb; ++j)
    1336        8712 :     gel(Fbi, j) = cgetg(dim+1, t_VECSMALL);
    1337        7008 :   for (i = 1; i <= n; ++i)
    1338             :   {
    1339        3504 :     GEN FAi = gel(F,i);
    1340        3504 :     gel(gram, i) = cgetg(nb+1, t_MAT);
    1341       12216 :     for (j = 1; j <= nb; ++j)
    1342      132798 :       for (k = 1; k <= dim; ++k)
    1343      124086 :         mael(Fbi,j,k) = zv_dotproduct(gel(FAi,k), gel(b,j));
    1344       12216 :     for (j = 1; j <= nb; ++j)
    1345             :     {
    1346        8712 :       GEN comFij = cgetg(nb+1, t_VECSMALL);
    1347       38760 :       for (k = 1; k <= nb; ++k)
    1348       30048 :         comFij[k] = zv_dotproduct(gel(b,j), gel(Fbi,k));
    1349        8712 :       gmael(gram,i,j) = comFij;
    1350             :     }
    1351             :   }
    1352        3504 :   return gram;
    1353             : }
    1354             : 
    1355             : static GEN
    1356         360 : gen_comb(long cdep, GEN A, GEN e, struct qfauto *qf, long lim)
    1357             : {
    1358         360 :   long i, dim = lg(A)-1;
    1359         360 :   GEN comb = cgetg(dim+1,t_VEC);
    1360        3864 :   for (i = 1; i <= dim; ++i)
    1361             :   {
    1362        3588 :     pari_sp av = avma;
    1363             :     GEN trans, ccoef, cF, B, BI;
    1364             :     GEN sumveclist, sumvecbase;
    1365        3588 :     GEN list = scpvecs(&sumveclist, i, e, cdep, qf);
    1366        3588 :     GEN M = zm_to_ZM(sumveclist);
    1367        3588 :     GEN T = lllgramint(qf_apply_ZM(A,M));
    1368        3588 :     if (lim && lg(T)-1>=lim) return NULL;
    1369        3504 :     B = ZM_mul(M,T);
    1370        3504 :     BI = RgM_solve(B,NULL);
    1371        3504 :     sumvecbase = ZM_trunc_to_zm(B);
    1372        3504 :     trans = ZM_trunc_to_zm(T);
    1373        3504 :     ccoef = ZM_trunc_to_zm(RgM_mul(BI,M));
    1374        3504 :     cF = scpforms(sumvecbase, qf);
    1375        3504 :     gel(comb,i) = gerepilecopy(av, mkvec4(list, trans, ccoef, cF));
    1376             :   }
    1377         276 :   return comb;
    1378             : }
    1379             : 
    1380             : static void
    1381          84 : init_comb(struct qfcand *cand, GEN A, GEN e, struct qfauto *qf)
    1382             : {
    1383          84 :   long dim = lg(A)-1;
    1384          84 :   GEN Am = zm_to_ZM(A);
    1385         240 :   for (cand->cdep = 1; ; cand->cdep++)
    1386             :   {
    1387         240 :     cand->comb = gen_comb(cand->cdep, Am, e, qf, (dim+1)>>1);
    1388         240 :     if (!cand->comb) break;
    1389         156 :   }
    1390          84 :   cand->cdep= maxss(1, cand->cdep-1);
    1391          84 :   cand->comb = gen_comb(cand->cdep, Am, e, qf, 0);
    1392          84 : }
    1393             : 
    1394             : static void
    1395         120 : init_flags(struct qfcand *cand, GEN A, struct fingerprint *fp,
    1396             :                                        struct qfauto *qf, GEN flags)
    1397             : {
    1398         120 :   if (!flags)
    1399             :   {
    1400          84 :     init_comb(cand, A, fp->e, qf);
    1401          84 :     cand->bacher_pol = init_bacher(0, fp, qf);
    1402             :   }
    1403             :   else
    1404             :   {
    1405             :     long cdep, bach;
    1406          36 :     if (typ(flags)!=t_VEC || lg(flags)!=3)
    1407           0 :       pari_err_TYPE("qfisominit",flags);
    1408          36 :     cdep = gtos(gel(flags,1));
    1409          36 :     bach = minss(gtos(gel(flags,2)),lg(fp->e)-1);
    1410          36 :     if (cdep<0 || bach<0) pari_err_FLAG("qfisom");
    1411          36 :     cand->cdep = cdep;
    1412          36 :     cand->comb = cdep ? gen_comb(cdep, zm_to_ZM(A), fp->e, qf, 0): NULL;
    1413          36 :     cand->bacher_pol = init_bacher(bach, fp, qf);
    1414             :   }
    1415         120 : }
    1416             : 
    1417             : static GEN
    1418          84 : gen_group(struct group *G, GEN U)
    1419             : {
    1420             :   GEN V;
    1421          84 :   long i, j, n=1, dim = lg(G->ord)-1;
    1422          84 :   GEN o = gen_1;
    1423        1062 :   for (i = 1; i <= dim; ++i)
    1424         978 :     o = muliu(o, G->ord[i]);
    1425        1062 :   for (i = 1; i <= dim; ++i)
    1426         978 :     n += G->ng[i]-G->nsg[i];
    1427          84 :   V = cgetg(n, t_VEC); n = 1;
    1428        1062 :   for (i = 1; i <= dim; ++i)
    1429        1494 :     for (j=G->nsg[i]+1; j<=G->ng[i]; j++)
    1430        1074 :       gel(V,n++) = U ? zm_mul(gel(U,1), zm_mul(gmael(G->g,i,j), gel(U,2)))
    1431         558 :                      : gmael(G->g,i,j);
    1432          84 :   return mkvec2(o, V);
    1433             : }
    1434             : 
    1435             : static long
    1436         300 : is_qfisom(GEN F)
    1437             : {
    1438         732 :   return (lg(F)==6 && typ(F)==t_VEC && typ(gel(F,1))==t_VEC
    1439         420 :                    && typ(gel(F,3))==t_VEC && typ(gel(F,4))==t_VEC);
    1440             : }
    1441             : 
    1442             : static GEN
    1443          60 : unpack_qfisominit(GEN F, GEN *norm, struct qfauto *qf,
    1444             :       struct fingerprint *fp, struct qfcand *cand)
    1445             : {
    1446          60 :   GEN QF = gel(F,3);
    1447          60 :   qf->F = gel(QF,1);
    1448          60 :   qf->V = gel(QF,2);
    1449          60 :   qf->W = gel(QF,3);
    1450          60 :   qf->v = gel(QF,4);
    1451          60 :   qf->p = itou(gel(QF,5));
    1452          60 :   qf->U = lg(QF)>6 ? gel(QF,6):NULL;
    1453          60 :   QF = gel(F,4);
    1454          60 :   fp->diag = gel(QF,1);
    1455          60 :   fp->per  = gel(QF,2);
    1456          60 :   fp->e    = gel(QF,3);
    1457          60 :   QF = gel(F,5);
    1458          60 :   cand->cdep =itos(gel(QF,1));
    1459          60 :   cand->comb = gel(QF,2);
    1460          60 :   cand->bacher_pol = gel(QF,3);
    1461          60 :   *norm = gel(F,2);
    1462          60 :   qf->dim = lg(gmael(F,1,1))-1;
    1463          60 :   return qf->F;
    1464             : }
    1465             : 
    1466             : static GEN
    1467         108 : qfisom_bestmat(GEN A, long *pt_max)
    1468             : {
    1469         108 :   long max = zm_maxdiag(A), max2;
    1470         108 :   GEN A1 = zm_to_ZM(A), A2;
    1471         108 :   GEN U = lllgramint(A1);
    1472         108 :   if (lg(U) != lg(A1))
    1473           0 :     pari_err_DOMAIN("qfisom","form","is not", strtoGENstr("positive definite"), A1);
    1474         108 :   A2 = ZM_to_zm(qf_apply_ZM(A1, U));
    1475         108 :   max2 = zm_maxdiag(A2);
    1476         108 :   if (max2 < max)
    1477             :   {
    1478          24 :     *pt_max = max2; return mkvec2(ZM_to_zm(U),ZM_to_zm(ZM_inv(U, gen_1)));
    1479             :   } else
    1480             :   {
    1481          84 :     *pt_max = max; return NULL;
    1482             :   }
    1483             : }
    1484             : 
    1485             : static GEN
    1486         180 : 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         180 :   if (is_qfisom(F))
    1491             :   {
    1492          60 :     F = unpack_qfisominit(F, &norm, qf, fp, cand);
    1493          60 :     A = gel(F,1);
    1494          60 :     *max = zm_maxdiag(A);
    1495          60 :     if (flags)
    1496           0 :       init_flags(cand, A, fp, qf, flags);
    1497             :   }
    1498             :   else
    1499             :   {
    1500         120 :     if (lg(F)<2) pari_err_TYPE("qfisom",F);
    1501         120 :     A = gel(F,1);
    1502         120 :     if (lg(A)<2) pari_err_TYPE("qfisom",A);
    1503         120 :     if (!minvec)
    1504             :     {
    1505         108 :       U = qfisom_bestmat(A, max);
    1506         108 :       if (DEBUGLEVEL) err_printf("QFIsom: max=%ld\n",*max);
    1507         108 :       if (U) F = zmV_apply_zm(F, gel(U,1));
    1508             :     } else
    1509             :     {
    1510          12 :       *max = zm_maxdiag(A); U = NULL;
    1511          12 :       if (typ(minvec)==t_VEC && lg(minvec)==4 && typ(gel(minvec,2))==t_INT)
    1512             :       {
    1513           6 :         long n = itos(gel(minvec,2));
    1514           6 :         if (n != *max)
    1515           0 :           pari_err_DOMAIN("qfisominit","m[2]","!=",stoi(*max),stoi(n));
    1516           6 :         minvec = gel(minvec, 3);
    1517             :       }
    1518          12 :       if (typ(minvec)!=t_MAT || lg(gel(minvec,1))!=lg(A))
    1519           0 :         pari_err_TYPE("qfisominit",minvec);
    1520             :     }
    1521         120 :     norm = init_qfauto(F, U, *max, qf, NULL, minvec);
    1522         120 :     fingerprint(fp, qf);
    1523         120 :     if (DEBUGLEVEL) err_printf("QFIsom: fp=%Ps\n",fp->diag);
    1524         120 :     init_flags(cand, A, fp, qf, flags);
    1525             :   }
    1526         180 :   return norm;
    1527             : }
    1528             : 
    1529             : GEN
    1530          84 : qfauto(GEN F, GEN flags)
    1531             : {
    1532          84 :   pari_sp av = avma;
    1533             :   struct fingerprint fp;
    1534             :   struct group G;
    1535             :   struct qfcand cand;
    1536             :   struct qfauto qf;
    1537             :   long max;
    1538          84 :   (void)init_qfisom(F, &fp, &cand, &qf, flags, &max, NULL);
    1539          84 :   init_qfgroup(&G, &fp, &qf);
    1540          84 :   autom(&G, &qf, &fp, &cand);
    1541          84 :   return gerepilecopy(av, gen_group(&G, qf.U));
    1542             : }
    1543             : 
    1544             : static GEN
    1545         174 : qf_to_zmV(GEN F)
    1546             : {
    1547         348 :   return typ(F)==t_MAT ?
    1548         156 :      (RgM_is_ZM(F) ? mkvec(ZM_to_zm(F)): NULL)
    1549         366 :      : typ(F)==t_VEC ? (RgV_is_ZMV(F) ? ZMV_to_zmV(F): NULL)
    1550          36 :      : NULL;
    1551             : }
    1552             : 
    1553             : GEN
    1554          84 : qfauto0(GEN x, GEN flags)
    1555             : {
    1556          84 :   pari_sp av = avma;
    1557             :   GEN F, G;
    1558          84 :   if (is_qfisom(x))
    1559          42 :     F = x;
    1560             :   else
    1561             :   {
    1562          42 :     F = qf_to_zmV(x);
    1563          42 :     if (!F) pari_err_TYPE("qfauto",x);
    1564             :   }
    1565          84 :   G = qfauto(F, flags);
    1566          84 :   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         408 : isostab(long pt, GEN G, GEN V, long Maxfail, ulong p)
    1574             : {
    1575         408 :   pari_sp av = avma;
    1576             :   long  len, cnd, orblen, tmplen, rpt;
    1577             :   GEN  w, flag, orb;
    1578             :   long  i, im, nH, fail;
    1579         408 :   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         408 :   H = cgetg(2,t_VEC);
    1592         408 :   nH = 0;
    1593             :   /* w[i+V.n] is a matrix that maps V.v[pt] on V.v[i] */
    1594         408 :   w = cgetg(2*n+2,t_MAT);
    1595         408 :   orb = zero_Flv(2*n);
    1596             :   /* orblen is the length of the orbit of a random vector in V.v */
    1597         408 :   orblen = 1;
    1598             :   /* if flag[i+V.n] = 1, then the point i is already in the orbit */
    1599         408 :   flag = zero_Flv(2*n+1);
    1600         408 :   orb[1] = pt;
    1601         408 :   flag[orb[1]+n+1] = 1;
    1602             :   /* w[pt+V.n] is the Identity */
    1603         408 :   gel(w,orb[1]+n+1) = matid_Flm(dim);
    1604         408 :   cnd = 1;
    1605         408 :   len = 1;
    1606             :   /* fail is the number of successive failures */
    1607         408 :   fail = 0;
    1608        1260 :   while (cnd <= len &&  fail < Maxfail)
    1609             :   {
    1610         516 :     for (i = 1; i <= nG  &&  fail < Maxfail; ++i)
    1611             :     {
    1612          72 :       im = operate(orb[cnd], gel(G,i), V);
    1613          72 :       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          36 :         orb[++len] = im;
    1618          36 :         flag[im+n+1] = 1;
    1619          36 :         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          36 :         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          36 :         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         444 :     ++cnd;
    1655             :   }
    1656         408 :   setlg(H,nH+1);
    1657         408 :   return gerepilecopy(av, H);
    1658             : }
    1659             : 
    1660             : /* the heart of the program: the recursion */
    1661             : 
    1662             : static long
    1663         444 : 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         444 :   long dim = qf->dim;
    1669         444 :   long found = 0;
    1670        1020 :   while (mael(C,step,1) != 0  &&  found == 0)
    1671             :   {
    1672         540 :     if (step < dim)
    1673             :     {
    1674             :       long nbc;
    1675             :       /* choose the image of the base-vector nr. step */
    1676         504 :       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         504 :       nbc = qfisom_candidates(gel(C,step+1), step+1, x, qf, qff, fp, cand);
    1680         504 :       if (nbc == fp->diag[step+1])
    1681             :       {
    1682             :         /* go deeper into the recursion */
    1683         408 :         Maxfail = 0;
    1684             :         /* determine the heuristic value of Maxfail for the break condition in
    1685             :            isostab */
    1686        3408 :         for (i = 1; i <= step; ++i)
    1687        3000 :           if (fp->diag[i] > 1)
    1688        2400 :             Maxfail += 1;
    1689        3408 :         for (i = step+1; i <= dim; ++i)
    1690        3000 :           if (fp->diag[i] > 1)
    1691        2520 :             Maxfail += 2;
    1692             :         /* compute the stabilizer H of x[step] in G */
    1693         408 :         H = isostab(x[step], G, qff->V, Maxfail,qff->p);
    1694         408 :         found = iso(step+1, x, C, qf, qff, fp, H, cand);
    1695             :       }
    1696         504 :       if (found == 1)
    1697         408 :         return 1;
    1698             :       /* delete the orbit of the chosen vector from the list of candidates */
    1699          96 :       orbsubtract(gel(C,step), x, step-1, 1, G, qff->V, NULL);
    1700             :     }
    1701             :     else
    1702             :     {
    1703             :       /* an isomorphism is found */
    1704          36 :       x[dim] = mael(C,dim,1);
    1705          36 :       found = 1;
    1706             :     }
    1707             :   }
    1708          36 :   return found;
    1709             : }
    1710             : 
    1711             : /* search for an isometry */
    1712             : 
    1713             : static GEN
    1714          36 : 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          36 :   long dim = qf->dim;
    1721          36 :   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         480 :   for (i = 1; i <= dim; ++i)
    1724         444 :     gel(C,i) = cgetg(fp->diag[i]+1, t_VECSMALL);
    1725          36 :   x = cgetg(dim+1, t_VECSMALL);
    1726             :   /* compute the candidates for x[1] */
    1727          36 :   qfisom_candidates(gel(C,1), 1, x, qf, qff, fp, cand);
    1728          36 :   found = iso(1, x, C, qf, qff, fp, G, cand);
    1729          36 :   return found ? matgen(x, fp->per, qff->V): NULL;
    1730             : }
    1731             : 
    1732             : GEN
    1733          60 : qfisominit(GEN F, GEN flags, GEN minvec)
    1734             : {
    1735          60 :   pari_sp av = avma;
    1736             :   struct fingerprint fp;
    1737             :   struct qfauto qf;
    1738             :   struct qfcand cand;
    1739             :   long max;
    1740          60 :   GEN norm = init_qfisom(F, &fp, &cand, &qf, flags, &max, minvec);
    1741         180 :   return gerepilecopy(av, mkvec5(F, norm,
    1742          60 :                           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          60 :                           mkvec3(stoi(cand.cdep),cand.comb?cand.comb:cgetg(1,t_VEC),
    1745             :                                  cand.bacher_pol)));
    1746             : }
    1747             : 
    1748             : GEN
    1749          60 : qfisominit0(GEN x, GEN flags, GEN minvec)
    1750             : {
    1751          60 :   pari_sp av = avma;
    1752          60 :   GEN F = qf_to_zmV(x);
    1753          60 :   if (!F) pari_err_TYPE("qfisom",x);
    1754          60 :   return gerepileupto(av, qfisominit(F, flags, minvec));
    1755             : }
    1756             : 
    1757             : GEN
    1758          36 : qfisom(GEN F, GEN FF, GEN flags)
    1759             : {
    1760          36 :   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          36 :   GEN norm = init_qfisom(F, &fp, &cand, &qf, flags, &max, NULL);
    1767          36 :   init_qfauto(FF, NULL, max, &qff, norm, NULL);
    1768          36 :   if (lg(qf.W)!=lg(qff.W)
    1769          36 :       || !zvV_equal(vecvecsmall_sort(qf.W), vecvecsmall_sort(qff.W)))
    1770           0 :     { avma=av; return gen_0; }
    1771          36 :   G = mkvec(scalar_Flm(-1, qff.dim));
    1772          36 :   res = isometry(&qf, &qff, &fp, G, &cand);
    1773          36 :   if (!res) { avma=av; return gen_0; }
    1774          36 :   return gerepilecopy(av, zm_to_ZM(qf.U? zm_mul(res,gel(qf.U, 2)):res));
    1775             : }
    1776             : 
    1777             : GEN
    1778          36 : qfisom0(GEN x, GEN y, GEN flags)
    1779             : {
    1780          36 :   pari_sp av = avma;
    1781             :   GEN F, FF;
    1782          36 :   if (is_qfisom(x))
    1783          18 :     F = x;
    1784             :   else
    1785             :   {
    1786          18 :     F = qf_to_zmV(x);
    1787          18 :     if (!F) pari_err_TYPE("qfisom",x);
    1788             :   }
    1789          36 :   FF = qf_to_zmV(y);
    1790          36 :   if (!FF) pari_err_TYPE("qfisom",y);
    1791          36 :   return gerepileupto(av, qfisom(F, FF, flags));
    1792             : }
    1793             : 
    1794             : static GEN
    1795          36 : ZM_to_GAP(GEN M)
    1796             : {
    1797          36 :   pari_sp ltop=avma;
    1798          36 :   long rows = nbrows(M), cols = lg(M)-1;
    1799             :   long i, j, c;
    1800          36 :   GEN comma = strtoGENstr(", ");
    1801          36 :   GEN bra = strtoGENstr("[");
    1802          36 :   GEN ket = strtoGENstr("]");
    1803          36 :   GEN s = cgetg(2*rows*cols+2*rows+2,t_VEC);
    1804          36 :   gel(s,1) = bra; c=2;
    1805         180 :   for (i = 1; i <= rows; ++i)
    1806             :   {
    1807         144 :     if (i > 1)
    1808         108 :       gel(s,c++) = comma;
    1809         144 :     gel(s,c++) = bra;
    1810         720 :     for (j = 1; j <= cols; ++j)
    1811             :     {
    1812         576 :       if (j > 1)
    1813         432 :         gel(s,c++) = comma;
    1814         576 :       gel(s,c++) = GENtoGENstr(gcoeff(M,i,j));
    1815             :     }
    1816         144 :     gel(s,c++) = ket;
    1817             :   }
    1818          36 :   gel(s,c++) = ket;
    1819          36 :   return gerepilecopy(ltop,shallowconcat1(s));
    1820             : }
    1821             : 
    1822             : GEN
    1823          12 : qfautoexport(GEN G, long flag)
    1824             : {
    1825          12 :   pari_sp av = avma;
    1826          12 :   long i, lgen,  c = 2;
    1827          12 :   GEN gen, str, comma = strtoGENstr(", ");
    1828          12 :   if (typ(G)!=t_VEC || lg(G)!=3) pari_err_TYPE("qfautoexport", G);
    1829          12 :   if (flag!=0 && flag!=1) pari_err_FLAG("qfautoexport");
    1830          12 :   gen = gel(G,2); lgen = lg(gen)-1;
    1831          12 :   str = cgetg(2+2*lgen,t_VEC);
    1832             :   /* in GAP or MAGMA the matrix group is called BG */
    1833          12 :   if (flag == 0)
    1834           6 :     gel(str,1) = strtoGENstr("Group(");
    1835             :   else
    1836             :   {
    1837           6 :     long dim = lg(gmael(gen,1,1))-1;
    1838           6 :     gel(str,1) = gsprintf("MatrixGroup<%d, Integers() |",dim);
    1839             :   }
    1840          48 :   for(i = 1; i <= lgen; i++)
    1841             :   {
    1842          36 :     if (i!=1)
    1843          24 :       gel(str,c++) = comma;
    1844          36 :     gel(str,c++) = ZM_to_GAP(gel(gen,i));
    1845             :   }
    1846          12 :   gel(str,c++) = strtoGENstr(flag ? ">":")");
    1847          12 :   return gerepilecopy(av, shallowconcat1(str));
    1848             : }
    1849             : 
    1850             : GEN
    1851          18 : qforbits(GEN G, GEN V)
    1852             : {
    1853          18 :   pari_sp av = avma;
    1854             :   GEN gen, w, W, p, v, orb, o;
    1855             :   long i, j, n, ng;
    1856          18 :   long nborbits = 0;
    1857          18 :   if (typ(G)==t_VEC && lg(G)==3 && typ(gel(G,1))==t_INT)
    1858          12 :     G = gel(G,2);
    1859          18 :   gen = qf_to_zmV(G);
    1860          18 :   if (!gen) pari_err_TYPE("qforbits", G);
    1861          18 :   if (typ(V)==t_VEC && lg(V)==4
    1862           6 :    && typ(gel(V,1))==t_INT && typ(gel(V,2))==t_INT)
    1863           6 :     V = gel(V,3);
    1864          18 :   if (typ(V)!=t_MAT || !RgM_is_ZM(V)) pari_err_TYPE("qforbits", V);
    1865          18 :   n = lg(V)-1; ng = lg(gen)-1;
    1866          18 :   W = ZM_to_zm_canon(V);
    1867          18 :   p = vecvecsmall_indexsort(W);
    1868          18 :   v = vecpermute(W, p);
    1869          18 :   w = zero_zv(n);
    1870          18 :   orb = cgetg(n+1, t_VEC);
    1871          18 :   o = cgetg(n+1, t_VECSMALL);
    1872          18 :   if (lg(v) != lg(V)) return gen_0;
    1873         306 :   for (i=1; i<=n; i++)
    1874             :   {
    1875         288 :     long cnd, no = 1;
    1876             :     GEN T;
    1877         288 :     if (w[i]) continue;
    1878          36 :     w[i] = ++nborbits;
    1879          36 :     o[1] = i;
    1880         324 :     for (cnd=1; cnd <= no; ++cnd)
    1881        3168 :       for(j=1; j <= ng; j++)
    1882             :       {
    1883             :         long k;
    1884        2880 :         GEN Vij = zm_zc_mul(gel(gen, j), gel(v, o[cnd]));
    1885        2880 :         (void) zv_canon(Vij);
    1886        2880 :         k = vecvecsmall_search(v, Vij, 0);
    1887        2880 :         if (k == 0) { avma = av; return gen_0; }
    1888        2880 :         if (w[k] == 0)
    1889             :         {
    1890         252 :           o[++no] = k;
    1891         252 :           w[k] = nborbits;
    1892             :         }
    1893             :       }
    1894          36 :     T = cgetg(no+1, t_VEC);
    1895          36 :     for (j=1; j<=no; j++) gel(T,j) = gel(V,p[o[j]]);
    1896          36 :     gel(orb, nborbits) = T;
    1897             :   }
    1898          18 :   setlg(orb, nborbits+1);
    1899          18 :   return gerepilecopy(av, orb);
    1900             : }

Generated by: LCOV version 1.11