Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - qfisom.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.0 lcov report (development 23036-b751c0af5) Lines: 997 1028 97.0 %
Date: 2018-09-26 05:46:06 Functions: 55 55 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2012  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. 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      544313 : Z_trunc(GEN z)
      21             : {
      22      544313 :   long s = signe(z);
      23      544313 :   return s==0 ? 0: (long)(s*mod2BIL(z));
      24             : }
      25             : 
      26             : static GEN
      27       75873 : ZV_trunc_to_zv(GEN z)
      28             : {
      29       75873 :   long i, l = lg(z);
      30       75873 :   GEN x = cgetg(l, t_VECSMALL);
      31       75873 :   for (i=1; i<l; i++) x[i] = Z_trunc(gel(z,i));
      32       75873 :   return x;
      33             : }
      34             : 
      35             : /* returns scalar product of vectors x and y with respect to Gram-matrix F */
      36             : static long
      37      245574 : scp(GEN x, GEN F, GEN y)
      38             : {
      39             :   long i, j;
      40      245574 :   ulong sum = 0;
      41      245574 :   long n = lg(F)-1;
      42     4107684 :   for (i = 1; i <= n; ++i)
      43             :   {
      44     3862110 :     ulong xi = uel(x,i);
      45     3862110 :     if (xi)
      46             :     {
      47     1570198 :       GEN Fi = gel(F, i);
      48    26311964 :       for (j = 1; j <= n; ++j)
      49    24741766 :         sum += xi * uel(Fi,j) * uel(y,j);
      50             :     }
      51             :   }
      52      245574 :   return sum;
      53             : }
      54             : 
      55             : static GEN
      56       15351 : ZM_trunc_to_zm(GEN z)
      57             : {
      58       15351 :   long i, l = lg(z);
      59       15351 :   GEN x = cgetg(l,t_MAT);
      60       15351 :   for (i=1; i<l; i++) gel(x,i) = ZV_trunc_to_zv(gel(z,i));
      61       15351 :   return x;
      62             : }
      63             : 
      64             : static GEN
      65        9947 : zm_divmod(GEN A, GEN B, ulong p)
      66             : {
      67        9947 :   pari_sp av = avma;
      68        9947 :   GEN Ap = zm_to_Flm(A,p), Bp = zm_to_Flm(B,p);
      69        9947 :   GEN C = Flm_center(Flm_mul(Flm_inv(Ap, p), Bp, p), p, p>>1);
      70        9947 :   return gerepileupto(av, C);
      71             : }
      72             : 
      73             : static int
      74    19679366 : zv_canon(GEN V)
      75             : {
      76    19679366 :   long l = lg(V), j, k;
      77    19679366 :   for (j = 1; j < l  &&  V[j] == 0; ++j);
      78    19679366 :   if (j < l  &&  V[j] < 0)
      79             :   {
      80    38267754 :     for (k = j; k < l; ++k)
      81    29513834 :       V[k] = -V[k];
      82     8753920 :     return -1;
      83             :   }
      84    10925446 :   return 1;
      85             : }
      86             : 
      87             : static GEN
      88         245 : ZM_to_zm_canon(GEN V)
      89             : {
      90         245 :   GEN W = ZM_to_zm(V);
      91         245 :   long i, l = lg(W);
      92      246155 :   for(i=1; i<l; i++)
      93      245910 :     (void)zv_canon(gel(W,i));
      94         245 :   return W;
      95             : }
      96             : 
      97             : static GEN
      98          42 : zm_apply_zm(GEN M, GEN U)
      99             : {
     100          42 :   return zm_mul(zm_transpose(U),zm_mul(M, U));
     101             : }
     102             : 
     103             : static GEN
     104          42 : zmV_apply_zm(GEN v, GEN U)
     105             : {
     106             :   long i, l;
     107          42 :   GEN V = cgetg_copy(v, &l);
     108          42 :   for(i=1; i<l; i++) gel(V,i) = zm_apply_zm(gel(v,i), U);
     109          42 :   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       12432 : possible(GEN F, GEN Ftr, GEN V, GEN W, GEN per, long I, long J)
     159             : {
     160             :   long i, j, k, count;
     161       12432 :   long n = lg(W)-1, f = lg(F)-1;
     162       12432 :   count = 0;
     163    19985868 :   for (j = 1; j <= n; ++j)
     164             :   {
     165    19973436 :     GEN Wj = gel(W,j), Vj = gel(V,j);
     166    19973436 :     i = I+1;
     167             :     /* check the length of the vector */
     168    37385880 :     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    25464866 :       for (i = 1; i <= I; ++i)
     171    24733226 :         if (zv_dotproduct(Vj,gmael(Ftr,k,per[i])) != coeff(gel(F,k),J,per[i]))
     172    16680804 :           break;
     173    19973436 :     if (k == f+1  &&  i > I)
     174      731640 :       ++count;
     175             :     /* the same for the negative vector */
     176    19973436 :     i = I+1;
     177    37385880 :     for (k = 1; k <= f  &&  i > I  &&  Wj[k] == mael3(F,k,J,J); ++k)
     178    25837413 :       for (i = 1; i <= I ; ++i)
     179    24990763 :         if (zv_dotproduct(Vj,gmael(Ftr,k,per[i])) != -coeff(gel(F,k),J,per[i]))
     180    16565794 :           break;
     181    19973436 :     if (k == f+1  &&  i > I)
     182      846650 :       ++count;
     183             :   }
     184       12432 :   return count;
     185             : }
     186             : 
     187             : static void
     188         168 : fingerprint(struct fingerprint *fp, struct qfauto *qf)
     189             : {
     190             :   pari_sp av;
     191         168 :   GEN V=qf->V, W=qf->W, F=qf->F;
     192             :   GEN Mf;
     193             :   long i, j, k, min;
     194         168 :   long dim = qf->dim, n = lg(V)-1, f = lg(F)-1;
     195             :   GEN Ftr;
     196         168 :   fp->per = identity_perm(dim);
     197         168 :   fp->e = cgetg(dim+1, t_VECSMALL);
     198         168 :   fp->diag = cgetg(dim+1, t_VECSMALL);
     199         168 :   av = avma;
     200         168 :   Ftr = cgetg(f+1,t_VEC);
     201         336 :   for (i = 1; i <= f; ++i)
     202         168 :     gel(Ftr,i) = zm_transpose(gel(F,i));
     203         168 :   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      173957 :   for (j = 1; j <= n; ++j)
     208             :   {
     209      173789 :     GEN Wj = gel(W,j);
     210     2888494 :     for (i = 1; i <= dim; ++i)
     211             :     {
     212     2714705 :       for (k = 1; k <= f  &&  Wj[k] == mael3(F,k,i,i); ++k);
     213     2714705 :       if (k == f+1)
     214     1749489 :         mael(Mf,1,i) += 2;
     215             :     }
     216             :   }
     217        1953 :   for (i = 1; i <= dim-1; ++i)
     218             :   {
     219             :     /* a minimal entry != 0 in the i-th row is chosen */
     220        1785 :     min = i;
     221       14217 :     for (j = i+1; j <= dim; ++j)
     222             :     {
     223       12432 :       if (mael(Mf,i,fp->per[j]) < mael(Mf,i,fp->per[min]))
     224         777 :         min = j;
     225             :     }
     226        1785 :     lswap(fp->per[i],fp->per[min]);
     227             :     /* the column below the minimal entry is set to 0 */
     228       14217 :     for (j = i+1; j <= dim; ++j)
     229       12432 :       mael(Mf,j,fp->per[i]) = 0;
     230             :     /* compute the row i+1 of the fingerprint */
     231       14217 :     for (j = i+1; j <= dim; ++j)
     232       12432 :       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        2121 :   for (i = 1; i <= dim; ++i)
     236        1953 :     fp->diag[i] = mael(Mf,i,fp->per[i]);
     237        2121 :   for (i = 1; i <= dim; ++i)
     238             :   {
     239        1953 :     fp->e[i] = vecvecsmall_search(V,vecsmall_ei(dim,fp->per[i]),0);
     240        1953 :     if (!fp->e[i])
     241           0 :       pari_err_BUG("qfisom, standard basis vector not found");
     242             :   }
     243         168 :   set_avma(av);
     244         168 : }
     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         168 : init_bacher(long bachdep, struct fingerprint *fp, struct qfauto *qf)
     337             : {
     338         168 :   GEN z = cgetg(bachdep+1,t_VEC);
     339             :   long i;
     340         210 :   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         168 :   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             :   /* the number of vectors with scalar product S is already different */
     383         154 :   if (nlist != sum) return gc_long(av,0);
     384         112 :   if (nlist == 0) return gc_long(av,1);
     385             :   /* listxy is the list of the nxy vectors from list that have scalar product S
     386             :      with v[list[i]] */
     387          56 :   listxy = cgetg(nlist+1,t_VECSMALL);
     388          56 :   co = zero_Flv(maxd - mind + 1);
     389        4648 :   for (i = 1; i <= nlist; ++i)
     390             :   {
     391        4592 :     long S1 = list[i] > 0 ? S : -S;
     392      381136 :     for (j = 1; j <= nlist; ++j)
     393      376544 :       listxy[j] = 0;
     394        4592 :     nxy = 0;
     395      381136 :     for (j = 1; j <= nlist; ++j)
     396             :     {
     397      376544 :       long S2 = list[j] > 0 ? S1 : -S1;
     398      376544 :       if (zv_dotproduct(gel(V,labs(list[i])), gel(Fv,labs(list[j]))) == S2)
     399      120960 :         listxy[++nxy] = list[j];
     400             :     }
     401             :     /* count is the number of pairs */
     402        4592 :     count = 0;
     403      125552 :     for (j = 1; j <= nxy  &&  count <= maxd; ++j)
     404             :     {
     405      120960 :       long S1 = listxy[j] > 0 ? S : -S;
     406     1693440 :       for (k = j+1; k <= nxy  &&  count <= maxd; ++k)
     407             :       {
     408     1572480 :         long S2 = listxy[k] > 0 ? S1 : -S1;
     409     1572480 :         long lj = labs(listxy[j]), lk = labs(listxy[k]);
     410     1572480 :         if (zv_dotproduct(gel(V,lj), gel(Fv,lk)) == S2)
     411      725760 :           count++;
     412             :       }
     413             :     }
     414        9184 :     if (count < mind  ||  count > maxd  ||
     415        4592 :         co[count-mind+1] >= coef[count-mind+1])
     416             :     /* if the number of pairs is smaller than pol.mind or larger than pol.maxd
     417             :        or if the coefficient of X^count becomes now larger than the one in pol,
     418             :        then the Bacher-polynomials can not be equal */
     419             :     {
     420           0 :       set_avma(av);
     421           0 :       return 0;
     422             :     }
     423             :     else
     424        4592 :       co[count-mind+1]++;
     425             :   }
     426             :   /* the Bacher-polynomials are equal */
     427          56 :   set_avma(av);
     428          56 :   return 1;
     429             : }
     430             : 
     431             : static GEN
     432         224 : checkvecs(GEN V, GEN F, GEN norm)
     433             : {
     434             :   long i, j, k;
     435         224 :   long n = lg(V)-1, f = lg(F)-1;
     436         224 :   GEN W = cgetg(n+1, t_MAT);
     437         224 :   j = 0;
     438      245798 :   for (i = 1; i <= n; ++i)
     439             :   {
     440      245574 :     GEN normvec = cgetg(f+1, t_VECSMALL);
     441      245574 :     GEN Vi = gel(V,i);
     442      491148 :     for (k = 1; k <= f; ++k)
     443      245574 :       normvec[k] = scp(Vi, gel(F,k), Vi);
     444      245574 :     if (!vecvecsmall_search(norm,normvec,0))
     445           0 :       ++j;
     446             :     else
     447             :     {
     448      245574 :       gel(V,i-j) = Vi;
     449      245574 :       gel(W,i-j) = normvec;
     450             :     }
     451             :   }
     452         224 :   setlg(V, n+1-j);
     453         224 :   setlg(W, n+1-j);
     454         224 :   return W;
     455             : }
     456             : 
     457             : static long
     458     4580121 : operate(long nr, GEN A, GEN V)
     459             : {
     460     4580121 :   pari_sp av = avma;
     461             :   long im,eps;
     462     4580121 :   GEN w = zm_zc_mul(A,gel(V,labs(nr)));
     463     4580121 :   eps = zv_canon(w);
     464     4580121 :   if (nr < 0) eps = -eps; /* -w */
     465     4580121 :   im = vecvecsmall_search(V,w,0);
     466     4580121 :   if (!im) pari_err_BUG("qfauto, image of vector not found");
     467     4580121 :   set_avma(av);
     468     4580121 :   return eps*im;
     469             : }
     470             : 
     471             : static GEN
     472        3185 : orbit(GEN pt, long ipt, long npt, GEN H, GEN V)
     473             : {
     474        3185 :   pari_sp av = avma;
     475             :   long i, cnd, im;
     476        3185 :   long n = lg(V)-1, nH = lg(H)-1, no = npt;
     477        3185 :   GEN flag = zero_Flv(2*n+1)+n+1; /*We need negative indices*/
     478        3185 :   GEN orb = cgetg(2*n+1,t_VECSMALL);
     479        5964 :   for (i = 1; i <= npt; ++i)
     480             :   {
     481        2779 :     orb[i] = pt[ipt+i];
     482        2779 :     flag[pt[ipt+i]] = 1;
     483             :   }
     484       58786 :   for (cnd=1; cnd <= no; ++cnd)
     485      347361 :     for (i = 1; i <= nH; ++i)
     486             :     {
     487      291760 :       im = operate(orb[cnd], gel(H,i), V);
     488      291760 :       if (flag[im] == 0)
     489             :         /* the image is a new point in the orbit */
     490             :       {
     491       52822 :         orb[++no] = im;
     492       52822 :         flag[im] = 1;
     493             :       }
     494             :     }
     495        3185 :   setlg(orb,no+1); return gerepileuptoleaf(av, orb);
     496             : }
     497             : 
     498             : /* return the length of the orbit of pt under the first nG matrices in G */
     499             : 
     500             : static long
     501       21518 : orbitlen(long pt, long orblen, GEN G, long nG, GEN V)
     502             : {
     503       21518 :   pari_sp av = avma;
     504             :   long i, len, cnd;
     505       21518 :   long n = lg(V)-1;
     506             :   GEN orb, flag;
     507             :   /* if flag[i + n+1] = 1, -n <= i <= n, then the point i is already in the
     508             :      orbit */
     509       21518 :   flag = zero_Flv(2*n + 1)+n+1;
     510       21518 :   orb  = zero_Flv(orblen);
     511       21518 :   orb[1] = pt;
     512       21518 :   flag[pt] = 1;
     513       21518 :   len = 1;
     514      711928 :   for(cnd = 1; cnd <= len  &&  len < orblen; ++cnd)
     515     4528776 :     for (i = 1; i <= nG  &&  len < orblen; ++i)
     516             :     {
     517     3838366 :       long im = operate(orb[cnd], gel(G,i), V);
     518     3838366 :       if (flag[im] == 0)
     519             :         /* the image is a new point in the orbit */
     520             :       {
     521      674219 :         orb[++len] = im;
     522      674219 :         flag[im] = 1;
     523             :       }
     524             :     }
     525       21518 :   set_avma(av);
     526       21518 :   return len;
     527             : }
     528             : 
     529             : /* delete the elements in orb2 from orb1, an entry 0 marks the end of the
     530             :  * list, returns the length of orb1 */
     531             : 
     532             : static long
     533        8253 : orbdelete(GEN orb1, GEN orb2)
     534             : {
     535             :   long i, j, len;
     536        8253 :   long l1 = lg(orb1)-1;
     537        8253 :   long l2 = lg(orb2)-1;
     538        8253 :   for (i = 1; i <= l1  &&  orb1[i] != 0; ++i);
     539        8253 :   len = i - 1;
     540       68922 :   for (i = 1; i <= l2  &&  orb2[i] != 0; ++i)
     541             :   {
     542       60669 :     long o2i = orb2[i];
     543       60669 :     for (j = 1; j <= len  &&  orb1[j] != o2i; ++j);
     544             :     /* orb1[j] = orb2[i], hence delete orb1[j] from orb1 */
     545       60669 :     if (j <= len)
     546             :     {
     547       48300 :       orb1[j] = orb1[len];
     548       48300 :       orb1[len--] = 0;
     549             :     }
     550             :   }
     551        8253 :   return len;
     552             : }
     553             : 
     554             : static long
     555        3185 : orbsubtract(GEN Cs, GEN pt, long ipt, long npt, GEN H, GEN V, long *len)
     556             : {
     557        3185 :   pari_sp av = avma;
     558        3185 :   GEN orb = orbit(pt, ipt, npt, H, V);
     559        3185 :   if (len) *len = lg(orb)-1;
     560        3185 :   return gc_long(av, orbdelete(Cs, orb));
     561             : }
     562             : 
     563             : /* Generates the matrix X which has as row per[i] the vector nr. x[i] from the
     564             :  * list V */
     565             : static GEN
     566       18984 : matgen(GEN x, GEN per, GEN V)
     567             : {
     568             :   long i, j;
     569       18984 :   long dim = lg(x)-1;
     570       18984 :   GEN X = cgetg(dim+1,t_MAT);
     571      288862 :   for (i = 1; i <= dim; ++i)
     572             :   {
     573      269878 :     long xi = x[i];
     574      269878 :     GEN Xp = cgetg(dim+1,t_VECSMALL);
     575     4234356 :     for (j = 1; j <= dim; ++j)
     576     3964478 :       Xp[j] = xi > 0? mael(V,xi,j): -mael(V,-xi,j);
     577      269878 :     gel(X,per[i]) = Xp;
     578             :   }
     579       18984 :   return X;
     580             : }
     581             : /* x1 corresponds to an element X1 mapping some vector e on p1, x2 to an
     582             :  * element X2 mapping e on p2 and G is a generator mapping p1 on p2, then
     583             :  * S = X1*G*X2^-1 stabilizes e
     584             :  */
     585             : 
     586             : static GEN
     587        9170 : stabil(GEN x1, GEN x2, GEN per, GEN G, GEN V, ulong p)
     588             : {
     589        9170 :   pari_sp av = avma;
     590             :   long i;
     591             :   GEN x, XG, X2;
     592        9170 :   long dim = lg(x1)-1;
     593        9170 :   x = cgetg(dim+1,t_VECSMALL);
     594      139594 :   for (i = 1; i <= dim; ++i)
     595      130424 :     x[i] = operate(x1[i], G, V);
     596             :   /* XG is the composite mapping of the matrix corresponding to x1 and G */
     597        9170 :   XG = matgen(x, per, V);
     598        9170 :   X2 = matgen(x2, per, V);
     599        9170 :   return gerepileupto(av,zm_divmod(X2,XG,p));
     600             : }
     601             : 
     602             : /* computes the orbit of fp.e[I] under the generators in G->g[I]...G->g[n-1]
     603             :  * and elements stabilizing fp.e[I], has some heuristic break conditions, the
     604             :  * generators in G->g[i] stabilize fp.e[0]...fp.e[i-1] but not fp.e[i],
     605             :  * G->ng[i] is the number of generators in G->g[i], the first G->nsg[i] of
     606             :  * which are elements which are obtained as stabilizer elements in
     607             :  * <G->g[0],...,G->g[i-1]>, G->ord[i] is the orbit length of fp.e[i] under
     608             :  * <G->g[i],...,G->g[n-1]> */
     609             : 
     610             : static void
     611         840 : stab(long I, struct group *G, struct fingerprint *fp, GEN V, ulong p)
     612             : {
     613             :   long len, cnd, tmplen;
     614             :   GEN orb, w, flag, H, Hj, S;
     615             :   long i, j, k, l, im, nH, nHj, fail;
     616             :   long Maxfail, Rest;
     617         840 :   long dim = lg(fp->diag)-1, n = lg(V)-1;
     618             :   /* Some heuristic break conditions for the computation of stabilizer elements:
     619             :      it would be too expensive to calculate all the stabilizer generators,
     620             :      which are obtained from the orbit, since this is highly redundant, on the
     621             :      other hand every new generator which enlarges the group is much cheaper
     622             :      than one obtained from the backtrack, after Maxfail subsequent stabilizer
     623             :      elements, that do not enlarge the group, Rest more elements are calculated
     624             :      even if they leave the group unchanged, since it turned out that this is
     625             :      often useful in the following steps, increasing the parameters will
     626             :      possibly decrease the number of generators for the group, but will
     627             :      increase the running time, there is no magic behind this heuristic, tuning
     628             :      might be appropriate */
     629             : 
     630        7658 :   for (Rest = 0, i = I; i <= dim; ++i)
     631        6818 :     if (fp->diag[i] > 1  &&  G->ord[i] < fp->diag[i])
     632        3633 :       ++Rest;
     633       12348 :   for (Maxfail = Rest, i = 1; i <= dim; ++i)
     634       11508 :     if (fp->diag[i] > 1)
     635        9324 :       ++Maxfail;
     636        7658 :   for (nH = 0, i = I; i <= dim; ++i)
     637        6818 :     nH += G->ng[i];
     638             :   /* H are the generators of the group in which the stabilizer is computed */
     639         840 :   H = cgetg(nH+1,t_MAT);
     640         840 :   Hj = cgetg(nH+2,t_MAT);
     641         840 :   k = 0;
     642        7658 :   for (i = I; i <= dim; ++i)
     643       16310 :     for (j = 1; j <= G->ng[i]; ++j)
     644        9492 :       gel(H,++k) = gmael(G->g,i,j);
     645             :   /* in w[V.n+i] an element is stored that maps fp.e[I] on v[i] */
     646         840 :   w = cgetg(2*n+2,t_VEC);
     647             :   /* orb contains the orbit of fp.e[I] */
     648         840 :   orb = zero_Flv(2*n);
     649             :   /* if flag[i + V.n] = 1, then the point i is already in the orbit */
     650         840 :   flag = zero_Flv(2*n+1);
     651         840 :   orb[1] = fp->e[I];
     652         840 :   flag[orb[1]+n+1] = 1;
     653         840 :   gel(w,orb[1]+n+1) = cgetg(dim+1,t_VECSMALL);
     654       12348 :   for (i = 1; i <= dim; ++i)
     655       11508 :     mael(w,orb[1]+n+1,i) = fp->e[i];
     656         840 :   cnd = 1;
     657         840 :   len = 1;
     658             :   /* fail is the number of successive failures */
     659         840 :   fail = 0;
     660        5845 :   while (cnd <= len && fail < Maxfail+Rest)
     661             :   {
     662       40348 :     for (i = 1; i <= nH  &&  fail < Maxfail+Rest; ++i)
     663             :     {
     664       36183 :       if (fail >= Maxfail)
     665             :         /* there have already been Maxfail successive failures, now a random
     666             :            generator is applied to a random point of the orbit to get Rest more
     667             :            stabilizer elements */
     668             :       {
     669        5362 :         cnd = 1+(long)random_Fl(len);
     670        5362 :         i = 1+(long)random_Fl(nH);
     671             :       }
     672       36183 :       im = operate(orb[cnd], gel(H,i), V);
     673       36183 :       if (flag[im+n+1] == 0)
     674             :         /* a new element is found, appended to the orbit and an element mapping
     675             :            fp.e[I] to im is stored in w[im+V.n] */
     676             :       {
     677             :         GEN wim;
     678       13321 :         orb[++len] = im;
     679       13321 :         flag[im+n+1] = 1;
     680       13321 :         wim = cgetg(dim+1,t_VECSMALL);
     681       13321 :         gel(w,im+n+1) = wim;
     682      178955 :         for (j = 1; j <= dim; ++j)
     683      165634 :           wim[j] = operate(mael(w,orb[cnd]+n+1,j), gel(H,i), V);
     684             :       }
     685             :       else
     686             :         /* the image was already in the orbit */
     687             :       {
     688             :         /* j is the first index where the images of the old and the new element
     689             :            mapping e[I] on im differ */
     690      121450 :         for (j = I; j <= dim; j++)
     691      116354 :           if (operate(mael(w,orb[cnd]+n+1,j), gel(H,i), V) != mael(w,im+n+1,j))
     692       17766 :             break;
     693       40628 :         if (j <= dim  &&
     694       28056 :             (G->ord[j] < fp->diag[j] || fail >= Maxfail))
     695        9170 :         {
     696        9170 :           GEN wo = gel(w,orb[cnd]+n+1);
     697             :       /* new stabilizer element S = w[orb[cnd]+V.n] * H[i] * (w[im+V.n])^-1 */
     698        9170 :           S = stabil(wo, gel(w,im+n+1), fp->per, gel(H,i), V, p);
     699        9170 :           gel(Hj,1) = S;
     700        9170 :           nHj = 1;
     701       96306 :           for (k = j; k <= dim; ++k)
     702      145551 :             for (l = 1; l <= G->ng[k]; ++l)
     703       58415 :               gel(Hj,++nHj) = gmael(G->g,k,l);
     704        9170 :           tmplen = orbitlen(fp->e[j], fp->diag[j], Hj, nHj, V);
     705        9170 :           if (tmplen > G->ord[j]  ||  fail >= Maxfail)
     706             :             /* the new stabilizer element S either enlarges the orbit of e[j]
     707             :                or it is one of the additional elements after MAXFAIL failures */
     708        4032 :           {
     709             :             GEN Ggj;
     710        4032 :             G->ord[j] = tmplen;
     711        4032 :             G->ng[j]++;
     712        4032 :             G->nsg[j]++;
     713             :             /* allocate memory for the new generator */
     714        4032 :             gel(G->g,j) = vec_lengthen(gel(G->g,j),G->ng[j]);
     715        4032 :             Ggj = gel(G->g,j);
     716             :             /* the new generator is inserted as stabilizer element
     717             :                nr. nsg[j]-1 */
     718        4032 :             for (k = G->ng[j]; k > G->nsg[j]; --k)
     719           0 :               gel(Ggj,k) = gel(Ggj,k-1);
     720        4032 :             gel(Ggj,G->nsg[j]) = S;
     721        4032 :             nH++;
     722        4032 :             H = vec_lengthen(H, nH);
     723        4032 :             Hj = vec_lengthen(Hj, nH+1);
     724             :             /* the new generator is appended to H */
     725        4032 :             gel(H,nH) = gel(Ggj,G->nsg[j]);
     726             :             /* the number of failures is reset to 0 */
     727        4032 :             if (fail < Maxfail)
     728         980 :               fail = 0;
     729             :             else
     730        3052 :               ++fail;
     731             :           }
     732             :           else
     733             :             /*the new stabilizer element S does not enlarge the orbit of e[j]*/
     734        5138 :             ++fail;
     735             :         }
     736       13692 :         else if ((j < dim  &&  fail < Maxfail)  ||
     737         126 :             (j == dim  &&  fail >= Maxfail))
     738        8470 :           ++fail;
     739             :         /* if S is the identity and fail < Maxfail, nothing is done */
     740             :       }
     741             :     }
     742        4165 :     if (fail < Maxfail)
     743        3374 :       ++cnd;
     744             :   }
     745         840 : }
     746             : 
     747             : /* tests, whether x[1],...,x[I-1] is a partial automorphism, using scalar
     748             :  * product combinations and Bacher-polynomials depending on the chosen options,
     749             :  * puts the candidates for x[I] (i.e. the vectors vec such that the scalar
     750             :  * products of x[1],...,x[I-1],vec are correct) on CI, returns their number
     751             :  * (should be fp.diag[I]) */
     752             : 
     753             : static long
     754         168 : qfisom_candidates_novec(GEN CI, long I, GEN x, struct qfauto *qf,
     755             :        struct qfauto *qff, struct fingerprint *fp)
     756             : {
     757         168 :   pari_sp av = avma;
     758             :   long i, j, k, okp, okm, nr, fail;
     759             :   GEN vec;
     760         168 :   GEN F =qf->F, V=qff->V, W=qff->W, v=qff->v;
     761         168 :   long n = lg(V)-1, f = lg(F)-1;
     762         168 :   vec = cgetg(I,t_VECSMALL);
     763             :   /* CI is the list for the candidates */
     764       34818 :   for (i = 1; i <= fp->diag[I]; ++i)
     765       34650 :     CI[i] = 0;
     766         168 :   nr = 0;
     767         168 :   fail = 0;
     768      173957 :   for (j = 1; j <= n  &&  fail == 0; ++j)
     769             :   {
     770      173789 :     GEN Vj = gel(V,j), Wj = gel(W, j);
     771      173789 :     okp = 0;
     772      173789 :     okm = 0;
     773      347578 :     for (i = 1; i <= f; ++i)
     774             :     {
     775      173789 :       GEN FAiI = gmael(F,i,fp->per[I]);
     776      173789 :       GEN FFvi = gel(v,i);
     777             :       /* vec is the vector of scalar products of V.v[j] with the first I base
     778             :          vectors x[0]...x[I-1] */
     779      173789 :       for (k = 1; k < I; ++k)
     780             :       {
     781           0 :         long xk = x[k];
     782           0 :         if (xk > 0)
     783           0 :           vec[k] = zv_dotproduct(Vj, gel(FFvi,xk));
     784             :         else
     785           0 :           vec[k] = -zv_dotproduct(Vj, gel(FFvi,-xk));
     786             :       }
     787      173789 :       for (k = 1; k < I  &&  vec[k] == FAiI[fp->per[k]]; ++k);
     788      173789 :       if (k == I  &&  Wj[i] == FAiI[fp->per[I]])
     789             :         /* V.v[j] is a candidate for x[I] with respect to the form F.A[i] */
     790       17325 :         ++okp;
     791      173789 :       for (k = 1; k < I  &&  vec[k] == -FAiI[fp->per[k]]; ++k);
     792      173789 :       if (k == I  &&  Wj[i] == FAiI[fp->per[I]])
     793             :         /* -V.v[j] is a candidate for x[I] with respect to the form F.A[i] */
     794       17325 :         ++okm;
     795             :     }
     796      173789 :     if (okp == f)
     797             :       /* V.v[j] is a candidate for x[I] */
     798             :     {
     799       17325 :       if (nr < fp->diag[I])
     800       17325 :         CI[++nr] = j;
     801             :       else
     802             :         /* there are too many candidates */
     803           0 :         fail = 1;
     804             :     }
     805      173789 :     if (okm == f)
     806             :       /* -V.v[j] is a candidate for x[I] */
     807             :     {
     808       17325 :       if (nr < fp->diag[I])
     809       17325 :         CI[++nr] = -j;
     810             :       else
     811             :         /* there are too many candidates */
     812           0 :         fail = 1;
     813             :     }
     814             :   }
     815         168 :   if (fail == 1)
     816           0 :     nr = 0;
     817         168 :   set_avma(av);
     818         168 :   return nr;
     819             : }
     820             : 
     821             : static long
     822       13132 : qfisom_candidates(GEN CI, long I, GEN x, struct qfauto *qf,
     823             :       struct qfauto *qff, struct fingerprint *fp, struct qfcand *qfcand)
     824             : {
     825       13132 :   pari_sp av = avma;
     826             :   long i, j, k, okp, okm, nr, fail;
     827             :   GEN vec;
     828             :   GEN xvec, xbase, Fxbase, scpvec;
     829             :   long vj, rank, num;
     830             :   GEN com, list, trans, ccoef, cF;
     831       13132 :   GEN F =qf->F, V=qff->V, W=qff->W, v=qff->v, FF= qff->F;
     832       13132 :   long dim = qf->dim, n = lg(V)-1, f = lg(F)-1;
     833             :   long nc;
     834       13132 :   long DEP = qfcand->cdep, len = f * DEP;
     835       13132 :   if (I >= 2  &&  I <= lg(qfcand->bacher_pol))
     836             :   {
     837         154 :     long t = labs(x[I-1]);
     838         154 :     GEN bpolI = gel(qfcand->bacher_pol,I-1);
     839         154 :     if (bachcomp(bpolI, t, V, W, gel(v,1)) == 0) return 0;
     840             :   }
     841       13090 :   if (I==1 || DEP ==0)
     842         168 :     return qfisom_candidates_novec(CI,I,x,qf,qff,fp);
     843       12922 :   vec = cgetg(I,t_VECSMALL);
     844       12922 :   scpvec = cgetg(len+1,t_VECSMALL);
     845       12922 :   com = gel(qfcand->comb,I-1);
     846       12922 :   list=gel(com,1); trans = gel(com,2); ccoef = gel(com,3); cF=gel(com,4);
     847       12922 :   rank = lg(trans)-1;
     848       12922 :   nc = lg(list)-2;
     849             :   /* xvec is the list of vector sums which are computed with respect to the
     850             :      partial basis in x */
     851       12922 :   xvec = zero_Flm_copy(dim,nc+1);
     852             :   /* xbase should be a basis for the lattice generated by the vectors in xvec,
     853             :      it is obtained via the transformation matrix comb[I-1].trans */
     854       12922 :   xbase = cgetg(rank+1,t_MAT);
     855       51716 :   for (i = 1; i <= rank; ++i)
     856       38794 :     gel(xbase,i) = cgetg(dim+1,t_VECSMALL);
     857             :   /* Fxbase is the product of a form F with the base xbase */
     858       12922 :   Fxbase = cgetg(rank+1,t_MAT);
     859       51716 :   for (i = 1; i <= rank; ++i)
     860       38794 :     gel(Fxbase,i) = cgetg(dim+1,t_VECSMALL);
     861             :   /* CI is the list for the candidates */
     862       87066 :   for (i = 1; i <= fp->diag[I]; ++i)
     863       74144 :     CI[i] = 0;
     864       12922 :   nr = 0;
     865       12922 :   fail = 0;
     866    14125776 :   for (j = 1; j <= n  &&  fail == 0; ++j)
     867             :   {
     868             :     long sign;
     869    14112854 :     GEN Vj = gel(V,j), Wj = gel(W, j);
     870    14112854 :     okp = 0;
     871    14112854 :     okm = 0;
     872    37653980 :     for (k = 1; k <= len; ++k)
     873    23541126 :       scpvec[k] = 0;
     874    28225708 :     for (i = 1; i <= f; ++i)
     875             :     {
     876    14112854 :       GEN FAiI = gmael(F,i,fp->per[I]);
     877    14112854 :       GEN FFvi = gel(v,i);
     878             :       /* vec is the vector of scalar products of V.v[j] with the first I base
     879             :          vectors x[0]...x[I-1] */
     880   127677354 :       for (k = 1; k < I; ++k)
     881             :       {
     882   113564500 :         long xk = x[k];
     883   113564500 :         if (xk > 0)
     884    70179998 :           vec[k] = zv_dotproduct(Vj, gel(FFvi,xk));
     885             :         else
     886    43384502 :           vec[k] = -zv_dotproduct(Vj, gel(FFvi,-xk));
     887             :       }
     888    14112854 :       for (k = 1; k < I  &&  vec[k] == FAiI[fp->per[k]]; ++k);
     889    14112854 :       if (k == I  &&  Wj[i] == FAiI[fp->per[I]])
     890             :         /* V.v[j] is a candidate for x[I] with respect to the form F.A[i] */
     891       32529 :         ++okp;
     892    14112854 :       for (k = 1; k < I  &&  vec[k] == -FAiI[fp->per[k]]; ++k);
     893    14112854 :       if (k == I  &&  Wj[i] == FAiI[fp->per[I]])
     894             :         /* -V.v[j] is a candidate for x[I] with respect to the form F.A[i] */
     895       22407 :         ++okm;
     896    37224152 :       for (k = I-1; k >= 1  &&  k > I-1-DEP; --k)
     897    23111298 :         scpvec[(i-1)*DEP+I-k] = vec[k];
     898             :     }
     899             :     /* check, whether the scalar product combination scpvec is contained in the
     900             :        list comb[I-1].list */
     901    14112854 :     if (!zv_equal0(scpvec))
     902             :     {
     903     9991996 :       sign = zv_canon(scpvec);
     904     9991996 :       num = vecvecsmall_search(list,scpvec,0);
     905     9991996 :       if (!num)
     906             :         /* scpvec is not found, hence x[0]...x[I-1] is not a partial
     907             :            automorphism */
     908           0 :         fail = 1;
     909             :       else
     910             :         /* scpvec is found and the vector is added to the corresponding
     911             :            vector sum */
     912             :       {
     913     9991996 :         GEN xnum = gel(xvec,num);
     914   161675360 :         for (k = 1; k <= dim; ++k)
     915   151683364 :           xnum[k] += sign * Vj[k];
     916             :       }
     917             :     }
     918    14112854 :     if (okp == f)
     919             :       /* V.v[j] is a candidate for x[I] */
     920             :     {
     921       32529 :       if (nr < fp->diag[I])
     922       32494 :         CI[++nr] = j;
     923             :       else
     924             :         /* there are too many candidates */
     925          35 :         fail = 1;
     926             :     }
     927    14112854 :     if (okm == f)
     928             :       /* -V.v[j] is a candidate for x[I] */
     929             :     {
     930       22407 :       if (nr < fp->diag[I])
     931       22358 :         CI[++nr] = -j;
     932             :       else
     933             :         /* there are too many candidates */
     934          49 :         fail = 1;
     935             :     }
     936             :   }
     937       12922 :   if (fail == 1)
     938          84 :     nr = 0;
     939       12922 :   if (nr == fp->diag[I])
     940             :   {
     941             :     /* compute the basis of the lattice generated by the vectors in xvec via
     942             :        the transformation matrix comb[I-1].trans */
     943       35861 :     for (i = 1; i <= rank; ++i)
     944             :     {
     945       27027 :       GEN comtri = gel(trans,i);
     946      417823 :       for (j = 1; j <= dim; ++j)
     947             :       {
     948      390796 :         long xbij = 0;
     949     5908140 :         for (k = 1; k <= nc+1; ++k)
     950     5517344 :           xbij += comtri[k] * mael(xvec,k,j);
     951      390796 :         mael(xbase,i,j) = xbij;
     952             :       }
     953             :     }
     954             :     /* check, whether the base xbase has the right scalar products */
     955       17668 :     for (i = 1; i <= f &&  nr > 0; ++i)
     956             :     {
     957       35861 :       for (j = 1; j <= rank; ++j)
     958      417823 :         for (k = 1; k <= dim; ++k)
     959      390796 :           mael(Fxbase,j,k) = zv_dotproduct(gmael(FF,i,k), gel(xbase,j));
     960       35861 :       for (j = 1; j <= rank  &&  nr > 0; ++j)
     961       94598 :         for (k = 1; k <= j  &&  nr > 0; ++k)
     962             :         {
     963       67571 :           if (zv_dotproduct(gel(xbase,j), gel(Fxbase,k)) != mael3(cF,i,j,k))
     964             :             /* a scalar product is wrong */
     965           0 :             nr = 0;
     966             :         }
     967             :     }
     968      103544 :     for (i = 1; i <= nc+1  &&  nr > 0; ++i)
     969             :     {
     970       94710 :       GEN comcoi = gel(ccoef,i);
     971     1462258 :       for (j = 1; j <= dim && nr > 0; ++j)
     972             :       {
     973     1367548 :         vj = 0;
     974     6884892 :         for (k = 1; k <= rank; ++k)
     975     5517344 :           vj += comcoi[k] * mael(xbase,k,j);
     976     1367548 :         if (vj != mael(xvec,i,j))
     977             :           /* an entry is wrong */
     978           0 :           nr = 0;
     979             :       }
     980             :     }
     981             :   }
     982       12922 :   set_avma(av);
     983       12922 :   return nr;
     984             : }
     985             : 
     986             : static long
     987        7294 : aut(long step, GEN x, GEN C, struct group *G, struct qfauto *qf,
     988             :                              struct fingerprint *fp, struct qfcand *cand)
     989             : {
     990        7294 :   long dim = qf->dim;
     991        7294 :   GEN orb = cgetg(2,t_VECSMALL);
     992        7294 :   while (mael(C,step,1) != 0)
     993             :   {
     994       10703 :     if (step < dim)
     995             :     {
     996             :       long nbc;
     997             :       /* choose the image of the base-vector nr. step */
     998       10164 :       gel(x,step) = gmael(C,step,1);
     999             :       /* check, whether x[0]...x[step] is a partial automorphism and compute
    1000             :          the candidates for x[step+1] */
    1001       10164 :       nbc = qfisom_candidates(gel(C,step+1), step+1, x, qf, qf, fp, cand);
    1002       10164 :       if (nbc == fp->diag[step+1])
    1003             :         /* go deeper into the recursion */
    1004        6622 :         if (aut(step+1, x, C, G, qf, fp, cand))
    1005        5096 :           return 1;
    1006        5068 :       orb[1] = x[step];
    1007             :       /* delete the chosen vector from the list of candidates */
    1008        5068 :       (void)orbdelete(gel(C,step), orb);
    1009             :     }
    1010             :     else
    1011             :     {
    1012             :       /* a new automorphism is found */
    1013         539 :       gel(x,dim) = gmael(C,dim,1);
    1014         539 :       return 1;
    1015             :     }
    1016             :   }
    1017        1659 :   return 0;
    1018             : }
    1019             : 
    1020             : /* search new automorphisms until on all levels representatives for all orbits
    1021             :  * have been tested */
    1022             : 
    1023             : static void
    1024         112 : autom(struct group *G, struct qfauto *qf, struct fingerprint *fp,
    1025             :                                           struct qfcand *cand)
    1026             : {
    1027             :   long i, j, step, im, nC, nH, found, tries;
    1028             :   GEN  x, bad, H;
    1029             :   long nbad;
    1030         112 :   GEN V = qf->V;
    1031         112 :   long dim = qf->dim, n = lg(V)-1;
    1032         112 :   long STAB = 1;
    1033         112 :   GEN C = cgetg(dim+1,t_VEC);
    1034             :   /* C[i] is the list of candidates for the image of the i-th base-vector */
    1035        1400 :   for (i = 1; i <= dim; ++i)
    1036        1288 :     gel(C,i) = cgetg(fp->diag[i]+1, t_VECSMALL);
    1037             :   /* x is the new base i.e. x[i] is the index in V.v of the i-th base-vector */
    1038         112 :   x = cgetg(dim+1, t_VECSMALL);
    1039        1400 :   for (step = STAB; step <= dim; ++step)
    1040             :   {
    1041        1288 :     if(DEBUGLEVEL) err_printf("QFAuto: Step %ld/%ld\n",step,dim);
    1042        1288 :     nH = 0;
    1043       10598 :     for (nH = 0, i = step; i <= dim; ++i)
    1044        9310 :       nH += G->ng[i];
    1045        1288 :     H = cgetg(nH+1,t_VEC);
    1046       10598 :     for (nH = 0, i = step; i <= dim; ++i)
    1047       21798 :       for (j = 1; j <= G->ng[i]; ++j)
    1048       12488 :         gel(H,++nH) = gmael(G->g,i,j);
    1049        1288 :     bad = zero_Flv(2*n);
    1050        1288 :     nbad = 0;
    1051             :     /* the first step base-vectors are fixed */
    1052        9310 :     for (i = 1; i < step; ++i)
    1053        8022 :       x[i] = fp->e[i];
    1054             :     /* compute the candidates for x[step] */
    1055        1288 :     if (fp->diag[step] > 1)
    1056             :       /* if fp.diag[step] > 1 compute the candidates for x[step] */
    1057        1043 :       nC = qfisom_candidates(gel(C,step), step, x, qf, qf, fp, cand);
    1058             :     else
    1059             :       /* if fp.diag[step] == 1, fp.e[step] is the only candidate */
    1060             :     {
    1061         245 :       mael(C,step,1) = fp->e[step];
    1062         245 :       nC = 1;
    1063             :     }
    1064             :     /* delete the orbit of the step-th base-vector from the candidates */
    1065        1288 :     nC = orbsubtract(gel(C,step), fp->e, step-1, 1, H, V, &(G->ord[step]));
    1066        3752 :     while (nC > 0  &&  (im = mael(C,step,1)) != 0)
    1067             :     {
    1068        1176 :       found = 0;
    1069             :       /* tries vector V.v[im] as image of the step-th base-vector */
    1070        1176 :       x[step] = im;
    1071        1176 :       if (step < dim)
    1072             :       {
    1073             :         long nbc;
    1074             :         /* check, whether x[0]...x[step] is a partial basis and compute the
    1075             :            candidates for x[step+1] */
    1076        1127 :         nbc = qfisom_candidates(gel(C,step+1), step+1, x, qf, qf, fp, cand);
    1077        1127 :         if (nbc == fp->diag[step+1])
    1078             :           /* go into the recursion */
    1079         672 :           found = aut(step+1, x, C, G, qf, fp, cand);
    1080             :         else
    1081         455 :           found = 0;
    1082             :       }
    1083             :       else
    1084          49 :         found = 1;
    1085        1176 :       if (found == 0)
    1086             :         /* x[0]...x[step] can not be continued to an automorphism */
    1087             :       {
    1088             :         /* delete the orbit of im from the candidates for x[step] */
    1089         588 :         nC = orbsubtract(gel(C,step),mkvecsmall(im), 0, 1, H, V, NULL);
    1090         588 :         bad[++nbad] = im;
    1091             :       }
    1092             :       else
    1093             :         /* a new generator has been found */
    1094             :       {
    1095             :         GEN Gstep;
    1096         588 :         ++G->ng[step];
    1097             :         /* append the new generator to G->g[step] */
    1098         588 :         Gstep = vec_lengthen(gel(G->g,step),G->ng[step]);
    1099         588 :         gel(Gstep,G->ng[step]) = matgen(x, fp->per, V);
    1100         588 :         gel(G->g,step) = Gstep;
    1101         588 :         ++nH;
    1102         588 :         H = cgetg(nH+1, t_VEC);
    1103        6811 :         for (nH = 0, i = step; i <= dim; ++i)
    1104       11116 :           for (j = 1; j <= G->ng[i]; ++j)
    1105        4893 :             gel(H,++nH) = gmael(G->g,i,j);
    1106         588 :         nC = orbsubtract(gel(C,step), fp->e, step-1, 1, H, V, &(G->ord[step]));
    1107         588 :         nC = orbsubtract(gel(C,step), bad, 0, nbad, H, V, NULL);
    1108             :       }
    1109             :     }
    1110             :     /* test, whether on step STAB some generators may be omitted */
    1111        1288 :     if (step == STAB)
    1112         595 :       for (tries = G->nsg[step]; tries <= G->ng[step]; ++tries)
    1113             :       {
    1114         483 :         nH = 0;
    1115         945 :         for (j = 1; j < tries; ++j)
    1116         462 :           gel(H,++nH) = gmael(G->g,step,j);
    1117         945 :         for (j = tries+1; j < G->ng[step]; ++j)
    1118         462 :           gel(H,++nH) = gmael(G->g,step,j);
    1119        5754 :         for (i = step+1; i <= dim; ++i)
    1120        5271 :           for (j = 1; j <= G->ng[i]; ++j)
    1121           0 :             gel(H,++nH) = gmael(G->g,i,j);
    1122         483 :         if (orbitlen(fp->e[step], G->ord[step], H, nH, V) == G->ord[step])
    1123             :           /* the generator g[step][tries] can be omitted */
    1124             :         {
    1125           0 :           G->ng[step]--;
    1126           0 :           for (i = tries; i < G->ng[step]; ++i)
    1127           0 :             gmael(G->g,step,i) = gmael(G->g,step,i+1);
    1128           0 :           tries--;
    1129             :         }
    1130             :       }
    1131        1288 :     if (step < dim  &&  G->ord[step] > 1)
    1132             :       /* calculate stabilizer elements fixing the basis-vectors
    1133             :          fp.e[0]...fp.e[step] */
    1134         840 :       stab(step, G, fp, V, qf->p);
    1135             :   }
    1136         112 : }
    1137             : 
    1138             : #define MAXENTRY (1L<<((BITS_IN_LONG-2)>>1))
    1139             : #define MAXNORM (1L<<(BITS_IN_LONG-2))
    1140             : 
    1141             : static long
    1142         406 : zm_maxdiag(GEN A)
    1143             : {
    1144         406 :   long dim = lg(A)-1;
    1145         406 :   long max = coeff(A,1,1);
    1146             :   long i;
    1147        4739 :   for (i = 2; i <= dim; ++i)
    1148        4333 :     if (coeff(A,i,i) > max)
    1149         189 :       max = coeff(A,i,i);
    1150         406 :   return max;
    1151             : }
    1152             : 
    1153             : static GEN
    1154         224 : init_qfauto(GEN F, GEN U, long max, struct qfauto *qf, GEN norm, GEN minvec)
    1155             : {
    1156             :   long i, j, k;
    1157             :   GEN W, v;
    1158         224 :   GEN M = minvec? minvec: gel(minim(zm_to_ZM(gel(F,1)), stoi(max), NULL), 3);
    1159         224 :   GEN V = ZM_to_zm_canon(M);
    1160         224 :   long n = lg(V)-1, f = lg(F)-1, dim = lg(gel(F,1))-1;
    1161      245798 :   for (i = 1; i <= n; ++i)
    1162             :   {
    1163      245574 :     GEN Vi = gel(V,i);
    1164     4107684 :     for (k = 1; k <= dim; ++k)
    1165             :     {
    1166     3862110 :       long l = labs(Vi[k]);
    1167     3862110 :       if (l > max)
    1168           0 :         max = l;
    1169             :     }
    1170             :   }
    1171         224 :   if (max > MAXENTRY) pari_err_OVERFLOW("qfisom [lattice too large]");
    1172         224 :   qf->p = unextprime(2*max+1);
    1173         224 :   V = vecvecsmall_sort_uniq(V);
    1174         224 :   if (!norm)
    1175             :   {
    1176         168 :     norm = cgetg(dim+1,t_VEC);
    1177        2121 :     for (i = 1; i <= dim; ++i)
    1178             :     {
    1179        1953 :       GEN Ni = cgetg(f+1,t_VECSMALL);
    1180        3906 :       for (k = 1; k <= f; ++k)
    1181        1953 :         Ni[k] = mael3(F,k,i,i);
    1182        1953 :       gel(norm,i) = Ni;
    1183             :     }
    1184         168 :     norm = vecvecsmall_sort_uniq(norm);
    1185             :   }
    1186         224 :   W = checkvecs(V, F, norm);
    1187         224 :   v = cgetg(f+1,t_VEC);
    1188             :   /* the product of the maximal entry in the short vectors with the maximal
    1189             :      entry in v[i] should not exceed MAXNORM to avoid overflow */
    1190         224 :   max = MAXNORM / max;
    1191         448 :   for (i = 1; i <= f; ++i)
    1192             :   {
    1193         224 :     GEN Fi = gel(F,i), vi;
    1194         224 :     vi = cgetg(n+1,t_MAT);
    1195         224 :     gel(v,i) = vi;
    1196      245798 :     for (j = 1; j <= n; ++j)
    1197             :     {
    1198      245574 :       GEN Vj = gel(V,j);
    1199      245574 :       GEN vij = cgetg(dim+1, t_VECSMALL);
    1200      245574 :       gel(vi,j) = vij;
    1201     4107684 :       for (k = 1; k <= dim; ++k)
    1202             :       {
    1203     3862110 :         vij[k] = zv_dotproduct(gel(Fi,k), Vj);
    1204     3862110 :         if (labs(vij[k]) > max) pari_err_OVERFLOW("qfisom [lattice too large]");
    1205             :       }
    1206             :     }
    1207             :   }
    1208         224 :   qf->dim = dim; qf->F = F; qf->V = V; qf->W = W; qf->v = v; qf->U = U;
    1209         224 :   return norm;
    1210             : }
    1211             : 
    1212             : static void
    1213         112 : init_qfgroup(struct group *G, struct fingerprint *fp, struct qfauto *qf)
    1214             : {
    1215         112 :   GEN H, M, V = qf->V;
    1216             :   long nH;
    1217             :   long i, j, k;
    1218         112 :   long dim = qf->dim;
    1219         112 :   G->ng  = zero_Flv(dim+1);
    1220         112 :   G->nsg = zero_Flv(dim+1);
    1221         112 :   G->ord = cgetg(dim+1,t_VECSMALL);
    1222         112 :   G->g = cgetg(dim+1,t_VEC);
    1223        1400 :   for (i = 1; i <= dim; ++i)
    1224        1288 :     gel(G->g,i) = mkvec(gen_0);
    1225         112 :   M = matid_Flm(dim);
    1226         112 :   gmael(G->g,1,1) = M;
    1227         112 :   G->ng[1] = 1;
    1228             :   /* -Id is always an automorphism */
    1229        1400 :   for (i = 1; i <= dim; ++i)
    1230        1288 :     mael(M,i,i) = -1;
    1231         112 :   nH = 0;
    1232        1400 :   for (i = 1; i <= dim; ++i)
    1233        1288 :     nH += G->ng[i];
    1234         112 :   H = cgetg(nH+1,t_MAT);
    1235             :   /* calculate the orbit lengths under the automorphisms known so far */
    1236        1400 :   for (i = 1; i <= dim; ++i)
    1237             :   {
    1238        1288 :     if (G->ng[i] > 0)
    1239             :     {
    1240         112 :       nH = 0;
    1241        1400 :       for (j = i; j <= dim; ++j)
    1242        1400 :         for (k = 1; k <= G->ng[j]; ++k)
    1243         112 :           gel(H,++nH) = gmael(G->g,j,k);
    1244         112 :       G->ord[i] = orbitlen(fp->e[i], fp->diag[i], H, nH, V);
    1245             :     }
    1246             :     else
    1247        1176 :       G->ord[i] = 1;
    1248             :   }
    1249         112 : }
    1250             : 
    1251             : /* calculates the scalar products of the vector w with the base vectors
    1252             :  * v[b[I]] down to v[b[I-dep+1]] with respect to all invariant forms and puts
    1253             :  * them on scpvec  */
    1254             : static GEN
    1255     7060298 : scpvector(GEN w, GEN b, long I, long dep, GEN v)
    1256             : {
    1257     7060298 :   long  i, j, n = lg(v)-1;
    1258     7060298 :   GEN scpvec = zero_Flv(dep*n);
    1259    18788658 :   for (i = I; i >= 1  &&  i > I-dep; --i)
    1260             :   {
    1261    11728360 :     long bi = b[i];
    1262    11728360 :     if (bi > 0)
    1263    23456720 :       for (j = 1; j <= n; ++j)
    1264    11728360 :         scpvec[1+(j-1)*dep + I-i] = zv_dotproduct(w, gmael(v,j,bi));
    1265             :     else
    1266           0 :       for (j = 1; j <= n; ++j)
    1267           0 :         scpvec[1+(j-1)*dep + I-i] = -zv_dotproduct(w, gmael(v,j,-bi));
    1268             :   }
    1269     7060298 :   return scpvec;
    1270             : }
    1271             : 
    1272             : /* computes the list of scalar product combinations of the vectors
    1273             :  * in V.v with the basis-vectors in b */
    1274             : static GEN
    1275        5243 : scpvecs(GEN *pt_vec, long I, GEN b, long dep, struct qfauto *qf)
    1276             : {
    1277             :   long  i, j, nr, sign;
    1278             :   GEN list, vec;
    1279             :   GEN vecnr;
    1280        5243 :   GEN V = qf->V, F = qf->F, v = qf->v;
    1281        5243 :   long n = lg(V)-1;
    1282        5243 :   long dim = lg(gel(F,1))-1;
    1283        5243 :   long len = (lg(F)-1)*dep;
    1284             :   /* the first vector in the list is the 0-vector and is not counted */
    1285        5243 :   list = mkmat(zero_Flv(len));
    1286        5243 :   vec  = mkmat(zero_Flv(dim));
    1287     7065541 :   for (j = 1; j <= n; ++j)
    1288             :   {
    1289     7060298 :     GEN Vvj = gel(V,j);
    1290     7060298 :     GEN scpvec = scpvector(Vvj, b, I, dep, v);
    1291     7060298 :     if (zv_equal0(scpvec))
    1292     2202655 :       nr = -1;
    1293             :     else
    1294             :     {
    1295     4857643 :       sign = zv_canon(scpvec);
    1296     4857643 :       nr = vecvecsmall_search(list,scpvec,0);
    1297             :     }
    1298             :     /* scpvec is already in list */
    1299     7060298 :     if (nr > 0)
    1300             :     {
    1301     4808874 :       vecnr = gel(vec,nr);
    1302    80345370 :       for (i = 1; i <= dim; ++i)
    1303    75536496 :         vecnr[i] += sign * Vvj[i];
    1304             :     }
    1305             :     /* scpvec is a new scalar product combination */
    1306     2251424 :     else if (nr==0)
    1307             :     {
    1308       48769 :       nr = vecvecsmall_search(list,scpvec,1);
    1309       48769 :       list=vec_insert(list,nr,scpvec);
    1310       48769 :       vec=vec_insert(vec,nr,sign < 0 ? zv_neg(Vvj) : zv_copy(Vvj));
    1311             :     }
    1312             :   }
    1313        5243 :   settyp(list,t_MAT);
    1314        5243 :   settyp(vec,t_MAT);
    1315        5243 :   *pt_vec = vec;
    1316        5243 :   return list;
    1317             : }
    1318             : 
    1319             : /* com->F[i] is the Gram-matrix of the basis b with respect to F.A[i] */
    1320             : static GEN
    1321        5117 : scpforms(GEN b, struct qfauto *qf)
    1322             : {
    1323             :   long i, j, k;
    1324        5117 :   GEN F = qf->F;
    1325        5117 :   long n = lg(F)-1, dim = lg(gel(F,1))-1;
    1326        5117 :   long nb = lg(b)-1;
    1327        5117 :   GEN gram = cgetg(n+1, t_VEC);
    1328             :   /* Fbi is the list of products of F.A[i] with the vectors in b */
    1329        5117 :   GEN Fbi = cgetg(nb+1, t_MAT);
    1330       18207 :   for (j = 1; j <= nb; ++j)
    1331       13090 :     gel(Fbi, j) = cgetg(dim+1, t_VECSMALL);
    1332       10234 :   for (i = 1; i <= n; ++i)
    1333             :   {
    1334        5117 :     GEN FAi = gel(F,i);
    1335        5117 :     gel(gram, i) = cgetg(nb+1, t_MAT);
    1336       18207 :     for (j = 1; j <= nb; ++j)
    1337      200669 :       for (k = 1; k <= dim; ++k)
    1338      187579 :         mael(Fbi,j,k) = zv_dotproduct(gel(FAi,k), gel(b,j));
    1339       18207 :     for (j = 1; j <= nb; ++j)
    1340             :     {
    1341       13090 :       GEN comFij = cgetg(nb+1, t_VECSMALL);
    1342       59318 :       for (k = 1; k <= nb; ++k)
    1343       46228 :         comFij[k] = zv_dotproduct(gel(b,j), gel(Fbi,k));
    1344       13090 :       gmael(gram,i,j) = comFij;
    1345             :     }
    1346             :   }
    1347        5117 :   return gram;
    1348             : }
    1349             : 
    1350             : static GEN
    1351         532 : gen_comb(long cdep, GEN A, GEN e, struct qfauto *qf, long lim)
    1352             : {
    1353         532 :   long i, dim = lg(A)-1;
    1354         532 :   GEN comb = cgetg(dim+1,t_VEC);
    1355        5649 :   for (i = 1; i <= dim; ++i)
    1356             :   {
    1357        5243 :     pari_sp av = avma;
    1358             :     GEN trans, ccoef, cF, B, BI;
    1359             :     GEN sumveclist, sumvecbase;
    1360        5243 :     GEN list = scpvecs(&sumveclist, i, e, cdep, qf);
    1361        5243 :     GEN M = zm_to_ZM(sumveclist);
    1362        5243 :     GEN T = lllgramint(qf_apply_ZM(A,M));
    1363        5243 :     if (lim && lg(T)-1>=lim) return NULL;
    1364        5117 :     B = ZM_mul(M,T);
    1365        5117 :     BI = RgM_inv(B);
    1366        5117 :     sumvecbase = ZM_trunc_to_zm(B);
    1367        5117 :     trans = ZM_trunc_to_zm(T);
    1368        5117 :     ccoef = ZM_trunc_to_zm(RgM_mul(BI,M));
    1369        5117 :     cF = scpforms(sumvecbase, qf);
    1370        5117 :     gel(comb,i) = gerepilecopy(av, mkvec4(list, trans, ccoef, cF));
    1371             :   }
    1372         406 :   return comb;
    1373             : }
    1374             : 
    1375             : static void
    1376         126 : init_comb(struct qfcand *cand, GEN A, GEN e, struct qfauto *qf)
    1377             : {
    1378         126 :   long dim = lg(A)-1;
    1379         126 :   GEN Am = zm_to_ZM(A);
    1380         364 :   for (cand->cdep = 1; ; cand->cdep++)
    1381             :   {
    1382         602 :     cand->comb = gen_comb(cand->cdep, Am, e, qf, (dim+1)>>1);
    1383         364 :     if (!cand->comb) break;
    1384             :   }
    1385         126 :   cand->cdep= maxss(1, cand->cdep-1);
    1386         126 :   cand->comb = gen_comb(cand->cdep, Am, e, qf, 0);
    1387         126 : }
    1388             : 
    1389             : static void
    1390         168 : init_flags(struct qfcand *cand, GEN A, struct fingerprint *fp,
    1391             :                                        struct qfauto *qf, GEN flags)
    1392             : {
    1393         168 :   if (!flags)
    1394             :   {
    1395         126 :     init_comb(cand, A, fp->e, qf);
    1396         126 :     cand->bacher_pol = init_bacher(0, fp, qf);
    1397             :   }
    1398             :   else
    1399             :   {
    1400             :     long cdep, bach;
    1401          42 :     if (typ(flags)!=t_VEC || lg(flags)!=3)
    1402           0 :       pari_err_TYPE("qfisominit",flags);
    1403          42 :     cdep = gtos(gel(flags,1));
    1404          42 :     bach = minss(gtos(gel(flags,2)),lg(fp->e)-1);
    1405          42 :     if (cdep<0 || bach<0) pari_err_FLAG("qfisom");
    1406          42 :     cand->cdep = cdep;
    1407          42 :     cand->comb = cdep ? gen_comb(cdep, zm_to_ZM(A), fp->e, qf, 0): NULL;
    1408          42 :     cand->bacher_pol = init_bacher(bach, fp, qf);
    1409             :   }
    1410         168 : }
    1411             : 
    1412             : static GEN
    1413         112 : gen_group(struct group *G, GEN U)
    1414             : {
    1415             :   GEN V;
    1416         112 :   long i, j, n=1, dim = lg(G->ord)-1;
    1417         112 :   GEN o = gen_1;
    1418        1400 :   for (i = 1; i <= dim; ++i)
    1419        1288 :     o = muliu(o, G->ord[i]);
    1420        1400 :   for (i = 1; i <= dim; ++i)
    1421        1288 :     n += G->ng[i]-G->nsg[i];
    1422         112 :   V = cgetg(n, t_VEC); n = 1;
    1423        1400 :   for (i = 1; i <= dim; ++i)
    1424        1988 :     for (j=G->nsg[i]+1; j<=G->ng[i]; j++)
    1425        1519 :       gel(V,n++) = U ? zm_mul(gel(U,1), zm_mul(gmael(G->g,i,j), gel(U,2)))
    1426         819 :                      : gmael(G->g,i,j);
    1427         112 :   return mkvec2(o, V);
    1428             : }
    1429             : 
    1430             : static long
    1431         420 : is_qfisom(GEN F)
    1432             : {
    1433         609 :   return (lg(F)==6 && typ(F)==t_VEC && typ(gel(F,1))==t_VEC
    1434         588 :                    && typ(gel(F,3))==t_VEC && typ(gel(F,4))==t_VEC);
    1435             : }
    1436             : 
    1437             : static GEN
    1438          84 : unpack_qfisominit(GEN F, GEN *norm, struct qfauto *qf,
    1439             :       struct fingerprint *fp, struct qfcand *cand)
    1440             : {
    1441          84 :   GEN QF = gel(F,3);
    1442          84 :   qf->F = gel(QF,1);
    1443          84 :   qf->V = gel(QF,2);
    1444          84 :   qf->W = gel(QF,3);
    1445          84 :   qf->v = gel(QF,4);
    1446          84 :   qf->p = itou(gel(QF,5));
    1447          84 :   qf->U = lg(QF)>6 ? gel(QF,6):NULL;
    1448          84 :   QF = gel(F,4);
    1449          84 :   fp->diag = gel(QF,1);
    1450          84 :   fp->per  = gel(QF,2);
    1451          84 :   fp->e    = gel(QF,3);
    1452          84 :   QF = gel(F,5);
    1453          84 :   cand->cdep =itos(gel(QF,1));
    1454          84 :   cand->comb = gel(QF,2);
    1455          84 :   cand->bacher_pol = gel(QF,3);
    1456          84 :   *norm = gel(F,2);
    1457          84 :   qf->dim = lg(gmael(F,1,1))-1;
    1458          84 :   return qf->F;
    1459             : }
    1460             : 
    1461             : static GEN
    1462         154 : qfisom_bestmat(GEN A, long *pt_max)
    1463             : {
    1464         154 :   long max = zm_maxdiag(A), max2;
    1465         154 :   GEN A1 = zm_to_ZM(A), A2;
    1466         154 :   GEN U = lllgramint(A1);
    1467         154 :   if (lg(U) != lg(A1))
    1468           0 :     pari_err_DOMAIN("qfisom","form","is not", strtoGENstr("positive definite"), A1);
    1469         154 :   A2 = ZM_to_zm(qf_apply_ZM(A1, U));
    1470         154 :   max2 = zm_maxdiag(A2);
    1471         154 :   if (max2 < max)
    1472             :   {
    1473          42 :     *pt_max = max2; return mkvec2(ZM_to_zm(U),ZM_to_zm(ZM_inv(U,NULL)));
    1474             :   } else
    1475             :   {
    1476         112 :     *pt_max = max; return NULL;
    1477             :   }
    1478             : }
    1479             : 
    1480             : static GEN
    1481         252 : init_qfisom(GEN F, struct fingerprint *fp, struct qfcand *cand,
    1482             :                    struct qfauto *qf, GEN flags, long *max, GEN minvec)
    1483             : {
    1484             :   GEN U, A, norm;
    1485         252 :   if (is_qfisom(F))
    1486             :   {
    1487          84 :     F = unpack_qfisominit(F, &norm, qf, fp, cand);
    1488          84 :     A = gel(F,1);
    1489          84 :     *max = zm_maxdiag(A);
    1490          84 :     if (flags)
    1491           0 :       init_flags(cand, A, fp, qf, flags);
    1492             :   }
    1493             :   else
    1494             :   {
    1495         168 :     if (lg(F)<2) pari_err_TYPE("qfisom",F);
    1496         168 :     A = gel(F,1);
    1497         168 :     if (lg(A)<2) pari_err_TYPE("qfisom",A);
    1498         168 :     if (!minvec)
    1499             :     {
    1500         154 :       U = qfisom_bestmat(A, max);
    1501         154 :       if (DEBUGLEVEL) err_printf("QFIsom: max=%ld\n",*max);
    1502         154 :       if (U) F = zmV_apply_zm(F, gel(U,1));
    1503             :     } else
    1504             :     {
    1505          14 :       *max = zm_maxdiag(A); U = NULL;
    1506          14 :       if (typ(minvec)==t_VEC && lg(minvec)==4 && typ(gel(minvec,2))==t_INT)
    1507             :       {
    1508           7 :         long n = itos(gel(minvec,2));
    1509           7 :         if (n != *max)
    1510           0 :           pari_err_DOMAIN("qfisominit","m[2]","!=",stoi(*max),stoi(n));
    1511           7 :         minvec = gel(minvec, 3);
    1512             :       }
    1513          14 :       if (typ(minvec)!=t_MAT || lg(gel(minvec,1))!=lg(A))
    1514           0 :         pari_err_TYPE("qfisominit",minvec);
    1515             :     }
    1516         168 :     norm = init_qfauto(F, U, *max, qf, NULL, minvec);
    1517         168 :     fingerprint(fp, qf);
    1518         168 :     if (DEBUGLEVEL) err_printf("QFIsom: fp=%Ps\n",fp->diag);
    1519         168 :     init_flags(cand, A, fp, qf, flags);
    1520             :   }
    1521         252 :   return norm;
    1522             : }
    1523             : 
    1524             : GEN
    1525         112 : qfauto(GEN F, GEN flags)
    1526             : {
    1527         112 :   pari_sp av = avma;
    1528             :   struct fingerprint fp;
    1529             :   struct group G;
    1530             :   struct qfcand cand;
    1531             :   struct qfauto qf;
    1532             :   long max;
    1533         112 :   (void)init_qfisom(F, &fp, &cand, &qf, flags, &max, NULL);
    1534         112 :   init_qfgroup(&G, &fp, &qf);
    1535         112 :   autom(&G, &qf, &fp, &cand);
    1536         112 :   return gerepilecopy(av, gen_group(&G, qf.U));
    1537             : }
    1538             : 
    1539             : static GEN
    1540         259 : qf_to_zmV(GEN F)
    1541             : {
    1542         259 :   return typ(F)==t_MAT ?
    1543         224 :      (RgM_is_ZM(F) ? mkvec(ZM_to_zm(F)): NULL)
    1544         553 :      : typ(F)==t_VEC ? (RgV_is_ZMV(F) ? ZMV_to_zmV(F): NULL)
    1545          70 :      : NULL;
    1546             : }
    1547             : 
    1548             : GEN
    1549         112 : qfauto0(GEN x, GEN flags)
    1550             : {
    1551         112 :   pari_sp av = avma;
    1552             :   GEN F, G;
    1553         112 :   if (is_qfisom(x))
    1554          49 :     F = x;
    1555             :   else
    1556             :   {
    1557          63 :     F = qf_to_zmV(x);
    1558          63 :     if (!F) pari_err_TYPE("qfauto",x);
    1559             :   }
    1560         112 :   G = qfauto(F, flags);
    1561         112 :   return gerepilecopy(av, mkvec2(gel(G,1), zmV_to_ZMV(gel(G,2))));
    1562             : }
    1563             : /* computes the orbit of V.v[pt] under the generators G[0],...,G[nG-1] and
    1564             :  * elements stabilizing V.v[pt], which are stored in H, returns the number of
    1565             :  * generators in H */
    1566             : 
    1567             : static GEN
    1568         609 : isostab(long pt, GEN G, GEN V, long Maxfail, ulong p)
    1569             : {
    1570         609 :   pari_sp av = avma;
    1571             :   long  len, cnd, orblen, tmplen, rpt;
    1572             :   GEN  w, flag, orb;
    1573             :   long  i, im, nH, fail;
    1574         609 :   long dim = lg(gel(V,1))-1, n = lg(V)-1, nG = lg(G)-1;
    1575             :   GEN H;
    1576             :   /* a heuristic break condition for the computation of stabilizer elements:
    1577             :      it would be too expensive to calculate all the stabilizer generators,
    1578             :      which are obtained from the orbit, since this is highly redundant,
    1579             :      on the other hand every new generator which enlarges the group reduces the
    1580             :      number of orbits and hence the number of candidates to be tested, after
    1581             :      Maxfail subsequent stabilizer elements, that do not enlarge the group, the
    1582             :      procedure stops, increasing Maxfail will possibly decrease the number of
    1583             :      tests, but will increase the running time of the stabilizer computation
    1584             :      there is no magic behind the heuristic, tuning might be appropriate */
    1585             :   /* H are the generators of the stabilizer of V.v[pt] */
    1586         609 :   H = cgetg(2,t_VEC);
    1587         609 :   nH = 0;
    1588             :   /* w[i+V.n] is a matrix that maps V.v[pt] on V.v[i] */
    1589         609 :   w = cgetg(2*n+2,t_MAT);
    1590         609 :   orb = zero_Flv(2*n);
    1591             :   /* orblen is the length of the orbit of a random vector in V.v */
    1592         609 :   orblen = 1;
    1593             :   /* if flag[i+V.n] = 1, then the point i is already in the orbit */
    1594         609 :   flag = zero_Flv(2*n+1);
    1595         609 :   orb[1] = pt;
    1596         609 :   flag[orb[1]+n+1] = 1;
    1597             :   /* w[pt+V.n] is the Identity */
    1598         609 :   gel(w,orb[1]+n+1) = matid_Flm(dim);
    1599         609 :   cnd = 1;
    1600         609 :   len = 1;
    1601             :   /* fail is the number of successive failures */
    1602         609 :   fail = 0;
    1603        2184 :   while (cnd <= len &&  fail < Maxfail)
    1604             :   {
    1605        2366 :     for (i = 1; i <= nG  &&  fail < Maxfail; ++i)
    1606             :     {
    1607        1400 :       im = operate(orb[cnd], gel(G,i), V);
    1608        1400 :       if (flag[im+n+1] == 0)
    1609             :         /* a new element is found, appended to the orbit and an element mapping
    1610             :            V.v[pt] to im is stored in w[im+V.n] */
    1611             :       {
    1612         420 :         orb[++len] = im;
    1613         420 :         flag[im+n+1] = 1;
    1614         420 :         gel(w,im+n+1)= zm_mul(gel(G,i), gel(w,orb[cnd]+n+1));
    1615             :       }
    1616             :       else /* the image was already in the orbit */
    1617             :       {
    1618         980 :         GEN B = zm_mul(gel(G,i), gel(w,orb[cnd]+n+1));
    1619             :         /* check whether the old and the new element mapping pt on im differ */
    1620         980 :         if (!zvV_equal(B, gel(w,im+n+1)))
    1621             :         {
    1622         777 :           gel(H,nH+1) = zm_divmod(gel(w,im+n+1),B,p);
    1623         777 :           rpt = 1+(long)random_Fl(n);
    1624         777 :           tmplen = orbitlen(rpt, 2*n, H, nH+1, V);
    1625       12530 :           while (tmplen < orblen)
    1626             :             /* the orbit of this vector is shorter than a previous one, hence
    1627             :                choose a new random vector */
    1628             :           {
    1629       10976 :             rpt = 1+(long)random_Fl(n);
    1630       10976 :             tmplen = orbitlen(rpt, 2*n, H, nH+1, V);
    1631             :           }
    1632         777 :           if (tmplen > orblen)
    1633             :             /* the new stabilizer element H[nH] enlarges the group generated by
    1634             :                H */
    1635             :           {
    1636         175 :             orblen = tmplen;
    1637             :             /* allocate memory for the new generator */
    1638         175 :             H = vec_lengthen(H, (++nH)+1);
    1639         175 :             fail = 0;
    1640             :           }
    1641             :           else
    1642             :             /* the new stabilizer element does not enlarge the orbit of a
    1643             :                random vector */
    1644         602 :             ++fail;
    1645             :         }
    1646             :         /* if H[nH] is the identity, nothing is done */
    1647             :       }
    1648             :     }
    1649         966 :     ++cnd;
    1650             :   }
    1651         609 :   setlg(H,nH+1);
    1652         609 :   return gerepilecopy(av, H);
    1653             : }
    1654             : 
    1655             : /* the heart of the program: the recursion */
    1656             : 
    1657             : static long
    1658         665 : iso(long step, GEN x, GEN C, struct qfauto *qf, struct qfauto *qff,
    1659             :     struct fingerprint *fp, GEN G, struct qfcand *cand)
    1660             : {
    1661             :   int  i, Maxfail;
    1662             :   GEN H;
    1663         665 :   long dim = qf->dim;
    1664         665 :   long found = 0;
    1665        1519 :   while (mael(C,step,1) != 0  &&  found == 0)
    1666             :   {
    1667         798 :     if (step < dim)
    1668             :     {
    1669             :       long nbc;
    1670             :       /* choose the image of the base-vector nr. step */
    1671         742 :       x[step] = mael(C,step,1);
    1672             :       /* check whether x[0]...x[step] is a partial automorphism and compute the
    1673             :          candidates for x[step+1] */
    1674         742 :       nbc = qfisom_candidates(gel(C,step+1), step+1, x, qf, qff, fp, cand);
    1675         742 :       if (nbc == fp->diag[step+1])
    1676             :       {
    1677             :         /* go deeper into the recursion */
    1678         609 :         Maxfail = 0;
    1679             :         /* determine the heuristic value of Maxfail for the break condition in
    1680             :            isostab */
    1681        5019 :         for (i = 1; i <= step; ++i)
    1682        4410 :           if (fp->diag[i] > 1)
    1683        3535 :             Maxfail += 1;
    1684        5019 :         for (i = step+1; i <= dim; ++i)
    1685        4410 :           if (fp->diag[i] > 1)
    1686        3710 :             Maxfail += 2;
    1687             :         /* compute the stabilizer H of x[step] in G */
    1688         609 :         H = isostab(x[step], G, qff->V, Maxfail,qff->p);
    1689         609 :         found = iso(step+1, x, C, qf, qff, fp, H, cand);
    1690             :       }
    1691         742 :       if (found == 1)
    1692         609 :         return 1;
    1693             :       /* delete the orbit of the chosen vector from the list of candidates */
    1694         133 :       orbsubtract(gel(C,step), x, step-1, 1, G, qff->V, NULL);
    1695             :     }
    1696             :     else
    1697             :     {
    1698             :       /* an isomorphism is found */
    1699          56 :       x[dim] = mael(C,dim,1);
    1700          56 :       found = 1;
    1701             :     }
    1702             :   }
    1703          56 :   return found;
    1704             : }
    1705             : 
    1706             : /* search for an isometry */
    1707             : 
    1708             : static GEN
    1709          56 : isometry(struct qfauto *qf, struct qfauto *qff, struct fingerprint *fp, GEN G,
    1710             :                                                 struct qfcand *cand)
    1711             : 
    1712             : {
    1713             :   long i, found;
    1714             :   GEN x;
    1715          56 :   long dim = qf->dim;
    1716          56 :   GEN C = cgetg(dim+1,t_VEC);
    1717             :   /* C[i] is the list of candidates for the image of the i-th base-vector */
    1718         721 :   for (i = 1; i <= dim; ++i)
    1719         665 :     gel(C,i) = cgetg(fp->diag[i]+1, t_VECSMALL);
    1720          56 :   x = cgetg(dim+1, t_VECSMALL);
    1721             :   /* compute the candidates for x[1] */
    1722          56 :   qfisom_candidates(gel(C,1), 1, x, qf, qff, fp, cand);
    1723          56 :   found = iso(1, x, C, qf, qff, fp, G, cand);
    1724          56 :   return found ? matgen(x, fp->per, qff->V): NULL;
    1725             : }
    1726             : 
    1727             : GEN
    1728          84 : qfisominit(GEN F, GEN flags, GEN minvec)
    1729             : {
    1730          84 :   pari_sp av = avma;
    1731             :   struct fingerprint fp;
    1732             :   struct qfauto qf;
    1733             :   struct qfcand cand;
    1734             :   long max;
    1735          84 :   GEN norm = init_qfisom(F, &fp, &cand, &qf, flags, &max, minvec);
    1736         252 :   return gerepilecopy(av, mkvec5(F, norm,
    1737          84 :                           mkvecn(qf.U?6:5, qf.F, qf.V, qf.W, qf.v, utoi(qf.p), qf.U),
    1738             :                           mkvec3(fp.diag, fp.per, fp.e),
    1739          84 :                           mkvec3(stoi(cand.cdep),cand.comb?cand.comb:cgetg(1,t_VEC),
    1740             :                                  cand.bacher_pol)));
    1741             : }
    1742             : 
    1743             : GEN
    1744          84 : qfisominit0(GEN x, GEN flags, GEN minvec)
    1745             : {
    1746          84 :   pari_sp av = avma;
    1747          84 :   GEN F = qf_to_zmV(x);
    1748          84 :   if (!F) pari_err_TYPE("qfisom",x);
    1749          84 :   return gerepileupto(av, qfisominit(F, flags, minvec));
    1750             : }
    1751             : 
    1752             : GEN
    1753          56 : qfisom(GEN F, GEN FF, GEN flags, GEN G)
    1754             : {
    1755          56 :   pari_sp av = avma;
    1756             :   struct fingerprint fp;
    1757             :   GEN res;
    1758             :   struct qfauto qf, qff;
    1759             :   struct qfcand cand;
    1760             :   long max;
    1761          56 :   GEN norm = init_qfisom(F, &fp, &cand, &qf, flags, &max, NULL);
    1762          56 :   init_qfauto(FF, NULL, max, &qff, norm, NULL);
    1763          56 :   if (lg(qf.W)!=lg(qff.W)
    1764          56 :       || !zvV_equal(vecvecsmall_sort(qf.W), vecvecsmall_sort(qff.W)))
    1765           0 :     { set_avma(av); return gen_0; }
    1766          56 :   if (!G) G = mkvec(scalar_Flm(-1, qff.dim));
    1767          56 :   res = isometry(&qf, &qff, &fp, G, &cand);
    1768          56 :   if (!res) { set_avma(av); return gen_0; }
    1769          56 :   return gerepilecopy(av, zm_to_ZM(qf.U? zm_mul(res,gel(qf.U, 2)):res));
    1770             : }
    1771             : 
    1772             : static GEN
    1773          35 : check_qfauto(GEN G)
    1774             : {
    1775          35 :   if (typ(G)==t_VEC && lg(G)==3 && typ(gel(G,1))==t_INT)
    1776          28 :     G = gel(G,2);
    1777          35 :   return qf_to_zmV(G);
    1778             : }
    1779             : 
    1780             : GEN
    1781          56 : qfisom0(GEN x, GEN y, GEN flags, GEN G)
    1782             : {
    1783          56 :   pari_sp av = avma;
    1784             :   GEN F, FF;
    1785          56 :   if (is_qfisom(x))
    1786          35 :     F = x;
    1787             :   else
    1788             :   {
    1789          21 :     F = qf_to_zmV(x);
    1790          21 :     if (!F) pari_err_TYPE("qfisom",x);
    1791             :   }
    1792          56 :   FF = qf_to_zmV(y);
    1793          56 :   if (!FF) pari_err_TYPE("qfisom",y);
    1794          56 :   if (G) G = check_qfauto(G);
    1795          56 :   return gerepileupto(av, qfisom(F, FF, flags, G));
    1796             : }
    1797             : 
    1798             : static GEN
    1799          42 : ZM_to_GAP(GEN M)
    1800             : {
    1801          42 :   pari_sp ltop=avma;
    1802          42 :   long rows = nbrows(M), cols = lg(M)-1;
    1803             :   long i, j, c;
    1804          42 :   GEN comma = strtoGENstr(", ");
    1805          42 :   GEN bra = strtoGENstr("[");
    1806          42 :   GEN ket = strtoGENstr("]");
    1807          42 :   GEN s = cgetg(2*rows*cols+2*rows+2,t_VEC);
    1808          42 :   gel(s,1) = bra; c=2;
    1809         210 :   for (i = 1; i <= rows; ++i)
    1810             :   {
    1811         168 :     if (i > 1)
    1812         126 :       gel(s,c++) = comma;
    1813         168 :     gel(s,c++) = bra;
    1814         840 :     for (j = 1; j <= cols; ++j)
    1815             :     {
    1816         672 :       if (j > 1)
    1817         504 :         gel(s,c++) = comma;
    1818         672 :       gel(s,c++) = GENtoGENstr(gcoeff(M,i,j));
    1819             :     }
    1820         168 :     gel(s,c++) = ket;
    1821             :   }
    1822          42 :   gel(s,c++) = ket;
    1823          42 :   return gerepilecopy(ltop,shallowconcat1(s));
    1824             : }
    1825             : 
    1826             : GEN
    1827          14 : qfautoexport(GEN G, long flag)
    1828             : {
    1829          14 :   pari_sp av = avma;
    1830          14 :   long i, lgen,  c = 2;
    1831          14 :   GEN gen, str, comma = strtoGENstr(", ");
    1832          14 :   if (typ(G)!=t_VEC || lg(G)!=3) pari_err_TYPE("qfautoexport", G);
    1833          14 :   if (flag!=0 && flag!=1) pari_err_FLAG("qfautoexport");
    1834          14 :   gen = gel(G,2); lgen = lg(gen)-1;
    1835          14 :   str = cgetg(2+2*lgen,t_VEC);
    1836             :   /* in GAP or MAGMA the matrix group is called BG */
    1837          14 :   if (flag == 0)
    1838           7 :     gel(str,1) = strtoGENstr("Group(");
    1839             :   else
    1840             :   {
    1841           7 :     long dim = lg(gmael(gen,1,1))-1;
    1842           7 :     gel(str,1) = gsprintf("MatrixGroup<%d, Integers() |",dim);
    1843             :   }
    1844          56 :   for(i = 1; i <= lgen; i++)
    1845             :   {
    1846          42 :     if (i!=1)
    1847          28 :       gel(str,c++) = comma;
    1848          42 :     gel(str,c++) = ZM_to_GAP(gel(gen,i));
    1849             :   }
    1850          14 :   gel(str,c++) = strtoGENstr(flag ? ">":")");
    1851          14 :   return gerepilecopy(av, shallowconcat1(str));
    1852             : }
    1853             : 
    1854             : GEN
    1855          21 : qforbits(GEN G, GEN V)
    1856             : {
    1857          21 :   pari_sp av = avma;
    1858             :   GEN gen, w, W, p, v, orb, o;
    1859             :   long i, j, n, ng;
    1860          21 :   long nborbits = 0;
    1861          21 :   gen = check_qfauto(G);
    1862          21 :   if (!gen) pari_err_TYPE("qforbits", G);
    1863          21 :   if (typ(V)==t_VEC && lg(V)==4
    1864           7 :    && typ(gel(V,1))==t_INT && typ(gel(V,2))==t_INT)
    1865           7 :     V = gel(V,3);
    1866          21 :   if (typ(V)!=t_MAT || !RgM_is_ZM(V)) pari_err_TYPE("qforbits", V);
    1867          21 :   n = lg(V)-1; ng = lg(gen)-1;
    1868          21 :   W = ZM_to_zm_canon(V);
    1869          21 :   p = vecvecsmall_indexsort(W);
    1870          21 :   v = vecpermute(W, p);
    1871          21 :   w = zero_zv(n);
    1872          21 :   orb = cgetg(n+1, t_VEC);
    1873          21 :   o = cgetg(n+1, t_VECSMALL);
    1874          21 :   if (lg(v) != lg(V)) return gen_0;
    1875         357 :   for (i=1; i<=n; i++)
    1876             :   {
    1877         336 :     long cnd, no = 1;
    1878             :     GEN T;
    1879         336 :     if (w[i]) continue;
    1880          42 :     w[i] = ++nborbits;
    1881          42 :     o[1] = i;
    1882         378 :     for (cnd=1; cnd <= no; ++cnd)
    1883        4032 :       for(j=1; j <= ng; j++)
    1884             :       {
    1885             :         long k;
    1886        3696 :         GEN Vij = zm_zc_mul(gel(gen, j), gel(v, o[cnd]));
    1887        3696 :         (void) zv_canon(Vij);
    1888        3696 :         k = vecvecsmall_search(v, Vij, 0);
    1889        3696 :         if (k == 0) { set_avma(av); return gen_0; }
    1890        3696 :         if (w[k] == 0)
    1891             :         {
    1892         294 :           o[++no] = k;
    1893         294 :           w[k] = nborbits;
    1894             :         }
    1895             :       }
    1896          42 :     T = cgetg(no+1, t_VEC);
    1897          42 :     for (j=1; j<=no; j++) gel(T,j) = gel(V,p[o[j]]);
    1898          42 :     gel(orb, nborbits) = T;
    1899             :   }
    1900          21 :   setlg(orb, nborbits+1);
    1901          21 :   return gerepilecopy(av, orb);
    1902             : }

Generated by: LCOV version 1.13