Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - qfsolve.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 17953-c39f2e6) Lines: 595 604 98.5 %
Date: 2015-08-31 Functions: 30 30 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 359 395 90.9 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2000-2004  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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13                 :            : 
      14                 :            : /* Copyright (C) 2014 Denis Simon
      15                 :            :  * Adapted from qfsolve.gp v. 09/01/2014
      16                 :            :  *   http://www.math.unicaen.fr/~simon/qfsolve.gp
      17                 :            :  *
      18                 :            :  * Author: Denis SIMON <simon@math.unicaen.fr> */
      19                 :            : 
      20                 :            : #include "pari.h"
      21                 :            : #include "paripriv.h"
      22                 :            : 
      23                 :            : /* LINEAR ALGEBRA */
      24                 :            : /* complete by 0s, assume l-1 <= n */
      25                 :            : static GEN
      26                 :      33685 : vecextend(GEN v, long n)
      27                 :            : {
      28                 :      33685 :   long i, l = lg(v);
      29                 :      33685 :   GEN w = cgetg(n+1, t_COL);
      30         [ +  + ]:     111074 :   for (i = 1; i < l; i++) gel(w,i) = gel(v,i);
      31         [ +  + ]:     220282 :   for (     ; i <=n; i++) gel(w,i) = gen_0;
      32                 :      33685 :   return w;
      33                 :            : }
      34                 :            : 
      35                 :            : /* Gives a unimodular matrix with the last column(s) equal to Mv.
      36                 :            :  * Mv can be a column vector or a rectangular matrix.
      37                 :            :  * redflag = 0 or 1. If redflag = 1, LLL-reduce the n-#v first columns. */
      38                 :            : static GEN
      39                 :      83497 : completebasis(GEN Mv, long redflag)
      40                 :            : {
      41                 :            :   GEN U;
      42                 :            :   long m, n;
      43                 :            : 
      44         [ +  + ]:      83497 :   if (typ(Mv) == t_COL) Mv = mkmat(Mv);
      45                 :      83497 :   n = lg(Mv)-1;
      46                 :      83497 :   m = nbrows(Mv); /* m x n */
      47         [ +  + ]:      83497 :   if (m == n) return Mv;
      48                 :      83427 :   (void)ZM_hnfall(shallowtrans(Mv), &U, 0);
      49                 :      83427 :   U = ZM_inv(shallowtrans(U), gen_1);
      50 [ +  - ][ +  + ]:      83427 :   if (m==1 || !redflag) return U;
      51                 :            :   /* LLL-reduce the m-n first columns */
      52                 :      83497 :   return shallowconcat(ZM_lll(vecslice(U,1,m-n), 0.99, LLL_INPLACE),
      53                 :      33488 :                        vecslice(U, m-n+1,m));
      54                 :            : }
      55                 :            : 
      56                 :            : /* Compute the kernel of M mod p.
      57                 :            :  * returns [d,U], where
      58                 :            :  * d = dim (ker M mod p)
      59                 :            :  * U in GLn(Z), and its first d columns span the kernel. */
      60                 :            : static GEN
      61                 :      49805 : kermodp(GEN M, GEN p, long *d)
      62                 :            : {
      63                 :            :   long j, l;
      64                 :            :   GEN K, B, U;
      65                 :            : 
      66                 :      49805 :   K = FpM_center(FpM_ker(M, p), p, shifti(p,-1));
      67                 :      49805 :   B = completebasis(K,0);
      68                 :      49805 :   l = lg(M); U = cgetg(l, t_MAT);
      69         [ +  + ]:     420560 :   for (j =  1; j < l; j++) gel(U,j) = gel(B,l-j);
      70                 :      49805 :   *d = lg(K)-1; return U;
      71                 :            : }
      72                 :            : 
      73                 :            : /* INVARIANTS COMPUTATIONS */
      74                 :            : 
      75                 :            : static GEN
      76                 :      28876 : principal_minor(GEN G, long  i) { return matslice(G,1,i,1,i); }
      77                 :            : static GEN
      78                 :        784 : det_minors(GEN G)
      79                 :            : {
      80                 :        784 :   long i, l = lg(G);
      81                 :        784 :   GEN v = cgetg(l+1, t_VEC);
      82                 :        784 :   gel(v,1) = gen_1;
      83         [ +  + ]:       3920 :   for (i = 2; i <= l; i++) gel(v,i) = ZM_det(principal_minor(G,i-1));
      84                 :        784 :   return v;
      85                 :            : }
      86                 :            : 
      87                 :            : /* Given a symmetric matrix G over Z, compute the Witt invariant
      88                 :            :  *  of G at the prime p (at real place if p = NULL)
      89                 :            :  * Assume that none of the determinant G[1..i,1..i] is 0. */
      90                 :            : static long
      91                 :        119 : qflocalinvariant(GEN G, GEN p)
      92                 :            : {
      93                 :        119 :   long i, j, c, l = lg(G);
      94                 :        119 :   GEN diag, v = det_minors(G);
      95                 :            :   /* Diagonalize G first. */
      96                 :        119 :   diag = cgetg(l, t_VEC);
      97         [ +  + ]:        595 :   for (i = 1; i < l; i++) gel(diag,i) = mulii(gel(v,i+1), gel(v,i));
      98                 :            : 
      99                 :            :   /* Then compute the product of the Hilbert symbols */
     100                 :            :   /* (diag[i],diag[j])_p for i < j */
     101                 :        119 :   c = 1;
     102         [ +  + ]:        476 :   for (i = 1; i < l-1; i++)
     103         [ +  + ]:       1071 :     for (j = i+1; j < l; j++)
     104         [ +  + ]:        714 :       if (hilbertii(gel(diag,i), gel(diag,j), p) < 0) c = -c;
     105                 :        119 :   return c;
     106                 :            : }
     107                 :            : 
     108                 :            : static GEN
     109                 :       7273 : hilberts(GEN a, GEN b, GEN P, long lP)
     110                 :            : {
     111                 :       7273 :   GEN v = cgetg(lP, t_VECSMALL);
     112                 :            :   long i;
     113         [ +  + ]:      42189 :   for (i = 1; i < lP; i++) v[i] = hilbertii(a, b, gel(P,i)) < 0;
     114                 :       7273 :   return v;
     115                 :            : }
     116                 :            : 
     117                 :            : /* G symmetrix matrix or qfb or list of quadratic forms with same discriminant.
     118                 :            :  * P must be equal to factor(-abs(2*matdet(G)))[,1]. */
     119                 :            : static GEN
     120                 :       4620 : qflocalinvariants(GEN G, GEN P)
     121                 :            : {
     122                 :            :   GEN sol;
     123                 :       4620 :   long i, j, l, lP = lg(P);
     124                 :            : 
     125                 :            :   /* convert G into a vector of symmetric matrices */
     126         [ +  + ]:       4620 :   G = (typ(G) == t_VEC)? shallowcopy(G): mkvec(G);
     127                 :       4620 :   l = lg(G);
     128         [ +  + ]:      12019 :   for (j = 1; j < l; j++)
     129                 :            :   {
     130                 :       7399 :     GEN g = gel(G,j);
     131 [ +  + ][ +  + ]:       7399 :     if (typ(g) == t_QFI || typ(g) == t_QFR) gel(G,j) = gtomat(g);
     132                 :            :   }
     133                 :       4620 :   sol = cgetg(l, t_MAT);
     134         [ +  + ]:       4620 :   if (lg(gel(G,1)) == 3)
     135                 :            :   { /* in dimension 2, each invariant is a single Hilbert symbol. */
     136                 :       3955 :     GEN d = negi(ZM_det(gel(G,1)));
     137         [ +  + ]:      10689 :     for (j = 1; j < l; j++)
     138                 :            :     {
     139                 :       6734 :       GEN a = gcoeff(gel(G,j),1,1);
     140                 :       6734 :       gel(sol,j) = hilberts(a, d, P, lP);
     141                 :            :     }
     142                 :            :   }
     143                 :            :   else /* in dimension n > 2, we compute a product of n Hilbert symbols. */
     144         [ +  + ]:       1330 :     for (j = 1; j <l; j++)
     145                 :            :     {
     146                 :        665 :       GEN g = gel(G,j), v = det_minors(g), w = cgetg(lP, t_VECSMALL);
     147                 :        665 :       long n = lg(v);
     148                 :        665 :       gel(sol,j) = w;
     149         [ +  + ]:       3997 :       for (i = 1; i < lP; i++)
     150                 :            :       {
     151                 :       3332 :         GEN p = gel(P,i);
     152                 :       3332 :         long k = n-2, h = hilbertii(gel(v,k), gel(v,k+1),p);
     153         [ +  + ]:      13328 :         for (k--; k >= 1; k--)
     154         [ +  + ]:       9996 :           if (hilbertii(negi(gel(v,k)), gel(v,k+1),p) < 0) h = -h;
     155                 :       3332 :         w[i] = h < 0;
     156                 :            :       }
     157                 :            :     }
     158                 :       4620 :   return sol;
     159                 :            : }
     160                 :            : 
     161                 :            : /* QUADRATIC FORM REDUCTION */
     162                 :            : static GEN
     163                 :       8456 : qfb(GEN D, GEN a, GEN b, GEN c)
     164                 :            : {
     165         [ +  + ]:       8456 :   if (signe(D) < 0) return mkqfi(a,b,c);
     166                 :       8456 :   retmkqfr(a,b,c,real_0(DEFAULTPREC));
     167                 :            : }
     168                 :            : 
     169                 :            : /* Gauss reduction of the binary quadratic form
     170                 :            :  * Q = a*X^2+2*b*X*Y+c*Y^2 of discriminant D (divisible by 4)
     171                 :            :  * returns the reduced form */
     172                 :            : static GEN
     173                 :       1722 : qfbreduce(GEN D, GEN Q)
     174                 :            : {
     175                 :       1722 :   GEN a = gel(Q,1), b = shifti(gel(Q,2),-1), c = gel(Q,3);
     176         [ +  - ]:       5215 :   while (signe(a))
     177                 :            :   {
     178                 :            :     GEN r, q, nexta, nextc;
     179                 :       5215 :     q = dvmdii(b,a, &r); /* FIXME: export as dvmdiiround ? */
     180 [ +  + ][ +  + ]:       5215 :     if (signe(r) > 0 && absi_cmp(shifti(r,1), a) > 0) {
     181                 :       1057 :       r = subii(r, absi(a)); q = addis(q, signe(a));
     182                 :            :     }
     183                 :       5215 :     nextc = a; nexta = subii(c, mulii(q, addii(r,b)));
     184         [ +  + ]:       5215 :     if (absi_cmp(nexta, a) >= 0) break;
     185                 :       5215 :     c = nextc; b = negi(r); a = nexta;
     186                 :            :   }
     187                 :       1722 :   return qfb(D,a,shifti(b,1),c);
     188                 :            : }
     189                 :            : 
     190                 :            : /* private version of qfgaussred:
     191                 :            :  * - early abort if k-th principal minor is singular, return stoi(k)
     192                 :            :  * - else return a matrix whose upper triangular part is qfgaussred(a) */
     193                 :            : static GEN
     194                 :      21203 : partialgaussred(GEN a)
     195                 :            : {
     196                 :      21203 :   long n = lg(a)-1, k;
     197                 :      21203 :   a = RgM_shallowcopy(a);
     198         [ +  + ]:      90231 :   for(k = 1; k < n; k++)
     199                 :            :   {
     200                 :      75167 :     GEN ak, p = gcoeff(a,k,k);
     201                 :            :     long i, j;
     202         [ +  + ]:      75167 :     if (isintzero(p)) return stoi(k);
     203                 :      69028 :     ak = row(a, k);
     204         [ +  + ]:     299277 :     for (i=k+1; i<=n; i++) gcoeff(a,k,i) = gdiv(gcoeff(a,k,i), p);
     205         [ +  + ]:     299277 :     for (i=k+1; i<=n; i++)
     206                 :            :     {
     207                 :     230249 :       GEN c = gel(ak,i);
     208         [ +  + ]:     230249 :       if (gequal0(c)) continue;
     209         [ +  + ]:     775766 :       for (j=i; j<=n; j++)
     210                 :     575765 :         gcoeff(a,i,j) = gsub(gcoeff(a,i,j), gmul(c,gcoeff(a,k,j)));
     211                 :            :     }
     212                 :            :   }
     213         [ +  + ]:      15064 :   if (isintzero(gcoeff(a,n,n))) return stoi(n);
     214                 :      21203 :   return a;
     215                 :            : }
     216                 :            : 
     217                 :            : /* LLL-reduce a positive definite qf QD bounding the indefinite G, dim G > 1.
     218                 :            :  * Then finishes by looking for trivial solution */
     219                 :            : static GEN qftriv(GEN G, GEN z, long base);
     220                 :            : static GEN
     221                 :      21203 : qflllgram_indef(GEN G, long base, int *fail)
     222                 :            : {
     223                 :            :   GEN M, R, g, DM, S, dR;
     224                 :      21203 :   long i, j, n = lg(G)-1;
     225                 :            : 
     226                 :      21203 :   *fail = 0;
     227                 :      21203 :   R = partialgaussred(G);
     228         [ +  + ]:      21203 :   if (typ(R) == t_INT) return qftriv(G, R, base);
     229                 :      15057 :   R = Q_remove_denom(R, &dR); /* avoid rational arithmetic */
     230                 :      15057 :   M = zeromatcopy(n,n);
     231                 :      15057 :   DM = zeromatcopy(n,n);
     232         [ +  + ]:      94750 :   for (i = 1; i <= n; i++)
     233                 :            :   {
     234                 :      79693 :     GEN d = absi_shallow(gcoeff(R,i,i));
     235         [ +  + ]:      79693 :     if (dR) {
     236                 :      70985 :       gcoeff(M,i,i) = dR;
     237                 :      70985 :       gcoeff(DM,i,i) = mulii(d,dR);
     238                 :            :     } else {
     239                 :       8708 :       gcoeff(M,i,i) = gen_1;
     240                 :       8708 :       gcoeff(DM,i,i) = d;
     241                 :            :     }
     242         [ +  + ]:     290021 :     for (j = i+1; j <= n; j++)
     243                 :            :     {
     244                 :     210328 :       gcoeff(M,i,j) = gcoeff(R,i,j);
     245                 :     210328 :       gcoeff(DM,i,j) = mulii(d, gcoeff(R,i,j));
     246                 :            :     }
     247                 :            :   }
     248                 :            :   /* G = M~*D*M, D diagonal, DM=|D|*M, g =  M~*|D|*M */
     249                 :      15057 :   g = ZM_transmultosym(M,DM);
     250                 :      15057 :   S = lllgramint(Q_primpart(g));
     251                 :      15057 :   R = qftriv(qf_apply_ZM(G,S), NULL, base);
     252      [ +  +  + ]:      15057 :   switch(typ(R))
     253                 :            :   {
     254                 :        224 :     case t_COL: return ZM_ZC_mul(S,R);
     255                 :       7210 :     case t_MAT: *fail = 1; return mkvec2(R, S);
     256                 :            :     default:
     257                 :       7623 :       gel(R,2) = ZM_mul(S, gel(R,2));
     258                 :      21203 :       return R;
     259                 :            :   }
     260                 :            : }
     261                 :            : 
     262                 :            : /* G symmetric, i < j, let E = E_{i,j}(a), G <- E~*G*E,  U <- U*E.
     263                 :            :  * Everybody integral */
     264                 :            : static void
     265                 :       4172 : qf_apply_transvect_Z(GEN G, GEN U, long i, long j, GEN a)
     266                 :            : {
     267                 :       4172 :   long k, n = lg(G)-1;
     268                 :       4172 :   gel(G, j) =  ZC_lincomb(gen_1, a, gel(G,j), gel(G,i));
     269         [ +  + ]:      23103 :   for (k = 1; k < n; k++) gcoeff(G, j, k) = gcoeff(G, k, j);
     270                 :       4172 :   gcoeff(G,j,j) = addmulii(gcoeff(G,j,j), a,
     271                 :       8344 :                            addmulii(gcoeff(G,i,j), a,gcoeff(G,i,i)));
     272                 :       4172 :   gel(U, j) =  ZC_lincomb(gen_1, a, gel(U,j), gel(U,i));
     273                 :       4172 : }
     274                 :            : 
     275                 :            : /* LLL reduction of the quadratic form G (Gram matrix)
     276                 :            :  * where we go on, even if an isotropic vector is found. */
     277                 :            : static GEN
     278                 :      10983 : qflllgram_indefgoon(GEN G)
     279                 :            : {
     280                 :            :   GEN red, U, A, U1,U2,U3,U5,U6, V, B, G2,G3,G4,G5, G6, a, g;
     281                 :      10983 :   long i, j, n = lg(G)-1;
     282                 :            :   int fail;
     283                 :            : 
     284                 :      10983 :   red = qflllgram_indef(G,1, &fail);
     285         [ +  + ]:      10983 :   if (fail) return red; /*no isotropic vector found: nothing to do*/
     286                 :            :   /* otherwise a solution is found: */
     287                 :      10976 :   U1 = gel(red,2);
     288                 :      10976 :   G2 = gel(red,1); /* G2[1,1] = 0 */
     289                 :      10976 :   U2 = gel(ZV_extgcd(row(G2,1)), 2);
     290                 :      10976 :   G3 = qf_apply_ZM(G2,U2);
     291                 :      10976 :   U = ZM_mul(U1,U2); /* qf_apply(G,U) = G3 */
     292                 :            :   /* G3[1,] = [0,...,0,g], g^2 | det G */
     293                 :      10976 :   g = gcoeff(G3,1,n);
     294                 :      10976 :   a = diviiround(negi(gcoeff(G3,n,n)), shifti(g,1));
     295         [ +  + ]:      10976 :   if (signe(a)) qf_apply_transvect_Z(G3,U,1,n,a);
     296                 :            :   /* G3[n,n] reduced mod 2g */
     297         [ +  + ]:      10976 :   if (n == 2) return mkvec2(G3,U);
     298                 :      10311 :   V = rowpermute(vecslice(G3, 2,n-1), mkvecsmall2(1,n));
     299                 :      10311 :   A = mkmat2(mkcol2(gcoeff(G3,1,1),gcoeff(G3,1,n)),
     300                 :      20622 :              mkcol2(gcoeff(G3,1,n),gcoeff(G3,2,2)));
     301                 :      10311 :   B = ground(RgM_neg(RgM_mul(RgM_inv(A), V)));
     302                 :      10311 :   U3 = matid(n);
     303         [ +  + ]:      51324 :   for (j = 2; j < n; j++)
     304                 :            :   {
     305                 :      41013 :     gcoeff(U3,1,j) = gcoeff(B,1,j-1);
     306                 :      41013 :     gcoeff(U3,n,j) = gcoeff(B,2,j-1);
     307                 :            :   }
     308                 :      10311 :   G4 = qf_apply_ZM(G3,U3); /* the last column of G4 is reduced */
     309                 :      10311 :   U = ZM_mul(U,U3);
     310         [ +  + ]:      10311 :   if (n == 3) return mkvec2(G4,U);
     311                 :            : 
     312                 :       8064 :   red = qflllgram_indefgoon(matslice(G4,2,n-1,2,n-1));
     313         [ -  + ]:       8064 :   if (typ(red) == t_MAT) return mkvec2(G4,U);
     314                 :            :   /* Let U5:=matconcat(diagonal[1,red[2],1])
     315                 :            :    * return [qf_apply_ZM(G5, U5), U*U5] */
     316                 :       8064 :   G5 = gel(red,1);
     317                 :       8064 :   U5 = gel(red,2);
     318                 :       8064 :   G6 = cgetg(n+1,t_MAT);
     319                 :       8064 :   gel(G6,1) = gel(G4,1);
     320                 :       8064 :   gel(G6,n) = gel(G4,n);
     321         [ +  + ]:      46830 :   for (j=2; j<n; j++)
     322                 :            :   {
     323                 :      38766 :     gel(G6,j) = cgetg(n+1,t_COL);
     324                 :      38766 :     gcoeff(G6,1,j) = gcoeff(G4,j,1);
     325                 :      38766 :     gcoeff(G6,n,j) = gcoeff(G4,j,n);
     326         [ +  + ]:     256368 :     for (i=2; i<n; i++) gcoeff(G6,i,j) = gcoeff(G5,i-1,j-1);
     327                 :            :   }
     328                 :       8064 :   U6 = mkvec3(mkmat(gel(U,1)), ZM_mul(vecslice(U,2,n-1),U5), mkmat(gel(U,n)));
     329                 :      10983 :   return mkvec2(G6, shallowconcat1(U6));
     330                 :            : }
     331                 :            : 
     332                 :            : /* qf_apply_ZM(G,H),  where H = matrix of \tau_{i,j}, i != j */
     333                 :            : static GEN
     334                 :      12375 : qf_apply_tau(GEN G, long i, long j)
     335                 :            : {
     336                 :      12375 :   long l = lg(G), k;
     337                 :      12375 :   G = RgM_shallowcopy(G);
     338                 :      12375 :   swap(gel(G,i), gel(G,j));
     339         [ +  + ]:      77683 :   for (k = 1; k < l; k++) swap(gcoeff(G,i,k), gcoeff(G,j,k));
     340                 :      12375 :   return G;
     341                 :            : }
     342                 :            : 
     343                 :            : /* LLL reduction of the quadratic form G (Gram matrix)
     344                 :            :  * in dim 3 only, with detG = -1 and sign(G) = [2,1]; */
     345                 :            : static GEN
     346                 :       2268 : qflllgram_indefgoon2(GEN G)
     347                 :            : {
     348                 :            :   GEN red, G2, a, b, c, d, e, f, u, v, r, r3, U2, G3;
     349                 :            :   int fail;
     350                 :            : 
     351                 :       2268 :   red = qflllgram_indef(G,1,&fail); /* always find an isotropic vector. */
     352                 :       2268 :   G2 = qf_apply_tau(gel(red,1),1,3); /* G2[3,3] = 0 */
     353                 :       2268 :   r = row(gel(red,2), 3);
     354                 :       2268 :   swap(gel(r,1), gel(r,3)); /* apply tau_{1,3} */
     355                 :       2268 :   a = gcoeff(G2,3,1);
     356                 :       2268 :   b = gcoeff(G2,3,2);
     357                 :       2268 :   d = bezout(a,b, &u,&v);
     358         [ -  + ]:       2268 :   if (!is_pm1(d))
     359                 :            :   {
     360                 :          0 :     a = diviiexact(a,d);
     361                 :          0 :     b = diviiexact(b,d);
     362                 :            :   }
     363                 :            :   /* for U2 = [-u,-b,0;-v,a,0;0,0,1]
     364                 :            :    * G3 = qf_apply_ZM(G2,U2) has known last row (-d, 0, 0),
     365                 :            :    * so apply to principal_minor(G3,2), instead */
     366                 :       2268 :   U2 = mkmat2(mkcol2(negi(u),negi(v)), mkcol2(negi(b),a));
     367                 :       2268 :   G3 = qf_apply_ZM(principal_minor(G2,2),U2);
     368                 :       2268 :   r3 = gel(r,3);
     369                 :       2268 :   r = ZV_ZM_mul(mkvec2(gel(r,1),gel(r,2)),U2);
     370                 :            : 
     371                 :       2268 :   a = gcoeff(G3,1,1);
     372                 :       2268 :   b = gcoeff(G3,1,2);
     373                 :       2268 :   c = negi(d); /* G3[1,3] */
     374                 :       2268 :   d = gcoeff(G3,2,2);
     375         [ +  + ]:       2268 :   if (mpodd(a))
     376                 :            :   {
     377                 :       1260 :     e = addii(b,d);
     378                 :       1260 :     a = addii(a, addii(b,e));
     379                 :       1260 :     e = diviiround(negi(e),c);
     380                 :       1260 :     f = diviiround(negi(a), shifti(c,1));
     381                 :       1260 :     a = addmulii(addii(gel(r,1),gel(r,2)), f,r3);
     382                 :            :   }
     383                 :            :   else
     384                 :            :   {
     385                 :       1008 :     e = diviiround(negi(b),c);
     386                 :       1008 :     f = diviiround(negi(shifti(a,-1)), c);
     387                 :       1008 :     a = addmulii(gel(r,1), f, r3);
     388                 :            :   }
     389                 :       2268 :   b = addmulii(gel(r,2), e, r3);
     390                 :       2268 :   return mkvec3(a,b, r3);
     391                 :            : }
     392                 :            : 
     393                 :            : /* QUADRATIC FORM MINIMIZATION */
     394                 :            : /* G symmetric, return ZM_Z_divexact(G,d) */
     395                 :            : static GEN
     396                 :      61033 : ZsymM_Z_divexact(GEN G, GEN d)
     397                 :            : {
     398                 :      61033 :   long i,j,l = lg(G);
     399                 :      61033 :   GEN H = cgetg(l, t_MAT);
     400         [ +  + ]:     441154 :   for(j=1; j<l; j++)
     401                 :            :   {
     402                 :     380121 :     GEN c = cgetg(l, t_COL), b = gel(G,j);
     403         [ +  + ]:    1600452 :     for(i=1; i<j; i++) gcoeff(H,j,i) = gel(c,i) = diviiexact(gel(b,i),d);
     404                 :     380121 :     gel(c,j) = diviiexact(gel(b,j),d);
     405                 :     380121 :     gel(H,j) = c;
     406                 :            :   }
     407                 :      61033 :   return H;
     408                 :            : }
     409                 :            : 
     410                 :            : /* write symmetric G as [A,B;B~,C], A dxd, C (n-d)x(n-d) */
     411                 :            : static void
     412                 :        497 : blocks4(GEN G, long d, long n, GEN *A, GEN *B, GEN *C)
     413                 :            : {
     414                 :        497 :   GEN G2 = vecslice(G,d+1,n);
     415                 :        497 :   *A = principal_minor(G, d);
     416                 :        497 :   *B = rowslice(G2, 1, d);
     417                 :        497 :   *C = rowslice(G2, d+1, n);
     418                 :        497 : }
     419                 :            : /* Minimization of the quadratic form G, deg G != 0, dim n >= 2
     420                 :            :  * G symmetric integral
     421                 :            :  * Returns [G',U,factd] with U in GLn(Q) such that G'=U~*G*U*constant
     422                 :            :  * is integral and has minimal determinant.
     423                 :            :  * In dimension 3 or 4, may return a prime p if the reduction at p is
     424                 :            :  * impossible because of local non-solvability.
     425                 :            :  * P,E = factor(+/- det(G)), "prime" -1 is ignored. Destroy E. */
     426                 :            : static GEN qfsolvemodp(GEN G, GEN p);
     427                 :            : static GEN
     428                 :       8211 : qfminimize(GEN G, GEN P, GEN E)
     429                 :            : {
     430                 :            :   GEN d, U, Ker, sol, aux, faE, faP;
     431                 :       8211 :   long n = lg(G)-1, lP = lg(P), i, dimKer, m;
     432                 :            : 
     433                 :       8211 :   faP = vectrunc_init(lP);
     434                 :       8211 :   faE = vecsmalltrunc_init(lP);
     435                 :       8211 :   U = NULL;
     436         [ +  + ]:      81592 :   for (i = 1; i < lP; i++)
     437                 :            :   {
     438                 :      73864 :     GEN p = gel(P,i);
     439                 :      73864 :     long vp = E[i];
     440 [ +  + ][ -  + ]:      73864 :     if (!vp || !p) continue;
     441                 :            : 
     442         [ -  + ]:      56868 :     if (DEBUGLEVEL >= 4) err_printf("    p^v = %Ps^%ld\n", p,vp);
     443                 :            :     /* The case vp = 1 can be minimized only if n is odd. */
     444 [ +  + ][ +  + ]:      56868 :     if (vp == 1 && n%2 == 0) {
     445                 :       7560 :       vectrunc_append(faP, p);
     446                 :       7560 :       vecsmalltrunc_append(faE, 1);
     447                 :       7560 :       continue;
     448                 :            :     }
     449                 :      49308 :     Ker = kermodp(G,p, &dimKer); /* dimKer <= vp */
     450         [ -  + ]:      49308 :     if (DEBUGLEVEL >= 4) err_printf("    dimKer = %ld\n",dimKer);
     451         [ -  + ]:      49308 :     if (dimKer == n)
     452                 :            :     { /* trivial case: dimKer = n */
     453         [ #  # ]:          0 :       if (DEBUGLEVEL >= 4) err_printf("     case 0: dimKer = n\n");
     454                 :          0 :       G = ZsymM_Z_divexact(G, p);
     455                 :          0 :       E[i] -= n;
     456                 :          0 :       i--; continue; /* same p */
     457                 :            :     }
     458                 :      49308 :     G = qf_apply_ZM(G, Ker);
     459         [ +  + ]:      49308 :     U = U? RgM_mul(U,Ker): Ker;
     460                 :            : 
     461                 :            :     /* 1st case: dimKer < vp */
     462                 :            :     /* then the kernel mod p contains a kernel mod p^2 */
     463         [ +  + ]:      49308 :     if (dimKer < vp)
     464                 :            :     {
     465         [ -  + ]:       2324 :       if (DEBUGLEVEL >= 4) err_printf("    case 1: dimker < vp\n");
     466         [ +  + ]:       2324 :       if (dimKer == 1)
     467                 :            :       {
     468                 :            :         long j;
     469                 :       1827 :         gel(G,1) = ZC_Z_divexact(gel(G,1), p);
     470         [ +  + ]:      11788 :         for (j = 1; j<=n; j++) gcoeff(G,1,j) = diviiexact(gcoeff(G,1,j), p);
     471                 :       1827 :         gel(U,1) = RgC_Rg_div(gel(U,1), p);
     472                 :       1827 :         E[i] -= 2;
     473                 :            :       }
     474                 :            :       else
     475                 :            :       {
     476                 :        497 :         GEN A,B,C, K2 = ZsymM_Z_divexact(principal_minor(G,dimKer),p);
     477                 :            :         long j, dimKer2;
     478                 :        497 :         K2 = kermodp(K2, p, &dimKer2);
     479         [ +  + ]:       1022 :         for (j = dimKer2+1; j <= dimKer; j++) gel(K2,j) = ZC_Z_mul(gel(K2,j),p);
     480                 :            :         /* Write G = [A,B;B~,C] and apply [K2,0;0,p*Id]/p by blocks */
     481                 :        497 :         blocks4(G, dimKer,n, &A,&B,&C);
     482                 :        497 :         A = ZsymM_Z_divexact(qf_apply_ZM(A,K2), sqri(p));
     483                 :        497 :         B = ZM_Z_divexact(ZM_transmul(B,K2), p);
     484                 :        497 :         G = shallowmatconcat(mkmat2(mkcol2(A,B),
     485                 :            :                                     mkcol2(shallowtrans(B), C)));
     486                 :            :         /* U *= [K2,0;0,Id] */
     487                 :        497 :         U = shallowconcat(RgM_Rg_div(RgM_mul(vecslice(U,1,dimKer),K2), p),
     488                 :            :                           vecslice(U,dimKer+1,n));
     489                 :        497 :         E[i] -= 2*dimKer2;
     490                 :            :       }
     491                 :       2324 :       i--; continue; /* same p */
     492                 :            :     }
     493                 :            : 
     494                 :            :    /* vp = dimKer
     495                 :            :     * 2nd case: kernel has dim >= 2 and contains an element of norm 0 mod p^2
     496                 :            :     * search for an element of norm p^2... in the kernel */
     497                 :      46984 :     sol = NULL;
     498         [ +  + ]:      46984 :     if (dimKer > 2) {
     499         [ -  + ]:      17969 :       if (DEBUGLEVEL >= 4) err_printf("    case 2.1\n");
     500                 :      17969 :       dimKer = 3;
     501                 :      17969 :       sol = qfsolvemodp(ZsymM_Z_divexact(principal_minor(G,3),p),  p);
     502                 :      17969 :       sol = FpC_red(sol, p);
     503                 :            :     }
     504         [ +  + ]:      29015 :     else if (dimKer == 2)
     505                 :            :     {
     506                 :      15526 :       GEN a = modii(diviiexact(gcoeff(G,1,1),p), p);
     507                 :      15526 :       GEN b = modii(diviiexact(gcoeff(G,1,2),p), p);
     508                 :      15526 :       GEN c = diviiexact(gcoeff(G,2,2),p);
     509                 :      15526 :       GEN di= modii(subii(sqri(b), mulii(a,c)), p);
     510         [ +  + ]:      15526 :       if (kronecker(di,p) >= 0)
     511                 :            :       {
     512         [ -  + ]:      15477 :         if (DEBUGLEVEL >= 4) err_printf("    case 2.2\n");
     513         [ +  + ]:      15477 :         sol = signe(a)? mkcol2(Fp_sub(Fp_sqrt(di,p), b, p), a): vec_ei(2,1);
     514                 :            :       }
     515                 :            :     }
     516         [ +  + ]:      46984 :     if (sol)
     517                 :            :     {
     518                 :            :       long j;
     519                 :      33446 :       sol = FpC_center(sol, p, shifti(p,-1));
     520                 :      33446 :       sol = Q_primpart(sol);
     521         [ -  + ]:      33446 :       if (DEBUGLEVEL >= 4) err_printf("    sol = %Ps\n", sol);
     522                 :      33446 :       Ker = completebasis(vecextend(sol,n), 1);
     523         [ +  + ]:     262367 :       for(j=1; j<n; j++) gel(Ker,j) = ZC_Z_mul(gel(Ker,j), p);
     524                 :      33446 :       G = ZsymM_Z_divexact(qf_apply_ZM(G, Ker), sqri(p));
     525                 :      33446 :       U = RgM_Rg_div(RgM_mul(U,Ker), p);
     526                 :      33446 :       E[i] -= 2;
     527                 :      33446 :       i--; continue; /* same p */
     528                 :            :     }
     529                 :            :     /* Now 1 <= vp = dimKer <= 2 and kernel contains no vector with norm p^2 */
     530                 :            :     /* exchanging kernel and image makes minimization easier ? */
     531                 :      13538 :     m = (n-3)/2;
     532         [ +  + ]:      13538 :     d = ZM_det(G); if (odd(m)) d = negi(d);
     533 [ +  + ][ +  + ]:      13538 :     if ((vp==1 && kronecker(gmod(gdiv(negi(d), gcoeff(G,1,1)),p), p) >= 0)
     534 [ +  + ][ +  + ]:       4949 :      || (vp==2 && odd(n) && n >= 5)
                 [ -  + ]
     535 [ +  + ][ +  - ]:       4935 :      || (vp==2 && !odd(n) && kronecker(modii(diviiexact(d,sqri(p)), p),p) < 0))
                 [ +  + ]
     536                 :            :     {
     537                 :            :       long j;
     538         [ -  + ]:       8624 :       if (DEBUGLEVEL >= 4) err_printf("    case 3\n");
     539                 :       8624 :       Ker = matid(n);
     540         [ +  + ]:      61572 :       for (j = dimKer+1; j <= n; j++) gcoeff(Ker,j,j) = p;
     541                 :       8624 :       G = ZsymM_Z_divexact(qf_apply_ZM(G, Ker), p);
     542                 :       8624 :       U = RgM_mul(U,Ker);
     543                 :       8624 :       E[i] -= 2*dimKer-n;
     544                 :       8624 :       i--; continue; /* same p */
     545                 :            :     }
     546                 :            : 
     547                 :            :     /* Minimization was not possible so far. */
     548                 :            :     /* If n == 3 or 4, this proves the local non-solubility at p. */
     549 [ +  + ][ -  + ]:       4914 :     if (n == 3 || n == 4)
     550                 :            :     {
     551         [ -  + ]:        483 :       if (DEBUGLEVEL >= 1) err_printf(" no local solution at %Ps\n",p);
     552                 :        483 :       return(p);
     553                 :            :     }
     554                 :       4431 :     vectrunc_append(faP, p);
     555                 :       4431 :     vecsmalltrunc_append(faE, vp);
     556                 :            :   }
     557         [ +  + ]:       7728 :   if (!U) U = matid(n);
     558                 :            :   else
     559                 :            :   { /* apply LLL to avoid coefficient explosion */
     560                 :       6671 :     aux = lllint(Q_primpart(U));
     561                 :       6671 :     G = qf_apply_ZM(G,aux);
     562                 :       6671 :     U = RgM_mul(U,aux);
     563                 :            :   }
     564                 :       8211 :   return mkvec4(G, U, faP, faE);
     565                 :            : }
     566                 :            : 
     567                 :            : /* CLASS GROUP COMPUTATIONS */
     568                 :            : 
     569                 :            : /* Compute the square root of the quadratic form q of discriminant D. Not
     570                 :            :  * fully implemented; it only works for detqfb squarefree except at 2, where
     571                 :            :  * the valuation is 2 or 3.
     572                 :            :  * mkmat2(P,zv_to_ZV(E)) = factor(2*abs(det q)) */
     573                 :            : static GEN
     574                 :       2268 : qfbsqrt(GEN D, GEN q, GEN P, GEN E)
     575                 :            : {
     576                 :       2268 :   GEN a = gel(q,1), b = shifti(gel(q,2),-1), c = gel(q,3), mb = negi(b);
     577                 :            :   GEN m,n, aux,Q1,M, A,B,C;
     578                 :       2268 :   GEN d = subii(mulii(a,c), sqri(b));
     579                 :            :   long i;
     580                 :            : 
     581                 :            :   /* 1st step: solve m^2 = a (d), m*n = -b (d), n^2 = c (d) */
     582                 :       2268 :   m = n = mkintmod(gen_0,gen_1);
     583                 :       2268 :   E[1] -= 3;
     584         [ +  + ]:       9800 :   for (i = 1; i < lg(P); i++)
     585                 :            :   {
     586                 :       7532 :     GEN p = gel(P,i), N, M;
     587         [ +  + ]:       7532 :     if (!E[i]) continue;
     588         [ +  + ]:       7462 :     if (dvdii(a,p)) {
     589                 :       1806 :       aux = Fp_sqrt(c, p);
     590                 :       1806 :       N = aux;
     591                 :       1806 :       M = Fp_div(mb, aux, p);
     592                 :            :     } else {
     593                 :       5656 :       aux = Fp_sqrt(a, p);
     594                 :       5656 :       M = aux;
     595                 :       5656 :       N = Fp_div(mb, aux, p);
     596                 :            :     }
     597                 :       7462 :     n = chinese(n, mkintmod(N,p));
     598                 :       7462 :     m = chinese(m, mkintmod(M,p));
     599                 :            :   }
     600                 :       2268 :   m = centerlift(m);
     601                 :       2268 :   n = centerlift(n);
     602         [ -  + ]:       2268 :   if (DEBUGLEVEL >=4) err_printf("    [m,n] = [%Ps, %Ps]\n",m,n);
     603                 :            : 
     604                 :            :   /* 2nd step: build Q1, with det=-1 such that Q1(x,y,0) = G(x,y) */
     605                 :       2268 :   A = diviiexact(subii(sqri(n),c), d);
     606                 :       2268 :   B = diviiexact(addii(b, mulii(m,n)), d);
     607                 :       2268 :   C = diviiexact(subii(sqri(m), a), d);
     608                 :       2268 :   Q1 = mkmat3(mkcol3(A,B,n), mkcol3(B,C,m), mkcol3(n,m,d));
     609                 :       2268 :   Q1 = gneg(adj(Q1));
     610                 :            : 
     611                 :            :   /* 3rd step: reduce Q1 to [0,0,-1;0,1,0;-1,0,0] */
     612                 :       2268 :   M = qflllgram_indefgoon2(Q1);
     613         [ +  + ]:       2268 :   if (signe(gel(M,1)) < 0) M = ZC_neg(M);
     614                 :       2268 :   a = gel(M,1);
     615                 :       2268 :   b = gel(M,2);
     616                 :       2268 :   c = gel(M,3);
     617         [ +  + ]:       2268 :   if (mpodd(a))
     618                 :       2212 :     return qfb(D, a, shifti(b,1), shifti(c,1));
     619                 :            :   else
     620                 :       2268 :     return qfb(D, c, shifti(negi(b),1), shifti(a,1));
     621                 :            : }
     622                 :            : 
     623                 :            : /* \prod gen[i]^e[i] as a Qfb, e in {0,1}^n non-zero */
     624                 :            : static GEN
     625                 :       3955 : qfb_factorback(GEN D, GEN gen, GEN e)
     626                 :            : {
     627                 :       3955 :   GEN q = NULL;
     628                 :       3955 :   long j, l = lg(gen), n = 0;
     629         [ +  + ]:      13482 :   for (j = 1; j < l; j++)
     630 [ +  + ][ +  + ]:       9527 :     if (e[j]) { n++; q = q? qfbcompraw(q, gel(gen,j)): gel(gen,j); }
     631         [ +  + ]:       3955 :   return (n <= 1)? q: qfbreduce(D, q);
     632                 :            : }
     633                 :            : 
     634                 :            : /* unit form, assuming 4 | D */
     635                 :            : static GEN
     636                 :        973 : id(GEN D)
     637                 :        973 : { return mkmat2(mkcol2(gen_1,gen_0),mkcol2(gen_0,shifti(negi(D),-2))); }
     638                 :            : 
     639                 :            : /* Shanks/Bosma-Stevenhagen algorithm to compute the 2-Sylow of the class
     640                 :            :  * group of discriminant D. Only works for D = fundamental discriminant.
     641                 :            :  * When D = 1(4), work with 4D.
     642                 :            :  * P2D,E2D = factor(abs(2*D))
     643                 :            :  * Pm2D = factor(-abs(2*D))[,1].
     644                 :            :  * Return a form having Witt invariants W at Pm2D */
     645                 :            : static GEN
     646                 :       2660 : quadclass2(GEN D, GEN P2D, GEN E2D, GEN Pm2D, GEN W, int n_is_4)
     647                 :            : {
     648                 :            :   GEN gen, Wgen, U2;
     649                 :            :   long i, n, r, m, vD;
     650                 :            : 
     651         [ +  + ]:       2660 :   if (mpodd(D))
     652                 :            :   {
     653                 :        210 :     D = shifti(D,2);
     654                 :        210 :     E2D = shallowcopy(E2D);
     655                 :        210 :     E2D[1] = 3;
     656                 :            :   }
     657         [ +  + ]:       2660 :   if (zv_equal0(W)) return id(D);
     658                 :            : 
     659                 :       1806 :   n = lg(Pm2D)-1; /* >= 3 since W != 0 */
     660                 :       1806 :   r = n-3;
     661         [ +  + ]:       1806 :   m = (signe(D)>0)? r+1: r;
     662                 :            :   /* n=4: look among forms of type q or 2*q, since Q can be imprimitive */
     663         [ +  + ]:       1806 :   U2 = n_is_4? mkmat(hilberts(gen_2, D, Pm2D, lg(Pm2D))): NULL;
     664 [ +  + ][ +  + ]:       1806 :   if (U2 && zv_equal(gel(U2,1),W)) return gmul2n(id(D),1);
     665                 :            : 
     666                 :       1687 :   gen = cgetg(m+1, t_VEC);
     667         [ +  + ]:       4571 :   for (i = 1; i <= m; i++) /* no need to look at Pm2D[1]=-1, nor Pm2D[2]=2 */
     668                 :            :   {
     669                 :       2884 :     GEN p = gel(Pm2D,i+2), d;
     670                 :       2884 :     long vp = Z_pvalrem(D,p, &d);
     671                 :       2884 :     gel(gen,i) = qfb(D, powiu(p,vp), gen_0, negi(shifti(d,-2)));
     672                 :            :   }
     673                 :       1687 :   vD = Z_lval(D,2);  /* = 2 or 3 */
     674 [ +  + ][ +  + ]:       1687 :   if (vD == 2 && smodis(D,16) != 4)
     675                 :            :   {
     676                 :        119 :     GEN q2 = qfb(D, gen_2,gen_2, shifti(subsi(4,D),-3));
     677                 :        119 :     m++; r++; gen = shallowconcat(gen, mkvec(q2));
     678                 :            :   }
     679         [ +  + ]:       1687 :   if (vD == 3)
     680                 :            :   {
     681                 :       1463 :     GEN q2 = qfb(D, gen_2,gen_0, negi(shifti(D,-3)));
     682                 :       1463 :     m++; r++; gen = shallowconcat(gen, mkvec(q2));
     683                 :            :   }
     684         [ -  + ]:       1687 :   if (!r) return id(D);
     685                 :       1687 :   Wgen = qflocalinvariants(gen,Pm2D);
     686                 :            :   for(;;)
     687                 :            :   {
     688                 :       3745 :     GEN Wgen2, gen2, Ker, indexim = gel(Flm_indexrank(Wgen,2), 2);
     689                 :            :     long dKer;
     690         [ +  + ]:       3745 :     if (lg(indexim)-1 >= r)
     691                 :            :     {
     692                 :       1687 :       GEN W2 = Wgen, V;
     693         [ +  + ]:       1687 :       if (lg(indexim) < lg(Wgen)) W2 = vecpermute(Wgen,indexim);
     694         [ +  + ]:       1687 :       if (U2) W2 = shallowconcat(W2,U2);
     695                 :       1687 :       V = Flm_Flc_invimage(W2, W,2);
     696         [ +  - ]:       1687 :       if (V) {
     697                 :       1687 :         GEN Q = qfb_factorback(D, vecpermute(gen,indexim), V);
     698                 :       1687 :         Q = gtomat(Q);
     699 [ +  + ][ +  + ]:       1687 :         if (U2 && V[lg(V)-1]) Q = gmul2n(Q,1);
     700                 :       1687 :         return Q;
     701                 :            :       }
     702                 :            :     }
     703                 :       2058 :     Ker = Flm_ker(Wgen,2); dKer = lg(Ker)-1;
     704                 :       2058 :     gen2 = cgetg(m+1, t_VEC);
     705                 :       2058 :     Wgen2 = cgetg(m+1, t_MAT);
     706         [ +  + ]:       4326 :     for (i = 1; i <= dKer; i++)
     707                 :            :     {
     708                 :       2268 :       GEN q = qfb_factorback(D, gen, gel(Ker,i));
     709                 :       2268 :       q = qfbsqrt(D,q,P2D,E2D);
     710                 :       2268 :       gel(gen2,i) = q;
     711                 :       2268 :       gel(Wgen2,i) = gel(qflocalinvariants(q,Pm2D), 1);
     712                 :            :     }
     713         [ +  + ]:       4438 :     for (; i <=m; i++)
     714                 :            :     {
     715                 :       2380 :       long j = indexim[i-dKer];
     716                 :       2380 :       gel(gen2,i) = gel(gen,j);
     717                 :       2380 :       gel(Wgen2,i) = gel(Wgen,j);
     718                 :            :     }
     719                 :       2058 :     gen = gen2; Wgen = Wgen2;
     720                 :       4718 :   }
     721                 :            : }
     722                 :            : 
     723                 :            : /* QUADRATIC EQUATIONS */
     724                 :            : /* is x*y = -1 ? */
     725                 :            : static int
     726                 :       4658 : both_pm1(GEN x, GEN y)
     727 [ +  + ][ +  + ]:       4658 : { return is_pm1(x) && is_pm1(y) && signe(x) == -signe(y); }
                 [ +  + ]
     728                 :            : 
     729                 :            : /* Try to solve G = 0 with small coefficients. This is proved to work if
     730                 :            :  * -  det(G) = 1, dim <= 6 and G is LLL reduced
     731                 :            :  * Returns G if no solution is found.
     732                 :            :  * Exit with a norm 0 vector if one such is found.
     733                 :            :  * If base == 1 and norm 0 is obtained, returns [H~*G*H,H] where
     734                 :            :  * the 1st column of H is a norm 0 vector */
     735                 :            : static GEN
     736                 :      21203 : qftriv(GEN G, GEN R, long base)
     737                 :            : {
     738                 :      21203 :   long n = lg(G)-1, i;
     739                 :            :   GEN s, H;
     740                 :            : 
     741                 :            :   /* case 1: A basis vector is isotropic */
     742         [ +  + ]:      78155 :   for (i = 1; i <= n; i++)
     743         [ +  + ]:      67675 :     if (!signe(gcoeff(G,i,i)))
     744                 :            :     {
     745         [ +  + ]:      10723 :       if (!base) return col_ei(n,i);
     746                 :      10107 :       H = matid(n); swap(gel(H,1), gel(H,i));
     747                 :      10107 :       return mkvec2(qf_apply_tau(G,1,i),H);
     748                 :            :     }
     749                 :            :   /* case 2: G has a block +- [1,0;0,-1] on the diagonal */
     750         [ +  + ]:      46512 :   for (i = 2; i <= n; i++)
     751 [ +  + ][ +  + ]:      39063 :     if (!signe(gcoeff(G,i-1,i)) && both_pm1(gcoeff(G,i-1,i-1),gcoeff(G,i,i)))
     752                 :            :     {
     753                 :       3031 :       s = col_ei(n,i); gel(s,i-1) = gen_m1;
     754         [ +  + ]:       3031 :       if (!base) return s;
     755                 :       2940 :       H = matid(n); gel(H,i) = gel(H,1); gel(H,1) = s;
     756                 :       2940 :       return mkvec2(qf_apply_ZM(G,H),H);
     757                 :            :     }
     758         [ +  + ]:       7449 :   if (!R) return G; /* fail */
     759                 :            :   /* case 3: a principal minor is 0 */
     760                 :        239 :   s = keri(principal_minor(G, itos(R)));
     761                 :        239 :   s = vecextend(Q_primpart(gel(s,1)), n);
     762         [ +  + ]:        239 :   if (!base) return s;
     763                 :        204 :   H = completebasis(s, 0);
     764                 :        204 :   gel(H,n) = ZC_neg(gel(H,1)); gel(H,1) = s;
     765                 :      21203 :   return mkvec2(qf_apply_ZM(G,H),H);
     766                 :            : }
     767                 :            : 
     768                 :            : /* p a prime number, G 3x3 symmetric. Finds X!=0 such that X^t G X = 0 mod p.
     769                 :            :  * Allow returning a shorter X: to be completed with 0s. */
     770                 :            : static GEN
     771                 :      17969 : qfsolvemodp(GEN G, GEN p)
     772                 :            : {
     773                 :            :   GEN a,b,c,d,e,f, v1,v2,v3,v4,v5, x1,x2,x3,N1,N2,N3,s,r;
     774                 :            : 
     775                 :            :   /* principal_minor(G,3) = [a,b,d; b,c,e; d,e,f] */
     776                 :      17969 :   a = modii(gcoeff(G,1,1), p);
     777         [ +  + ]:      17969 :   if (!signe(a)) return mkcol(gen_1);
     778                 :      15275 :   v1 = a;
     779                 :      15275 :   b = modii(gcoeff(G,1,2), p);
     780                 :      15275 :   c = modii(gcoeff(G,2,2), p);
     781                 :      15275 :   v2 = modii(subii(mulii(a,c), sqri(b)), p);
     782         [ +  + ]:      15275 :   if (!signe(v2)) return mkcol2(Fp_neg(b,p), a);
     783                 :      12194 :   d = modii(gcoeff(G,1,3), p);
     784                 :      12194 :   e = modii(gcoeff(G,2,3), p);
     785                 :      12194 :   f = modii(gcoeff(G,3,3), p);
     786                 :      12194 :   v4 = modii(subii(mulii(c,d), mulii(e,b)), p);
     787                 :      12194 :   v5 = modii(subii(mulii(a,e), mulii(d,b)), p);
     788                 :      12194 :   v3 = subii(mulii(v2,f), addii(mulii(v4,d), mulii(v5,e))); /* det(G) */
     789                 :      12194 :   v3 = modii(v3, p);
     790                 :      12194 :   N1 =  Fp_neg(v2,  p);
     791                 :      12194 :   x3 = mkcol3(v4, v5, N1);
     792         [ +  + ]:      12194 :   if (!signe(v3)) return x3;
     793                 :            : 
     794                 :            :   /* now, solve in dimension 3... reduction to the diagonal case: */
     795                 :      10577 :   x1 = mkcol3(gen_1, gen_0, gen_0);
     796                 :      10577 :   x2 = mkcol3(negi(b), a, gen_0);
     797         [ +  + ]:      10577 :   if (kronecker(N1,p) == 1) return ZC_lincomb(Fp_sqrt(N1,p),gen_1,x1,x2);
     798                 :       4212 :   N2 = Fp_div(Fp_neg(v3,p), v1, p);
     799         [ +  + ]:       4212 :   if (kronecker(N2,p) == 1) return ZC_lincomb(Fp_sqrt(N2,p),gen_1,x2,x3);
     800                 :       2093 :   N3 = Fp_mul(v2, N2, p);
     801         [ +  + ]:       2093 :   if (kronecker(N3,p) == 1) return ZC_lincomb(Fp_sqrt(N3,p),gen_1,x1,x3);
     802                 :       1106 :   r = gen_1;
     803                 :            :   for(;;)
     804                 :            :   {
     805                 :       1925 :     s = Fp_sub(gen_1, Fp_mul(N1,Fp_sqr(r,p),p), p);
     806         [ +  + ]:       1925 :     if (kronecker(s, p) <= 0) break;
     807                 :        819 :     r = randomi(p);
     808                 :        819 :   }
     809                 :       1106 :   s = Fp_sqrt(Fp_div(s,N3,p), p);
     810                 :      17969 :   return ZC_add(x1, ZC_lincomb(r,s,x2,x3));
     811                 :            : }
     812                 :            : 
     813                 :            : /* assume G square integral */
     814                 :            : static void
     815                 :       4333 : check_symmetric(GEN G)
     816                 :            : {
     817                 :       4333 :   long i,j, l = lg(G);
     818         [ +  + ]:      27804 :   for (i = 1; i < l; i++)
     819         [ +  + ]:      82012 :     for(j = 1; j < i; j++)
     820         [ +  + ]:      58541 :       if (!equalii(gcoeff(G,i,j), gcoeff(G,j,i)))
     821                 :          7 :         pari_err_TYPE("qfsolve [not symmetric]",G);
     822                 :       4326 : }
     823                 :            : 
     824                 :            : /* Given a square matrix G of dimension n >= 1, */
     825                 :            : /* solves over Z the quadratic equation X^tGX = 0. */
     826                 :            : /* G is assumed to have integral coprime coefficients. */
     827                 :            : /* The solution might be a vectorv or a matrix. */
     828                 :            : /* If no solution exists, returns an integer, that can */
     829                 :            : /* be a prime p such that there is no local solution at p, */
     830                 :            : /* or -1 if there is no real solution, */
     831                 :            : /* or 0 in some rare cases. */
     832                 :            : static  GEN
     833                 :       4291 : qfsolve_i(GEN G)
     834                 :            : {
     835                 :            :   GEN M, signG, Min, U, G1, M1, G2, M2, solG2, P, E;
     836                 :            :   GEN solG1, sol, Q, d, dQ, detG2, fam2detG;
     837                 :            :   long n, np, codim, dim;
     838                 :            :   int fail;
     839                 :            : 
     840         [ -  + ]:       4291 :   if (typ(G) != t_MAT) pari_err_TYPE("qfsolve", G);
     841                 :       4291 :   n = lg(G)-1;
     842         [ -  + ]:       4291 :   if (n == 0) pari_err_DOMAIN("qfsolve", "dimension" , "=", gen_0, G);
     843         [ -  + ]:       4291 :   if (n != nbrows(G)) pari_err_DIM("qfsolve");
     844                 :       4291 :   G = Q_primpart(G); RgM_check_ZM(G, "qfsolve");
     845                 :       4291 :   check_symmetric(G);
     846                 :            : 
     847                 :            :   /* Trivial case: det = 0 */
     848                 :       4284 :   d = ZM_det(G);
     849         [ +  + ]:       4284 :   if (!signe(d))
     850                 :            :   {
     851         [ +  - ]:          7 :     if (n == 1) return mkcol(gen_1);
     852                 :          0 :     sol = keri(G);
     853         [ #  # ]:          0 :     if (lg(sol) == 2) sol = gel(sol,1);
     854                 :          0 :     return sol;
     855                 :            :   }
     856                 :            : 
     857                 :            :   /* Small dimension: n <= 2 */
     858         [ +  + ]:       4277 :   if (n == 1) return gen_m1;
     859         [ +  + ]:       4270 :   if (n == 2)
     860                 :            :   {
     861                 :         21 :     GEN t, a =  gcoeff(G,1,1);
     862         [ +  + ]:         21 :     if (!signe(a)) return mkcol2(gen_1, gen_0);
     863         [ -  + ]:         14 :     if (signe(d) > 0) return gen_m1; /* no real solution */
     864         [ +  + ]:         14 :     if (!Z_issquareall(negi(d), &t)) return gen_m2;
     865                 :         21 :     return mkcol2(subii(t,gcoeff(G,1,2)), a);
     866                 :            :   }
     867                 :            : 
     868                 :            :   /* 1st reduction of the coefficients of G */
     869                 :       4249 :   M = qflllgram_indef(G,0,&fail);
     870         [ +  + ]:       4249 :   if (typ(M) == t_COL) return M;
     871                 :       4235 :   G = gel(M,1);
     872                 :       4235 :   M = gel(M,2);
     873                 :            : 
     874                 :            :   /* real solubility */
     875                 :       4235 :   signG = ZV_to_zv(qfsign(G));
     876                 :            :   {
     877                 :       4235 :     long r =  signG[1], s = signG[2];
     878 [ +  + ][ +  + ]:       4235 :     if (!r || !s) return gen_m1;
     879         [ +  + ]:       4165 :     if (r < s) { G = ZM_neg(G); signG = mkvecsmall2(s,r);  }
     880                 :            :   }
     881                 :            : 
     882                 :            :   /* factorization of the determinant */
     883                 :       4165 :   fam2detG = Z_factor( absi(d) );
     884                 :       4165 :   P = gel(fam2detG,1);
     885                 :       4165 :   E = ZV_to_zv(gel(fam2detG,2));
     886                 :            :   /* P,E = factor(|det(G)|) */
     887                 :            : 
     888                 :            :   /* Minimization and local solubility */
     889                 :       4165 :   Min = qfminimize(G, P, E);
     890         [ +  + ]:       4165 :   if (typ(Min) == t_INT) return Min;
     891                 :            : 
     892                 :       3682 :   M = RgM_mul(M, gel(Min,2));
     893                 :       3682 :   G = gel(Min,1);
     894                 :       3682 :   P = gel(Min,3);
     895                 :       3682 :   E = gel(Min,4);
     896                 :            :   /* P,E = factor(|det(G))| */
     897                 :            : 
     898                 :            :   /* Now, we know that local solutions exist (except maybe at 2 if n==4)
     899                 :            :    * if n==3, det(G) = +-1
     900                 :            :    * if n==4, or n is odd, det(G) is squarefree.
     901                 :            :    * if n>=6, det(G) has all its valuations <=2. */
     902                 :            : 
     903                 :            :   /* Reduction of G and search for trivial solutions. */
     904                 :            :   /* When |det G|=1, such trivial solutions always exist. */
     905                 :       3682 :   U = qflllgram_indef(G,0,&fail);
     906         [ +  + ]:       3682 :   if(typ(U) == t_COL) return Q_primpart(RgM_RgC_mul(M,U));
     907                 :       2954 :   G = gel(U,1);
     908                 :       2954 :   M = RgM_mul(M, gel(U,2));
     909                 :            :   /* P,E = factor(|det(G))| */
     910                 :            : 
     911                 :            :   /* If n >= 6 is even, need to increment the dimension by 1 to suppress all
     912                 :            :    * squares from det(G) */
     913                 :       2954 :   np = lg(P)-1;
     914 [ +  + ][ +  + ]:       2954 :   if (n < 6 || odd(n) || !np)
                 [ -  + ]
     915                 :            :   {
     916                 :       1568 :     codim = 0;
     917                 :       1568 :     G1 = G;
     918                 :       1568 :     M1 = NULL;
     919                 :            :   }
     920                 :            :   else
     921                 :            :   {
     922                 :            :     GEN aux;
     923                 :            :     long i;
     924                 :       1386 :     codim = 1; n++;
     925                 :            :     /* largest square divisor of d */
     926                 :       1386 :     aux = gen_1;
     927         [ +  + ]:       6524 :     for (i = 1; i <= np; i++)
     928         [ +  + ]:       5138 :       if (E[i] == 2) { aux = mulii(aux, gel(P,i)); E[i] = 3; }
     929                 :            :     /* Choose sign(aux) so as to balance the signature of G1 */
     930         [ +  + ]:       1386 :     if (signG[1] > signG[2])
     931                 :            :     {
     932                 :        546 :       signG[2]++;
     933                 :        546 :       aux = negi(aux);
     934                 :            :     }
     935                 :            :     else
     936                 :        840 :       signG[1]++;
     937                 :       1386 :     G1 = shallowmatconcat(diagonal_shallow(mkvec2(G,aux)));
     938                 :            :     /* P,E = factor(|det G1|) */
     939                 :       1386 :     Min = qfminimize(G1, P, E);
     940                 :       1386 :     G1 = gel(Min,1);
     941                 :       1386 :     M1 = gel(Min,2);
     942                 :       1386 :     P = gel(Min,3);
     943                 :       1386 :     E = gel(Min,4);
     944                 :       1386 :     np = lg(P)-1;
     945                 :            :   }
     946                 :            : 
     947                 :            :   /* now, d is squarefree */
     948         [ +  + ]:       2954 :   if (!np)
     949                 :            :   { /* |d| = 1 */
     950                 :        259 :      G2 = G1;
     951                 :        259 :      M2 = NULL;
     952                 :            :   }
     953                 :            :   else
     954                 :            :   { /* |d| > 1: increment dimension by 2 */
     955                 :            :     GEN factdP, factdE, W;
     956                 :            :     long i, lfactdP;
     957                 :       2695 :     codim += 2;
     958                 :       2695 :     d = ZV_prod(P); /* d = abs(matdet(G1)); */
     959         [ +  + ]:       2695 :     if (odd(signG[2])) togglesign_safe(&d); /* d = matdet(G1); */
     960                 :            :     /* solubility at 2 (this is the only remaining bad prime). */
     961 [ +  + ][ +  + ]:       2695 :     if (n == 4 && smodis(d,8) == 1 && qflocalinvariant(G,gen_2) == 1)
                 [ +  + ]
     962                 :         35 :       return gen_2;
     963                 :            : 
     964         [ +  + ]:       2660 :     P = shallowconcat(mpodd(d)? mkvec2(NULL,gen_2): mkvec(NULL), P);
     965                 :            :     /* build a binary quadratic form with given Witt invariants */
     966                 :       2660 :     W = const_vecsmall(lg(P)-1, 0);
     967                 :            :     /* choose signature of Q (real invariant and sign of the discriminant) */
     968                 :       2660 :     dQ = absi(d);
     969         [ +  + ]:       2660 :     if (signG[1] > signG[2]) togglesign_safe(&dQ); /* signQ = [2,0]; */
     970 [ +  + ][ +  + ]:       2660 :     if (n == 4 && smodis(dQ,4) != 1) dQ = shifti(dQ,2);
     971         [ +  + ]:       2660 :     if (n >= 5) dQ = shifti(dQ,3);
     972                 :            : 
     973                 :            :     /* p-adic invariants */
     974         [ +  + ]:       2660 :     if (n == 4)
     975                 :            :     {
     976                 :        665 :       GEN t = qflocalinvariants(ZM_neg(G1),P);
     977         [ +  + ]:       2667 :       for (i = 3; i < lg(P); i++) W[i] = ucoeff(t,i,1);
     978                 :            :     }
     979                 :            :     else
     980                 :            :     {
     981         [ +  + ]:       1995 :       long s = signe(dQ) == signe(d)? 1: -1;
     982                 :            :       GEN t;
     983         [ +  + ]:       1995 :       if (odd((n-3)/2)) s = -s;
     984         [ +  + ]:       1995 :       t = s > 0? utoipos(8): utoineg(8);
     985         [ +  + ]:       6013 :       for (i = 3; i < lg(P); i++)
     986                 :       4018 :         W[i] = hilbertii(t, gel(P,i), gel(P,i)) > 0;
     987                 :            :     }
     988                 :            :     /* for p = 2, the choice is fixed from the product formula */
     989                 :       2660 :     W[2] = Flv_sum(W, 2);
     990                 :            : 
     991                 :            :     /* Construction of the 2-class group of discriminant dQ until some product
     992                 :            :      * of the generators gives the desired invariants. */
     993                 :       2660 :     factdP = vecsplice(P, 1); lfactdP =  lg(factdP);
     994                 :       2660 :     factdE = cgetg(lfactdP, t_VECSMALL);
     995         [ +  + ]:      11340 :     for (i = 1; i < lfactdP; i++) factdE[i] = Z_pval(dQ, gel(factdP,i));
     996                 :       2660 :     factdE[1]++;
     997                 :            :     /* factdP,factdE = factor(2|dQ|), P = factor(-2|dQ|)[,1] */
     998                 :       2660 :     Q = quadclass2(dQ, factdP,factdE, P, W, n == 4);
     999                 :            :     /* Build a form of dim=n+2 potentially unimodular */
    1000                 :       2660 :     G2 = shallowmatconcat(diagonal_shallow(mkvec2(G1,ZM_neg(Q))));
    1001                 :            :     /* Minimization of G2 */
    1002                 :       2660 :     detG2 = mulii(d, ZM_det(Q));
    1003         [ +  + ]:      11340 :     for (i = 1; i < lfactdP; i++) factdE[i] = Z_pval(detG2, gel(factdP,i));
    1004                 :            :     /* factdP,factdE = factor(|det G2|) */
    1005                 :       2660 :     Min = qfminimize(G2, factdP,factdE);
    1006                 :       2660 :     M2 = gel(Min,2);
    1007                 :       2660 :     G2 = gel(Min,1);
    1008                 :            :   }
    1009                 :            :   /* |det(G2)| = 1, find a totally isotropic subspace for G2 */
    1010                 :       2919 :   solG2 = qflllgram_indefgoon(G2);
    1011                 :            :   /* G2 must have a subspace of solutions of dimension > codim */
    1012                 :       2919 :   dim = codim+2;
    1013         [ +  + ]:       4270 :   while(gequal0(principal_minor(gel(solG2,1), dim))) dim ++;
    1014                 :       2919 :   solG2 = vecslice(gel(solG2,2), 1, dim-1);
    1015                 :            : 
    1016         [ +  + ]:       2919 :   if (!M2)
    1017                 :        259 :     solG1 = solG2;
    1018                 :            :   else
    1019                 :            :   { /* solution of G1 is simultaneously in solG2 and x[n+1] = x[n+2] = 0*/
    1020                 :            :     GEN K;
    1021                 :       2660 :     solG1 = RgM_mul(M2,solG2);
    1022                 :       2660 :     K = ker(rowslice(solG1,n+1,n+2));
    1023                 :       2660 :     solG1 = RgM_mul(rowslice(solG1,1,n), K);
    1024                 :            :   }
    1025         [ +  + ]:       2919 :   if (!M1)
    1026                 :       1533 :     sol = solG1;
    1027                 :            :   else
    1028                 :            :   { /* solution of G1 is simultaneously in solG2 and x[n] = 0 */
    1029                 :            :     GEN K;
    1030                 :       1386 :     sol = RgM_mul(M1,solG1);
    1031                 :       1386 :     K = ker(rowslice(sol,n,n));
    1032                 :       1386 :     sol = RgM_mul(rowslice(sol,1,n-1), K);
    1033                 :            :   }
    1034                 :       2919 :   sol = Q_primpart(RgM_mul(M, sol));
    1035         [ +  + ]:       2919 :   if (lg(sol) == 2) sol = gel(sol,1);
    1036                 :       4284 :   return sol;
    1037                 :            : }
    1038                 :            : GEN
    1039                 :       4291 : qfsolve(GEN G)
    1040                 :            : {
    1041                 :       4291 :   pari_sp av = avma;
    1042                 :       4291 :   return gerepilecopy(av, qfsolve_i(G));
    1043                 :            : }
    1044                 :            : 
    1045                 :            : /* G is a symmetric 3x3 matrix, and sol a solution of sol~*G*sol=0.
    1046                 :            :  * Returns a parametrization of the solutions with the good invariants,
    1047                 :            :  * as a matrix 3x3, where each line contains
    1048                 :            :  * the coefficients of each of the 3 quadratic forms.
    1049                 :            :  * If fl!=0, the fl-th form is reduced. */
    1050                 :            : GEN
    1051                 :         42 : qfparam(GEN G, GEN sol, long fl)
    1052                 :            : {
    1053                 :         42 :   pari_sp av = avma;
    1054                 :            :   GEN U, G1, G2, a, b, c, d, e;
    1055                 :         42 :   long n, tx = typ(sol);
    1056                 :            : 
    1057         [ -  + ]:         42 :   if (typ(G) != t_MAT) pari_err_TYPE("qfsolve", G);
    1058         [ -  + ]:         42 :   if (!is_vec_t(tx)) pari_err_TYPE("qfsolve", G);
    1059         [ +  + ]:         42 :   if (tx == t_VEC) sol = shallowtrans(sol);
    1060                 :         42 :   n = lg(G)-1;
    1061         [ -  + ]:         42 :   if (n == 0) pari_err_DOMAIN("qfsolve", "dimension" , "=", gen_0, G);
    1062 [ +  - ][ +  - ]:         42 :   if (n != nbrows(G) || n != 3 || lg(sol) != 4) pari_err_DIM("qfsolve");
                 [ -  + ]
    1063                 :         42 :   G = Q_primpart(G); RgM_check_ZM(G,"qfsolve");
    1064                 :         42 :   check_symmetric(G);
    1065                 :         42 :   sol = Q_primpart(sol); RgV_check_ZV(sol,"qfsolve");
    1066                 :            :   /* build U such that U[,3] = sol, and |det(U)| = 1 */
    1067                 :         42 :   U = completebasis(sol,1);
    1068                 :         42 :   G1 = qf_apply_ZM(G,U); /* G1 has a 0 at the bottom right corner */
    1069                 :         42 :   a = shifti(gcoeff(G1,1,2),1);
    1070                 :         42 :   b = shifti(negi(gcoeff(G1,1,3)),1);
    1071                 :         42 :   c = shifti(negi(gcoeff(G1,2,3)),1);
    1072                 :         42 :   d = gcoeff(G1,1,1);
    1073                 :         42 :   e = gcoeff(G1,2,2);
    1074                 :         42 :   G2 = mkmat3(mkcol3(b,gen_0,d), mkcol3(c,b,a), mkcol3(gen_0,c,e));
    1075                 :         42 :   sol = ZM_mul(U,G2);
    1076         [ +  + ]:         42 :   if (fl)
    1077                 :            :   {
    1078                 :         21 :     GEN v = row(sol,fl);
    1079                 :            :     int fail;
    1080                 :         21 :     a = gel(v,1);
    1081                 :         21 :     b = gmul2n(gel(v,2),-1);
    1082                 :         21 :     c = gel(v,3);
    1083                 :         21 :     U = qflllgram_indef(mkmat2(mkcol2(a,b),mkcol2(b,c)), 1, &fail);
    1084                 :         21 :     U = gel(U,2);
    1085                 :         21 :     a = gcoeff(U,1,1); b = gcoeff(U,1,2);
    1086                 :         21 :     c = gcoeff(U,2,1); d = gcoeff(U,2,2);
    1087                 :         21 :     U = mkmat3(mkcol3(sqri(a),mulii(a,c),sqri(c)),
    1088                 :            :                mkcol3(shifti(mulii(a,b),1), addii(mulii(a,d),mulii(b,c)),
    1089                 :            :                       shifti(mulii(c,d),1)),
    1090                 :            :                mkcol3(sqri(b),mulii(b,d),sqri(d)));
    1091                 :         21 :     sol = ZM_mul(sol,U);
    1092                 :            :   }
    1093                 :         42 :   return gerepileupto(av, sol);
    1094                 :            : }

Generated by: LCOV version 1.9