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 - qfsolve.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29375-be3334fb81) Lines: 621 632 98.3 %
Date: 2024-05-26 08:08:42 Functions: 34 34 100.0 %
Legend: Lines: hit not hit

          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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /* Copyright (C) 2014 Denis Simon
      16             :  * Adapted from qfsolve.gp v. 09/01/2014
      17             :  *   http://www.math.unicaen.fr/~simon/qfsolve.gp
      18             :  *
      19             :  * Author: Denis SIMON <simon@math.unicaen.fr> */
      20             : 
      21             : #include "pari.h"
      22             : #include "paripriv.h"
      23             : 
      24             : #define DEBUGLEVEL DEBUGLEVEL_qfsolve
      25             : 
      26             : /* LINEAR ALGEBRA */
      27             : /* complete by 0s, assume l-1 <= n */
      28             : static GEN
      29      115883 : vecextend(GEN v, long n)
      30             : {
      31      115883 :   long i, l = lg(v);
      32      115883 :   GEN w = cgetg(n+1, t_COL);
      33      357839 :   for (i = 1; i < l; i++) gel(w,i) = gel(v,i);
      34      385465 :   for (     ; i <=n; i++) gel(w,i) = gen_0;
      35      115883 :   return w;
      36             : }
      37             : 
      38             : /* Gives a unimodular matrix with the last column(s) equal to Mv.
      39             :  * Mv can be a column vector or a rectangular matrix.
      40             :  * redflag = 0 or 1. If redflag = 1, LLL-reduce the n-#v first columns. */
      41             : static GEN
      42      407482 : completebasis(GEN Mv, long redflag)
      43             : {
      44             :   GEN U;
      45             :   long m, n;
      46             : 
      47      407482 :   if (typ(Mv) == t_COL) Mv = mkmat(Mv);
      48      407482 :   n = lg(Mv)-1;
      49      407482 :   m = nbrows(Mv); /* m x n */
      50      407482 :   if (m == n) return Mv;
      51      406957 :   (void)ZM_hnfall_i(shallowtrans(Mv), &U, 0);
      52      406957 :   U = ZM_inv(shallowtrans(U),NULL);
      53      406957 :   if (m==1 || !redflag) return U;
      54             :   /* LLL-reduce the m-n first columns */
      55      122480 :   return shallowconcat(ZM_lll(vecslice(U,1,m-n), 0.99, LLL_INPLACE),
      56      122480 :                        vecslice(U, m-n+1,m));
      57             : }
      58             : 
      59             : /* Return U in GLn(Z) whose first d columns span Ker (M mod p). */
      60             : static GEN
      61      284739 : kermodp(GEN M, GEN p, long *d)
      62             : {
      63             :   long j, l;
      64             :   GEN K, B, U;
      65             : 
      66      284739 :   K = FpM_center(FpM_ker(M, p), p, shifti(p,-1));
      67      284739 :   B = completebasis(K,0);
      68      284739 :   l = lg(M); U = cgetg(l, t_MAT);
      69     1359421 :   for (j =  1; j < l; j++) gel(U,j) = gel(B,l-j);
      70      284739 :   *d = lg(K)-1; return U;
      71             : }
      72             : 
      73             : /* INVARIANTS COMPUTATIONS */
      74             : 
      75             : static GEN
      76       33246 : principal_minor(GEN G, long  i) { return matslice(G,1,i,1,i); }
      77             : static GEN
      78         707 : det_minors(GEN G)
      79             : {
      80         707 :   long i, l = lg(G);
      81         707 :   GEN v = cgetg(l+1, t_VEC);
      82         707 :   gel(v,1) = gen_1;
      83        3535 :   for (i = 2; i <= l; i++) gel(v,i) = ZM_det(principal_minor(G,i-1));
      84         707 :   return v;
      85             : }
      86             : 
      87             : static GEN
      88        7280 : hilberts(GEN a, GEN b, GEN P)
      89             : {
      90        7280 :   long i, lP = lg(P);
      91        7280 :   GEN v = cgetg(lP, t_VECSMALL);
      92       42217 :   for (i = 1; i < lP; i++) v[i] = hilbertii(a, b, gel(P,i)) < 0;
      93        7280 :   return v;
      94             : }
      95             : 
      96             : /* 4 | disc(q); special case of gtomat */
      97             : static GEN
      98        1666 : qfbmat(GEN q)
      99             : {
     100        1666 :   GEN a = gel(q,1), b = shifti(gel(q,2), -1), c = gel(q,3);
     101        1666 :   retmkmat22(a, b, b, c);
     102             : }
     103             : /* 2*qfbmat(q) */
     104             : static GEN
     105          21 : qfbmat2(GEN q)
     106             : {
     107          21 :   GEN a = shifti(gel(q,1), 1), b = gel(q,2), c = shifti(gel(q,3), 1);
     108          21 :   retmkmat22(a, b, b, c);
     109             : }
     110             : /* Given a symmetric matrix G over Z, compute the Witt invariant of G at p
     111             :  * v = det_minors(G) [G diagonalized]; assume that none of the v[i] is 0. */
     112             : static long
     113        2128 : witt(GEN v, GEN p)
     114             : {
     115        2128 :   long k = lg(v)-2, h = hilbertii(gel(v,k), gel(v,k+1), p);
     116        8512 :   for (k--; k >= 1; k--)
     117        6384 :     if (hilbertii(negi(gel(v,k)), gel(v,k+1),p) < 0) h = -h;
     118        2128 :   return h;
     119             : }
     120             : 
     121             : /* QUADRATIC FORM REDUCTION */
     122             : /* private version of qfgaussred:
     123             :  * - early abort if k-th principal minor is singular, return stoi(k)
     124             :  * - else return a matrix whose upper triangular part is qfgaussred(a) */
     125             : static GEN
     126       44370 : partialgaussred(GEN a)
     127             : {
     128       44370 :   long n = lg(a)-1, k;
     129       44370 :   a = RgM_shallowcopy(a);
     130      147689 :   for(k = 1; k < n; k++)
     131             :   {
     132      112733 :     GEN ak, p = gcoeff(a,k,k);
     133             :     long i, j;
     134      112733 :     if (isintzero(p)) return stoi(k);
     135      103319 :     ak = row(a, k);
     136      383421 :     for (i=k+1; i<=n; i++) gcoeff(a,k,i) = gdiv(gcoeff(a,k,i), p);
     137      383421 :     for (i=k+1; i<=n; i++)
     138             :     {
     139      280102 :       GEN c = gel(ak,i);
     140      280102 :       if (gequal0(c)) continue;
     141      842974 :       for (j=i; j<=n; j++)
     142      615182 :         gcoeff(a,i,j) = gsub(gcoeff(a,i,j), gmul(c,gcoeff(a,k,j)));
     143             :     }
     144             :   }
     145       34956 :   if (isintzero(gcoeff(a,n,n))) return stoi(n);
     146       34935 :   return a;
     147             : }
     148             : 
     149             : /* LLL-reduce a positive definite qf QD bounding the indefinite G, dim G > 1.
     150             :  * Then finishes by looking for trivial solution */
     151             : static GEN qftriv(GEN G, GEN z, long base);
     152             : static GEN
     153       44370 : qflllgram_indef(GEN G, long base, int *fail)
     154             : {
     155             :   GEN M, R, g, DM, S, dR;
     156       44370 :   long i, j, n = lg(G)-1;
     157             : 
     158       44370 :   *fail = 0;
     159       44370 :   R = partialgaussred(G);
     160       44370 :   if (typ(R) == t_INT) return qftriv(G, R, base);
     161       34935 :   R = Q_remove_denom(R, &dR); /* avoid rational arithmetic */
     162       34935 :   M = zeromatcopy(n,n);
     163       34935 :   DM = zeromatcopy(n,n);
     164      168478 :   for (i = 1; i <= n; i++)
     165             :   {
     166      133543 :     GEN d = absi_shallow(gcoeff(R,i,i));
     167      133543 :     if (dR) {
     168      109392 :       gcoeff(M,i,i) = dR;
     169      109392 :       gcoeff(DM,i,i) = mulii(d,dR);
     170             :     } else {
     171       24151 :       gcoeff(M,i,i) = gen_1;
     172       24151 :       gcoeff(DM,i,i) = d;
     173             :     }
     174      393529 :     for (j = i+1; j <= n; j++)
     175             :     {
     176      259986 :       gcoeff(M,i,j) = gcoeff(R,i,j);
     177      259986 :       gcoeff(DM,i,j) = mulii(d, gcoeff(R,i,j));
     178             :     }
     179             :   }
     180             :   /* G = M~*D*M, D diagonal, DM=|D|*M, g =  M~*|D|*M */
     181       34935 :   g = ZM_transmultosym(M,DM);
     182       34935 :   S = lllgramint(Q_primpart(g));
     183       34935 :   R = qftriv(qf_ZM_apply(G,S), NULL, base);
     184       34935 :   switch(typ(R))
     185             :   {
     186        6083 :     case t_COL: return ZM_ZC_mul(S,R);
     187       21153 :     case t_MAT: *fail = 1; return mkvec2(R, S);
     188        7699 :     default:
     189        7699 :       gel(R,2) = ZM_mul(S, gel(R,2));
     190        7699 :       return R;
     191             :   }
     192             : }
     193             : 
     194             : /* G symmetric, i < j, let E = E_{i,j}(a), G <- E~*G*E,  U <- U*E.
     195             :  * Everybody integral */
     196             : static void
     197        4205 : qf_apply_transvect_Z(GEN G, GEN U, long i, long j, GEN a)
     198             : {
     199        4205 :   long k, n = lg(G)-1;
     200        4205 :   gel(G, j) =  ZC_lincomb(gen_1, a, gel(G,j), gel(G,i));
     201       23089 :   for (k = 1; k < n; k++) gcoeff(G, j, k) = gcoeff(G, k, j);
     202        4205 :   gcoeff(G,j,j) = addmulii(gcoeff(G,j,j), a,
     203        4205 :                            addmulii(gcoeff(G,i,j), a,gcoeff(G,i,i)));
     204        4205 :   gel(U, j) =  ZC_lincomb(gen_1, a, gel(U,j), gel(U,i));
     205        4205 : }
     206             : 
     207             : /* LLL reduction of the quadratic form G (Gram matrix)
     208             :  * where we go on, even if an isotropic vector is found. */
     209             : static GEN
     210       11095 : qflllgram_indefgoon(GEN G)
     211             : {
     212             :   GEN red, U, A, U1,U2,U3,U5,U6, V, B, G2,G3,G4,G5, G6, a, g;
     213       11095 :   long i, j, n = lg(G)-1;
     214             :   int fail;
     215             : 
     216       11095 :   red = qflllgram_indef(G,1, &fail);
     217       11095 :   if (fail) return red; /*no isotropic vector found: nothing to do*/
     218             :   /* otherwise a solution is found: */
     219       11074 :   U1 = gel(red,2);
     220       11074 :   G2 = gel(red,1); /* G2[1,1] = 0 */
     221       11074 :   U2 = gel(ZV_extgcd(row(G2,1)), 2);
     222       11074 :   G3 = qf_ZM_apply(G2,U2);
     223       11074 :   U = ZM_mul(U1,U2); /* qf_apply(G,U) = G3 */
     224             :   /* G3[1,] = [0,...,0,g], g^2 | det G */
     225       11074 :   g = gcoeff(G3,1,n);
     226       11074 :   a = diviiround(negi(gcoeff(G3,n,n)), shifti(g,1));
     227       11074 :   if (signe(a)) qf_apply_transvect_Z(G3,U,1,n,a);
     228             :   /* G3[n,n] reduced mod 2g */
     229       11074 :   if (n == 2) return mkvec2(G3,U);
     230       10409 :   V = rowpermute(vecslice(G3, 2,n-1), mkvecsmall2(1,n));
     231       10409 :   A = mkmat22(gcoeff(G3,1,1),gcoeff(G3,1,n),gcoeff(G3,1,n),gcoeff(G3,2,2));
     232       10409 :   B = ground(RgM_neg(QM_mul(QM_inv(A), V)));
     233       10409 :   U3 = matid(n);
     234       51835 :   for (j = 2; j < n; j++)
     235             :   {
     236       41426 :     gcoeff(U3,1,j) = gcoeff(B,1,j-1);
     237       41426 :     gcoeff(U3,n,j) = gcoeff(B,2,j-1);
     238             :   }
     239       10409 :   G4 = qf_ZM_apply(G3,U3); /* the last column of G4 is reduced */
     240       10409 :   U = ZM_mul(U,U3);
     241       10409 :   if (n == 3) return mkvec2(G4,U);
     242             : 
     243        8141 :   red = qflllgram_indefgoon(matslice(G4,2,n-1,2,n-1));
     244        8141 :   if (typ(red) == t_MAT) return mkvec2(G4,U);
     245             :   /* Let U5:=matconcat(diagonal[1,red[2],1])
     246             :    * return [qf_ZM_apply(G5, U5), U*U5] */
     247        8141 :   G5 = gel(red,1);
     248        8141 :   U5 = gel(red,2);
     249        8141 :   G6 = cgetg(n+1,t_MAT);
     250        8141 :   gel(G6,1) = gel(G4,1);
     251        8141 :   gel(G6,n) = gel(G4,n);
     252       47299 :   for (j=2; j<n; j++)
     253             :   {
     254       39158 :     gel(G6,j) = cgetg(n+1,t_COL);
     255       39158 :     gcoeff(G6,1,j) = gcoeff(G4,j,1);
     256       39158 :     gcoeff(G6,n,j) = gcoeff(G4,j,n);
     257      258958 :     for (i=2; i<n; i++) gcoeff(G6,i,j) = gcoeff(G5,i-1,j-1);
     258             :   }
     259        8141 :   U6 = mkvec3(mkmat(gel(U,1)), ZM_mul(vecslice(U,2,n-1),U5), mkmat(gel(U,n)));
     260        8141 :   return mkvec2(G6, shallowconcat1(U6));
     261             : }
     262             : 
     263             : /* qf_ZM_apply(G,H),  where H = matrix of \tau_{i,j}, i != j */
     264             : static GEN
     265       13303 : qf_apply_tau(GEN G, long i, long j)
     266             : {
     267       13303 :   long l = lg(G), k;
     268       13303 :   G = RgM_shallowcopy(G);
     269       13303 :   swap(gel(G,i), gel(G,j));
     270       80493 :   for (k = 1; k < l; k++) swap(gcoeff(G,i,k), gcoeff(G,j,k));
     271       13303 :   return G;
     272             : }
     273             : 
     274             : /* LLL reduction of the quadratic form G (Gram matrix)
     275             :  * in dim 3 only, with detG = -1 and sign(G) = [2,1]; */
     276             : static GEN
     277        2275 : qflllgram_indefgoon2(GEN G)
     278             : {
     279             :   GEN red, G2, a, b, c, d, e, f, u, v, r, r3, U2, G3;
     280             :   int fail;
     281             : 
     282        2275 :   red = qflllgram_indef(G,1,&fail); /* always find an isotropic vector. */
     283        2275 :   G2 = qf_apply_tau(gel(red,1),1,3); /* G2[3,3] = 0 */
     284        2275 :   r = row(gel(red,2), 3);
     285        2275 :   swap(gel(r,1), gel(r,3)); /* apply tau_{1,3} */
     286        2275 :   a = gcoeff(G2,3,1);
     287        2275 :   b = gcoeff(G2,3,2);
     288        2275 :   d = bezout(a,b, &u,&v);
     289        2275 :   if (!equali1(d))
     290             :   {
     291           0 :     a = diviiexact(a,d);
     292           0 :     b = diviiexact(b,d);
     293             :   }
     294             :   /* for U2 = [-u,-b,0;-v,a,0;0,0,1]
     295             :    * G3 = qf_ZM_apply(G2,U2) has known last row (-d, 0, 0),
     296             :    * so apply to principal_minor(G3,2), instead */
     297        2275 :   U2 = mkmat22(negi(u),negi(b),negi(v),a);
     298        2275 :   G3 = qf_ZM_apply(principal_minor(G2,2),U2);
     299        2275 :   r3 = gel(r,3);
     300        2275 :   r = ZV_ZM_mul(mkvec2(gel(r,1),gel(r,2)),U2);
     301             : 
     302        2275 :   a = gcoeff(G3,1,1);
     303        2275 :   b = gcoeff(G3,1,2);
     304        2275 :   c = negi(d); /* G3[1,3] */
     305        2275 :   d = gcoeff(G3,2,2);
     306        2275 :   if (mpodd(a))
     307             :   {
     308        1204 :     e = addii(b,d);
     309        1204 :     a = addii(a, addii(b,e));
     310        1204 :     e = diviiround(negi(e),c);
     311        1204 :     f = diviiround(negi(a), shifti(c,1));
     312        1204 :     a = addmulii(addii(gel(r,1),gel(r,2)), f,r3);
     313             :   }
     314             :   else
     315             :   {
     316        1071 :     e = diviiround(negi(b),c);
     317        1071 :     f = diviiround(negi(shifti(a,-1)), c);
     318        1071 :     a = addmulii(gel(r,1), f, r3);
     319             :   }
     320        2275 :   b = addmulii(gel(r,2), e, r3);
     321        2275 :   return mkvec3(a,b, r3);
     322             : }
     323             : 
     324             : /* QUADRATIC FORM MINIMIZATION */
     325             : /* G symmetric, return ZM_Z_divexact(G,d) */
     326             : static GEN
     327       20328 : ZsymM_Z_divexact(GEN G, GEN d)
     328             : {
     329       20328 :   long i,j,l = lg(G);
     330       20328 :   GEN H = cgetg(l, t_MAT);
     331       79149 :   for (j = 1; j < l; j++)
     332             :   {
     333       58821 :     GEN c = cgetg(l, t_COL), b = gel(G,j);
     334      115486 :     for (i=1; i < j; i++) gcoeff(H,j,i) = gel(c,i) = diviiexact(gel(b,i),d);
     335       58821 :     gel(c,j) = diviiexact(gel(b,j),d);
     336       58821 :     gel(H,j) = c;
     337             :   }
     338       20328 :   return H;
     339             : }
     340             : /* by blocks, in place: G[1,1] /= d, G[2,2] *= d */
     341             : static void
     342       92032 : ZsymM_Z_divexact_partial(GEN G, long n,  GEN d)
     343             : {
     344       92032 :   long i, j, l = lg(G);
     345      184652 :   for(j = 1; j <= n; j++)
     346             :   {
     347       93208 :     for (i = 1; i < j; i++)
     348         588 :       gcoeff(G,i,j) = gcoeff(G,j,i) = diviiexact(gcoeff(G,i,j), d);
     349       92620 :     gcoeff(G,j,j) = diviiexact(gcoeff(G,j,j), d);
     350             :   }
     351      309220 :   for (; j < l; j++)
     352             :   {
     353      454113 :     for (i = n+1; i < j; i++)
     354      236925 :       gcoeff(G,i,j) = gcoeff(G,j,i) = mulii(gcoeff(G,i,j), d);
     355      217188 :     gcoeff(G,j,j) = mulii(gcoeff(G,j,j), d);
     356             :   }
     357       92032 : }
     358             : 
     359             : /* write symmetric G as [A,B;B~,C], A dxd, C (n-d)x(n-d) */
     360             : static void
     361        2289 : blocks4(GEN G, long d, long n, GEN *A, GEN *B, GEN *C)
     362             : {
     363        2289 :   GEN G2 = vecslice(G,d+1,n);
     364        2289 :   *A = principal_minor(G, d);
     365        2289 :   *B = rowslice(G2, 1, d);
     366        2289 :   *C = rowslice(G2, d+1, n);
     367        2289 : }
     368             : static GEN qfsolvemodp(GEN G, GEN p);
     369             : static void
     370       72452 : update_fm(GEN f, GEN a, long i)
     371             : {
     372       72452 :   if (!odd(i))
     373       51050 :     gel(f,1) = a;
     374             :   else
     375             :   {
     376       21402 :     long v = vals(i+1), k;
     377       21402 :     GEN b = gel(f,1), u = QM_mul(b, a);
     378       21402 :     gel(f,1) = gen_0;
     379       21402 :     if (v+2 >= lg(f)) pari_err_BUG("update_fm");
     380       26047 :     for (k = 1; k < v; k++)
     381             :     {
     382        4645 :       u = QM_mul(gel(f, k+1), u);
     383        4645 :       gel(f,k+1) = gen_0; /* for gerepileall */
     384             :     }
     385       21402 :     gel(f,v+1) = u;
     386             :   }
     387       72452 : }
     388             : static GEN
     389       41167 : prod_fm(GEN f, long i)
     390             : {
     391       41167 :   long v = vals(i), k;
     392       41167 :   GEN u = gel(f, ++v);
     393       47700 :   for (i >>= v, k = v+1; i; i >>= 1, k++)
     394        6533 :     if (odd(i)) u = QM_mul(gel(f,k), u);
     395       41167 :   return u;
     396             : }
     397             : 
     398             : /* Minimization of the quadratic form G, deg G != 0, dim n >= 2
     399             :  * G symmetric integral
     400             :  * Returns [G',U,factd] with U in GLn(Q) such that G'=U~*G*U*constant
     401             :  * is integral and has minimal determinant.
     402             :  * In dimension 3 or 4, may return a prime p if the reduction at p is
     403             :  * impossible because of local nonsolvability.
     404             :  * P,E = factor(+/- det(G)), d = det(G) "prime" -1 is ignored,
     405             :  * Either E or d should be NULL, but not both */
     406             : static GEN
     407       15980 : qfminimize_fact(GEN G, GEN P, GEN E, GEN d, long loc)
     408             : {
     409       15980 :   GEN U = NULL, Ker = NULL, faE, faP;
     410       15980 :   long n = lg(G)-1, lP = lg(P), i;
     411       15980 :   faP = vectrunc_init(lP);
     412       15980 :   faE = vecsmalltrunc_init(lP);
     413      155819 :   for (i = 1; i < lP; i++)
     414             :   {
     415             :     long Ei, vp, wp;
     416      140322 :     GEN p = gel(P,i);
     417      140322 :     if (is_pm1(p)) continue;
     418      139720 :     Ei = E ? E[i]: Z_pval(d, p); vp = Ei; wp = vp;
     419      139720 :     if (!vp) continue;
     420      139258 :     if (DEBUGLEVEL >= 3) err_printf("qfminimize: for %Ps^%ld:", p,vp);
     421      344342 :     while (vp)
     422             :     {
     423      247040 :       long idx = 0, dimKer = 0; /* -Wall */
     424      247040 :       GEN d, sol = NULL, FU = zerovec(2*expu(vp)+100);
     425      247040 :       pari_sp av = avma;
     426      319492 :       while (vp) /* loop until vp <= n */
     427             :       {
     428      290017 :         if (DEBUGLEVEL>=3 && vp <= wp)
     429           0 :         { err_printf(" %ld%%", (Ei-vp)*100/Ei); wp -= Ei/100; }
     430             :         /* The case vp = 1 can be minimized only if n is odd. */
     431      290017 :         if (vp == 1 && !odd(n)) break;
     432      282450 :         Ker = kermodp(G,p, &dimKer); /* dimKer <= vp */
     433      282450 :         if (DEBUGLEVEL >= 4) err_printf("    dimKer = %ld\n",dimKer);
     434      282450 :         if (dimKer == n) break;
     435      282450 :         G = qf_ZM_apply(G, Ker);
     436             :         /* 1st case: dimKer < vp */
     437             :         /* then the kernel mod p contains a kernel mod p^2 */
     438      282450 :         if (dimKer >= vp) break;
     439             : 
     440       72452 :         if (DEBUGLEVEL >= 4) err_printf("    case 1: dimker < vp\n");
     441       72452 :         if (dimKer == 1)
     442             :         {
     443             :           long j;
     444       70163 :           gcoeff(G,1,1) = diviiexact(gcoeff(G,1,1), sqri(p));
     445      214969 :           for (j = 2; j <= n; j++)
     446      144806 :             gcoeff(G,1,j) = gcoeff(G,j,1) = diviiexact(gcoeff(G,j,1), p);
     447       70163 :           gel(Ker,1) = RgC_Rg_div(gel(Ker,1), p);
     448       70163 :           vp -= 2;
     449             :         }
     450             :         else
     451             :         {
     452             :           GEN A, B, C, K2;
     453             :           long j, dimKer2;
     454        2289 :           blocks4(G, dimKer,n, &A,&B,&C); A = ZsymM_Z_divexact(A, p);
     455        2289 :           K2 = kermodp(A, p, &dimKer2);
     456             :           /* Write G = [pA,B;B~,C] and apply [K2/p,0;0,Id] by blocks */
     457        2289 :           A = qf_ZM_apply(A,K2); ZsymM_Z_divexact_partial(A, dimKer2, p);
     458        2289 :           B = ZM_transmul(B,K2);
     459        5131 :           for (j = 1; j <= dimKer2; j++) gel(B,j) = ZC_Z_divexact(gel(B,j), p);
     460        2289 :           G = shallowmatconcat(mkmat22(A,shallowtrans(B),B,C));
     461             :           /* Ker *= [K2,0;0,Id] */
     462        2289 :           B = ZM_mul(vecslice(Ker,1,dimKer),K2);
     463        5131 :           for (j = 1; j <= dimKer2; j++) gel(B,j) = RgC_Rg_div(gel(B,j), p);
     464        2289 :           Ker = shallowconcat(B, vecslice(Ker,dimKer+1,n));
     465        2289 :           vp -= 2*dimKer2;
     466             :         }
     467       72452 :         update_fm(FU, Ker, idx++);
     468       72452 :         if (gc_needed(av, 1))
     469             :         {
     470           0 :           if (DEBUGMEM >= 2) pari_warn(warnmem,"qfminimize");
     471           0 :           gerepileall(av, 2, &G, &FU);
     472             :         }
     473             :       }
     474      247040 :       if (idx)
     475             :       {
     476       41167 :         GEN PU = prod_fm(FU, idx);
     477       41167 :         U = U ? QM_mul(U, PU): PU;
     478             :       }
     479      259038 :       if (vp == 0) break;
     480      217565 :       if (vp == 1 && !odd(n))
     481             :       {
     482        7567 :         vectrunc_append(faP, p);
     483        7567 :         vecsmalltrunc_append(faE, 1);
     484        7567 :         break;
     485             :       }
     486      209998 :       if (dimKer == n)
     487             :       { /* trivial case: dimKer = n */
     488           0 :         if (DEBUGLEVEL >= 4) err_printf("     case 0: dimKer = n\n");
     489           0 :         G = ZsymM_Z_divexact(G, p);
     490      205084 :         vp -= n; continue;
     491             :       }
     492      209998 :       U = U ? QM_mul(U, Ker): Ker;
     493             :       /* vp = dimKer
     494             :        * 2nd case: kernel has dim >= 2 and contains an elt of norm 0 mod p^2,
     495             :        * find it */
     496      209998 :       if (dimKer > 2) {
     497       18039 :         if (DEBUGLEVEL >= 4) err_printf("    case 2.1\n");
     498       18039 :         dimKer = 3;
     499       18039 :         sol = qfsolvemodp(ZsymM_Z_divexact(principal_minor(G,3),p),  p);
     500       18039 :         sol = FpC_red(sol, p);
     501             :       }
     502      191959 :       else if (dimKer == 2)
     503             :       {
     504       97351 :         GEN a = modii(diviiexact(gcoeff(G,1,1),p), p);
     505       97351 :         GEN b = modii(diviiexact(gcoeff(G,1,2),p), p);
     506       97351 :         GEN c = modii(diviiexact(gcoeff(G,2,2),p), p);
     507       97351 :         GEN D = modii(subii(sqri(b), mulii(a,c)), p);
     508       97351 :         if (kronecker(D,p) >= 0)
     509             :         {
     510       97302 :           if (DEBUGLEVEL >= 4) err_printf("    case 2.2\n");
     511       97302 :           sol = signe(a)? mkcol2(Fp_sub(Fp_sqrt(D,p), b, p), a): vec_ei(2,1);
     512             :         }
     513             :       }
     514      209998 :       if (sol)
     515             :       {
     516             :         long j;
     517      115341 :         sol = FpC_center(sol, p, shifti(p,-1));
     518      115341 :         sol = Q_primpart(sol);
     519      115341 :         if (DEBUGLEVEL >= 4) err_printf("    sol = %Ps\n", sol);
     520      115341 :         Ker = completebasis(vecextend(sol,n), 1);
     521      115341 :         G = qf_ZM_apply(G, Ker);
     522      508829 :         for (j = 1; j < n; j++)
     523      393488 :           gcoeff(G,n,j) = gcoeff(G,j,n) = diviiexact(gcoeff(G,j,n), p);
     524      115341 :         gcoeff(G,n,n) = diviiexact(gcoeff(G,n,n), sqri(p));
     525      115341 :         U = QM_mul(U,Ker); gel(U,n) = RgC_Rg_div(gel(U,n), p);
     526      115341 :         vp -= 2; continue;
     527             :       }
     528             :       /* Now 0 < vp = dimKer < 3 and kernel contains no vector with norm p^2 */
     529             :       /* exchanging kernel and image makes minimization easier ? */
     530       94657 :       d = ZM_det(G); if (odd((n-3) / 2)) d = negi(d);
     531       94657 :       if ((vp==1 && kronecker(gmod(gdiv(negi(d), gcoeff(G,1,1)),p), p) >= 0)
     532        4949 :           || (vp==2 && odd(n) && n >= 5)
     533        4935 :           || (vp==2 && !odd(n) && kronecker(modii(diviiexact(d,sqri(p)), p),p) < 0))
     534             :       {
     535             :         long j;
     536       89743 :         if (DEBUGLEVEL >= 4) err_printf("    case 3\n");
     537       89743 :         ZsymM_Z_divexact_partial(G, dimKer, p);
     538      305069 :         for (j = dimKer+1; j <= n; j++) gel(U,j) = RgC_Rg_mul(gel(U,j), p);
     539       89743 :         vp -= 2*dimKer-n; continue;
     540             :       }
     541             : 
     542             :       /* Minimization was not possible so far. */
     543             :       /* If n == 3 or 4, this proves the local nonsolubility at p. */
     544        4914 :       if (loc && (n == 3 || n == 4))
     545             :       {
     546         483 :         if (DEBUGLEVEL >= 1) err_printf(" no local solution at %Ps\n",p);
     547         483 :         return p;
     548             :       }
     549        4431 :       vectrunc_append(faP, p);
     550        4431 :       vecsmalltrunc_append(faE, vp);
     551        4431 :       break;
     552             :     }
     553      138775 :     if (DEBUGLEVEL >= 3) err_printf("\n");
     554             :   }
     555       15497 :   if (!U) U = matid(n);
     556             :   else
     557             :   { /* apply LLL to avoid coefficient explosion */
     558       14426 :     GEN u = lllint(Q_primpart(U));
     559       14426 :     G = qf_ZM_apply(G, u);
     560       14426 :     U = QM_mul(U, u);
     561             :   }
     562       15497 :   return mkvec4(G, U, faP, faE);
     563             : }
     564             : 
     565             : /* assume G square integral */
     566             : static void
     567       19668 : check_symmetric(GEN G)
     568             : {
     569       19668 :   long i,j, l = lg(G);
     570       89200 :   for (i = 1; i < l; i++)
     571      174393 :     for(j = 1; j < i; j++)
     572      104861 :       if (!equalii(gcoeff(G,i,j), gcoeff(G,j,i)))
     573           7 :         pari_err_TYPE("qfsolve [not symmetric]",G);
     574       19661 : }
     575             : /* assume G symmetric and integral */
     576             : static void
     577          14 : symmetric_non0_coeff(GEN G, long *pi, long *pj)
     578             : {
     579          14 :   long i, j, l = lg(G);
     580          14 :   *pi = *pj = 0;
     581          14 :   for (i = 1; i < l; i++)
     582          14 :     for (j = 1; j <= i; j++)
     583          14 :       if (signe(gcoeff(G,i,j))) { *pi = i; *pj = j; return; }
     584             : }
     585             : 
     586             : GEN
     587          14 : qfminimize(GEN G)
     588             : {
     589          14 :   pari_sp av = avma;
     590             :   GEN c, d, F, H, U;
     591          14 :   long i, j, n = lg(G)-1;
     592          14 :   if (typ(G) != t_MAT) pari_err_TYPE("qfminimize", G);
     593          14 :   if (n == 0) pari_err_DOMAIN("qfminimize", "dimension" , "=", gen_0, G);
     594          14 :   if (n != nbrows(G)) pari_err_DIM("qfminimize");
     595          14 :   G = Q_primpart(G); RgM_check_ZM(G, "qfminimize");
     596          14 :   check_symmetric(G);
     597          14 :   d = ZM_det(G);
     598          14 :   if (!signe(d)) pari_err_DOMAIN("qfminimize", "det" , "=", gen_0, gen_0);
     599          14 :   F = absZ_factor(d);
     600          14 :   H = qfminimize_fact(G, gel(F,1), ZV_to_zv(gel(F,2)), NULL, 0);
     601          14 :   symmetric_non0_coeff(G, &i, &j);
     602          14 :   U = gel(H,2); H = gel(H,1);
     603          14 :   c = gdiv(gcoeff(H,i,j), RgV_dotproduct(gel(U,i), RgM_RgC_mul(G, gel(U,j))));
     604          14 :   return gerepilecopy(av, mkvec3(H, U, c));
     605             : }
     606             : 
     607             : /* CLASS GROUP COMPUTATIONS */
     608             : 
     609             : /* Compute the square root of the quadratic form q of discriminant D = 4 * md
     610             :  * Not fully implemented; only works for D squarefree except at 2, where the
     611             :  * valuation is 2 or 3. Finally, [P,E] = factor(2*abs(D)) if valuation is 3 and
     612             :  * factor(abs(D / 4)) otherwise */
     613             : static GEN
     614        2275 : qfbsqrt(GEN D, GEN md, GEN q, GEN P)
     615             : {
     616        2275 :   GEN a = gel(q,1), b = shifti(gel(q,2),-1), c = gel(q,3), B = negi(b);
     617        2275 :   GEN m, n, Q, M, N, d = negi(md); /* ac - b^2 */
     618        2275 :   long i, lP = lg(P);
     619             : 
     620             :   /* 1) solve m^2 = a, m*n = -b, n^2 = c in Z/dZ => q(n,m) = 0 mod d */
     621        2275 :   M = cgetg(lP, t_VEC);
     622        2275 :   N = cgetg(lP, t_VEC);
     623        9716 :   for (i = 1; i < lg(P); i++)
     624             :   {
     625        7441 :     GEN p = gel(P,i);
     626        7441 :     if (dvdii(a,p)) { n = Fp_sqrt(c, p); m = Fp_div(B, n, p); }
     627        5628 :     else            { m = Fp_sqrt(a, p); n = Fp_div(B, m, p); }
     628        7441 :     gel(M, i) = m;
     629        7441 :     gel(N, i) = n;
     630             :   }
     631        2275 :   m = ZV_chinese_center(M, P, NULL);
     632        2275 :   n = ZV_chinese_center(N, P, NULL);
     633             : 
     634             :   /* 2) build Q, with det=-1 such that Q(x,y,0) = G(x,y) */
     635        2275 :   N = diviiexact(addii(mulii(a,n), mulii(b,m)), d);
     636        2275 :   M = diviiexact(addii(mulii(b,n), mulii(c,m)), d);
     637        2275 :   Q = diviiexact(subiu(addii(mulii(m,M), mulii(n,N)), 1), d); /*(q(n,m)-d)/d^2 */
     638        2275 :   Q = mkmat3(mkcol3(a,b,N), mkcol3(b,c,M), mkcol3(N,M,Q)); /* det = -1 */
     639             : 
     640             :   /* 3) reduce Q to [0,0,-1; 0,1,0; -1,0,0] */
     641        2275 :   M = qflllgram_indefgoon2(Q);
     642        2275 :   if (signe(gel(M,1)) < 0) M = ZC_neg(M);
     643        2275 :   a = gel(M,1);
     644        2275 :   b = gel(M,2);
     645        2275 :   c = gel(M,3);
     646        2275 :   if (!mpodd(a)) { swap(a, c); togglesign_safe(&b); }
     647        2275 :   return mkqfb(a, shifti(b,1), shifti(c,1), D);
     648             : }
     649             : 
     650             : /* \prod gen[i]^e[i] as a Qfb, e in {0,1}^n nonzero */
     651             : static GEN
     652        3962 : qfb_factorback(GEN gen, GEN e, GEN isqrtD)
     653             : {
     654        3962 :   GEN q = NULL;
     655        3962 :   long j, l = lg(gen), n = 0;
     656       13496 :   for (j = 1; j < l; j++)
     657        9534 :     if (e[j]) { n++; q = q? qfbcompraw(q, gel(gen,j)): gel(gen,j); }
     658        3962 :   return (n <= 1)? q: qfbred0(q, 0, isqrtD, NULL);
     659             : }
     660             : 
     661             : /* unit form discriminant 4d */
     662             : static GEN
     663        1001 : id(GEN d)
     664        1001 : { retmkmat22(gen_1, gen_0, gen_0, negi(d)); }
     665             : 
     666             : /* Shanks/Bosma-Stevenhagen algorithm to compute the 2-Sylow of the class
     667             :  * group of discriminant D. Only works for D = fundamental discriminant.
     668             :  * When D = 1(4), work with 4D.
     669             :  * P2D,E2D = factor(abs(2*D))
     670             :  * Pm2D = factor(-abs(2*D))[,1].
     671             :  * Return a form having Witt invariants W at Pm2D */
     672             : static GEN
     673        2688 : quadclass2(GEN D, GEN P2D, GEN E2D, GEN Pm2D, GEN W, int n_is_4)
     674             : {
     675        2688 :   GEN U2 = NULL, gen, Wgen, isqrtD, d;
     676             :   long i, r, m;
     677        2688 :   int splice2 = mpodd(D);
     678             : 
     679        2688 :   if (!splice2) d = shifti(D,-2); else { d = D; D = shifti(D,2); }
     680             :   /* D = 4d */
     681        2688 :   if (zv_equal0(W)) return id(d);
     682        1806 :   r = lg(Pm2D) - 4; /* >= 0 since W != 0 */
     683        1806 :   m = (signe(D) > 0)? r+1: r;
     684        1806 :   if (n_is_4)
     685             :   { /* n = 4: look among forms of type q or 2*q, since Q can be imprimitive */
     686         539 :     U2 = hilberts(gen_2, d, Pm2D);
     687         539 :     if (zv_equal(U2,W)) return gmul2n(id(d),1);
     688             :   }
     689             : 
     690        1687 :   gen = cgetg(m+2, t_VEC);
     691        4571 :   for (i = 1; i <= m; i++)
     692             :   { /* no need to look at P2D[1]=2*/
     693        2884 :     GEN q = powiu(gel(P2D,i+1), E2D[i+1]);
     694        2884 :     gel(gen,i) = mkqfb(q, gen_0, negi(diviiexact(d,q)), D);
     695             :   }
     696        1687 :   if (!mpodd(d))
     697             :   {
     698        1463 :     gel(gen, ++m) = mkqfb(gen_2, gen_0, negi(shifti(d,-1)), D);
     699        1463 :     r++;
     700             :   }
     701         224 :   else if (Mod4(d) != 1)
     702             :   {
     703         119 :     gel(gen, ++m) = mkqfb(gen_2, gen_2, shifti(subsi(1,d),-1), D);
     704         119 :     r++;
     705             :   }
     706         105 :   else setlg(gen, m+1);
     707        1687 :   if (!r) return id(d);
     708             :   /* remove 2^3; leave alone 2^4 */
     709        1687 :   if (splice2) P2D = vecsplice(P2D, 1);
     710        1687 :   Wgen = cgetg(m+1, t_MAT);
     711        6153 :   for (i = 1; i <= m; i++) gel(Wgen,i) = hilberts(gmael(gen,i,1), d, Pm2D);
     712        1687 :   isqrtD = signe(D) > 0? sqrti(D) : NULL;
     713             :   for(;;)
     714        2065 :   {
     715        3752 :     GEN Wgen2, gen2, Ker, indexim = gel(Flm_indexrank(Wgen,2), 2);
     716             :     long dKer;
     717        3752 :     if (lg(indexim)-1 >= r)
     718             :     {
     719        1687 :       GEN W2 = Wgen, V;
     720        1687 :       if (lg(indexim) < lg(Wgen)) W2 = vecpermute(Wgen,indexim);
     721        1687 :       if (U2) W2 = vec_append(W2,U2);
     722        1687 :       V = Flm_Flc_invimage(W2, W, 2);
     723        1687 :       if (V)
     724             :       {
     725        1687 :         GEN Q = qfb_factorback(vecpermute(gen,indexim), V, isqrtD);
     726        1687 :         return (U2 && V[lg(V)-1])? qfbmat2(Q): qfbmat(Q);
     727             :       }
     728             :     }
     729        2065 :     Ker = Flm_ker(Wgen,2); dKer = lg(Ker)-1;
     730        2065 :     gen2 = cgetg(m+1, t_VEC);
     731        2065 :     Wgen2 = cgetg(m+1, t_MAT);
     732        4340 :     for (i = 1; i <= dKer; i++)
     733             :     {
     734        2275 :       GEN q = qfb_factorback(gen, gel(Ker,i), isqrtD);
     735        2275 :       gel(gen2,i) = q = qfbsqrt(D, d, q, P2D);
     736        2275 :       gel(Wgen2,i) = hilberts(gel(q,1), d, Pm2D);
     737             :     }
     738        4445 :     for (; i <=m; i++)
     739             :     {
     740        2380 :       long j = indexim[i-dKer];
     741        2380 :       gel(gen2,i) = gel(gen,j);
     742        2380 :       gel(Wgen2,i) = gel(Wgen,j);
     743             :     }
     744        2065 :     gen = gen2; Wgen = Wgen2;
     745             :   }
     746             : }
     747             : 
     748             : /* QUADRATIC EQUATIONS */
     749             : /* is x*y = -1 ? */
     750             : static int
     751       19443 : both_pm1(GEN x, GEN y)
     752       19443 : { return is_pm1(x) && is_pm1(y) && signe(x) == -signe(y); }
     753             : 
     754             : /* Try to solve G = 0 with small coefficients. This is proved to work if
     755             :  * -  det(G) = 1, dim <= 6 and G is LLL reduced
     756             :  * Returns G if no solution is found.
     757             :  * Exit with a norm 0 vector if one such is found.
     758             :  * If base == 1 and norm 0 is obtained, returns [H~*G*H,H] where
     759             :  * the 1st column of H is a norm 0 vector */
     760             : static GEN
     761       44370 : qftriv(GEN G, GEN R, long base)
     762             : {
     763       44370 :   long n = lg(G)-1, i;
     764             :   GEN s, H;
     765             : 
     766             :   /* case 1: A basis vector is isotropic */
     767      140949 :   for (i = 1; i <= n; i++)
     768      115464 :     if (!signe(gcoeff(G,i,i)))
     769             :     {
     770       18885 :       if (!base) return col_ei(n,i);
     771       11028 :       H = matid(n); swap(gel(H,1), gel(H,i));
     772       11028 :       return mkvec2(qf_apply_tau(G,1,i),H);
     773             :     }
     774             :   /* case 2: G has a block +- [1,0;0,-1] on the diagonal */
     775       84580 :   for (i = 2; i <= n; i++)
     776       62885 :     if (!signe(gcoeff(G,i-1,i)) && both_pm1(gcoeff(G,i-1,i-1),gcoeff(G,i,i)))
     777             :     {
     778        3790 :       s = col_ei(n,i); gel(s,i-1) = gen_m1;
     779        3790 :       if (!base) return s;
     780        2995 :       H = matid(n); gel(H,i) = gel(H,1); gel(H,1) = s;
     781        2995 :       return mkvec2(qf_ZM_apply(G,H),H);
     782             :     }
     783       21695 :   if (!R) return G; /* fail */
     784             :   /* case 3: a principal minor is 0 */
     785         542 :   s = ZM_ker(principal_minor(G, itos(R)));
     786         542 :   s = vecextend(Q_primpart(gel(s,1)), n);
     787         542 :   if (!base) return s;
     788         263 :   H = completebasis(s, 0);
     789         263 :   gel(H,n) = ZC_neg(gel(H,1)); gel(H,1) = s;
     790         263 :   return mkvec2(qf_ZM_apply(G,H),H);
     791             : }
     792             : 
     793             : /* p a prime number, G 3x3 symmetric. Finds X!=0 such that X^t G X = 0 mod p.
     794             :  * Allow returning a shorter X: to be completed with 0s. */
     795             : static GEN
     796       18039 : qfsolvemodp(GEN G, GEN p)
     797             : {
     798             :   GEN a,b,c,d,e,f, v1,v2,v3,v4,v5, x1,x2,x3,N1,N2,N3,s,r;
     799             : 
     800             :   /* principal_minor(G,3) = [a,b,d; b,c,e; d,e,f] */
     801       18039 :   a = modii(gcoeff(G,1,1), p);
     802       18039 :   if (!signe(a)) return mkcol(gen_1);
     803       15234 :   v1 = a;
     804       15234 :   b = modii(gcoeff(G,1,2), p);
     805       15234 :   c = modii(gcoeff(G,2,2), p);
     806       15234 :   v2 = modii(subii(mulii(a,c), sqri(b)), p);
     807       15234 :   if (!signe(v2)) return mkcol2(Fp_neg(b,p), a);
     808       12319 :   d = modii(gcoeff(G,1,3), p);
     809       12319 :   e = modii(gcoeff(G,2,3), p);
     810       12319 :   f = modii(gcoeff(G,3,3), p);
     811       12319 :   v4 = modii(subii(mulii(c,d), mulii(e,b)), p);
     812       12319 :   v5 = modii(subii(mulii(a,e), mulii(d,b)), p);
     813       12319 :   v3 = subii(mulii(v2,f), addii(mulii(v4,d), mulii(v5,e))); /* det(G) */
     814       12319 :   v3 = modii(v3, p);
     815       12319 :   N1 =  Fp_neg(v2,  p);
     816       12319 :   x3 = mkcol3(v4, v5, N1);
     817       12319 :   if (!signe(v3)) return x3;
     818             : 
     819             :   /* now, solve in dimension 3... reduction to the diagonal case: */
     820       10620 :   x1 = mkcol3(gen_1, gen_0, gen_0);
     821       10620 :   x2 = mkcol3(negi(b), a, gen_0);
     822       10620 :   if (kronecker(N1,p) == 1) return ZC_lincomb(Fp_sqrt(N1,p),gen_1,x1,x2);
     823        4256 :   N2 = Fp_div(Fp_neg(v3,p), v1, p);
     824        4256 :   if (kronecker(N2,p) == 1) return ZC_lincomb(Fp_sqrt(N2,p),gen_1,x2,x3);
     825        2129 :   N3 = Fp_mul(v2, N2, p);
     826        2129 :   if (kronecker(N3,p) == 1) return ZC_lincomb(Fp_sqrt(N3,p),gen_1,x1,x3);
     827        1065 :   r = gen_1;
     828             :   for(;;)
     829             :   {
     830        2046 :     s = Fp_sub(gen_1, Fp_mul(N1,Fp_sqr(r,p),p), p);
     831        2046 :     if (kronecker(s, p) <= 0) break;
     832         981 :     r = randomi(p);
     833             :   }
     834        1065 :   s = Fp_sqrt(Fp_div(s,N3,p), p);
     835        1065 :   return ZC_add(x1, ZC_lincomb(r,s,x2,x3));
     836             : }
     837             : 
     838             : /* Given a square rational matrix G of dimension n >= 1, solves over Z the
     839             :  * quadratic equation X^tGX = 0. The solution is a t_VEC (a solution) or a
     840             :  * t_MAT (totally isotropic subspace). If no solution exists, returns an
     841             :  * integer: a prime p (no local solution at p), or -1 (no real solution), or
     842             :  * -2 (n = 2 and -deg(G) not a square). */
     843             : static  GEN
     844       12515 : qfsolve_i(GEN G)
     845             : {
     846       12515 :   GEN M, signG, Min, U, G1, M1, G2, M2, solG2, P = NULL, E = NULL;
     847             :   GEN solG1, sol, Q, d, dQ, detG2;
     848             :   long n, np, codim, dim;
     849             :   int fail;
     850             : 
     851       12515 :   if (typ(G) == t_VEC && lg(G)==3)
     852             :   {
     853        8217 :     P = gel(G,2);
     854        8217 :     G = gel(G,1);
     855        8217 :     if (typ(P)==t_MAT)
     856             :     {
     857           7 :       if (!is_Z_factornon0(P)) pari_err_TYPE("qfsolve", P);
     858           7 :       P = gel(P,1);
     859             :     } else
     860        8210 :       if (!is_vec_t(typ(P)) || !RgV_is_ZVnon0(P))
     861           0 :         pari_err_TYPE("qfsolve", P);
     862             :   }
     863       12515 :   if (typ(G) != t_MAT) pari_err_TYPE("qfsolve", G);
     864       12515 :   n = lg(G)-1;
     865       12515 :   if (n == 0) pari_err_DOMAIN("qfsolve", "dimension" , "=", gen_0, G);
     866       12515 :   if (n != nbrows(G)) pari_err_DIM("qfsolve");
     867       12515 :   G = Q_primpart(G); RgM_check_ZM(G, "qfsolve");
     868       12515 :   check_symmetric(G);
     869             : 
     870             :   /* Trivial case: det = 0 */
     871       12508 :   d = ZM_det(G);
     872       12508 :   if (!signe(d))
     873             :   {
     874           7 :     if (n == 1) return mkcol(gen_1);
     875           0 :     sol = ZM_ker(G);
     876           0 :     if (lg(sol) == 2) sol = gel(sol,1);
     877           0 :     return sol;
     878             :   }
     879             : 
     880             :   /* Small dimension: n <= 2 */
     881       12501 :   if (n == 1) return gen_m1;
     882       12494 :   if (n == 2)
     883             :   {
     884          21 :     GEN t, a =  gcoeff(G,1,1);
     885          21 :     if (!signe(a)) return mkcol2(gen_1, gen_0);
     886          14 :     if (signe(d) > 0) return gen_m1; /* no real solution */
     887          14 :     if (!Z_issquareall(negi(d), &t)) return gen_m2;
     888           7 :     return mkcol2(subii(t,gcoeff(G,1,2)), a);
     889             :   }
     890             : 
     891             :   /* 1st reduction of the coefficients of G */
     892       12473 :   M = qflllgram_indef(G,0,&fail);
     893       12473 :   if (typ(M) == t_COL) return M;
     894       11962 :   G = gel(M,1);
     895       11962 :   M = gel(M,2);
     896             : 
     897             :   /* real solubility */
     898       11962 :   signG = ZV_to_zv(qfsign(G));
     899             :   {
     900       11962 :     long r =  signG[1], s = signG[2];
     901       11962 :     if (!r || !s) return gen_m1;
     902       11892 :     if (r < s) { G = ZM_neg(G); signG = mkvecsmall2(s,r);  }
     903             :   }
     904             : 
     905             :   /* factorization of the determinant */
     906       11892 :   if (!P)
     907             :   {
     908        4172 :     GEN F = absZ_factor(d);
     909        4172 :     P = gel(F,1); E = ZV_to_zv(gel(F,2));
     910             :   }
     911             : 
     912             :   /* Minimization and local solubility */
     913       11892 :   Min = qfminimize_fact(G, P, E, d, 1);
     914       11892 :   if (typ(Min) == t_INT) return Min;
     915             : 
     916       11409 :   M = QM_mul(M, gel(Min,2));
     917       11409 :   G = gel(Min,1);
     918       11409 :   P = gel(Min,3);
     919       11409 :   E = gel(Min,4);
     920             :   /* P,E = factor(|det(G))| */
     921             : 
     922             :   /* Now, we know that local solutions exist (except maybe at 2 if n==4)
     923             :    * if n==3, det(G) = +-1
     924             :    * if n==4, or n is odd, det(G) is squarefree.
     925             :    * if n>=6, det(G) has all its valuations <=2. */
     926             : 
     927             :   /* Reduction of G and search for trivial solutions. */
     928             :   /* When |det G|=1, such trivial solutions always exist. */
     929       11409 :   U = qflllgram_indef(G,0,&fail);
     930       11409 :   if (typ(U) == t_COL) return Q_primpart(RgM_RgC_mul(M,U));
     931        2989 :   G = gel(U,1);
     932        2989 :   M = QM_mul(M, gel(U,2));
     933             :   /* P,E = factor(|det(G))| */
     934             : 
     935             :   /* If n >= 6 is even, need to increment the dimension by 1 to suppress all
     936             :    * squares from det(G) */
     937        2989 :   np = lg(P)-1;
     938        2989 :   if (n < 6 || odd(n) || !np)
     939             :   {
     940        1603 :     codim = 0;
     941        1603 :     G1 = G;
     942        1603 :     M1 = NULL;
     943             :   }
     944             :   else
     945             :   {
     946             :     GEN aux;
     947             :     long i;
     948        1386 :     codim = 1; n++;
     949        1386 :     aux = gen_1;
     950        6524 :     for (i = 1; i <= np; i++)
     951        5138 :       if (E[i] == 2) { aux = mulii(aux, gel(P,i)); E[i] = 3; }
     952             :     /* aux = largest square divisor of d; choose sign(aux) so as to balance
     953             :      * the signature of G1 */
     954        1386 :     if (signG[1] > signG[2])
     955             :     {
     956         546 :       signG[2]++;
     957         546 :       aux = negi(aux);
     958             :     }
     959             :     else
     960         840 :       signG[1]++;
     961        1386 :     G1 = shallowmatconcat(diagonal_shallow(mkvec2(G,aux)));
     962             :     /* P,E = factor(|det G1|) */
     963        1386 :     Min = qfminimize_fact(G1, P, E, NULL, 1);
     964        1386 :     G1 = gel(Min,1);
     965        1386 :     M1 = gel(Min,2);
     966        1386 :     P = gel(Min,3);
     967        1386 :     E = gel(Min,4);
     968        1386 :     np = lg(P)-1;
     969             :   }
     970             : 
     971             :   /* now, d is squarefree */
     972        2989 :   if (!np) { G2 = G1; M2 = NULL; } /* |d| = 1 */
     973             :   else
     974             :   { /* |d| > 1: increment dimension by 2 */
     975        2723 :     GEN factdP, factdE, W, v = NULL;
     976             :     long i, lfactdP, lP;
     977        2723 :     codim += 2;
     978        2723 :     d = ZV_prod(P); /* d = abs(matdet(G1)); */
     979        2723 :     if (odd(signG[2])) togglesign_safe(&d); /* d = matdet(G1); */
     980        2723 :     if (n == 4)
     981             :     { /* solubility at 2, the only remaining bad prime */
     982         707 :       v = det_minors(G1);
     983         707 :       if (Mod8(d) == 1 && witt(v, gen_2) == 1) return gen_2;
     984             :     }
     985        2688 :     P = shallowconcat(mpodd(d)? mkvec2(NULL,gen_2): mkvec(NULL), P);
     986             :     /* build a binary quadratic form with given Witt invariants */
     987        2688 :     lP = lg(P); W = const_vecsmall(lP-1, 0);
     988             :     /* choose signature of Q (real invariant and sign of the discriminant) */
     989        2688 :     dQ = absi(d);
     990        2688 :     if (signG[1] > signG[2]) togglesign_safe(&dQ); /* signQ = [2,0]; */
     991        2688 :     if (n == 4 && Mod4(dQ) != 1) dQ = shifti(dQ,2);
     992        2688 :     if (n >= 5) dQ = shifti(dQ,3);
     993             : 
     994             :     /* p-adic invariants */
     995        2688 :     if (n == 4)
     996             :     {
     997         672 :       togglesign(gel(v,2));
     998         672 :       togglesign(gel(v,4)); /* v = det_minors(-G1) */
     999        2681 :       for (i = 3; i < lP; i++) W[i] = witt(v, gel(P,i)) < 0;
    1000             :     }
    1001             :     else
    1002             :     {
    1003        2016 :       long s = signe(dQ) == signe(d)? 1: -1;
    1004             :       GEN t;
    1005        2016 :       if (odd((n-3)/2)) s = -s;
    1006        2016 :       t = s > 0? utoipos(8): utoineg(8);
    1007        6055 :       for (i = 3; i < lP; i++) W[i] = hilbertii(t, gel(P,i), gel(P,i)) > 0;
    1008             :     }
    1009        2688 :     W[2] = Flv_sum(W, 2); /* for p = 2, use product formula */
    1010             : 
    1011             :     /* Construction of the 2-class group of discriminant dQ until some product
    1012             :      * of the generators gives the desired invariants. */
    1013        2688 :     factdP = vecsplice(P, 1); lfactdP =  lg(factdP);
    1014        2688 :     factdE = cgetg(lfactdP, t_VECSMALL);
    1015       11424 :     for (i = 1; i < lfactdP; i++) factdE[i] = Z_pval(dQ, gel(factdP,i));
    1016        2688 :     factdE[1]++;
    1017             :     /* factdP,factdE = factor(2|dQ|), P = factor(-2|dQ|)[,1] */
    1018        2688 :     Q = quadclass2(dQ, factdP,factdE, P, W, n == 4);
    1019             :     /* Build a form of dim=n+2 potentially unimodular */
    1020        2688 :     G2 = shallowmatconcat(diagonal_shallow(mkvec2(G1,ZM_neg(Q))));
    1021             :     /* Minimization of G2 */
    1022        2688 :     detG2 = mulii(d, ZM_det(Q));
    1023             :     /* factdP,factdE = factor(|det G2|) */
    1024        2688 :     Min = qfminimize_fact(G2, factdP, NULL, detG2, 1);
    1025        2688 :     M2 = gel(Min,2);
    1026        2688 :     G2 = gel(Min,1);
    1027             :   }
    1028             :   /* |det(G2)| = 1, find a totally isotropic subspace for G2 */
    1029        2954 :   solG2 = qflllgram_indefgoon(G2);
    1030             :   /* G2 must have a subspace of solutions of dimension > codim */
    1031        2954 :   dim = codim+1;
    1032        7273 :   while(gequal0(principal_minor(gel(solG2,1), dim))) dim ++;
    1033        2954 :   if (dim == codim+1)
    1034          14 :     pari_err_IMPL("qfsolve, dim >= 10");
    1035        2940 :   solG2 = vecslice(gel(solG2,2), 1, dim-1);
    1036             : 
    1037        2940 :   if (!M2) solG1 = solG2;
    1038             :   else
    1039             :   { /* solution of G1 is simultaneously in solG2 and x[n+1] = x[n+2] = 0*/
    1040             :     GEN K;
    1041        2681 :     solG1 = QM_mul(M2,solG2);
    1042        2681 :     K = QM_ker(rowslice(solG1,n+1,n+2));
    1043        2681 :     solG1 = QM_mul(rowslice(solG1,1,n), K);
    1044             :   }
    1045        2940 :   if (!M1) sol = solG1;
    1046             :   else
    1047             :   { /* solution of G1 is simultaneously in solG2 and x[n] = 0 */
    1048             :     GEN K;
    1049        1386 :     sol = QM_mul(M1,solG1);
    1050        1386 :     K = QM_ker(rowslice(sol,n,n));
    1051        1386 :     sol = QM_mul(rowslice(sol,1,n-1), K);
    1052             :   }
    1053        2940 :   sol = Q_primpart(QM_mul(M, sol));
    1054        2940 :   if (lg(sol) == 2) sol = gel(sol,1);
    1055        2940 :   return sol;
    1056             : }
    1057             : GEN
    1058       12515 : qfsolve(GEN G)
    1059       12515 : { pari_sp av = avma; return gerepilecopy(av, qfsolve_i(G)); }
    1060             : 
    1061             : /* G is a symmetric 3x3 matrix, and sol a solution of sol~*G*sol=0.
    1062             :  * Returns a parametrization of the solutions with the good invariants,
    1063             :  * as a 3x3 matrix, where each line contains the coefficients of each of the 3
    1064             :  * quadratic forms. If fl!=0, the fl-th form is reduced. */
    1065             : GEN
    1066        7139 : qfparam(GEN G, GEN sol, long fl)
    1067             : {
    1068        7139 :   pari_sp av = avma;
    1069             :   GEN U, G1, G2, a, b, c, d, e;
    1070        7139 :   long n, tx = typ(sol);
    1071             : 
    1072        7139 :   if (typ(G) != t_MAT) pari_err_TYPE("qfsolve", G);
    1073        7139 :   if (!is_vec_t(tx)) pari_err_TYPE("qfsolve", G);
    1074        7139 :   if (tx == t_VEC) sol = shallowtrans(sol);
    1075        7139 :   n = lg(G)-1;
    1076        7139 :   if (n == 0) pari_err_DOMAIN("qfsolve", "dimension" , "=", gen_0, G);
    1077        7139 :   if (n != nbrows(G) || n != 3 || lg(sol) != 4) pari_err_DIM("qfsolve");
    1078        7139 :   G = Q_primpart(G); RgM_check_ZM(G,"qfsolve");
    1079        7139 :   check_symmetric(G);
    1080        7139 :   sol = Q_primpart(sol); RgV_check_ZV(sol,"qfsolve");
    1081             :   /* build U such that U[,3] = sol, and |det(U)| = 1 */
    1082        7139 :   U = completebasis(sol,1);
    1083        7139 :   G1 = qf_ZM_apply(G,U); /* G1 has a 0 at the bottom right corner */
    1084        7139 :   a = shifti(gcoeff(G1,1,2),1);
    1085        7139 :   b = shifti(negi(gcoeff(G1,1,3)),1);
    1086        7139 :   c = shifti(negi(gcoeff(G1,2,3)),1);
    1087        7139 :   d = gcoeff(G1,1,1);
    1088        7139 :   e = gcoeff(G1,2,2);
    1089        7139 :   G2 = mkmat3(mkcol3(b,gen_0,d), mkcol3(c,b,a), mkcol3(gen_0,c,e));
    1090        7139 :   sol = ZM_mul(U,G2);
    1091        7139 :   if (fl)
    1092             :   {
    1093        7118 :     GEN v = row(sol,fl);
    1094             :     int fail;
    1095        7118 :     a = gel(v,1);
    1096        7118 :     b = gmul2n(gel(v,2),-1);
    1097        7118 :     c = gel(v,3);
    1098        7118 :     U = qflllgram_indef(mkmat22(a, b, b, c), 1, &fail);
    1099        7118 :     U = gel(U,2);
    1100        7118 :     a = gcoeff(U,1,1); b = gcoeff(U,1,2);
    1101        7118 :     c = gcoeff(U,2,1); d = gcoeff(U,2,2);
    1102        7118 :     U = mkmat3(mkcol3(sqri(a),mulii(a,c),sqri(c)),
    1103             :                mkcol3(shifti(mulii(a,b),1), addii(mulii(a,d),mulii(b,c)),
    1104             :                       shifti(mulii(c,d),1)),
    1105             :                mkcol3(sqri(b),mulii(b,d),sqri(d)));
    1106        7118 :     sol = ZM_mul(sol,U);
    1107             :   }
    1108        7139 :   return gerepileupto(av, sol);
    1109             : }

Generated by: LCOV version 1.14