source: git/Singular/LIB/paraplanecurves.lib @ 1e1ec4

///////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Algebraic Geometry";
info="
LIBRARY:  paraplanecurves.lib Rational parametrization of rational plane curves

AUTHORS:  J. Boehm, j.boehm at mx.uni-saarland.de @*
          W. Decker, decker at mathematik.uni-kl.de> @*
          S. Laplagne, slaplagn at dm.uba.ar @*
          F. Seelisch, seelisch at mathematik.uni-kl.de

OVERVIEW:

Suppose C = {f(x,y,z)=0} is a rational plane curve, where f is homogeneous
of degree n with coefficients in Q and absolutely irreducible (these
conditions are checked automatically.) @*
After a first step, realized by a projective automorphism in the procedure
adjointIdeal, C satisfies: @*
- C does not have singularities at infinity z=0. @*
- C does not contain the point (0:1:0) (that is, the dehomogenization of f
  with respect to z is monic as a polynomial in y). @*
Considering C in the chart z<>0, the algorithm regards x as transcendental
and y as algebraic and computes an integral basis in C(x)[y] of the integral
closure of C[x] in C(x,y) using the  normalization algorithm from
@ref{normal_lib}: see @ref{integralbasis_lib}. In a future edition of the
library, also van Hoeij's algorithm for computing the integral basis will
be available. @*
From the integral basis, the adjoint ideal is obtained by linear algebra.
Alternatively, the algorithm starts with a local analysis of the singular
locus of C. Then, for each  primary component of the singular locus which
does not correspond to ordinary multiple points or cusps, the integral
basis algorithm is applied separately. The ordinary multiple points and
cusps, in turn, are addressed by a straightforward direct algorithm. The
adjoint ideal is obtained by intersecting all ideals obtained locally.
The local variant of the algorithm is used by default. @*
The linear system corresponding to the adjoint ideal maps the curve
birationally to a rational normal curve in P^(n-2). @*
Iterating the anticanonical map, the algorithm projects the rational normal
curve to PP1 for n odd resp. to a conic C2 in PP2 for n even. @*
In case n is even, the algorithm tests whether there is a rational point on
C2 and if so gives a parametrization of C2 which is defined over Q. Otherwise,
the parametrization given is defined over a quadratic field extension of Q. @*
By inverting the birational map of C to PP1 resp. to C2, a parametrization
of C is obtained (defined over Q or the quadratic field extension).

REFERENCES:

Janko Boehm: Parametrisierung rationaler Kurven, Diploma Thesis,
http://www.math.uni-sb.de/ag/schreyer/jb/diplom%20janko%20boehm.pdf

Theo de Jong: An algorithm for computing the integral closure,
Journal of Symbolic Computation 26 (3) (1998), p. 273-277

Gert-Martin Greuel, Santiago Laplagne, Frank Seelisch: Normalization of Rings,
Journal of Symbolic Computation 9 (2010), p. 887-901

Mark van Hoeij: An Algorithm for Computing an Integral Basis in an Algebraic
Function Field, Journal of Symbolic Computation 18 (1994), p. 353-363,
http://www.math.fsu.edu/~hoeij/papers/comments/jsc1994.html

KEYWORDS:
Curves; Parametrization; Rational curves; Adjoint ideal; Geometric genus


PROCEDURES:

adjointIdeal(poly, [...]);        Adjoint ideal of a plane curve
invertBirMap(ideal,ideal);        Invert a birational map of algebraic
                                  varieties
paraPlaneCurve(poly, [...]);      Compute a rational parametrization of a
                                  rational plane curve
rncAntiCanonicalMap(ideal);       Anticanonical map of a rational normal curve
rationalPointConic(poly);         Finds a point on the conic. This point has
                                  either coefficients in Q or in a quadratic
                                  extension field of Q
mapToRatNormCurve(poly,ideal);    Map a plane rational curve to a rational
                                  normal curve (RNC)
rncItProjOdd(ideal);              Map a RNC via successive anticanonical maps
                                  to PP1
rncItProjEven(ideal);             Map a RNC via successive anticanonical maps
                                  to a conic in PP2
paraConic(poly);                  Compute a rational parametrization of a conic
testParametrization(poly,ring);   Checks whether a given curve is parametrized
                                  by a given rational map (defined in the
                                  given ring)
testPointConic(poly,ring);        Checks whether a given point (defined in the
                                  given ring) lies on the given conic.
";

LIB "elim.lib";
LIB "general.lib";
LIB "primdec.lib";
LIB "absfact.lib";
LIB "matrix.lib";
LIB "random.lib";
LIB "homolog.lib";
LIB "integralbasis.lib";
LIB "normal.lib";


///////////////////////////////////////////////////////////////////////////////
proc invertBirMap(ideal phi, ideal I)
"USAGE: invertBirMap(phi, I); phi ideal, I ideal
ASSUME: The ideal phi in the basering R represents a birational map of the
        variety given by the ideal I in R to its image in projective space
        P = PP^(size(phi)-1).
NOTE:   The procedure might fail or give a wrong output if phi does
        not define a birational map.
RETURN: ring, the coordinate ring of P, with an ideal named J and an ideal
        named psi.@*
        The ideal J defines the image of phi.@*
        The ideal psi gives the inverse of phi.@*
        Note that the entries of psi should be considered as representatives
        of classes in the quotient ring R/J.@*
THEORY: We compute the ideal I(G) in R**S of the graph G of phi.@*
        The ideal J is given by the intersection of I(G) with S.@*
        The map psi is given by a relation mod J of those relations
        in I(G) which are linear in the variables of R.@*
KEYWORDS: birational map, image, inverse.
EXAMPLE: example invertBirMap; shows an example
"
{
   def Roriginal = basering;
   int n = nvars(Roriginal);
   int m = size(phi);
   /*phi: P^(n-1) --> P^(m-1)*/
   list rl = ringlist(Roriginal);
   int k;
   for(k = 1; k <= n; k++)
     {
        rl[2][k] = "x("+string(k)+")";
     }
   for(k = 1; k <= m; k++)
     {
        rl[2][k+n] = "y("+string(k)+")";
     }
   rl[3]= list(list("dp",1:(n+m)),list("C",0));
   /*Use Hilbert driven Buchberger*/
   def Rbig0 = ring(rl);
   setring Rbig0;
   ideal I = fetch(Roriginal,I);
   ideal phi = fetch(Roriginal,phi);
   ideal mi = maxideal(1);
   ideal xv = mi[1..n];
   ideal yv  = mi[n+1..n+m];
   matrix HM[2][m] = concat(transpose(yv),transpose(phi));
   ideal graph = sat(I+minor(HM,2),phi)[1];
   graph = sat(graph,xv)[1];
   intvec Hgraph = hilb(graph,1);
   setring Roriginal;
   rl[3]= list(list("dp",1:n),list("dp",1:m),list("C",0));
   def Rbig = ring(rl);
   setring Rbig;
   ideal graph = imap(Rbig0,graph);
   graph = std(graph,Hgraph);
   ideal xv = imap(Rbig0,xv);
   /*The ideal J defines the image of phi*/
   ideal J = graph;
   for(k = 1; k <= n; k++)
     {
        J = subst(J,xv[k],0);
     }
   J = compress(J);
   /*now we start inverting phi to psi*/
   matrix relpsi = diff(xv,graph);
   for(k = 1; k <= n; k++)
     {
        relpsi = subst(relpsi,xv[k],0);
     }
    relpsi = compress(relpsi);
    list rl = ringlist(Rbig);
    list rl2 = rl[2];
    rl[2] = list(rl2[n+1..n+m]);
    rl[3]= list(list("dp",1:m),list("C",0));
    def Rtarget = ring(rl);
    setring Rtarget;
    ideal J = imap(Rbig,J);
    qring RtargetmodJ = std(J);
    matrix relpsi = imap(Rbig,relpsi);
    relpsi = syz(transpose(relpsi));
    ideal psi = submat(relpsi,1..nrows(relpsi),1);
    setring Rtarget;
    ideal psi = imap(RtargetmodJ,psi);
    export(J,psi);
    int p = printlevel - voice + 3;
    dbprint(p,"// 'invertBirMap' created a ring together with two ideals J and psi.");
    dbprint(p,"// Supposing you typed, say,  def RPn = invertBirMap(phi,I);");
    dbprint(p,"// you may access the ideals by typing");
    dbprint(p,"//      setring RPn; J; psi;");
    return(Rtarget);
 }

example
{ "EXAMPLE:"; echo=2;
  ring R = 0,(x,y,z),dp;
  poly f = y^8-x^3*(z+x)^5;
  ideal adj = adjointIdeal(f);
  def Rn = invertBirMap(adj,ideal(f));
  setring(Rn);
  J;
  psi;
}

///////////////////////////////////////////////////////////////////////////////
static proc checkAssumptions(poly f)
"USAGE:  checkAssumptions(f); f poly
RETURN:  1 if assumptions are satisfied, 0 otherwise.@*
         Assumptions checked are: basering is polynomial ring in 3 variables
         with coefficients in Q, f is homogeneous and absolutely irreducible
"
{
  def Roriginal = basering;
  list rl = ringlist(Roriginal);
  rl[3] = list(list("dp",1:3),list("C",0));
  if(size(rl[1])>1){ERROR("ground field is not Q");}
  if(rl[1]!=0){ERROR("ground field is not Q");}
  if(nvars(Roriginal)!=3)
  { ERROR("not a projective plane curve: wrong number of variables"); }
  if(homog(f)==0)
  { ERROR("not a projective plane curve: polynomial is not homogeneous"); }
  def RP2 = ring(rl);
  setring RP2;
  poly f = fetch(Roriginal,f);
  if(isIrreducible(f)==0){ERROR("curve is not absolutely irreducible");}
  setring Roriginal;
}

///////////////////////////////////////////////////////////////////////////////
proc paraPlaneCurve(poly f, list #)
"USAGE:  paraPlaneCurve(f [, s]); f poly , s optional string@*
         optional string s can be: @*
         'normal': compute integral basis via normalization. @*
         'local':  make local analysis of singularities first and apply
                   normalization separately. @*
         The default is 2. @*
ASSUME:  The basering must be a polynomial ring in three variables, say x,y,z,
         with coefficients in Q. @*
         The polynomial f must be homogeneous and absolutely irreducible. @*
         The curve C = {f = 0} must be rational, i.e., have geometric genus 0
         (see @ref{genus}). @*
         These conditions will be checked automatically.
RETURN:  ring with an ideal PARA which contains a rational parametrization of
         the rational plane curve given by f; the ground field of the returned
         polynomial ring is either Q or some algebraic extension Q(a); PARA
         consists of three generators that parametrize the three coordinates
         of the rational curve
THEORY:  After a first step, realized by a projective automorphism in the
         procedure adjointIdeal, C satisfies: @*
- C does not have singularities at infinity z=0. @*
- C does not contain the point (0:1:0) (that is, the dehomogenization of f
  with respect to z is monic as a polynomial in y). @*
Considering C in the chart z<>0, the algorithm regards x as transcendental
and y as algebraic and computes an integral basis in C(x)[y] of the integral
closure of C[x] in C(x,y) using the normalization algorithm from @ref{normal_lib}:
see @ref{integralbasis_lib}. In a future edition of the library, also van Hoeij's
algorithm for computing the integral basis will be available. @*
From the integral basis, the adjoint ideal is obtained by linear algebra.
Alternatively, the algorithm starts with a local analysis of the singular
locus of C. Then, for each  primary component of the singular locus which
does not correspond to ordinary multiple points or cusps, the integral
basis algorithm is applied separately. The ordinary multiple points and
cusps, in turn, are addressed by a straightforward direct algorithm. The
adjoint ideal is obtained by intersecting all ideals obtained locally.
The local variant of the algorithm is used by default. @*
The linear system corresponding to the adjoint ideal maps the curve
birationally to a rational normal curve in P^(n-2). @*
Iterating the anticanonical map, the algorithm projects the rational normal
curve to PP1 for n odd resp. to a conic C2 in PP2 for n even. @*
In case n is even, the algorithm tests whether there is a rational point on C2
and if so gives a parametrization of C2 which is defined over Q. Otherwise the
parametrization is defined over a quadratic field extension of Q. @*
By inverting the birational map of C to PP1 resp. to C2, a parametrization of
C is obtained (defined over Q or the quadratic field extension).
KEYWORDS: rational curves, rational parametrization of rational curves.
EXAMPLE:  example paraPlaneCurve; shows an example
"
{
  int choice = 2;
  if (size(#) != 0)
  {
    if (typeof(#[1]) == "string")
    {
      string s = string(#[1]);
      if (s == "normal") { choice = 1; }
      else
      {
        if (s == "local") { choice = 2; }
        else { ERROR("expected optional argument to be either"
                   + " 'local' or 'normal'"); }
      }
    }
    else { ERROR("expected optional argument to be a string"); }
  }
  def Roriginal = basering;
  /*checking assumptions and handling the conic case*/
  checkAssumptions(f);
  list rl = ringlist(Roriginal);
  rl[2] = list("x","y","z");
  rl[3] = list(list("dp",1:3),list("C",0));
  def RP2 = ring(rl);
  setring RP2;
  poly f = fetch(Roriginal,f);
  int d = deg(f);
  if(d==2)
    {
      def RP1 = paraConic(f); // (ring, PARACONIC)
      setring RP1;
      ideal PARA = PARACONIC;
      export(PARA);
      "// 'paraPlaneCurve' created a ring together with an ideal PARA.";
      "// Supposing you typed, say,  def RP1 = paraPlaneCurve(f);";
      "// you may access the ideal by typing";
      "//      setring RP1; PARA;";
      return(RP1);
    }
  int k;
  /*the adjoint ideal*/
  ideal AI = adjointIdeal(f,list(choice,"rattestyes/firstchecksdone"));
  /*rattestyes -> causes error message if curve is not rational*/
  /*firstchecksdone -> prevents that check of assumptions will be done again*/
  /*mapping the curve to a rational normal curve V(RNC) in Proj(Rrnc) via AI*/
  def Rrnc = mapToRatNormCurve(f, AI);  // ring containing ideal RNC
  setring Rrnc;
  int m = d-1;  // the size of AI
  /*the odd dimensional case*/
  if((d mod 2) == 1)
    {
      /*mapping the rational normal curve to P^1 creating PHI*/
      ideal PHI = rncItProjOdd(RNC);
      /*composing the maps AI and PHI*/
      def Rbig = Rrnc + RP2;
      setring Rbig;
      ideal AI = imap(RP2,AI);
      ideal PROJ = imap(Rrnc,PHI);
      for(k = 1; k <= m; k++)
        {
          PROJ = subst(PROJ,var(k),AI[k]);
        }
      setring RP2;
      ideal If = ideal(fetch(Roriginal,f));
      ideal PROJ = imap(Rbig,PROJ);
      /*inverting the composed map to psi*/
      def rp1 = invertBirMap(PROJ,If); // ring containing ideal psi
      setring rp1;
      list rl1 = ringlist(rp1);
      rl1[2] = list("s","t");
      def RP1 = ring(rl1);
      setring RP1;
      ideal PARA = fetch(rp1,psi);
      export(PARA);
      "// 'paraPlaneCurve' created a ring together with an ideal PARA.";
      "// Supposing you typed, say,  def RP1 = paraPlaneCurve(f);";
      "// you may access the ideal by typing";
      "//      setring RP1; PARA;";
      return(RP1);
    }
  /*the even dimensional case*/
  /*mapping the rational normal curve to a CONIC in P^2* creating PHI*/
  def RP2conic = rncItProjEven(RNC);  // exports PHI, returns ring
                                      // containing CONIC
  setring RP2conic;
  /*mapping the conic to P^1 via pencil defined by Ipoint*/
  def RP2conicNew = projConic(CONIC);  // ring containing ideal Ipoint
                               // possibly defined over algebraic number field
                               // variables u,v,w
  /*composing the maps AI and PHI and Ipoint in two steps*/
  def Rbig = RP2conicNew + Rrnc;
  setring Rbig;
  ideal PHI = imap(Rrnc,PHI);
  ideal PROJ = imap(RP2conicNew,Ipoint);
  for(k = 1; k <= 3; k++)
    {
      PROJ = subst(PROJ,var(k),PHI[k]);
    }
  ideal AI = fetch(RP2,AI);
  for(k = 1; k <= m; k++)
    {
       PROJ = subst(PROJ,var(k+3),AI[k]);
    }
  setring RP2conicNew;
  ideal If = ideal(fetch(Roriginal,f));
  ideal PROJ = imap(Rbig,PROJ);
  /*inverting the composed map to psi*/
  def rp1 = invertBirMap(PROJ,If); // (ring, (J,psi))
  setring rp1;
  list rl1 = ringlist(rp1);
  rl1[2] = list("s","t");
  def RP1 = ring(rl1);
  setring RP1;
  ideal PARA = fetch(rp1,psi);
  export(PARA);
  "// 'paraPlaneCurve' created a ring together with an ideal PARA.";
  "// Supposing you typed, say,  def RP1 = paraPlaneCurve(f);";
  "// you may access the ideal by typing";
  "//      setring RP1; PARA;";
  return(RP1);
}

example
{ "EXAMPLE:"; echo=2;
  ring R = 0,(x,y,z),dp;
  poly f1 = 1/2*x^5+x^2*y*z^2+x^3*y*z+1/2*x*y^2*z^2-2*x*y^3*z+y^5;
  def Rp1 = paraPlaneCurve(f1);
  setring Rp1;
  PARA;
  setring R;
  poly f2 = x6+3x4y2+3x2y4+y6-4x4z2-34x3yz2-7x2y2z2+12xy3z2+6y4z2;
  f2 = f2+36x2z4+36xyz4+9y2z4;
  def Rp2 = paraPlaneCurve(f2);
  setring Rp2;
  PARA;
}

///////////////////////////////////////////////////////////////////////////////
// compute the weighted degree of p;
// this code is an exact copy of the proc in paraplanecurves.lib
// (since we do not want to make it non-static)
static proc w_deg(poly p, intvec v)
{
   if(p==0){return(-1);}
   int d=0;
   while(jet(p,d,v)==0){d++;}
   d=(transpose(leadexp(jet(p,d,v)))*v)[1];
   return(d);
}

///////////////////////////////////////////////////////////////////////////////
static proc findCoordChange(poly f, ideal JAC)
"USAGE:  findCoordChange(f, JAC); f poly, JAC ideal.
ASSUME:  The polynomial f is homogeneous in three variables, JAC is
         the Jacobi ideal of f.
RETURN:  intvec, say a,b,c. After the coordinate change
         var(3) --> a*var(1)+b*var(2)+c*var(3), the curve {f=0}
         has no singularities at infinity {var(3)=0}.
"
{
  int h = 2;
  int a,b,c;
  ideal Jfinfty;
  while(h)
    {
      c = 1;
      while(c<=h)
         {
           b = 0;
           while(b<=(h-c))
              {
                a = h-b-c;
                Jfinfty = JAC,a*var(1)+b*var(2)+c*var(3);
                if(dim(std(Jfinfty)) == 0)
                  {
                    return(a,b,c);
                  }
                b = b+1;
              }
            c = c+1;
          }
       h = h+1;
    }
}

///////////////////////////////////////////////////////////////////////////////
proc adjointIdeal(poly f, list #)
"USAGE:  adjointIdeal(f [, choices]); f polynomial in three variables, choices
         optional list consisting of one integer or of one string or of one
         integer followed by one string. @*
         Optional integer can be: @*
         1: compute integral basis via normalization. @*
         2: make local analysis of singularities first and apply normalization
            separately. @*
         3: normalization via ideal quotient. @*
         4: normalization via local ideal quotient. @*
         The default is 2. @*
         Optional string may contain substrings: @*
         - rattestyes -> causes error message if curve is not rational. @*
         - firstchecksdone -> prevents that check of assumptions will be done
           more than once.
ASSUME:  The basering must be a polynomial ring in three variables, say x,y,z,
         with coefficients in Q. @*
         The polynomial f must be homogeneous and absolutely irreducible.@*
         All these conditions will be checked automatically.@*
RETURN:  ideal, the adjoint ideal of the curve defined by f.
THEORY:  Considering C in the chart z<>0, the algorithm regards x as
transcendental and y as algebraic and computes an integral basis in C(x)[y] of
the integral closure of C[x] in C(x,y) using the normalization algorithm
from @ref{normal_lib}: see @ref{integralbasis_lib}. In a future edition of the library,
also van Hoeij's algorithm for computing the integral basis will be available.@*
From the integral basis, the adjoint ideal is obtained by linear algebra.
Alternatively, the algorithm starts with a local analysis of the singular
locus of C. Then, for each  primary component of the singular locus which
does not correspond to ordinary multiple points or cusps, the integral
basis algorithm is applied separately. The ordinary multiple points and
cusps, in turn, are addressed by a straightforward direct algorithm. The
adjoint ideal is obtained by intersecting all ideals obtained locally.
The local variant of the algorithm is used by default. @*
KEYWORDS: integral basis; normalization.
EXAMPLE: example adjointIdeal; shows an example
"
{
  list choices = #;
  if(size(#)==0)
    {
      checkAssumptions(f);
      choices = list(2, "rattestno");
    }
  if(size(#)==1)
    {
      if(typeof(choices[1])=="int")
        {
          checkAssumptions(f);
          choices = list(choices[1], "rattestno");
        }
      else
         {
           if(not(find(choices[1], "firstchecksdone")))
             {
               checkAssumptions(f);
             }
           else
             {
               choices = list(2, choices[1]);
             }
        }
    }
  if(size(#) == 2)
    {
      if(not(find(choices[2],"firstchecksdone")))
        {
          checkAssumptions(f);
        }
    }
  ideal JAC = diff(maxideal(1),ideal(f));
  ideal Jfinfty = JAC,var(3);
  /*applying a change of coordinates if (f=0) has singularities at infinity*/
  int bb1;
  if(dim(std(Jfinfty)) >= 1)
    {
       bb1 = 1;
       int a,b,c  = findCoordChange(f,JAC);
       f = subst(f,var(3),var(3)-number(a)/c*var(1)-number(b)/c*var(2));
    }
  /*applying a change of coordinates if the point (0:1:0) lies on the curve*/
  matrix co = coeffs(f,var(2));
  int bb2 = ((size(co)-1) != deg(f));
  if(bb2)
    {
       co = coeffs(f,var(1));
       int bb2x = ((size(co)-1) == deg(f));
       if(bb2x)
         {
           map perm = basering, var(2), var(1), var(3);
           f = perm(f);
         }
       else
         {
           f = subst(f,var(1),var(1)+var(2));
         }
    }
  co = coeffs(f,var(2));
  f = f/co[size(co),1];
  /*the actual computation*/
  ideal AI = adjointIdealAtWork(f,choices);
  /*reversing the changes of coordinates if needed*/
  if(bb2)
    {
       if(bb2x)
         {
           map perm = basering, var(2), var(1), var(3);
           AI = mstd(perm(AI))[2];
         }
       else
         {
           AI = mstd(substitute(AI,var(1),var(1)-var(2)))[2];
         }

    }
  if(bb1==1)
   {
      AI = mstd(substitute(AI,var(3),var(3)+number(a)/c*var(1)+number(b)/c*var(2)))[2];
    }
  return(AI);
}

example
{ "EXAMPLE:"; echo=2;
  ring R = 0,(x,y,z),dp;
  poly f = y^8-x^3*(z+x)^5;
  adjointIdeal(f);
}

///////////////////////////////////////////////////////////////////////////////
static proc adjointIdealAtWork(poly f, list choices)
"USAGE:  adjointIdealAtWork(f, choices); f polynomial in three variables,
         choices list consisting of one integer followed by one string. @*
         integer can be: @*
         1: compute integral basis via normalization. @*
         2: make local analysis of singularities first and apply normalization
            separately. @*
         3: normalization via ideal quotient. @*
         4: normalization via local ideal quotient. @*
         The default is 2. @*
         string  may contain substring: @*
         - rattestyes -> causes error message if curve is not rational. @*
ASSUME:  The basering must be a polynomial ring in three variables, say x,y,z,
         with coefficients in Q. @*
         The polynomial f must be homogeneous and absolutely irreducible. @*
         Its dehomogenization with respect to the third variable must be monic
         as a polynomial in the second variable (that is, C does not contain
         the point (0:1:0)).@*
         The curve C is not allowed to have singularities
         at infinity (z = 0). @*
RETURN:  ideal, the adjoint ideal of the curve defined by f.
"
{
  def Roriginal = basering;
  list rl = ringlist(Roriginal);
  rl[3] = list(list("dp",1:nvars(Roriginal)),list("C",0));
  def RP2 = ring(rl);
  setring RP2;
  poly f = imap(Roriginal,f);
  poly dhf = subst(f,var(3),1);
  int n = deg(f);
  if((choices[1]==1) || (choices[1]==3)) // no local analysis of singularities
    {
      ideal AI;
      if (choices[1]==1)
      { AI = adjointIdealIB(f,insert(choices,ideal(0),size(choices))); }
      else
      { AI = adjointIdealIQ(f,insert(choices,ideal(0),size(choices))); }
      AI = homog(std(AI),var(3));
      AI = sat(AI, maxideal(1))[1];
      AI = minbase(AI);
      setring Roriginal;
      return(imap(RP2,AI));
    }
  list LL = geomGenusLA(f);  // local analysis of singularities
  int sizeLL2 = size(LL[2]);
  int sizeLL3 = size(LL[3]);
  list LL3 = LL[3];
  ideal LL4 = LL[4];
  if((LL[1]!=0) && (find(choices[2],"rattestyes")))
  { ERROR("not a rational curve"); }
  if((LL[2]==0) && (sizeLL3==0)  && (LL4==1))  // smooth case
    {
      setring Roriginal;
      return(ideal(1));
    }
  int j,k;
  list rl = ringlist(RP2);
  rl[2] = list(var(1), var(2));
  rl[3] = list(list("dp",1:2),list("C",0));
  def Rdummy = ring(rl);
  ideal B;
  if(sizeLL3==0){B = 1;} // no ordinary multiple points
                               // (except possibly nodes)
  else                         // there are ordinary multiple points
                               // (other than nodes)
    {
      setring Rdummy;
      list OMP = imap(RP2,LL3);
      int ub;
      for(k=1;k<=size(OMP);k++)
        {
          if(OMP[k][1]>ub)
            {
              ub = OMP[k][1];
            }
        }
     int lb = ub;
     for(k=1;k<=size(OMP);k++)
        {
          if(OMP[k][1]<lb)
            {
              lb = OMP[k][1];
            }
        }
     for(k=lb;k<=ub;k++)
       {
         ideal A(k) = 1;
       }
     for(k=1;k<=size(OMP);k++)
       {
         A(OMP[k][1]) = intersect(A(OMP[k][1]), OMP[k][2]);
       }
     int i = ub;
     setring RP2;
     for(k=lb;k<=ub;k++)
       {
         ideal A(k) = homog(std(fetch(Rdummy,A(k))),var(3));
       }
     B = maxideal(n-i);
     ideal A;
     while(i>=lb)
        {
          A = A(i)**(i-1);
          j=1;
          while(j<=ncols(A))
            {
              if(deg(A[j]>(n-2)))
                {
                  A = sat(A, maxideal(1))[1];
                  break;
                 }
              j = j+1;
            }
          B = intersect(B,A);
          i = i-1;
       }
    }  //end else
  B = intersect(B,homog(std(LL4),var(3)));  // add nodes and cusps
  if(sizeLL2==0)  // ordinary multiple points plus cusps only
    {
      ideal AI = sat(B, maxideal(1))[1];
      AI = minbase(AI);
      setring Roriginal;
      return(imap(RP2,AI));
    }
  setring Rdummy;
  poly f = imap(RP2,dhf);
  ideal SL = jacob(f),f;
  SL = sat(SL, fetch(RP2,LL4))[1];
  if(sizeLL3!=0)
    {
      for(k=lb;k<=ub;k++)
        {
          SL = sat(SL, A(k))[1];
        }
    }
  list PD = primdecGTZ(SL);  // could be made faster -- see minAssGTZ
                             // in deltaLocMod -- only one PD needed
  int pd = size(PD);
  setring RP2;
  list PD = imap(Rdummy,PD);
  ideal AI = 1;
  for(k=1;k<=pd;k++)
    {
      if (choices[1]==2)
      { AI = intersect(AI,adjointIdealIB(f,insert(choices,PD[k][1],
                                                 size(choices)))); }
      else
      { AI = intersect(AI,adjointIdealIQ(f,insert(choices,PD[k][1],
                                                 size(choices)))); }
    }
  AI = homog(std(AI),var(3));
  AI = intersect(AI,B);
  AI = sat(AI, maxideal(1))[1];
  AI = minbase(AI);
  setring Roriginal;
  return(imap(RP2,AI));
}

///////////////////////////////////////////////////////////////////////////////
static proc adjointIdealIB(poly f, list choices)
"USAGE:  adjointIdealIB(f, choices); f polynomial in three variables, choices
         list consisting of one integer followed by one string followed by one
         ideal. @*
         integer can be: @*
         1, 2 : compute integral basis via normalization @*
         The default is 2. @*
         string  may contain substring: @*
         - rattestyes -> causes error message if curve is not rational. @*
         ideal serves as input for @ref integralBasis.
ASSUME:  The basering must be a polynomial ring in three variables, say x,y,z,
         with coefficients in Q. @*
         The polynomial f must be homogeneous and absolutely irreducible.@*
         Its dehomogenization with respect to the third variable must be monic
         as a polynomial in the second variable (that is, C does not contain
         the point (0:1:0)).@*
         The curve C is not allowed to have singularities
         at infinity (z = 0). @*
RETURN:  ideal containing the adjoint ideal of the curve defined by f. @*
"
{
  poly dhf = subst(f,var(3),1);
  def Roriginal = basering;
  list rl = ringlist(Roriginal);
  rl[2] = list(var(1), var(2));
  rl[3] = list(list("dp",1:2),list("C",0));
  def Rdummy = ring(rl);
  setring Rdummy;
  poly f = imap(Roriginal,dhf);
  poly d2f = diff(f,var(2));
  list DATA = imap(Roriginal,choices);
  /* Creating rings for later use */
  list rl = ringlist(Rdummy);
  rl[2] = list(var(2), var(1));
  rl[3] = list(list("lp",1:2),list("C",0));
  def Rred = ring(rl);   // make var(2) > var(1)
  rl = ringlist(Rdummy);
  rl[1] = list(0,list(var(1)),list(list("dp",1)),ideal(0));
  rl[2] = list(var(2));
  rl[3] = list(list("dp",1),list("C",0));
  def QF = ring(rl);   // make var(1) transcendental
  list LIntB;
  if(DATA[1] <= 4)  // use normalization algorithm
    {
      LIntB = integralBasis(f, 2, list(list("inputC", DATA[3]),"isIrred"));
    }
  else                                 // use van Hoeij's algorithm
    {
      LIntB = integralBasisVH(f,DATA[3],2);  // van Hoeij in future version
                                             // used when DATA[1] = 5
    }
  if(find(DATA[2],"rattestyes") && (DATA[3]==0))
    {
      setring Roriginal;
      int gg = geomGenusIB(f,imap(Rdummy, LIntB));
      if(gg!=0){ERROR("not a rational curve");}
      setring Rdummy;
    }
  int i,j,k,l;
  ideal IB = LIntB[1];
  poly d = LIntB[2];
  int sL=size(IB);
  setring Rred;
  ideal IB = imap(Rdummy,IB);
  ideal fred = std(imap(Rdummy,f));
  IB = reduce(IB,fred);
  matrix M = coeffs(IB,var(1));
  setring QF;
  matrix M = imap(Rred,M);
  poly d = imap(Rdummy,d);
  M=1/d*M;
  list LUM = ludecomp(M);
  list LS;
  matrix dummyvector[sL][1];
  matrix Gij[sL][sL];
  matrix Tr[sL][sL];
  setring Rred;
  poly eiej;
  list Iij, empty;
  matrix Gij[sL][sL];
  for(i = 1; i <= sL; i++)
     {
       for(j = i; j <= sL; j++)
          {
            setring Rred;
            Gij = 0;
            eiej = IB[i]*IB[j];
            Iij=empty;
            for(k = 1; k <= sL; k++)
               {
                  Iij[k] = reduce(eiej*IB[k],fred);
               }
            Gij = coeffs(ideal(Iij[1..sL]),var(1));
            setring QF;
            Gij = imap (Rred, Gij);
            for(k = 1; k <= sL; k++)
               {
                  dummyvector = Gij[1..sL,k];
                  LS = lusolve(LUM[1], LUM[2], LUM[3], dummyvector);
                  Tr[i,j] = Tr[i,j] + 1/d^3*LS[2][k,1];
               }
          }
     }
  for(i = 1; i <= sL; i++)
     {
       for(j = 1; j < i; j++)
          {
             Tr[i,j] = Tr[j,i];
          }
     }
  LUM = ludecomp(Tr);
  setring Rred;
  poly d2f = imap(Rdummy,d2f);
  IB = d2f*IB;
  IB = reduce(IB,fred);
  setring QF;
  matrix IB = transpose(matrix(imap(Rred,IB)));
  IB = 1/d*IB;
  LS = lusolve(LUM[1], LUM[2], LUM[3], IB);
  ideal LL = ideal(LS[2]);
  setring Roriginal;
  ideal AI = imap(QF,LL);
  return(AI);
}

///////////////////////////////////////////////////////////////////////////////
static proc adjointIdealIQ(poly f, list choices)
{
  poly dhf = subst(f,var(3),1);
  def Roriginal = basering;
  list rl = ringlist(Roriginal);
  rl[2] = list(var(1), var(2));
  rl[3] = list(list("dp",1:2),list("C",0));
  def Rdummy = ring(rl);
  setring Rdummy;
  list DATA = imap(Roriginal,choices);
  poly f = imap(Roriginal,dhf);
  list LIntB;
  if(DATA[1] <= 4)  // use normalization algorithm
    {
      LIntB = integralBasis(f, 2, list(list("inputC", DATA[3]),"isIrred"));
    }
  else                                 // use van Hoeij's algorithm
    {
      LIntB = integralBasisVH(f,DATA[3],2);  // van Hoeij in future version
                                             // used when DATA[1] = 5
    }
  if(find(DATA[2],"rattestyes") && (DATA[3]==0))
    {
      setring Roriginal;
      int gg = geomGenusIB(f,imap(Rdummy, LIntB));
      if(gg!=0){ERROR("not a rational curve");}
      setring Rdummy;
    }
  int i,j,k,l;
  ideal IB = LIntB[1];
  poly d = LIntB[2];
  ideal fd = f, d;
  ideal IBf = IB, f;
  ideal AI = quotient(fd, IBf);
//"#### IB:"; IB;
//"#### d:", d;
  setring Roriginal;
  ideal AI = imap(Rdummy,AI);
//"#### AI:"; AI;
  return(AI);
}

///////////////////////////////////////////////////////////////////////////////
proc mapToRatNormCurve(poly f, ideal AI)
"USAGE:  mapToRatNormCurve(f, AI); f polynomial, AI ideal
ASSUME:  The polynomial f is homogeneous in three variables and absolutely
         irreducible.
         The plane curve C defined by f is rational.
         The ideal AI is the adjoint ideal of C.
RETURN:  ring with an ideal RNC.
EXAMPLE: example mapToRatNormCurve; shows an example
"
{
   int n = size(AI);
   int k;
   //if(n!=deg(f)-1){ERROR("not a rational curve");}
   def Roriginal = basering;
   ideal IC = f;
   list rl = ringlist(Roriginal);
   /* begin workaround elimination*/
   for(k = 1; k <= 3; k++)
     {
        rl[2][k] = "x("+string(k)+")";
     }
   for(k = 1; k <= n; k++)
     {
        rl[2][k+3] = "y("+string(k)+")";
     }
   rl[3]= list(list("dp",1:(3+n)),list("C",0));
   def Relim = ring(rl);
   setring Relim;
   ideal IC = fetch(Roriginal,IC);
   ideal AI = fetch(Roriginal,AI);
   ideal J;
   J = IC;
   for(k=1;k<=n;k++)
     {
       J=J,var(k+3)-AI[k];
     }
   ideal SJ = std(J);
   intvec HJ = hilb(SJ,1);
   ideal RNC = eliminate(J,x(1)*x(2)*x(3),HJ);
   list rl = ringlist(Relim);
   list rl2 = rl[2];
   rl[2] = list(rl2[4..n+3]);
   rl[3]= list(list("dp",1:n),list("C",0));
   def Rtarget = ring(rl);
   setring Rtarget;
   ideal RNC = imap(Relim,RNC);
   /* end workaround elimination*/
   export(RNC);
   int p = printlevel - voice + 3;
   dbprint(p,"//'mapToRatNorm' created a ring together with an ideal RNC.");
   dbprint(p,"// Supposing you typed, say,  def RPn = mapToRatNorm(f,AI);");
   dbprint(p,"// you may access the ideal by typing");
   dbprint(p,"//      setring RPn; RNC;");
   return(Rtarget);
}

example
{ "EXAMPLE:"; echo=2;
  ring R = 0,(x,y,z),dp;
  poly f = y^8-x^3*(z+x)^5;
  ideal adj = adjointIdeal(f);
  def Rn = mapToRatNormCurve(f,adj);
  setring(Rn);
  RNC;
}

///////////////////////////////////////////////////////////////////////////////
proc rncAntiCanonicalMap(ideal I)
"USAGE:  rncAntiCanonicalMap(I); I ideal
ASSUME:  I is a homogeneous ideal in the basering
         defining a rational normal curve C in PP^n.
NOTE:   The procedure will fail or give a wrong output if I is not the
        ideal of a rational normal curve.
RETURN:  ideal defining the anticanonical map  C --> PP^(n-2). @*
         Note that the entries of the ideal should be considered as
         representatives of elements in R/I, where R is the basering.
THEORY:  The anti-canonical map of a rational normal curve
         maps C isomorpically to a rational normal curve in PP^(n-2).
KEYWORDS: rational normal curve, projection.
EXAMPLE: example rncAntiCanonicalMap; shows an example
"
{
  def Roriginal = basering;
  list rl = ringlist(Roriginal);
  rl[3] = list(list("dp",1:nvars(Roriginal)),list("C",0));
  def RoriginalDP = ring(rl);
  setring RoriginalDP;
  ideal I = imap(Roriginal,I);
  int cc = nvars(RoriginalDP)-2;
  module AKD = Ext_R(cc,I);
  qring qI = std(I);
  matrix AKD = imap(RoriginalDP,AKD);
  AKD = syz(transpose(AKD));
  ideal PR = submat(AKD,1..nrows(AKD),1);
  setring Roriginal;
  return(imap(qI,PR));
}

example
{ "EXAMPLE:"; echo=2;
  ring R = 0,(x,y,z),dp;
  poly f = y^8-x^3*(z+x)^5;
  ideal adj = adjointIdeal(f);
  def Rn = mapToRatNormCurve(f,adj);
  setring(Rn);
  RNC;
  rncAntiCanonicalMap(RNC);
}


///////////////////////////////////////////////////////////////////////////////
proc rncItProjOdd(ideal I)
"USAGE:  rncItProjOdd(I); I ideal
ASSUME:  I is a homogeneous ideal in the basering with n+1 variables
         defining a rational normal curve C in PP^n with n odd.
NOTE:    The procedure will fail or give a wrong output if I is not the
         ideal of a rational normal curve. It will test whether n is odd.
RETURN:  ideal PHI defining an isomorphic projection of C to PP^1.@*
         Note that the entries of PHI should be considered as
         representatives of elements in R/I, where R is the basering.
THEORY:  We iterate the procedure @ref{rncAntiCanonicalMap} to obtain PHI.
KEYWORDS: rational normal curve, projection.
SEE ALSO: rncItProjEven.
EXAMPLE: example rncItProjOdd; shows an example
"
{
  int n = nvars(basering);
  if((n mod 2) == 1){ERROR("Pn has even dimension");}
  def Roriginal = basering;
  list rlo = ringlist(Roriginal);
  rlo[3]= list(list("dp",1:n),list("C",0));
  int k;
  for(k = 1; k <= n; k++)
    {
      rlo[2][k] = "z("+string(k)+")";
    }
  def RoriginalCopy = ring(rlo);
  for(k = 1; k <= n; k++)
    {
      rlo[2][k] = "y("+string(k)+")";
    }
  def Rold = ring(rlo);
  setring RoriginalCopy;
  ideal PHI  = maxideal(1);
  setring Rold;
  ideal J = fetch(Roriginal,I);
  list rl2;
  def Rnew;
  def Rbig;
  def Relim;
  intvec HJJ;
  while(n>2)
     {
        ideal PR = rncAntiCanonicalMap(J);
        list rl = ringlist(Rold);
        Rbig = Rold + RoriginalCopy;
        setring Rbig;
        ideal PHI = imap(RoriginalCopy,PHI);
        ideal dummy = imap(Rold,PR);
        for(k = 1; k <= n; k++)
          {
             dummy = subst(dummy,var(k),PHI[k]);
          }
        setring RoriginalCopy;
        PHI = imap(Rbig,dummy);
        /* begin workaround elimination*/
        setring Rold;
        for(k = 1; k <= n; k++)
          {
            rl[2][k] = "x("+string(k)+")";
          }
        for(k = 1; k <= n-2; k++)
          {
            rl[2][k+n] = "y("+string(k)+")";
          }
        rl[3]= list(list("dp",1:(2*n-2)),list("C",0));
        Relim = ring(rl);
        setring Relim;
        ideal J = fetch(Rold,J);
        ideal PR = fetch(Rold,PR);
        ideal JJ = J;
        poly pvar=1;
        for(k = 1; k <= n; k++)
          {
            pvar = pvar*var(k);
          }
        for(k=1;k<=n-2;k++)
          {
            JJ=JJ,var(k+n)-PR[k];
          }
        ideal SJJ = std(JJ);
        HJJ = hilb(SJJ,1);
        J = eliminate(JJ,pvar,HJJ);
        list rl = ringlist(Relim);
        rl2 = rl[2];
        rl[2] = list(rl2[n+1..2*n-2]);
        rl[3]= list(list("dp",1:(n-2)),list("C",0));
        Rnew = ring(rl);
        setring Rnew;
        ideal J = imap(Relim,J);
        /* end workaround elimination*/
        Rold = Rnew;
        setring Rold;
        n = n-2;
     }
  setring Roriginal;
  return(fetch(RoriginalCopy,PHI));
}

example
{ "EXAMPLE:"; echo=2;
  ring R = 0,(x,y,z),dp;
  poly f = -x7-10x5y2-10x4y3-3x3y4+8x2y5+7xy6+11y7+3x6+10x5y +30x4y2
           +26x3y3-13x2y4-29xy5-33y6-3x5-20x4y-33x3y2-8x2y3+37xy4+33y5
           +x4+10x3y+13x2y2-15xy3-11y4;
  f = homog(f,z);
  ideal adj = adjointIdeal(f);
  def Rn = mapToRatNormCurve(f,adj);
  setring(Rn);
  RNC;
  rncItProjOdd(RNC);
}


///////////////////////////////////////////////////////////////////////////////
proc rncItProjEven(ideal I)
"USAGE:  rncItProjEven(I); I ideal
ASSUME:  I is a homogeneous ideal in the basering with n+1 variables
         defining a rational normal curve C in PP^n with n even.
NOTE:    The procedure will fail or give a wrong output if I is not the
         ideal of a rational normal curve. It will test whether n is odd.
RETURN:  ring with an ideal CONIC defining a conic C2 in PP^2.@*
         In addition, an ideal PHI in the basering defining an isomorphic
         projection of C to C2 will be exported.@*
         Note that the entries of PHI should be considered as
         representatives of elements in R/I, where R is the basering.
THEORY:  We iterate the procedure @ref{rncAntiCanonicalMap} to obtain PHI.
KEYWORDS: rational normal curve, projection.
SEE ALSO: rncItProjOdd.
EXAMPLE: example rncItProjEven; shows an example
"
{
  int n = nvars(basering);
  if((n mod 2) == 0){ERROR("Pn has odd dimension");}
  def Roriginal = basering;
  list rlo = ringlist(Roriginal);
  rlo[3]= list(list("dp",1:n),list("C",0));
  int k;
  for(k = 1; k <= n; k++)
    {
      rlo[2][k] = "z("+string(k)+")";
    }
  def RoriginalCopy = ring(rlo);
  for(k = 1; k <= n; k++)
    {
      rlo[2][k] = "y("+string(k)+")";
    }
  def Rold = ring(rlo);
  setring RoriginalCopy;
  ideal PHI  = maxideal(1);
  setring Rold;
  ideal J = fetch(Roriginal,I);
  list rl2;
  def Rnew;
  def Rbig;
  def Relim;
  intvec HJJ;
  while(n>3)
     {
        ideal PR = rncAntiCanonicalMap(J);
        list rl = ringlist(Rold);
        Rbig = Rold + RoriginalCopy;
        setring Rbig;
        ideal PHI = imap(RoriginalCopy,PHI);
        ideal dummy = imap(Rold,PR);
        for(k = 1; k <= n; k++)
          {
             dummy = subst(dummy,var(k),PHI[k]);
          }
        setring RoriginalCopy;
        PHI = imap(Rbig,dummy);
        /* begin workaround elimination*/
        setring Rold;
        for(k = 1; k <= n; k++)
          {
            rl[2][k] = "x("+string(k)+")";
          }
        for(k = 1; k <= n-2; k++)
          {
            rl[2][k+n] = "y("+string(k)+")";
          }
        rl[3]= list(list("dp",1:(2*n-2)),list("C",0));
        Relim = ring(rl);
        setring Relim;
        ideal J = fetch(Rold,J);
        ideal PR = fetch(Rold,PR);
        ideal JJ = J;
        poly pvar=1;
        for(k = 1; k <= n; k++)
          {
            pvar = pvar*var(k);
          }
        for(k=1;k<=n-2;k++)
          {
            JJ=JJ,var(k+n)-PR[k];
          }
        ideal SJJ = std(JJ);
        HJJ = hilb(SJJ,1);
        J = eliminate(JJ,pvar,HJJ);
        list rl = ringlist(Relim);
        rl2 = rl[2];
        rl[2] = list(rl2[n+1..2*n-2]);
        rl[3]= list(list("dp",1:(n-2)),list("C",0));
        Rnew = ring(rl);
        setring Rnew;
        ideal J = imap(Relim,J);
        /* end workaround elimination*/
        Rold = Rnew;
        setring Rold;
        n = n-2;
     }
  poly CONIC = J[1];
  export(CONIC);
  setring Roriginal;
  ideal PHI = fetch(RoriginalCopy,PHI);
  export(PHI);
  return(Rold);
}

example
{ "EXAMPLE:"; echo=2;
  ring R = 0,(x,y,z),dp;
  poly f = y^8-x^3*(z+x)^5;
  ideal adj = adjointIdeal(f);
  def Rn = mapToRatNormCurve(f,adj);
  setring(Rn);
  RNC;
  def Rc = rncItProjEven(RNC);
  PHI;
  setring Rc;
  CONIC;
}

///////////////////////////////////////////////////////////////////////////////
static proc geomGenusIB(poly f, list #)
"USAGE:  geomGenusIB(f [, L]); f poly, L optional list representing the
         integral basis as returned by @ref integralBasisJ.
ASSUME:  The basering must be a polynomial ring in three variables, say x,y,z,
         with coefficients in Q. @*
         The polynomial f must be homogeneous and absolutely irreducible.@*
         Its dehomogenization with respect to the third variable must be monic
         as a polynomial in the second variable (that is, the curve C = {f = 0}
         does not contain the point (0:1:0)).@*
         The curve C is not allowed to have singularities
         at infinity (z = 0). @*
NOTE:    The last two conditions can be met by a suitable change of coordinates in PGL(3)
         as applied in the procedure @ref adjointIdeal. The other conditions
         can be tested using @ref checkAssumptions.@*
RETURN:  int, the geometric genus of C.
THEORY:  We compute an integral basis of the integral closure of the coordinate
         ring of C and from that the geometric genus.@*
KEYWORDS: geometric genus, plane curves.
SEE ALSO: genus.
"
{
  int bb = size(#);
  poly dhf = subst(f,var(3),1);
  def Roriginal = basering;
  list rl = ringlist(Roriginal);
  rl[2] = list(var(1), var(2));
  rl[3] = list(list("dp",1:2),list("C",0));
  def Rdummy = ring(rl);
  setring Rdummy;
  poly f = imap(Roriginal,dhf);
  list LIntB;
  if(bb == 0)
    {
      LIntB = integralBasis(f,2,"isIrred");
    }
  else
    {
      LIntB = imap(Roriginal,#);
    }
  ideal IB = LIntB[1];
  poly d = LIntB[2];
  int ud = deg(d);
  int sL = size(IB);
  int k;
  int gg = (sL-1)*(sL-2) div 2-sL*ud;
  for(k = 1; k <= sL; k++)
     {
        gg = gg + deg(gcd(d,IB[k]));
     }
  setring Roriginal;
  return(gg);
}


///////////////////////////////////////////////////////////////////////////////
static proc geomGenusLA(poly F)
"USAGE:  geomGenusLA(F); F polynomial
ASSUME:  The basering must be a polynomial ring in three variables. @*
         The polynomial F must be homogeneous.@*
RETURN:  list L:
@texinfo
@table @asis
@item @code{L[1]}; int:
         the geometric genus p_g = p_a - delta of the projective
         curve C defined by F, where p_a is the arithmetic genus.
@item @code{L[2]}; int:
         is positive if C has singularities other
         than ordinary multiple points.@*
@item @code{L[3]}; list:
         consists of one list for each primary component
         of the singular locus of C which correponds to a set of conjugated
         ordinary multiple points. Each list consists of an int, the
         multiplicity of the points, and an ideal, the primary component.
@end table
@end texinfo
NOTE:    delta is the sum of all local delta-invariants of the singularities,
         i.e. dim(R'/R), R' the normalization of the local ring R of the
         singularity. @*
SEE ALSO: genus
"
{
   int w = printlevel-voice+2;  // w=printlevel (default: w=0)
   int d = deg(F);
   def R = basering;
   execute("ring S=("+charstr(R)+"),(x,y,t),dp;");
   execute("ring C=("+charstr(R)+"),(x,y),ds;");
   int genus=(d-1)*(d-2) div 2;
   if(w>=1){"the arithmetic genus of the plane curve:";genus;pause();}

   int delt,deltaloc,deltainf,tau,tauinf,cusps,iloc,iglob,l,nsing,
       tauloc,tausing,k,rat,nbranchinf,nbranch,nodes,cuspsinf,nodesinf;
   list inv;
   execute("ring newR=("+charstr(R)+"),(x,y),dp;");
   //the singularities at the affine part
   map sigma=R,var(1),var(2),1;
   ideal I=sigma(F);

   list OMPButNodes;
   int sizeOMPButNodes;
   int NotOnlyOMPPlusCusps;

   ideal I1=jacob(I);
   matrix Hess[2][2]=jacob(I1);
   ideal ID=I+I1+ideal(det(Hess));//singular locus of I+I1
   ideal radID=std(radical(ID));//the non-nodal locus
   if(w>=1){"the non-nodal locus:";"";radID;pause();"";}
   if(deg(radID[1])==0)
   {
     ideal IDsing=1;
   }
   else
   {
     ideal IDsing=minor(jacob(ID),2)+radID;//singular locus of ID
   }

   iglob=vdim(std(IDsing));

   ideal radIDsing = 1;

   if(iglob!=0)//computation of the radical of IDsing
   {
      radIDsing=reduce(IDsing,radID);
      if(size(radIDsing)==0)
      {
         radIDsing=radID;
         attrib(radIDsing,"isSB",1);
      }
      else
      {
         radIDsing=std(radical(IDsing));
      }
      iglob=vdim(radIDsing);
      if((w>=1)&&(iglob))
          {"the non-nodal-cuspidal locus:";radIDsing;pause();"";}
   }
   cusps=vdim(radID)-iglob;

   ideal NodesPlusCusps  = radical(sat(I+I1, radIDsing)[1]);

   nsing=nsing+cusps;

   if(iglob==0)
   {
      if(w>=1){"             there are only cusps and nodes";"";}
      tau=vdim(std(I+jacob(I)));
      tauinf=tauinf+tau;
      nodes=tau-2*cusps;
      delt=nodes+cusps;
      nbranch=2*tau-3*cusps;
      nsing=nsing+nodes;
   }
   else
   {
       if(w>=1){"the non-nodal-cuspidal singularities";"";}
       setring C;
       ideal I1=imap(newR,radIDsing);
       iloc=vdim(std(I1));
       if(iglob==iloc)
       {
          if(w>=1){"only cusps and nodes outside (0,0,1)";}
          setring newR;
          tau=vdim(std(I+jacob(I)));
          tauinf=tauinf+tau;
          inv=deltaLocMod(I[1],maxideal(1));
          delt=inv[1];
          tauloc=inv[2];
          nodes=tau-tauloc-2*cusps;
          nsing=nsing+nodes;
          if(inv[2]!=0)
            {
              nsing=nsing+1;
            }
          nbranch=inv[3]+ 2*nodes+cusps;
          delt=delt+nodes+cusps;
          if((w>=1)&&(inv[2]==0)){"smooth at (0,0,1)";}
          if(inv[4]!=0)
            {
              OMPButNodes = insert(OMPButNodes,list(inv[4],maxideal(1)),
                                   sizeOMPButNodes);
              sizeOMPButNodes = size(OMPButNodes); // new
            }
          else
            {
              NotOnlyOMPPlusCusps = NotOnlyOMPPlusCusps + 1;
            }
        }
        else
        {
           setring newR;
           list pr=minAssGTZ(radIDsing);
           if(w>=1){pr;}

           for(k=1;k<=size(pr);k++)
           {
              if(w>=1){nsing=nsing+vdim(std(pr[k]));}
              inv=deltaLocMod(I[1],pr[k]);
              delt=delt+inv[1];
              tausing=tausing+inv[2];
              nbranch=nbranch+inv[3];
              if(inv[4]!=0)
                {
                  OMPButNodes = insert(OMPButNodes,list(inv[4],pr[k]),
                                       sizeOMPButNodes);
                  sizeOMPButNodes = size(OMPButNodes);
                }
              else
                {
                  NotOnlyOMPPlusCusps = NotOnlyOMPPlusCusps + 1;
                }
           }
           tau=vdim(std(I+jacob(I)));
           tauinf=tauinf+tau;
           nodes=tau-tausing-2*cusps;
           nsing=nsing+nodes;
           delt=delt+nodes+cusps;
           nbranch=nbranch+2*nodes+cusps;
        }
   }
   genus=genus-delt-deltainf;
   if(w>=1)
   {
      "The projected plane curve has locally:";"";
      "singularities:";nsing;
      "branches:";nbranch+nbranchinf;
      "nodes:"; nodes+nodesinf;
      "cusps:";cusps+cuspsinf;
      "Tjurina number:";tauinf;
      "Milnor number:";2*(delt+deltainf)-nbranch-nbranchinf+nsing;
      "delta of the projected curve:";delt+deltainf;
      //"delta of the curve:";p_a-genus;
      "genus:";genus;
      "====================================================";
      "";
   }
   setring R;
   if(sizeOMPButNodes>0)
     {
       list OMPButNodes = fetch(newR,OMPButNodes);
     }
   return(list(genus,NotOnlyOMPPlusCusps,OMPButNodes,
               fetch(newR,NodesPlusCusps)));
}


///////////////////////////////////////////////////////////////////////////////
static proc deltaLocMod(poly f,ideal singL)
"USAGE:  deltaLoc(f,J);  f poly, J ideal
ASSUME: f is reduced bivariate polynomial; basering has exactly two variables;
        J is irreducible prime component of the singular locus of f (e.g., one
        entry of the output of @code{minAssGTZ(I);}, I = <f,jacob(f)>).
RETURN:  list L:
@texinfo
@table @asis
@item @code{L[1]}; int:
         the sum of (local) delta invariants of f at the (conjugated) singular
         points given by J.
@item @code{L[2]}; int:
         the sum of (local) Tjurina numbers of f at the (conjugated) singular
         points given by J.
@item @code{L[3]}; int:
         the sum of (local) number of branches of f at the (conjugated)
         singular points given by J.
@item @code{L[3]}; int:
         the multiplicity of f at the (conjugated) singular points given by J,
         if these are ordinary multiple points, and 0 otherwise.
@end table
@end texinfo
NOTE:    procedure makes use of @code{execute}; increasing printlevel displays
         more comments (default: printlevel=0).
SEE ALSO: deltaLoc, delta, tjurina
KEYWORDS: delta invariant; Tjurina number
"
{
   option(redSB);
   def R=basering;
   execute("ring S=("+charstr(R)+"),(x,y),lp;");
   map phi=R,x,y;
   ideal singL=phi(singL);
   singL=simplify(std(singL),1);
   attrib(singL,"isSB",1);
   int d=vdim(singL);
   poly f=phi(f);
   int i;
   int w = printlevel-voice+2;  // w=printlevel (default: w=0)
   if(d==1)
   {
      map alpha=S,var(1)-singL[2][2],var(2)-singL[1][2];
      f=alpha(f);
      execute("ring C=("+charstr(S)+"),("+varstr(S)+"),ds;");
      poly f=imap(S,f);
      ideal singL=imap(S,singL);
      if((w>=1)&&(ord(f)>=2))
      {
        "local analysis of the singularities";"";
        basering;
        singL;
        f;
        pause();
      }
   }
   else
   {
      poly p;
      poly c;
      map psi;
      number co;

      while((deg(lead(singL[1]))>1)&&(deg(lead(singL[2]))>1))
      {
         psi=S,x,y+random(-100,100)*x;
         singL=psi(singL);
         singL=std(singL);
          f=psi(f);
      }

      if(deg(lead(singL[2]))==1)
      {
         p=singL[1];
         c=singL[2]-lead(singL[2]);
         co=leadcoef(singL[2]);
      }
      if(deg(lead(singL[1]))==1)
      {
         psi=S,y,x;
         f=psi(f);
         singL=psi(singL);
         p=singL[2];
         c=singL[1]-lead(singL[1]);
         co=leadcoef(singL[1]);
      }

      execute("ring B=("+charstr(S)+"),a,dp;");
      map beta=S,a,a;
      poly p=beta(p);

      execute("ring C=("+charstr(S)+",a),("+varstr(S)+"),ds;");
      number p=number(imap(B,p));
      minpoly=p;

      map iota=S,a,a;
      number c=number(iota(c));
      number co=iota(co);

      map alpha=S,x-c/co,y+a;
      poly f=alpha(f);
      f=cleardenom(f);
      if((w>=1)&&(ord(f)>=2))
      {
        "local analysis of the singularities";"";
        basering;
        alpha;
        f;
        pause();
        "";
      }
   }
   int intMult = deg(lead(f));
   poly fdummy = f;
   poly gdummy = lead(f);
   int ivr = 1;
   while(ivr)
      {
        fdummy = fdummy - lead(fdummy);
        if((fdummy ==0) || (deg(lead(fdummy))>intMult)){break;}
        gdummy = gdummy + lead(fdummy);
      }
   poly SQRpart = sqrfree(gdummy, 3);
   int IntForRet;
   if(deg(SQRpart)==intMult)
     {
        IntForRet = intMult;
     }
   option(noredSB);
   ideal fstd=std(ideal(f)+jacob(f));
   poly hc=highcorner(fstd);
   int tau=vdim(fstd);
   int o=ord(f);
   int delt,nb;

   if(tau==0)                 //smooth case
   {
      setring R;
      return(list(0,0,1,0));
   }
   if((char(basering)>=181)||(char(basering)==0))
   {
      if(o==2)                //A_k-singularity
      {
        if(w>=1){"A_k-singularity";"";}
         setring R;
         delt=(tau+1) div 2;
         return(list(d*delt,d*tau,d*(2*delt-tau+1),IntForRet));
      }
      if((lead(f)==var(1)*var(2)^2)||(lead(f)==var(1)^2*var(2)))
      {
        if(w>=1){"D_k- singularity";"";}

         setring R;
         delt=(tau+2) div 2;
         return(list(d*delt,d*tau,d*(2*delt-tau+1),IntForRet));
      }

      int mu=vdim(std(jacob(f)));
      poly g=f+var(1)^mu+var(2)^mu;  //to obtain a convenient Newton-polygon

      list NP=newtonpoly(g);
      if(w>=1){"Newton-Polygon:";NP;"";}
      int s=size(NP);

      if(is_NND(f,mu,NP))
      { // the Newton-polygon is non-degenerate
        // compute nb, the number of branches
        for(i=1;i<=s-1;i++)
        {
          nb=nb+gcd(NP[i][2]-NP[i+1][2],NP[i][1]-NP[i+1][1]);
        }
        if(w>=1){"Newton-Polygon is non-degenerated";"";}
        return(list(d*(mu+nb-1) div 2,d*tau,d*nb,IntForRet));
      }

      if(w>=1){"Newton-Polygon is degenerated";"";}

      // the following can certainly be made more efficient when replacing
      // 'hnexpansion' (used only for computing number of branches) by
      // successive blowing-up + test if Newton polygon degenerate:
      if(s>2)    //  splitting of f
      {
         if(w>=1){"Newton polygon can be used for splitting";"";}
         intvec v=NP[1][2]-NP[2][2],NP[2][1];
         int de=w_deg(g,v);
         int st=w_deg(hc,v)+v[1]+v[2];
         poly f1=var(2)^NP[2][2];
         poly f2=jet(g,de,v)/var(2)^NP[2][2];
         poly h=g-f1*f2;
         de=w_deg(h,v);
         poly k;
         ideal wi=var(2)^NP[2][2],f2;
         matrix li;
         while(de<st)
         {
           k=jet(h,de,v);
           li=lift(wi,k);
           f1=f1+li[2,1];
           f2=f2+li[1,1];
           h=g-f1*f2;
           de=w_deg(h,v);
         }
         nb=deltaLocMod(f1,maxideal(1))[3]+deltaLocMod(f2,maxideal(1))[3];
         setring R;
         return(list(d*(mu+nb-1) div 2,d*tau,d*nb,IntForRet));
      }

      f=jet(f,deg(hc)+2);
      if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";}
      ideal fac=factorize(f,1);
      if(size(fac)>1)
      {
         nb=0;
         for(i=1;i<=size(fac);i++)
         {
            nb=nb+deltaLocMod(fac[i],maxideal(1))[3];
         }
         setring R;
         return(list(d*(mu+nb-1) div 2,d*tau,d*nb,IntForRet));
      }
      list HNEXP=hnexpansion(f);
      if (typeof(HNEXP[1])=="ring") {
        def altring = basering;
        def HNEring = HNEXP[1]; setring HNEring;
        nb=size(hne);
        setring R;
        kill HNEring;
      }
      else
      {
        nb=size(HNEXP);
      }
      return(list(d*(mu+nb-1) div 2,d*tau,d*nb,IntForRet));
   }
   else             //the case of small characteristic
   {
      f=jet(f,deg(hc)+2);
      if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";}
      delt=delta(f);
      return(list(d*delt,d*tau,d,IntForRet));
   }
}

///////////////////////////////////////////////////////////////////////////////
proc paraConic(poly q)
"USAGE:  paraConic(q); q poly
ASSUME:  The basering must be a polynomial ring in three variables with
         coefficients in Q. @*
         The polynomial q must be homogeneous of degree 2 and absolutely
         irreducible. @*
NOTE:    The procedure might fail or give a wrong output if the assumptions
         do not hold.

RETURN:  ring with an ideal PARACONIC. The ring should be considered as the
         homogeneous coordinate ring of PP^1, the ideal defines a rational
         parametrization PP^1 --> C2 = {q=0}.

THEORY:  We compute a point on C2 via @ref{rationalPointConic}. The pencil of
         lines through this point projects C2 birationally to PP^1. Inverting
         the projection gives the result.
KEYWORDS: conic, parametrization, rational point.
SEE ALSO: rationalPointConic.
EXAMPLE: example paraConic; shows an example
"
{
  def Roriginal = basering;
  def RP2 = projConic(q);  // ring with ideal Ipoint
                           // possibly defined over algebraic number field
  setring RP2;
  def rp1 = invertBirMap(Ipoint, ideal(fetch(Roriginal,q)));
  setring rp1;
  list rl = ringlist(rp1);
  rl[2] = list("s","t");
  def RP1 = ring(rl);
  setring RP1;
  ideal PARACONIC = fetch(rp1,psi);
  export(PARACONIC);
  "// 'paraConic' created a ring together with an ideal RNC.";
  "// Supposing you typed, say,  def RP1 = paraConic(q);";
  "// you may access the ideal by typing";
  "//      setring RP1; PARACONIC;";
  return(RP1);
}
example
{ "EXAMPLE:"; echo=2;
  ring R = 0,(x,y,z),dp;
  poly f = y^8-x^3*(z+x)^5;
  ideal adj = adjointIdeal(f);
  def Rn = invertBirMap(adj,ideal(f));
  setring(Rn);
  J;
  def Rc = rncItProjEven(J);
  PHI;
  setring Rc;
  CONIC;
  def RPc = paraConic(CONIC);
  setring RPc;
  PARACONIC;
}

///////////////////////////////////////////////////////////////////////////////
static proc projConic(poly q)
"USAGE:  projConic(q); q poly
ASSUME:  The basering must be a polynomial ring in three variables with
         coefficients in Q. @*
         The polynomial q must be homogeneous of degree 2 and absolutely
         irreducible. @*
NOTE:    The procedure might fail or give a wrong output if the assumptions
         do not hold.
RETURN:  ring with an ideal Ipoint defining a pencil of lines through a point
         on the conic C2 = {q=0}. This point has either coefficients in Q or
         in a quadratic extension field of Q.
THEORY:  We compute the point on C2 via @ref rationalPointConic.
KEYWORDS: conic, parametrization, rational point.
SEE ALSO: rationalPointConic.
"
{
  def Roriginal = basering;
  list rl = ringlist(Roriginal);
  rl[3] = list(list("dp",1:3),list("C",0));
  def RP20 = ring(rl);
  setring RP20;
  poly q = imap(Roriginal,q);
  def RP21 = rationalPointConic(q);  //  ring with ideal point representing
                                     //  point on conic
                                     //  possibly defined over algebraic number
                                     //  field
  setring RP21;
  list rl1 = ringlist(RP21);
  rl1[2] = list("u","v","w");
  rl1[3] = list(list("dp",1:3),list("C",0));
  def RP2 = ring(rl1);
  setring RP2;
  ideal point = fetch(RP21,point);
  matrix bP = syz(point);
  ideal Ipoint = matrix(maxideal(1))*bP;  // defines pencil of lines through
                                          // point
  export(Ipoint);
  return(RP2);
}


///////////////////////////////////////////////////////////////////////////////
static proc isIrreducible(poly f)
"USAGE:  isIrreducible(f); f poly
RETURN:  1 iff the given polynomial f is irreducible; 0 otherwise.
THEORY:  This test is performed by computing the absolute factorization of f.
KEYWORDS: irreducible.
"
{
  def r = basering;
  def s = absFactorize(f);
  setring s;
  list L = absolute_factors;
  int result = 0;
  if (L[4] == 1){result = 1;}
  setring r;
  kill s;
  return (result);
}


///////////////////////////////////////////////////////////////////////////////
static proc isQuadratic(poly p)
"USAGE:  isQuadratic(p); p poly
RETURN:  checks whether p is a homogeneous, quadratic polynomial in the
         first three ring variables, x, y, z say;
         If so, the method extracs the coefficients a, b, ..., f such that
         p = a*x2 + b*xy + c * y2 + d * xz + e * yz + f * z2
         and returns them as a list of seven entries, [1, a, b, c, d, e, f];
         otherwise, a list with the single entry [0] is returned
"
{
  bigint a = bigint(leadcoef(subst(p, var(2), 0, var(3), 0)));
  bigint c = bigint(leadcoef(subst(p, var(1), 0, var(3), 0)));
  bigint f = bigint(leadcoef(subst(p, var(1), 0, var(2), 0)));
  poly h = p - a * var(1)^2 - c * var(2)^2 - f * var(3)^2;
  bigint b = bigint(leadcoef(subst(h, var(3), 0)));
  bigint d = bigint(leadcoef(subst(h, var(2), 0)));
  bigint e = bigint(leadcoef(subst(h, var(1), 0)));
  list L = 0;
  if (h - b * var(1) * var(2) - d * var(1) * var(3)
        - e * var(2) * var(3) != 0) { return (L); }
  L = 1, a, b, c, d, e, f;
  return (L);
}


///////////////////////////////////////////////////////////////////////////////
static proc largestSquare(bigint n)
"USAGE:  largestSquare(n); n bigint
ASSUME:  n <> 0
RETURN:  returns the largest positive number m (as bigint) such that m^2
         divides n.
THEORY:  This computation is done by prime factorization of n.
KEYWORDS: prime factorization.
"
{
  if (n == 0) { "ERROR: largestSquare(0) had been invoked"; }

  bigint nn = n; if (nn < 0) { nn = -n; }
  list L = primefactors(nn);
  if (L[3] != 1)
  { "WARNING: command 'primefactors(.)' did not find all prime factors"; }
  int i; bigint m = bigint(1); int e; int j;
  for (i = 1; i <= size(L[1]); i++)
  {
    e = L[2][i] div 2;
    for (j = 1; j <= e; j++) { m = m * bigint(L[1][i]); }
  }
  return (m);
}


///////////////////////////////////////////////////////////////////////////////
static proc jIndex(bigint a, bigint b, bigint c)
"USAGE:  jIndex(a, b, c); a, b, c bigint's
RETURN:  returns the middle of the three numbers |ab|, |bc|, and |ca|.
"
{
  bigint n1 = a*b; if (n1 < 0) { n1 = -n1; }
  bigint n2 = b*c; if (n2 < 0) { n2 = -n2; }
  bigint n3 = c*a; if (n3 < 0) { n3 = -n3; }
  if ((n1 <= n2) && (n2 <= n3)) { return (n2); }
  if ((n1 >= n2) && (n2 >= n3)) { return (n2); }
  if ((n2 <= n1) && (n1 <= n3)) { return (n1); }
  if ((n2 >= n1) && (n1 >= n3)) { return (n1); }
  return (n3);
}


///////////////////////////////////////////////////////////////////////////////
static proc aMod(bigint a, bigint b)
"USAGE:  aMod(a,b); a, b bigint
RETURN:  r bigint
THEORY:  The asymmetric residue r of the division with remainder a mod b.
"
{
  bigint c = a mod b;
  if (c<0)
  {
    return(c+b);
  }
return(c);
}

///////////////////////////////////////////////////////////////////////////////
static proc aDiv(bigint a, bigint b)
"USAGE:  aDiv(a,b); a, b bigint
RETURN:  q bigint
THEORY:  Quotient with remainder q = a div b with asymmetric residue.
"
{
  bigint q = a div b;
  if ((a mod b)<0)
  {
    return(q-1);
  }
return(q);
}



///////////////////////////////////////////////////////////////////////////////
static proc polyModP(poly q, bigint p)
"USAGE:  polyModP(q, p); q poly, p bigint
RETURN:  takes each coefficient of q modulo p and returns the resulting poly
"
{
  poly qq = q; poly res = 0;
  bigint c;
  while (qq != 0)
  {
    c = bigint(leadcoef(qq) mod p);
    res = res + c * leadmonom(qq);
    qq = qq - lead(qq);
  }
  return (res);
}


///////////////////////////////////////////////////////////////////////////////
static proc rootModP(bigint r, bigint p)
"USAGE:  rootModP(r, p); r, p bigint's
ASSUME:  0 <= r < p, and p prime;
         Furthermore it is assumes that there is some x in {0, 1, ..., p-1}
         such that x^2 = r mod p;
RETURN:  an x in {0, 1, ..., p-1} such that x^2 = r mod p;
THEORY:  For p larger than 32003, this computation is done using Cantor-
         Zassenhaus' algorithm. Otherwise a brute force approach is used.
KEYWORDS: Cantor-Zassenhaus algorithm.
"
{
  if (r == 0) { return (0); }
  if (r == 1) { return (1); }
  if (p <= 32003)
  {
    /* For small p, we use a brute force approach: */
    int i;
    for (i = 2; i < p; i++)
    {
      if (((i*i) mod p) == r) { return (i); }
    }
    /* should never be reached: */
    return (-1);
  }

  /* For p > 32003, we use Cantor-Zassenhaus' algorithm: */
  def br = basering;
  ring rTemp = 0, x, dp;
  bigint b; bigint exponent; poly factor;
  poly h = x^2 - r;
  ideal redI = h; redI = std(redI);
  poly q = x^2; bigint root = 0;
  while (root == 0)
  {
    b = bigint(random(1, 2^30));
    exponent = bigint((p - 1) div 2);
    /* We need to compute q^exponent mod (x^2 - a) and mod p: */
    factor = x + b; q = 1;
    while (exponent > 0)
    {
      if ((exponent mod 2) == 1)
      {
        q = q * factor;
        q = reduce(q, redI);
        q = polyModP(q, p);
      }
      exponent = bigint(exponent div 2);
      factor = factor * factor;
      factor = reduce(factor, redI);
      factor = polyModP(factor, p);
    }
    if (deg(q) == 1)
    {
      q = q - 1;
      b = inverseModP(bigint(leadcoef(q)), p);
      q = q - lead(q);
      root = aMod((bigint(q) * b),p);
      if (((root * root - r) mod p) != 0) { root = 0; }
    }
  }
  setring br; kill rTemp;
  return (root);
}


///////////////////////////////////////////////////////////////////////////////
static proc inverseModP(bigint r, bigint p)
"USAGE:  inverseModP(r, p); r, p bigint's
ASSUME:  0 <= r < p, and r and p coprime;
RETURN:  returns the inverse of r in Z/p represented by an element in
         {1, 2, ..., p-1}
THEORY:  This uses Euclid's extended gcd algorithm.
"
{
  list L = extgcd(r, p);
  if (L[1] != 1) { ERROR("GCD of", r, "and", p, "should be 1."); }
  L[2] = aMod(L[2],p);
  return (L[2]);
}


///////////////////////////////////////////////////////////////////////////////
static proc squareRoot(bigint r, bigint m, int justCheck)
"USAGE:  squareRoot(r, m, j); r, m bigint's, j int
RETURN:  checks whether r is a square modulo m, i.e., checks whether there is
         some x such that x^2 = r mod m;
         If justCheck is 1, then the method will terminate after the check
         and return 1 if r is a square and -1 otherwise.
         If justCheck is 0 and r is a square, then the method continues and
         computes a solution x in {0, 1, m-1} with x^2 = r mod m, which will
         then be returned
THEORY:  This algorithm checks solvability by computing the Legendre symbols
         modulo all primes in m. Then, individual roots will be computed and
         lifted to the desired square root modulo m using Chinese
         remaindering.
"
{
  if (m == 0) { "ERROR: squareRoot had been invoked with m = 0"; }

  list L = primefactors(m);
  if ((L[3] != 1) && (L[3] != -1))
  { "WARNING: command 'primefactors(.)' did not find all prime factors"; }
  int i;
  for (i = 1; i <= size(L[2]); i++)
  {
    if (legendreSymbol(r, L[1][i]) == -1) { return (-1); }
  }
  /* now we know that there is some x in {0, 1, m-1} with
     x^2 = r mod m */
  if (justCheck == 1) { return (1); }
  else
  {
    // now we need to compute x; this works in two stages:
    // 1) write m = p1^e1 * ... * pk^ek (we already have that),
    // 2) for each i in {1, 2, ..., k}
    //    2.1) compute a yi such that yi^2 = r mod pi,
    //    2.2) lift yi to an xi such that xi^2 = r mod (pi^ei),
    // 3) lift (x1, x2, ..., xk) in Z/p1^e1 * ... * Z/pk^ek
    //    to x in Z/m via Chinese remainder theorem

    list roots;
    // 2.1):
    for (i = 1; i <= size(L[1]); i++)
    {
      roots = insert(roots, rootModP(aMod(r,L[1][i]), L[1][i]), size(roots));
    }

    // 2.2):
    bigint c; bigint l; bigint temp; bigint pPower; int e;
    for (i = 1; i <= size(roots); i++)
    {
      pPower = bigint(L[1][i]);
      for (e = 2; e <= L[2][i]; e++)
      {
        c = bigint(roots[i]); l = pPower;
        temp = r - c * c; l = bigint(2) * c * l; c = temp;
        c = aDiv(c,pPower); l = aDiv(l,pPower);
        c = aMod(c,L[1][i]); l = aMod(l,L[1][i]);
        c = aMod((c * bigint(inverseModP(l, L[1][i]))), L[1][i]);
        c = bigint(roots[i]) + c * pPower;
        pPower = pPower * L[1][i]; roots[i] = c;
      }
    }

    // 2.3):
    list mm; bigint z; int j;
    for (i = 1; i <= size(L[1]); i++)
    {
      z = bigint(L[1][i]);
      for (j = 2; j <= L[2][i]; j++)
      {
        z = z * bigint(L[1][i]);
      }
      mm = insert(mm, z, size(mm));
    }
    return (aMod(chinrem(roots, mm) , m));
  }
}


///////////////////////////////////////////////////////////////////////////////
static proc chineseRemainder(list rr, list mm)
"USAGE:  chineseRemainder(rr, mm); rr, mm lists of bigint's
ASSUME:  lists rr and mm must have same sizes;
         Furthermore the entries of mm must be mutually coprime.
RETURN:  an x which fulfills the simultaneous remainder conditions
         x = rr[i] mod mm[i], 1 <= i <= size(rr)
KEYWORDS: Chinese remainder.
"
{
  bigint x = bigint(0); int i; bigint N; list l;
  bigint M = bigint(mm[1]);
  for (i = 2; i <= size(mm); i++) { M = M * bigint(mm[i]); }
  for (i = 1; i <= size(mm); i++)
  {
    N = aDiv(M,mm[i]);
    l = extgcd(mm[i], N);
    x = x + rr[i]*l[3]*N;
  }
  return (x);
}


///////////////////////////////////////////////////////////////////////////////
static proc rationalPointSpecial(bigint b1, bigint c1)
"USAGE:  rationalPointSpecial(b1, c1); b1, c1 bigint's
ASSUME:  b1 <> 0 and c1 <> 0;
RETURN:  with poly p = var(1)^2 + b1 * var(2)^2 + c1 * var(3)^2, the method
         returns a list L with either one entry or four entries:
         case 'three entries':
           L[1] = 0 signaling that there is no rational point on V(p),
           L[2] the largest number b such that b^2 divides b1
                (for subsequent use by the caller of this method),
           L[3] the largest number c such that c^2 divides c1
                (for subsequent use by the caller of this method);
         case 'four entries':
           L[1] = 1 signaling that there is a rational point on V(p),
           L[2], L[3], L[4] rational numbers such that the tuple
                (L[2], L[3], L[4]) is on V(p)
"
{
  if (b1 == 0) { "ERROR: rationalPointSpecial(0, c1) had been invoked"; }
  if (c1 == 0) { "ERROR: rationalPointSpecial(b1, 0) had been invoked"; }

  bigint b_s = largestSquare(b1); bigint b_r = b1/b_s/b_s;
  bigint c_s = largestSquare(c1); bigint c_r = c1/c_s/c_s;
  bigint g = gcd(b_r, c_r);
  def S=basering;
  ideal mi = maxideal(1);
  map mm = basering, mi; map mTemp;
  mm[1] = var(1); mm[2] = var(2)/b_s/g; mm[3] = var(3)/c_s/g;
  bigint a = g;     bigint aa = a; if (aa <= 0) { aa = -aa; }
  bigint b = b_r/g; bigint bb = b; if (bb <= 0) { bb = -bb; }
  bigint c = c_r/g; bigint cc = c; if (cc <= 0) { cc = -cc; }
  bigint R1 = squareRoot(-a*b, cc, 1);
  if (R1 == -1) { list L = 0, b_s, c_s; return (L); }
  bigint R2 = squareRoot(-a*c, bb, 1);
  if (R2 == -1) { list L = 0, b_s, c_s; return (L); }
  bigint R3 = squareRoot(-b*c, aa, 1);
  if (R3 == -1) { list L = 0, b_s, c_s; return (L); }
  bigint t; bigint r1; bigint Q; bigint A; bigint B; bigint C;
  bigint alpha; bigint beta; bigint gamma;
  while (jIndex(a, b, c) > 1)
  {
    mTemp = basering, mi;
    if (aa > cc)
    {
      t = a; a = c; c = t;
      t = aa; aa = cc; cc = t;
      mTemp = basering, mi;
      mTemp[1] = var(3); mTemp[3] = var(1); mm = mTemp(mm);
    }
    if (bb > cc)
    {
      t = b; b = c; c = t;
      t = bb; bb = cc; cc = t;
      mTemp = basering, mi;
      mTemp[2] = var(3); mTemp[3] = var(2); mm = mTemp(mm);
    }
    if (bb < aa)
    {
      t = b; b = a; a = t;
      t = bb; bb = aa; aa = t;
      mTemp = basering, mi;
      mTemp[1] = var(2); mTemp[2] = var(1); mm = mTemp(mm);
    }
    /* now, we have established |a| <= |b| <= |c|; and permuted
       the map mm, accordingly */
    cc = c; if (cc <= 0) { cc = -cc; }
    R1 = squareRoot(-a*b, cc, 0);
    r1 = aMod((R1 * inverseModP(a, cc)), cc);
    if (r1*bigint(2) > cc) { r1 = r1 - cc; }
    Q = (a*r1*r1 + b)/c;
    if (Q == 0)
    {
      list L = 1, subst(mm[1], var(1), 1, var(2), 1, var(3), 0),
                  subst(mm[2], var(1), 1, var(2), 1, var(3), 0),
                  subst(mm[3], var(1), 1, var(2), 1, var(3), 0);
      return (L);
    }
    A = gcd(gcd(a*r1*r1, b), c*Q);
    alpha = r1/A; beta = b/A;
    B = a*beta;
    gamma = largestSquare(Q/A);
    C = Q/A/gamma/gamma;
    mTemp = basering, mi;
    mTemp[1] = A*alpha*var(1) - beta*var(2);
    mTemp[2] = var(1) + a*alpha*var(2);
    mTemp[3] = C*gamma*var(3);
    mm = mTemp(mm);
    a = A; b = B; c = C;
    aa = a; if (aa <= 0) { aa = -aa; }
    bb = b; if (bb <= 0) { bb = -bb; }
    cc = c; if (cc <= 0) { cc = -cc; }
  }
  if (a*b < 0)
  {
    list L = 1, subst(mm[1], var(1), 1, var(2), 1, var(3), 0),
                subst(mm[2], var(1), 1, var(2), 1, var(3), 0),
                subst(mm[3], var(1), 1, var(2), 1, var(3), 0);
    return (L);
  }
  if (a*c < 0)
  {
    list L = 1, subst(mm[1], var(1), 1, var(2), 0, var(3), 1),
                subst(mm[2], var(1), 1, var(2), 0, var(3), 1),
                subst(mm[3], var(1), 1, var(2), 0, var(3), 1);
    return (L);
  }
  if (b*c < 0)
  {
    list L = 1, subst(mm[1], var(1), 0, var(2), 1, var(3), 1),
                subst(mm[2], var(1), 0, var(2), 1, var(3), 1),
                subst(mm[3], var(1), 0, var(2), 1, var(3), 1);
    return (L);
  }
  list L = 0, b_s, c_s; return (L);
}


///////////////////////////////////////////////////////////////////////////////
static proc extendedEuclid(bigint a, bigint b)
"USAGE:  extendedEuclid(a, b); a, b bigint's
ASSUME:  a <> 0 or b <> 0;
RETURN:  returns a list with three entries:
         _[1]: gcd(a,b) > 0,
         _[2], _[3]: s, t, such that s*a + t*b = gcd(a,b)
KEYWORDS: extended Euclidean algorithm.
"
{
  list l = 0; bigint temp;
  if (a == 0)         { l = b, 0, 1; if (b < 0) { l = -b, 0, -1; } }
  if (b == 0)         { l = a, 1, 0; if (a < 0) { l = -a, -1, 0; } }
  if (aMod(a , b) == 0) { l = b, 0, 1; if (b < 0) { l = -b, 0, -1; } }
  if (aMod(b , a) == 0) { l = a, 1, 0; if (a < 0) { l = -a, -1, 0; } }
  if (size(l) > 1) { return (l); }

  temp = aMod(a , b);
  l = extendedEuclid(b, temp);
  temp = (a - temp) / b;
  temp = bigint(l[2]) - temp * bigint(l[3]);
  l = l[1], l[3], temp;
  return (l);
}

static proc legendreSymbol(bigint r, bigint p)
"assumes p prime;
returns the Legendre symbol (r/p), that is
 1 if r appears as residue modulo p of a square,
-1 if not,
 0 if r is a multiple of p
"
{
  bigint rr = aMod(r , p);
  if (rr == 0) { return (0) }
  if (rr == 1) { return (1) }
  /* now, p must be at least 3 */
  bigint e = (p - 1) / bigint(2);
  bigint result = 1;
  bigint power = rr;
  while (e > 0)
  {
    if ((e mod 2) == 1) { result = aMod((result * power), p); }
    e = e / bigint(2);
    power = aMod((power * power), p);
  }
  if (result > 1) { result = result - p; /* should be -1 */ }
  return (result);
}


///////////////////////////////////////////////////////////////////////////////
static proc buildExtension(bigint b, bigint c, bigint bs, bigint cs)
"USAGE:  buildExtension(b, c, bs, cs); b, c, bs, cs bigint's
ASSUME:  assumes that bs is the largest positive number such that bs^2
         divides b; analogously for cs regarding c;
         Assumes furthermore that there is no rational point on the conic
         X^2 + b*Y^2 + c*Z^2 = 0.
         Assumes that the ground field of the basering is Q.
RETURN:  builds an appropriate quadratic field extension Q(a) in which a
         point exists that lies on the given conic. This point is stored in
         a newly defined and exported (1x3) matrix named 'point'.
         The method returns the resulting polynomial ring over Q(a).
"
{
  bigint br = b/bs/bs;
  bigint cr = c/cs/cs;
  /* X^2 + br*bs^2*Y^2 + cr*cs^2*Z^2 = 0 */
  def bRing = basering;
  list L = ringlist(bRing);

  if (b != 0)
  {
    L[1] = list(0, list("a"), list(list("lp", 1)), ideal(0));
    def RTemp = ring(L);
    setring RTemp; list L = ringlist(RTemp);
    L[1][4] = ideal(a^2 + br);
    def R = ring(L);
    setring R; kill RTemp;
    matrix point[1][3];
    point[1, 1] = a * bs; point[1, 2] = 1; point[1, 3] = 0;
    export point;
    setring bRing;
    return (R);
  }
  if (c != 0)
  {
    L[1] = list(0, list("a"), list(list("lp", 1)), ideal(0));
    def RTemp = ring(L);
    setring RTemp; list L = ringlist(RTemp);
    L[1][4] = ideal(a^2 + cr);
    def R = ring(L);
    setring R; kill RTemp;
    matrix point[1][3];
    point[1, 1] = a * cs; point[1, 2] = 0; point[1, 3] = 1;
    export point;
    setring bRing;
    return (R);
  }

  "ERROR: unexpectedly encountered conic X^2 + 0*Y^2 + 0*Z^2 = 0";
  return (bRing);
}


///////////////////////////////////////////////////////////////////////////////
static proc testRationalPointConic(poly pp)
"USAGE:  testRationalPointConic(pp); pp poly
RETURN:  returns 0 in case of unexpected input (e.g. non-quadratic,
         reducible); 1 otherwise
NOTE:    This method calles rationalPointConic, measures time consumption
         and checks whether the computed point lies indeed on the conic pp.
         The results will be printed to standard output.
"
{
  "testing rationalPointConic(poly) for poly p:";
  "p =", pp;
  if (isQuadratic(pp)[1] == 1) { "p is quadratic."; }
  else                         { "p is not quadratic.";   return (0); }
  if (isIrreducible(pp) == 1)  { "p is irreducible."; }
  else                         { "p is not irreducible."; return (0); }
  def rOrig = basering;
  int t = rtimer;
  def rNew = rationalPointConic(pp);
  t = rtimer - t;
  "time for finding a point on the conic [sec] =", t;
  setring rNew;
  poly ff = fetch(rOrig, pp);
  if (minpoly == 0)
  { "there is a rational point on the conic p";
    "x =", point[1,1], "  y =", point[1,2], "  z =", point[1,3];
    "check (should be zero):", subst(ff, var(1), point[1,1],
                                         var(2), point[1,2],
                                         var(3), point[1,3]);
  }
  else
  {
    "there is no rational point on the conic p";
    "but there is a point on the conic in the field extension Q(a),";
    "with minpoly =", minpoly;
    "x =", point[1,1], "  y =", point[1,2], "  z =", point[1,3];
    "check (should be zero):", subst(ff, var(1), point[1,1],
                                         var(2), point[1,2],
                                         var(3), point[1,3]);
  }
  setring rOrig;
}

example
{ "EXAMPLE:"; echo=2;
  ring r = 0, (x,y,z, u, v, w), dp;
  poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z;
  testRationalPointConic(p);
}

///////////////////////////////////////////////////////////////////////////////
proc rationalPointConic(poly p)
"USAGE:  rationalPointConic(p); p poly
ASSUME:  assumes that p is an irreducible quadratic polynomial in the first
         three ring variables;
         ground field is expected to be Q.
RETURN:  The method finds a point on the given conic. There are two
         possibilities:
         1) There is a rational point on the curve.
         2) There is no rational point on the curve.
         In the second case, the method creates a modification of the current
         basering which is a polynomial ring over some quadratic field
         extension Q(a) of Q. Apart from the replacement of Q by Q(a), the
         new polynomial ring, R say, is the same as the original basering.
         (In the first case, R is identical with the basering.)
         In both cases, the method will then define a (1x3) matrix named
         'point' which lives in R and which contains the coordinates of the
         desired point on q.
         Finally, the method returns the ring R (which will in the 1st case
         be the original base ring).
EXAMPLE: example rationalPointConic; shows an example
"
{
  list L = isQuadratic(p);
  bigint a = bigint(L[2]); bigint b = bigint(L[3]); bigint c = bigint(L[4]);
  bigint d = bigint(L[5]); bigint e = bigint(L[6]); bigint f = bigint(L[7]);
  bigint x; bigint y; bigint z; bigint nn;
  def R = basering;

  if (b^2 == 4*a*c)
  {
    if (c == 0)
    {
      x = -2*d*e; y = d^2-4*a*f; z = e*4*a;
      nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
      matrix point[1][3];
      point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
      export point; return (R);
    }
    else
    {
      bigint fs = 4*c*f - e^2;
      bigint ds = 4*c*d - 2*b*e;
      x = -fs*2*c; y = b*fs-e*ds; z = ds*2*c;
      nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
      matrix point[1][3];
      point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
      export point; return (R);
    }
  }

  if (d^2 == 4*a*f)
  {
    if (f == 0)
    {
      x = -b*e*2; y = e*4*a; z = b^2-4*a*c;
      nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
      matrix point[1][3];
      point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
      export point; return (R);
    }
    else
    {
      bigint c_s = 4*c*f - e^2;
      bigint b_s = 4*f*b - 2*d*e;
      x = -c_s*2*f; y = b_s*2*f; z = d*c_s-e*b_s;
      nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
      matrix point[1][3];
      point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
      export point; return (R);
    }
  }

  if (e^2 == 4*c*f)
  {
    if (c == 0)
    {
      x = b*4*f; y = d^2-4*a*f; z = -b*d*2;
      nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
      matrix point[1][3];
      point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
      export point; return (R);
    }
    else
    {
      bigint as = 4*c*a - b^2;
      bigint ds = 4*c*d - 2*b*e;
      x = ds*2*c; y = e*as-b*ds; z = -as*2*c;
      nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
      matrix point[1][3];
      point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
      export point; return (R);
    }
  }

  ideal mi = maxideal(1);
  map mm = R, mi;
  bigint B; bigint C; bigint D;

  if ((a == 0) && (c == 0))
  {
    B = -1; C = 4*b*f - 4*d*e;
    /* now, b <> 0 since otherwise p would have the factor z,
       and hence not be irreducible */
    mm[1] = (var(1)+var(2)-2*e*var(3))/(2*b);
    mm[2] = (var(1)-var(2)-2*d*var(3))/(2*b);
  }
  if ((a != 0) && (c == 0))
  {
    mm[1] = var(2);
    mm[2] = var(1);
    bigint t = a; a = c; c = t;
    t = e; e = d; d = t;
  }
  if (c != 0)
  {
    D = 4*a*c-b^2;
    mm[2] = (var(2)-e*var(3)-b*var(1))/(2*c);
    map mTemp = basering, mi;
    mTemp[1] = (var(1)-2*d*c*var(3)+b*e*var(3))/D;
    mm = mTemp(mm);
    B = D;
    C = 16*a*c^2*f-4*a*c*e^2-4*b^2*c*f+4*b*c*d*e-4*c^2*d^2;
  }
  list K;
  if ((B > 0) && (C >= 0)) { K = 0; }
  if ((B >= 0) && (C > 0)) { K = 0; }
  if (B == 0)
  {
    /* looking for a point on X^2 = |C| * Z^2 */
    bigint root = largestSquare(absValue(C));
    if (absValue(C)/root/root == 1) { K = 1, root, 0, 1; }
    else                            { K = 0; }
  }
  if (C == 0)
  {
    /* looking for a point on X^2 = |B| * Y^2 */
    bigint root = largestSquare(absValue(B));
    if (absValue(B)/root/root == 1) { K = 1, root, 1, 0; }
    else                            { K = 0; }
  }
  else { K = rationalPointSpecial(B, C); }
  if (K[1] == 0)
  {
    /* no rational point on conic;
       we need to move to an appropriate field extension Q(a) */
    poly h1 = mm[1]; poly h2 = mm[2]; poly h3 = mm[3];
    def extendedR = buildExtension(B, C, K[2], K[3]);
    setring extendedR;
    poly g1 = fetch(R, h1);
    poly g2 = fetch(R, h2);
    poly g3 = fetch(R, h3);
    matrix temp[1][3];
    temp[1, 1] = subst(g1, var(1), point[1, 1], var(2), point[1, 2],
                           var(3), point[1, 3]);
    temp[1, 2] = subst(g2, var(1), point[1, 1], var(2), point[1, 2],
                           var(3), point[1, 3]);
    temp[1, 3] = subst(g3, var(1), point[1, 1], var(2), point[1, 2],
                           var(3), point[1, 3]);
    point[1, 1] = temp[1, 1]; point[1, 2] = temp[1, 2];
    point[1, 3] = temp[1, 3];
    setring R;
    return (extendedR);
  }
  else
  {
    string dummyString = string(K); // without this useless line, we
                                    // sometimes get a seg fault because
                                    // mm is corrupted; strange!?!?!?!?
    number nx = number(subst(mm[1], var(1), K[2], var(2), K[3], var(3), K[4]));
    number ny = number(subst(mm[2], var(1), K[2], var(2), K[3], var(3), K[4]));
    number nz = number(subst(mm[3], var(1), K[2], var(2), K[3], var(3), K[4]));
    /* the point (nx, ny, nz) is already a solution;
       the following lines will just remove denominators and reduce
       numerators in order to return a nice tuple from Z^3 */
    bigint nxd = bigint(denominator(absValue(nx)));
    bigint nyd = bigint(denominator(absValue(ny)));
    bigint nzd = bigint(denominator(absValue(nz)));
    nn = nxd * nyd / gcd(nxd, nyd);
    nn =  nn * nzd / gcd(nn, nzd);
    x = bigint(nx*nn); y = bigint(ny*nn); z = bigint(nz*nn);
    nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
    matrix point[1][3];
    point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
    export point;
    return (R);
  }
}

example
{ "EXAMPLE:"; echo=2;
ring R = 0, (x,y,z), dp;
system("random", 4711);
poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z;
def S = rationalPointConic(p); // quadratic field extension,
                               // minpoly = a^2 - 2
testPointConic(p, S);
setring R;
p = x^2 - 1857669520 * y^2 + 86709575222179747132487270400 * z^2;
S = rationalPointConic(p); // same as current basering,
                           // no extension needed
testPointConic(p, S);
}
///////////////////////////////////////////////////////////////////////////////
proc testParametrization(poly f, def rTT)
"USAGE:  testParametrization(f, rTT); f poly, rTT ring
ASSUME:  The assumptions on the basering and the polynomial f are as required
         by @ref{paraPlaneCurve}. The ring rTT has two variables and contains
         an ideal PARA (such as the ring obtained by applying
         @ref{paraPlaneCurve} to f).
RETURN: int which is 1 if PARA defines a parametrization of the curve
        {f=0} and 0, otherwise.
THEORY: We compute the polynomial defining the image of PARA
        and compare it with f.
KEYWORDS: Parametrization, image.
EXAMPLE: example testParametrization; shows an example
"
{
  def Roriginal = basering;
  setring rTT;
  /* begin workaround elimination*/
  int k;
  list rl = ringlist(rTT);
  rl[2] = list("s","t","x","y","z");
  rl[3]= list(list("dp",1:5),list("C",0));
  def Relim = ring(rl);
  setring Relim;
  ideal PARA = fetch(rTT,PARA);
  ideal JJ;
  for(k=1;k<=3;k++)
     {
       JJ=JJ,var(k+2)-PARA[k];
     }
  ideal SJJ = std(JJ);
  intvec HJJ = hilb(SJJ,1);
  ideal J = eliminate(JJ,var(1)*var(2),HJJ);
  setring rTT;
  /*end workaround elimination*/
  rl[2] = list("x","y","z");
  rl[3] = list(list("dp",1:3),list("C",0));
  def RP2 = ring(rl);
  setring RP2;
  ideal f = fetch(Roriginal,f);
  ideal ftest = imap(Relim,J);
  poly g = reduce(f[1],std(ftest));
  if(g!=0){return(0)}
  g = reduce(ftest[1],std(ideal(f)));
  if(g!=0){return(0)}
  return (1);
}

example
{ "EXAMPLE:"; echo=2;
  ring R = 0,(x,y,z),dp;
  poly f = y^8-x^3*(z+x)^5;
  def RP1 = paraPlaneCurve(f);
  testParametrization(f, RP1);
}

///////////////////////////////////////////////////////////////////////////////
proc testPointConic(poly p, def r)
"USAGE:  testPointConic(p, r); p poly, r ring
ASSUME:  assumes that p is a homogeneous quadratic polynomial in the
        first three ring variables of the current basering;
        Assumes that there is a (1x3) matrix named 'point' in r with
        entries from the ground field of r.
RETURN:  returns 1 iff the point named 'point', residing in r, lies on
        the conic given by p; 0 otherwise
NOTE:    This method temporarily changes the basering to r. Afterwards,
        the basering will be the same as before.
EXAMPLE: example testPointConic; shows an example
"
{
 def rOrig = basering;
 "conic:", p;
 setring r;
 string s = "point: " + string(point[1,1]) + ", " + string(point[1,2]);
 s = s + ", " + string(point[1,3]);
 s;
 if (minpoly != 0) { "minpoly:", minpoly; }
 poly f = fetch(rOrig, p);
 poly g = subst(f, var(1), point[1,1],
                   var(2), point[1,2],
                   var(3), point[1,3]);
 int result = 0; if (g == 0) { result = 1; }
 setring rOrig;
 return (result);
}

example
{ "EXAMPLE:"; echo=2;
 ring R = 0, (x,y,z), dp;
 system("random", 4711);
 poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z;
 def S = rationalPointConic(p);
 if (testPointConic(p, S) == 1)
 { "point lies on conic"; }
 else
 { "point does not lie on conic"; }
}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
/*
/////////////////////////////////////////////////////////////////////////////
/// Further examples for testing the main procedures
/// Timings on wawa Sept 29
/////////////////////////////////////////////////////////////////////////////
LIB"paraplanecurves.lib";
// -------------------------------------------------------
// Example 1
// -------------------------------------------------------
ring RR = 0, (x,y,z), dp;
poly f = 7*z2+11*y2+13*z*y+17*x2+19*x*y; // conic
def RP1 = paraConic(f);
setring RP1; PARACONIC;
setring RR;
RP1 = paraPlaneCurve(f);
testParametrization(f,RP1);
setring RP1; PARA;
kill RR;kill RP1;
// -------------------------------------------------------
// Example 2
// -------------------------------------------------------
ring RR = 0, (x,y,z), dp;
poly f = y3-x2z;  // cusp at origin
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 0
testParametrization(f,RP1);
setring RP1; PARA;
kill RR;kill RP1;
// -------------------------------------------------------
// Example 3
// -------------------------------------------------------
ring RR = 0, (x,y,z), dp;
poly f=(xz-y^2)^2-x*y^3; // 1 sing at origin, 1 cusp, no OMPs
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f); // time 0
testParametrization(f,RP1);
setring RP1; PARA;
kill RR;kill RP1;
// -------------------------------------------------------
// Example 4
// -------------------------------------------------------
ring RR = 0, (x,y,z), dp;
poly f = y5-y4x+4y2x2z-x4z;  // 1 sing at origin, no OMPs, no cusps
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 0
testParametrization(f,RP1);
setring RP1; PARA;
kill RR;kill RP1;
// -------------------------------------------------------
// Example 5
// -------------------------------------------------------
ring RR = 0, (x,y,z), dp;
poly f = 259x5-31913x4y+939551x3y2+2871542x2y3+2845801xy4;
f = f+914489y5+32068x4z-1884547x3yz-8472623x2y2z-11118524xy3z;
f = f-4589347y4z+944585x3z2+8563304x2yz2+16549772xy2z2+9033035y3z2;
f = f-2962425x2z3-11214315xyz3-8951744y2z3+2937420xz4+4547571yz4-953955z5;
// 6 nodes
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 0
testParametrization(f,RP1);
setring RP1; PARA;
kill RR;kill RP1;
// -------------------------------------------------------
// Example 7
// -------------------------------------------------------
ring RR = 0, (x,y,z), dp;
poly f = y^8-x^3*(z+x)^5;  // 1 sing at origin, 1 further sing, no OMPs,
                           // no cusps
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 0
testParametrization(f,RP1);
setring RP1; PARA;
kill RR;kill RP1;
// -------------------------------------------------------
// Example 8
// -------------------------------------------------------
ring RR = 0, (x,y,z), dp;
poly f = 11y7+7y6x+8y5x2-3y4x3-10y3x4-10y2x5-x7-33y6-29y5x-13y4x2+26y3x3;
f = f+30y2x4+10yx5+3x6+33y5+37y4x-8y3x2-33y2x3-20yx4-3x5-11y4-15y3x;
f = f+13y2x2+10yx3+x4; // 3 OMPs of mult 3, 1 OMP of mult 4
f = homog(f,z);
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 0
testParametrization(f,RP1);
setring RP1; PARA;
kill RR;kill RP1;
// -------------------------------------------------------
// Example 9
// -------------------------------------------------------
ring RR = 0, (x,y,z), dp;
poly f = y^8-x^3*(z+x)^5;  // 1 sing at origin, 1 further sing, no OMPs,
                           // no cusps
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 0
testParametrization(f,RP1);
setring RP1; PARA;
kill RR;kill RP1;
// -------------------------------------------------------
// Example 10
// -------------------------------------------------------
ring SS = 0, (u,v,z), dp;
poly f = u^4-14*u^2*v^2+v^4+8*u^2*v*z+8*v^3*z; // 1 OMP of mult 3 at orgin
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 0
testParametrization(f,RP1);
setring RP1; PARA;
kill SS;kill RP1;
// -------------------------------------------------------
// Example 11
// -------------------------------------------------------
ring SS = 0, (u,v,z), dp;
poly f = 14440*u^5-16227*u^4*v+10812*u^3*v^2-13533*u^2*v^3+3610*u*v^4;
f = f+1805*v^5+14440*u^4*z-18032*u^3*v*z+16218*u^2*v^2*z-12626*u*v^3*z;
f = f+3610*v^4*z+3610*u^3*z^2-4508*u^2*v*z^2+5406*u*v^2*z^2-2703*v^3*z^2;
// 1 OMP of mult 3 at origin, 2 nodes
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 0
testParametrization(f,RP1);
setring RP1; PARA;
kill SS;kill RP1;
// -------------------------------------------------------
// Example 12
// -------------------------------------------------------
ring SS = 0, (u,v,z), dp;
poly f = u^6+3*u^4*v^2+3*u^2*v^4+v^6-4*u^4*z^2-34*u^3*v*z^2-7*u^2*v^2*z^2;
f = f+12*u*v^3*z^2+6*v^4*z^2+36*u^2*z^4+36*u*v*z^4+9*v^2*z^4;
// needs field extension *** 6 nodes, 2 cusps, 1 sing at 0
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 0
testParametrization(f,RP1);
setring RP1; PARA;
kill SS;kill RP1;
// -------------------------------------------------------
// Example 13
// -------------------------------------------------------
ring SS = 0, (u,v,z), dp;
poly f = -24135/322*u^6-532037/6440*u^5*v+139459/560*u^4*v^2;
f = f-1464887/12880*u^3*v^3+72187/25760*u^2*v^4+9/8*u*v^5+1/8*v^6;
f = f-403511/3220*u^5*z-40817/920*u^4*v*z+10059/80*u^3*v^2*z;
f = f-35445/1288*u^2*v^3*z+19/4*u*v^4*z+3/4*v^5*z-20743/805*u^4*z^2;
f = f+126379/3220*u^3*v*z^2-423417/6440*u^2*v^2*z^2+11/2*u*v^3*z^2;
f = f+3/2*v^4*z^2+3443/140*u^3*z^3+u^2*v*z^3+u*v^2*z^3+v^3*z^3;
// 2 OMPs of mult 3 (1 at origin), 4 nodes
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 14
testParametrization(f,RP1);
setring RP1; PARA;
kill SS;kill RP1;
// -------------------------------------------------------
// Example 14
// -------------------------------------------------------
ring SS = 0, (u,v,z), dp;
poly f =
2*u^7+u^6*v+3*u^5*v^2+u^4*v^3+2*u^3*v^4+u^2*v^5+2*u*v^6+v^7
-7780247/995328*u^6*z-78641/9216*u^5*v*z-10892131/995328*u^4*v^2*z
-329821/31104*u^3*v^3*z-953807/331776*u^2*v^4*z-712429/248832*u*v^5*z
+1537741/331776*v^6*z+2340431/248832*u^5*z^2+5154337/248832*u^4*v*z^2
+658981/41472*u^3*v^2*z^2+1737757/124416*u^2*v^3*z^2
-1234733/248832*u*v^4*z^2-1328329/82944*v^5*z^2-818747/248832*u^4*z^3
-1822879/124416*u^3*v*z^3-415337/31104*u^2*v^2*z^3
+1002655/124416*u*v^3*z^3+849025/82944*v^4*z^3;
// 3 OMPs of mult 3, 1 OMP of mult 4 at origin
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 1
testParametrization(f,RP1);
setring RP1; PARA;
kill SS;kill RP1;
// -------------------------------------------------------
// Example 15
// -------------------------------------------------------
ring SS = 0, (u,v,z), dp;
poly f = 590819418867856650536224u7-147693905508217596067968u6v;
f = f+229117518934972047619978u5v2-174050799674982973889542u4v3;
f = f-92645796479789150855110u3v4-65477418713685583062704u2v5;
f = f+4529961835917468460168uv6+7715404057796585983136v7;
f = f-413640780091141905428104u6z+571836835577486968144618u5vz;
f = f-551807810327826605739444u4v2z-488556410340789283359926u3v3z;
f = f-473466023008413178155962u2v4z+48556741573432247323608uv5z;
f = f+77647371229172269259528v6z+340450118906560552282893u5z2;
f = f-433598825064368371610344u4vz2-937281070591684636591672u3v2z2;
f = f-1388949843915129934647751u2v3z2+204081793110898617103998uv4z2;
f = f+335789953068251652554308v5z2+6485661002496681852577u4z3;
f = f-772700266516318390630202u3vz3-2068348417248100329533330u2v2z3;
f = f+440320154612359641806108uv3z3+808932515589210854581618v4z3;
f = f-229384307132237615286548u3z4-1564303565658228216055227u2vz4;
f = f+520778334468674798322974uv2z4+1172483905704993294097655v3z4;
f = f-480789741398016816562100u2z5+322662751598958620410786uvz5;
f = f+1022525576391791616258310v2z5+82293493608853837667471uz6;
f = f+496839109904761426785889vz6+103766136235628614937587z7; // 15 nodes
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 72
testParametrization(f,RP1);
setring RP1; PARA;
kill SS;kill RP1;

// -------------------------------------------------------
// Example 16
// -------------------------------------------------------
ring SS = 0, (u,v,z), dp;
poly f = 25*u^8+184*u^7*v+518*u^6*v^2+720*u^5*v^3+576*u^4*v^4+282*u^3*v^5;
f = f+84*u^2*v^6+14*u*v^7+v^8+244*u^7*z+1326*u^6*v*z+2646*u^5*v^2*z;
f = f+2706*u^4*v^3*z+1590*u^3*v^4*z+546*u^2*v^5*z+102*u*v^6*z+8*v^7*z;
f = f+854*u^6*z^2+3252*u^5*v*z^2+4770*u^4*v^2*z^2+3582*u^3*v^3*z^2;
f = f+1476*u^2*v^4*z^2+318*u*v^5*z^2+28*v^6*z^2+1338*u^5*z^3+3740*u^4*v*z^3;
f = f+4030*u^3*v^2*z^3+2124*u^2*v^3*z^3+550*u*v^4*z^3+56*v^5*z^3+1101*u^4*z^4;
f = f+2264*u^3*v*z^4+1716*u^2*v^2*z^4+570*u*v^3*z^4+70*v^4*z^4+508*u^3*z^5;
f = f+738*u^2*v*z^5+354*u*v^2*z^5+56*v^3*z^5+132*u^2*z^6+122*u*v*z^6;
f = f+28*v^2*z^6+18*u*z^7+8*v*z^7+z^8; // 3 nodes, 1 sing
adjointIdeal(f,1);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 20
testParametrization(f,RP1);
setring RP1; PARA;
kill SS;kill RP1;
// -------------------------------------------------------
// Example 17
// -------------------------------------------------------
ring SS = 0, (u,v,z), dp;
poly f = -2*u*v^4*z^4+u^4*v^5+12*u^4*v^3*z^2+12*u^2*v^4*z^3-u^3*v*z^5;
f = f+11*u^3*v^2*z^4-21*u^3*v^3*z^3-4*u^4*v*z^4+2*u^4*v^2*z^3-6*u^4*v^4*z;
f = f+u^5*z^4-3*u^5*v^2*z^2+u^5*v^3*z-3*u*v^5*z^3-2*u^2*v^3*z^4+u^3*v^4*z^2;
f = f+v^5*z^4; // 2 OMPs of mult 4, 1 OMP of mult 5, 1 sing at origin
f = subst(f,z,u+v+z);
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 5
testParametrization(f,RP1);
setring RP1; PARA;
kill SS;kill RP1;
// -------------------------------------------------------
// Example 18
// -------------------------------------------------------
ring SS = 0, (u,v,z), dp;
poly f = u^5*v^5+21*u^5*v^4*z-36*u^4*v^5*z-19*u^5*v^3*z^2+12*u^4*v^4*z^2;
f = f+57*u^3*v^5*z^2+u^5*v^2*z^3+u^4*v^3*z^3-53*u^3*v^4*z^3-19*u^2*v^5*z^3;
f = f+u^5*v*z^4+43*u^3*v^3*z^4+u*v^5*z^4+u^5*z^5-15*u^3*v^2*z^5+u^2*v^3*z^5;
f = f+u*v^4*z^5+v^5*z^5; // 1 OMP of mult 4, 3 OMPs of mult 5 (1 at origin)
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 8
testParametrization(f,RP1);
setring RP1; PARA;
kill SS;kill RP1;
// -------------------------------------------------------
// Example 19
// -------------------------------------------------------
ring SS = 0, (u,v,z), dp;
poly f = u^10+6*u^9*v-30*u^7*v^3-15*u^6*v^4+u^5*v^5+u^4*v^6+6*u^3*v^7;
f = f+u^2*v^8+7*u*v^9+v^10+5*u^9*z+24*u^8*v*z-30*u^7*v^2*z-120*u^6*v^3*z;
f = f-43*u^5*v^4*z+5*u^4*v^5*z+20*u^3*v^6*z+10*u^2*v^7*z+29*u*v^8*z+5*v^9*z;
f = f+10*u^8*z^2+36*u^7*v*z^2-105*u^6*v^2*z^2-179*u^5*v^3*z^2-38*u^4*v^4*z^2;
f = f+25*u^3*v^5*z^2+25*u^2*v^6*z^2+46*u*v^7*z^2+10*v^8*z^2+10*u^7*z^3;
f = f+24*u^6*v*z^3-135*u^5*v^2*z^3-117*u^4*v^3*z^3-u^3*v^4*z^3+25*u^2*v^5*z^3;
f = f+34*u*v^6*z^3+10*v^7*z^3+5*u^6*z^4+6*u^5*v*z^4-75*u^4*v^2*z^4;
f = f-27*u^3*v^3*z^4+10*u^2*v^4*z^4+11*u*v^5*z^4+5*v^6*z^4+u^5*z^5;
f = f-15*u^3*v^2*z^5+u^2*v^3*z^5+u*v^4*z^5+v^5*z^5;
// 1 OMP of mult 4, 3 OMPs of mult 5 (1 at origin)
adjointIdeal(f,2);
def RP1 = paraPlaneCurve(f);  // time 2 // see Ex. 18
testParametrization(f,RP1);
setring RP1; PARA;
kill SS;kill RP1;
// -------------------------------------------------------
// Example 20
// -------------------------------------------------------
ring R = 0, (x,y,z), dp;
system("random", 4711);
poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z;
def S = rationalPointConic(p); // quadratic field extension,
                              // minpoly = a^2 - 2
if (testPointConic(p, S) == 1)
{ "point lies on conic"; }
else
{ "point does not lie on conic"; }
kill R;kill S;
// -------------------------------------------------------
// Example 21
// -------------------------------------------------------
ring R = 0, (x,y,z), dp;
system("random", 4711);
poly p = x^2 - 1857669520 * y^2 + 86709575222179747132487270400 * z^2;
def S = rationalPointConic(p); // same as current basering,
                              // no extension needed
if (testPointConic(p, S) == 1)
{ "point lies on conic"; }
else
{ "point does not lie on conic"; }
kill R;kill S;
// -------------------------------------------------------
// Example 21
// -------------------------------------------------------
ring RR = 0, (x,y,z), dp;
poly f = -1965466244509920x5y+34871245546721380061760x4y2;
f = f+104613747941595046117320x3y3+113331564241941002407560x2y4;
f = f+52306876673313609259800xy5+8717812860780028397880y6;
f = f+1040297748510024x5z+4468147845634872x4yz;
f = f-22398508728211453743258x3y2z-33223996581074443306854x2y3z;
f = f-10638598235041298082366xy4z+186886189971594356382y5z;
f = f-1385078844909312x4z2-34893092731637052532683x3yz2;
f = f-98591463214095439056609x2y2z2-92339459334829609336485xy3z2;
f = f-24923289542522905755711y4z2+472440640471377x3z3;
f = f+33821511925664516716011x2yz3+49745237303968344397437xy2z3;
f = f+11040465960074786720475y3z3+8728735735878837099404x2z4;
f = f+17676785754519678518537xyz4+17935885079051421934609y2z4;
f = f-11314701999743172607075xz5-16164284825803158969425yz5;
f = f+3666695988537425618750z6;
// 4 nodes, 1 OMP of mult 4
adjointIdeal(f,2);
kill RR;
*/