source: git/Singular/LIB/enumpoints.lib @ 61fbaf

///////////////////////////////////////////////////////////////////////////////
version="version enumpoints.lib 4.3.0.0 Jan_2022 "; // $Id$

category="General purpose";
info="
LIBRARY: enumpoints.lib   enumerating rational points
AUTHOR: Jieao Song (8d1h)
KEYWORDS: multivariate equations; finite field

PROCEDURES:
 points(ideal I)         points in the zero set of I
 projPoints(ideal I)     projective points in the zero set of the homogeneous I
";
///////////////////////////////////////////////////////////////////////////////
// utility function that lists all elements in the ground field
static proc listFieldElements()
{
  int k;
  int ch = char(basering);
  list els = list();
  if (npars(basering)==0) { // Fp for p prime
      for (k=ch-1;k>=0;k--) { els[k+1]=k; }
  }
  else
  { // Fq defined by a minpoly
    ideal fct = factorize(var(1)^(ch^pardeg(minpoly))-var(1),1);
    for (k=size(fct);k>0;k--)
    {
      els[k]= -coeffs(fct[k],var(1))[1][1];
    }
  }
  return(els);
}
///////////////////////////////////////////////////////////////////////////////
// enumerate points coordinate by coordinate using elimination
// i: the current coordinate to enumerate
// part: partially determined coordinates
// els: elements of the field
static proc enumPoints(ideal I,int i,list part,list els)
{
  int k; ideal I1; list part1;
  list ans = list();
  int n = nvars(basering);
  poly rest = 1;
  for (k=i+1;k<=n;k++) { rest = rest*var(k); }
  poly g = eliminate(I,rest)[1];
  if (g!=0)
  { // when the eliminant is non-trivial, only enumerate its zeros
    vector cc; poly v;
    ideal fct = factorize(g,1);
    for (k=1;k<=ncols(fct);k++)
    {
      if (deg(fct[k])==1)
      {
        cc = coeffs(fct[k],var(i))[1];
        part1 = part + list(-cc[1]/cc[2]);
        if (i==n) { ans[size(ans)+1]=part1; }
        else
        {
          I1 = std(reduce(I,fct[k]));
          ans = ans + enumPoints(I1,i+1,part1,els);
        }
      }
      // this can help find the splitting field
      // else { printf("  // %p", fct[k]); }
    }
  }
  else
  { // otherwise enumerate over the entire field
    for (k=1;k<=size(els);k++)
    {
      I1 = std(subst(I,var(i),els[k]));
      if (dim(I1)>-1)
      {
        part1 = part + list(els[k]);
        if (i==n) { ans[size(ans)+1]=part1; }
        else { ans = ans + enumPoints(I1,i+1,part1,els); }
      }
    }
  }
  return(ans);
}
///////////////////////////////////////////////////////////////////////////////
// enumerate points on a homogeneous ideal
static proc enumProjPoints(ideal I,int i,list part,list els)
{
  int j; int k; ideal I1; list part1;
  list ans = list();
  int n = nvars(basering);
  for (k=i;k<=n;k++)
  { // find points of the form (1,..), (0,1,..), etc.
    I1 = I;
    part1 = part;
    for (j=i;j<k;j++)
    {
      part1 = part1 + list(0);
      I1 = subst(I1,var(j),0);
    }
    part1 = part1 + list(1);
    I1 = std(subst(I1,var(k),1));
    if (dim(I1)>-1)
    {
      if (k==n) { ans = ans + list(part1); }
      else { ans = ans + enumPoints(I1,k+1,part1,els); }
    }
  }
  return(ans);
}
///////////////////////////////////////////////////////////////////////////////
// main interface
proc points(ideal I)
"USAGE: points(I);    I homogeneous ideal
RETURN: list of coordinates of points on the zero set of I
ASSUMES: dim(I)==0
"
{
    if (char(basering)==0 and dim(I) > 0) { ERROR("over number fields the ideal must be 0-dimensional")};
    return(enumPoints(I,1,list(),listFieldElements()));
}
///////////////////////////////////////////////////////////////////////////////
proc projPoints(ideal I)
// enumerate points on a homogeneous ideal
"USAGE: projPoints(I);    I homogeneous ideal
RETURN: list of coordinates of points on the zero set of I
ASSUMES: dim(I)==0
"
{
    if (!homog(I)) { ERROR("the ideal is not homogeneous"); }
    if (char(basering)==0 and dim(I) > 1) { ERROR("over number fields the ideal must be 0-dimensional")};
    return(enumProjPoints(I,1,list(),listFieldElements()));
}
example
{
  // The set of nodes is a 0-dimensional variety over a number field.
  ring R2 = (0,q),(x,y,z,w),dp;
  minpoly = q16-q12+q8-q4+1;
  poly s = 2q13-q9-q7+q5-q3-q;
  ideal Togliatti = 64*(x-w)*(x^4 -4*x^3*w -10*x^2*y^2 -4*x^2*w^2 +16*x*w^3 -20*x*y^2*w+5*y^4 +16*w^4 -20*y^2*w^2) -5*s*(2*z -s*w)*(4*(x^2+y^2-z^2) +(1+3*(5-s^2))*w^2)^2;
  matrix Jac = jacob(Togliatti);
  ideal I2 = Togliatti+Jac;
  list L=projPoints(std(I2));
  L[1];
  size(L);
}
///////////////////////////////////////////////////////////////////////////////