Re: size of the coefficients returned by bnfisnorm()

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Sat, Oct 30, 2021 at 10:12:07PM -0400, Max Alekseyev wrote:
> Dear Bill,
> 
> I did not have a chance to thank you for your suggestion on reducing
> coefficients of bnfisnorm() based on qfparam(), but now I have a similar
> question about qfparam() itself.
> Consider an example:
> 
> ? G = matdiagonal([650, -104329, -104329]);
> ? M = qfparam(G, qfsolve(G))
> %1 =
> [-33698267 -161709950 -194002198]
> [  -521645   -2487100   -2964370]
> [ -2608225  -12519480  -15023350]
> 
> I claim that the following matrix works equally well (i.e. it could have
> been returned by qfparam), but it has much smaller entries:
> 
> ? M2 = [323, 0, 323; 5, 50, -5; 25,- 10, - 25]
> %2 =
> [323   0 323]
> [  5  50  -5]
> [ 25 -10 -25]
> 
> So, I wonder about an approach to reduce the coefficients of the matrix
> returned by qfparam(), at least in the case of diagonal matrix G.

I join a GP script (mostly by Denis Simon) that might do what you want:

? A=matdiagonal([650, -104329, -104329]);
? qfparammin(A)
%2 = [-323,0,-323;-5,-50,5;-25,10,25]

Cheers,
Bill

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\       Copyright (C) 2015 Denis Simon
\\
\\ Distributed under the terms of the GNU General Public License (GPL)
\\
\\    This code is distributed in the hope that it will be useful,
\\    but WITHOUT ANY WARRANTY; without even the implied warranty of
\\    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
\\    General Public License for more details.
\\
\\ The full text of the GPL is available at:
\\
\\                 http://www.gnu.org/licenses/
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

/*
  Author:
  Denis SIMON -> denis.simon@unicaen.fr
  address of the file:
  www.math.unicaen.fr/~simond/qfsolve.gp

  *********************************************
  *          VERSION 14/01/2015               *
  *********************************************

  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
  \\                English help                \\
  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

  This package provides functions to solve quadratic equations over Q.
  language: GP (version 2.7.2)

  It can be run under GP by the command
  gp > \r qfsolve.gp

  This package contains 4 main functions:

  - qfsolve(G,{factD}): Solve over Q the quadratic equation X~*G*X = 0.
  G must be a symmetric matrix n*n, with coefficients in Z.
  The solution might be a single vector (vectorv)
  or a matrix (whose columns generate a totally isotropic subspace).
  If no solution exists, the output is a prime number
  indicating that there is no solution in the local field Q_p
  (-1 for the reals, p for Q_p).
  factD is an optional parameter. If present, it must be equal to
  the factorization of -abs(2*matdet(G)). This saves a lot of time.

  Example:
  gp > G = [1,0,0;0,1,0;0,0,-34];
  gp > qfsolve(G)
  %1 = [-3, -5, 1]~

  - qfparam(G,sol,fl): Coefficients of quadratic forms that parametrize the
  solutions of the ternary quadratic form G, using the particular
  solution sol.
  fl is optional and can be 1, 2, or 3, in which case the 'fl'th form is
  reduced. The default is fl=3.

  Example:
  gp > qfparam(G,[-3,-5,1]~)
  %2 = 
  [ 3 -10 -3]

  [-5  -6  5]

  [ 1   0  1]
  Indeed, the solutions can be parametrized as
  [3*x^2 - 10*y*x - 3*y^2, -5*x^2 - 6*y*x + 5*y^2, x^2 + y^2]~
  
  - qflllgram_indef(G,c): Solve or reduce the quadratic form G with
  integral coefficients. G might be definite or indefinite.
  This is an lll-type algorithm with a constant 1/4<c<=1.
  c is optional and the default is c=1.

  - quadclass2(d,factd): Compute the 2-Sylow of the (narrow) class group
  of discriminant d. d must be a fondamental discriminant.
  factD is an optional parameter. If present, it must be equal to
  the factorization of abs(2*d). In this case, the 
  algorithm runs in polynomial time.

  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
  \\         Description des fonctions          \\
  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

  Programme de resolution des equations quadratiques
  langage: GP (version 2.7.2)

  pour l'utiliser, lancer gp, puis taper
  \r qfsolve.gp

  Ce fichier contient 4 principales fonctions:

  - qfsolve(G,{factD}): pour resoudre l'equation quadratique X^t*G*X = 0
  G doit etre une matrice symetrique n*n, a coefficients dans Z.
  S'il n'existe pas de solution, la reponse est un entier
  indiquant un corps local dans lequel aucune solution n'existe
  (-1 pour les reels, p pour Q_p).
  factD est optionnel. Il doit etre egal a -abs(2*matdet(G)),
  Cela entraine un gain de temps.

  Exemple:
  gp > G = [1,0,0;0,1,0;0,0,-34];
  gp > qfsolve(G)
  %1 = [-3, -5, 1]~

  - qfparam(G,sol,fl): pour parametrer les solutions de la forme
  quadratique ternaire G, en utilisant la solution particuliere sol.
  fl est optionnel et peut prendre les valeurs 1,2 ou 3.
  Il indique que la 'fl'eme forme quadratique est reduite.
  Par defaut, fl=3

  Exemple:
  gp > qfparam(G,[-3,-5,1]~)
  %2 = 
  [ 3 -10 -3]

  [-5  -6  5]

  [ 1   0  1]
  Ici, les solutions sont parametrees par
  [3*x^2 - 10*y*x - 3*y^2, -5*x^2 - 6*y*x + 5*y^2, x^2 + y^2]~

  - qflllgram_indef(G,c): pour resoudre ou reduire la forme quadratique
  G a coefficients entiers. Il s'agit d'un algorithme type LLL, avec la
  constante 1/4<c<=1.
  c est optionnel et par defaut c=1.

  - quadclass2(d,factd): determine le 2-Sylow du (narrow) groupe de classes de
  discriminant d, ou d est un discriminant fondamental.
  factd est optionnel. Il doit etre egal a -abs(2*d),
  et dans ce cas le reste de l'algorithme est polynomial.

*/

\\
\\ Usual global variables
\\ 

global(DEBUGLEVEL_qfsolve);

  DEBUGLEVEL_qfsolve = 0;

\\ use the variable DEBUGLEVEL_qfsolve :
\\ From 0 to 5 : choose a higher value to have
\\ more details printed.


\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\          SCRIPT                             \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{default_qfsolve(DEBUGLEVEL_qfsolve_val = 0) =

  DEBUGLEVEL_qfsolve = DEBUGLEVEL_qfsolve_val;
  print("  DEBUGLEVEL_qfsolve = ",DEBUGLEVEL_qfsolve);
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\          TYPE CONVERSIONS                   \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\\ build the matrix whose diagonal blocks are listed in the vector v.
{matdiagonalblock(v) = matconcat(matdiagonal(v));}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\          LINEAR ALGEBRA                     \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\\ Gives a unimodular matrix with the last column equal to v.
\\ redflag = 0 or 1. If redflag = 1, then the n-#v first columns are reduced.
{completebasis(v,redflag=0) =
my(Mv,U,re);
my(n);

  if( type(v) != "t_COL" && type(v) != "t_MAT",
    error("wrong type in completebasis"));

  if( type(v) == "t_COL", Mv = Mat(v), Mv = v);
  n = length(Mv[,1]);
  if( n == length(Mv), return(Mv));
  U = (mathnf(Mv~,1)[2]~)^-1;
  if( n==1 || !redflag, return(U));
\\ extract the n-#v columns and LLL-reduce them
  re = qflll(vecextract(U,1<<n-1,1<<(n-#Mv)-1));
\\ apply the reduction
  re = U*matdiagonalblock([re,matid(#Mv)]);

return(re);
}

\\ Compute the kernel of M mod p.
\\ returns [d,U], where
\\ d = dim (ker M mod p)
\\ U in GLn(Z), and its first d columns span the kernel.
{kermodp(M,p) =
my(U);
my(n,d);

  n = length(M);
  U = centerlift(matker(M*Mod(1,p)));
  d = length(U);
  U = completebasis(U,0);
  U = matrix(n,n,i,j,U[i,n+1-j]);

return([d,U]);
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\          INVARIANTS COMPUTATIONS            \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\\ Compute the Hilbert symbol at p
\\ where p = -1 means real place and not p = 0 as in gp
{myhilbert(a,b,p) =
  if( sign(p) < 0, return(hilbert(a,b,0)));
  return(hilbert(a,b,p));
}

\\ Given a symmetric matrix G over Z, compute the local invariant
\\ (=Witt invariant) of G at the prime p (at real place if p = -1)
\\ Assume that none of the determinant G[1..i,1..i] is 0.
{qflocalinvariant(G,p) =
my(vdet,diag);
my(n,c);

  n = length(G);
\\ Diagonalize G first.
  vdet  = vector( n+1, i, matdet(matrix(i-1,i-1,j,k,G[j,k])));
  diag = vector( n, i, vdet[i+1]*vdet[i]);

\\ Then compute the product of the Hilbert symbols
\\ (diag[i],diag[j])_p for i < j
  c = 1;
  for( i = 1, n,
    for( j = i+1, n,
          c *= myhilbert( diag[i], diag[j], p)));
return(c);
}

\\ G is a quadratic form, or a symmetrix matrix,
\\ or a list of quadratic forms with the same discriminant.
\\ If given, fa must be equal to factor(-abs(2*matdet(G)))[,1].
{qflocalinvariants(G,fa=[]) =
my(vG,sol,vdet);
my(lG);

if( DEBUGLEVEL_qfsolve >= 4, print("    starting qflocalinvariants ",G));

\\ convert G into a vector of symmetric matrices
  if( type(G) == "t_VEC", vG = G , vG = [G]);
  lG = length(vG);
  for( j = 1, lG,
    if( type(vG[j]) == "t_QFI" || type(vG[j]) == "t_QFR",
      vG[j] = Mat(vG[j])));

\\ compute the factorization of -2*abs(det(G))
  if( !length(fa), 
    fa = (factor(-abs(2*matdet(vG[1])))[,1]));

\\ in dimension 2, each invariant is a single Hilbert symbol.
  if( length(vG[1]) == 2,
    sol = matrix(length(fa),lG,i,j,
      myhilbert(vG[j][1,1],-matdet(vG[1]),fa[i]) < 0);
if( DEBUGLEVEL_qfsolve >= 4, print("    end of qflocalinvariants"));
    return([fa,sol])
  );

  sol = matrix(length(fa),lG);
  for( j = 1, lG,
    my( n = length(vG[j]));
\\ in dimension n, we need to compute a product of n Hilbert symbols.
    vdet = vector(n+1, i, matdet(matrix(i-1,i-1,k,m,vG[j][k,m])));
    for( i = 1, length(fa),
      my( p = fa[i]);
      my( h = 1);
      for( k = 1, n-1, h *= myhilbert(-vdet[k],vdet[k+1],p));
      h *= myhilbert(vdet[n],vdet[n+1],p);
      sol[i,j] = h < 0;
    )
  );
if( DEBUGLEVEL_qfsolve >= 4, print("    end of qflocalinvariants"));
return([fa,sol]);
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\        QUADRATIC FORMS REDUCTION            \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\\ M = [a,b;b;c] has integral coefficients.
\\ Gauss reduction of the binary quadratic form
\\   qf = (a,b,c)=a*X^2+2*b*X*Y+c*Y^2
\\ Returns the reduction matrix with det = +1.
{qfbreduce(M) =
my(H,di,aux);
my(a,b,c,q,r,nexta,nextb,nextc);
my(test);

if( DEBUGLEVEL_qfsolve >= 5, print("     starting qfbreduce with ",M));

  a = M[1,1]; b = M[1,2]; c = M[2,2];

  H = matid(2); test = 1;
  while( test && a,
    di = divrem(b,a); q = di[1]; r = di[2];
    if( 2*r > abs(a), r -= abs(a); q += sign(a));
    H[,2] -= q*H[,1];
    nextc = a; nextb = -r; nexta= (nextb-b)*q+c;

    test = abs(nexta) < abs(a);
    if( test,
      c = nextc; b = nextb; a = nexta;
      aux = H[,1]; H[,1] = -H[,2]; H[,2] = aux
    )
  );

if( DEBUGLEVEL_qfsolve >= 5, print("     end of qfbreduce with ",H));
return(H);
}

\\ Performs first a LLL reduction on a positive definite
\\ quadratic form QD bounding the indefinite G.
\\ Then finishes the reduction with qfsolvetriv().
{qflllgram_indef(G,c=1,base=0) =
my(M,QD,M1,S,red);
my(n);

if( DEBUGLEVEL_qfsolve >= 4, print("    qflllgram_indef with G = ",G));
  n = length(G);
  M = matid(n);
  QD = G;
  for( i = 1, n-1,
    if( !QD[i,i],
      return(qfsolvetriv(G,base))
    );
    M1 = matid(n);
    M1[i,] = -QD[i,]/QD[i,i];
    M1[i,i] = 1;
    M = (M*M1);
    QD = (M1~*QD*M1)
  );
  M = (M^(-1));
  QD = (M~*abs(QD)*M);
  S = qflllgram(QD/content(QD));
  if( #S < n, S = completebasis(S));
  red = qfsolvetriv(S~*G*S,base);
  if( type(red) == "t_COL",
    red = (S*red);
    return(red));
  red[2] = S*red[2];
  if( length(red) == 3,
    red[3] = S*red[3]);
return(red);
}

\\ LLL reduction of the quadratic form G (Gram matrix)
\\ where we go on, even if an isotropic vector is found.
{qflllgram_indefgoon(G,c=1) =
my(red,U1,G2,U2,G3,U3,G4,U,V,B,U4,G5,U5,G6);
my(n);

  red = qflllgram_indef(G,c,1);
\\ If no isotropic vector is found, nothing to do.
  if( length(red) == 2, return(red));
\\ otherwise a solution is found:
  U1 = red[2];
  G2 = red[1]; \\ On a G2[1,1] = 0
  U2 = mathnf(Mat(G2[1,]),4)[2];
  G3 = U2~*G2*U2;
\\ The first line of the matrix G3 only contains 0,
\\ except some 'g' on the right, where g^2| det G.
  n = length(G);
  U3 = matid(n); U3[1,n] = round(-G3[n,n]/G3[1,n]/2);
  G4 = U3~*G3*U3;
\\ The coeff G4[n,n] is reduced modulo 2g
  U = vecextract(G4,[1,n],[1,n]);
  if( n == 2,
    V = matrix(2,0)
  , V = vecextract(G4,[1,n],1<<(n-1)-2));
  B = round(-U^-1*V);
  U4 = matid(n);
  for( j = 2, n-1,
    U4[1,j] = B[1,j-1];
    U4[n,j] = B[2,j-1]
  );
  G5 = U4~*G4*U4;
\\ The last column of G5 is reduced
  if( n < 4, return([G5,U1*U2*U3*U4]));

  red = qflllgram_indefgoon(matrix(n-2,n-2,i,j,G5[i+1,j+1]),c);
  U5 = matdiagonalblock([Mat(1),red[2],Mat(1)]);
  G6 = U5~*G5*U5;
return([G6,U1*U2*U3*U4*U5]);
}

\\ LLL reduction of the quadratic form G (Gram matrix)
\\ in dim 3 only, with detG = -1 and sign(G) = [2,1];
{qflllgram_indefgoon2(G,c=1) =
my(red,U1,G2,bez,U2,G3,U3);
my(cc);

  red = qflllgram_indef(G,c,1);
\\ We always find an isotropic vector.
  U1 = [0,0,1;0,1,0;1,0,0];
  G2 = U1~*red[1]*U1;
\\ G2 has a 0 at the bottom right corner.
  bez = bezout(G2[3,1],G2[3,2]);
  U2 = [bez[1],G2[3,2]/bez[3],0;bez[2],-G2[3,1]/bez[3],0;0,0,-1];
  G3 = U2~*G2*U2;
\\ G3 has 0 under the co-diagonal.
  cc = (G3[1,1])%2;
  U3 = [1,0,0;  cc,1,0;
        round(-(G3[1,1]+cc*(2*G3[1,2]+G3[2,2]*cc))/2/G3[1,3]),
        round(-(G3[1,2]+cc*G3[2,2])/G3[1,3]),1];

return([U3~*G3*U3,red[2]*U1*U2*U3]);
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\        QUADRATIC FORMS MINIMIZATION         \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\\ Minimization of the quadratic form G, with nonzero determinant.
\\ of dimension n>=2.
\\ G must by symmetric and have integral coefficients.
\\ Returns [G',U,factd] with U in GLn(Q) such that G'=U~*G*U*constant
\\ is integral and has minimal determinant.
\\ In dimension 3 or 4, may return a prime p
\\ if the reduction at p is impossible because of the local non solvability.
\\ If given, factdetG must be equal to factor(abs(det(G))).
{qfminimize(G,factdetG) =
my(factd,U,Ker,Ker2,sol,aux,di);
my(p);
my(n,lf,i,vp,dimKer,dimKer2,m);

  n = length(G);
  factd = matrix(0,2);
  if( !factdetG, factdetG = factor(matdet(G)));

  lf = length(factdetG[,1]);
  i = 1; U = matid(n);

  while(i <= lf,
    vp = factdetG[i,2];
    if( vp == 0, i++; next);
    p = factdetG[i,1];
    if( p == -1, i++; next);
if( DEBUGLEVEL_qfsolve >= 4, print("    p = ",p,"^",vp));

\\ The case vp = 1 can be minimized only if n is odd.
    if( vp == 1 && n%2 == 0,
      factd = concat(factd~, Mat([p,1])~)~;
      i++; next
    );
    Ker = kermodp(G,p); dimKer = Ker[1]; Ker = Ker[2];

\\ Rem: we must have dimKer <= vp
if( DEBUGLEVEL_qfsolve >= 4, print("    dimKer = ",dimKer));
\\ trivial case: dimKer = n
    if( dimKer == n, 
if( DEBUGLEVEL_qfsolve >= 4, print("     case 0: dimKer = n"));
      G /= p;
      factdetG[i,2] -= n;
      next
    );
    G = Ker~*G*Ker;
    U = U*Ker;

\\ 1st case: dimKer < vp
\\ then the kernel mod p contains a kernel mod p^2
    if( dimKer < vp,
if( DEBUGLEVEL_qfsolve >= 4, print("    case 1: dimker < vp"));
      if( dimKer == 1,
\\        G[,1] /= p; G[1,] /= p;
        G[,1] /= p; G[1,] = G[1,]/p;
        U[,1] /= p;
        factdetG[i,2] -= 2
      , 
        Ker2 = kermodp(matrix(dimKer,dimKer,j,k,G[j,k]/p),p);
        dimKer2 = Ker2[1]; Ker2 = Ker2[2];
        for( j = 1, dimKer2, Ker2[,j] /= p);
        Ker2 = matdiagonalblock([Ker2,matid(n-dimKer)]);
        G = Ker2~*G*Ker2;
        U = U*Ker2;
        factdetG[i,2] -= 2*dimKer2
);

if( DEBUGLEVEL_qfsolve >= 4, print("    end of case 1"));
      next
    );

\\ Now, we have vp = dimKer 
\\ 2nd case: the dimension of the kernel is >=2
\\ and contains an element of norm 0 mod p^2

\\ search for an element of norm p^2... in the kernel
    if( dimKer > 2 || 
       (dimKer == 2 && issquare( di = Mod((G[1,2]^2-G[1,1]*G[2,2])/p^2,p))),
      if( dimKer > 2,
if( DEBUGLEVEL_qfsolve >= 4, print("    case 2.1"));
        dimKer = 3;
        sol = qfsolvemodp(matrix(3,3,j,k,G[j,k]/p),p)
      ,
if( DEBUGLEVEL_qfsolve >= 4, print("    case 2.2"));
        if( G[1,1]%p^2 == 0, 
          sol = [1,0]~
        , sol = [-G[1,2]/p+sqrt(di),Mod(G[1,1]/p,p)]~
        )
      );
      sol = centerlift(sol);
      sol /= content(sol);
if( DEBUGLEVEL_qfsolve >= 4, print("    sol = ",sol));
      Ker = vectorv(n, j, if( j<= dimKer, sol[j], 0)); \\ fill with 0's
      Ker = completebasis(Ker,1);
      Ker[,n] /= p;
      G = Ker~*G*Ker;
      U = U*Ker;
      factdetG[i,2] -= 2;
if( DEBUGLEVEL_qfsolve >= 4, print("    end of case 2"));
      next
    );

\\ Now, we have vp = dimKer <= 2 
\\   and the kernel contains no vector with norm p^2...

\\ In some cases, exchanging the kernel and the image
\\ makes the minimization easy.

    m = (n-1)\2-1;
    if( ( vp == 1 && issquare(Mod(-(-1)^m*matdet(G)/G[1,1],p)))
     || ( vp == 2 && n%2 == 1 && n >= 5)
     || ( vp == 2 && n%2 == 0 && !issquare(Mod((-1)^m*matdet(G)/p^2,p)))
    , 
if( DEBUGLEVEL_qfsolve >= 4, print("    case 3"));
      Ker = matid(n);
      for( j = dimKer+1, n, Ker[j,j] = p);
      G = Ker~*G*Ker/p;
      U = U*Ker;
      factdetG[i,2] -= 2*dimKer-n;
if( DEBUGLEVEL_qfsolve >= 4, print("    end of case 3"));
      next
    );

\\ Minimization was not possible se far.
\\ If n == 3 or 4, this proves the local non-solubility at p.
    if( n == 3 || n == 4, 
if( DEBUGLEVEL_qfsolve >= 1, print(" no local solution at ",p));
      return(p));

if( DEBUGLEVEL_qfsolve >= 4, print("    prime ",p," finished"));
    factd = concat(factd~,Mat([p,vp])~)~;
    i++
  );
\\ apply LLL to avoid coefficients explosion
  aux = qflll(U/content(U));
return([aux~*G*aux,U*aux,factd]);
}

/*Added by Bill Allombert */

qfparammin(M)=
{
  my(F=factor(abs(matdet(M))));
  my([H,U]=qfminimize(M,F));
  my(S=qfsolve(H));
  my(P=U*qfparam(H,S));
  P/content(P);
}