Code coverage tests

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

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

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

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

Generated by: LCOV version 1.11