source: git/Singular/LIB/brnoeth.lib @ bee06d

//
version="$Id$";
category="Coding theory";
info="
LIBRARY:   brnoeth.lib  Brill-Noether Algorithm, Weierstrass-SG and AG-codes
AUTHORS:   Jose Ignacio Farran Martin, ignfar@eis.uva.es
           Christoph Lossen,           lossen@mathematik.uni-kl.de

OVERVIEW:
 Implementation of the Brill-Noether algorithm for solving the
 Riemann-Roch problem and applications to Algebraic Geometry codes.
 The computation of Weierstrass semigroups is also implemented.@*
 The procedures are intended only for plane (singular) curves defined over
 a prime field of positive characteristic.@*
 For more information about the library see the end of the file brnoeth.lib.

PROCEDURES:
 Adj_div(f)            computes the conductor of a curve
 NSplaces(h,A)         computes non-singular places with given degrees
 BrillNoether(D,C)     computes a vector space basis of the linear system L(D)
 Weierstrass(P,m,C)    computes the Weierstrass semigroup of C at P up to m
 extcurve(d,C)         extends the curve C to an extension of degree d
 AGcode_L(G,D,E)       computes the evaluation AG code with divisors G and D
 AGcode_Omega(G,D,E)   computes the residual AG code with divisors G and D
 prepSV(G,D,F,E)       preprocessing for the basic decoding algorithm
 decodeSV(y,K)         decoding of a word with the basic decoding algorithm
 closed_points(I)      computes the zero-set of a zero-dim. ideal in 2 vars
 dual_code(C)          computes the dual code
 sys_code(C)           computes an equivalent systematic code
 permute_L(L,P)        applies a permutation to a list

KEYWORDS:  Weierstrass semigroup; Algebraic Geometry codes;
           Brill-Noether algorithm
";

LIB "matrix.lib";
LIB "triang.lib";
LIB "hnoether.lib";
LIB "inout.lib";

///////////////////////////////////////////////////////////////////////////////

// **********************************************************
// * POINTS, PLACES AND DIVISORS OF (SINGULAR) PLANE CURVES *
// **********************************************************

proc closed_points (ideal I)
"USAGE:     closed_points(I); I an ideal
RETURN:     list of prime ideals (each a Groebner basis), corresponding to
            the (distinct affine closed) points of V(I)
NOTE:       The ideal must have dimension 0, the basering must have 2
            variables, the ordering must be lp, and the base field must
            be finite and prime.@*
            It might be convenient to set the option(redSB) in advance.
SEE ALSO:   triang_lib
EXAMPLE:    example closed_points; shows an example
"
{
  ideal II=std(I);
  if (II==1)
  {
    return(list());
  }
  list TL=triangMH(II);
  int s=size(TL);
  list L=list();
  int i,j,k;
  ideal Facts;
  poly G2;
  for (i=1;i<=s;i=i+1)
  {
    Facts=factorize(TL[i][1],1);
    k=size(Facts);
    G2=TL[i][2];
    for (j=1;j<=k;j=j+1)
    {
      L=L+pd2(Facts[j],G2);
    }
  }
  // eliminate possible repetitions
  s=size(L);
  list LP=list();
  LP[1]=std(L[1]);
  int counter=1;
  for (i=2;i<=s;i=i+1)
  {
    if (isPinL(L[i],LP)==0)
    {
      counter=counter+1;
      LP[counter]=std(L[i]);
    }
  }
  return(LP);
}
example
{
  "EXAMPLE:";  echo = 2;
  ring s=2,(x,y),lp;
  // this is just the affine plane over F_4 :
  ideal I=x4+x,y4+y;
  list L=closed_points(I);
  // and here you have all the points :
  L;
}
///////////////////////////////////////////////////////////////////////////////

static proc pd2 (poly g1,poly g2)
{
  // If g1,g2 is a std. resp. lex. in (x,y) then the procedure
  //    factorizes g2 in the "extension given by g1"
  //    (then g1 must be irreducible) and returns a list of
  //    ideals with always g1 as first component and the
  //    distinct factors of g2 as second components
  list L=list();
  ideal J=g1;
  int i,s;
  if (deg(g1)==1)
  {
    poly A=-subst(g1,var(2),0);
    poly B=subst(g2,var(2),A);
    ideal facts=factorize(B,1);
    s=size(facts);
    for (i=1;i<=s;i=i+1)
    {
      J[2]=facts[i];
      L[i]=J;
    }
  }
  else
  {
    def BR=basering;
    poly A=g1;
    poly B=g2;
    ring raux1=char(basering),(x,y,@a),lp;
    poly G;
    ring raux2=(char(basering),@a),(x,y),lp;
    map psi=BR,x,@a;
    minpoly=number(psi(A));
    poly f=psi(B);
    ideal facts=factorize(f,1);
    s=size(facts);
    poly g;
    string sg;
    for (i=1;i<=s;i=i+1)
    {
      g=facts[i];
      sg=string(g);
      setring raux1;
      execute("G="+sg+";");
      G=subst(G,@a,y);
      setring BR;
      map ppssii=raux1,var(1),var(2),0;
      J[2]=ppssii(G);
      L[i]=J;
      kill ppssii;
      setring raux2;
    }
    setring BR;
    kill raux1,raux2;
  }
  return(L);
}
///////////////////////////////////////////////////////////////////////////////
static proc isPinL (ideal P,list L)
{
  // checks if a (plane) point P is in a list of (plane) points L
  //    by just comparing generators
  // it works only if all (prime) ideals are given in a "canonical way",
  // namely:
  //    the first generator is monic and irreducible,
  //        and depends only on the second variable,
  //    and the second one is monic in the first variable
  //        and irreducible over the field extension determined by
  //        the second variable and the first generator as minpoly
  int s=size(L);
  int i;
  for (i=1;i<=s;i=i+1)
  {
    if ( P[1]==L[i][1] && P[2]==L[i][2] )
    {
      return(1);
    }
  }
  return(0);
}
///////////////////////////////////////////////////////////////////////////////
static proc s_locus (poly f)
{
  // computes : ideal of affine singular locus
  // the equation f must be affine
  // warning : if there is an error message then the output is "none"
  // option(redSB) is convenient to be set in advance
  ideal I=f,jacob(f);
  I=std(I);
  if (dim(I)>0)
  {
    // dimension check (has to be 0)
    ERROR("something was wrong; possibly non-reduced curve");
  }
  else
  {
    return(I);
  }
}
///////////////////////////////////////////////////////////////////////////////
static proc curve (poly f)
"USAGE:      curve(f), where f is a polynomial (affine or projective)
CREATE:     poly CHI in both rings aff_r=p,(x,y),lp and Proj_R=p,(x,y,z),lp
            also ideal (std) Aff_SLocus of affine singular locus in the ring
            aff_r
RETURN:     list (size 3) with two rings aff_r,Proj_R and an integer deg(f)
NOTE:       f must be absolutely irreducible, but this is not checked
            it is not implemented yet for extensions of prime fields
"
{
  def base_r=basering;
  ring aff_r=char(basering),(x,y),lp;
  ring Proj_R=char(basering),(x,y,z),lp;
  setring base_r;
  int degX=deg(f);
  if (nvars(basering)==2)
  {
    setring aff_r;
    map embpol=base_r,x,y;
    poly CHI=embpol(f);
    export(CHI);
    kill embpol;
    ideal Aff_SLocus=s_locus(CHI);
    export(Aff_SLocus);
    setring Proj_R;
    poly CHI=homog(imap(aff_r,CHI),z);
    export(CHI);
    setring base_r;
    list L=list();
    L[1]=aff_r;
    L[2]=Proj_R;
    L[3]=degX;
    kill aff_r,Proj_R;
    return(L);
  }
  if (nvars(basering)==3)
  {
    setring Proj_R;
    map embpol=base_r,x,y,z;
    poly CHI=embpol(f);
    export(CHI);
    kill embpol;
    string s=string(subst(CHI,z,1));
    setring aff_r;
    execute("poly CHI="+s+";");
    export(CHI);
    ideal Aff_SLocus=s_locus(CHI);
    export(Aff_SLocus);
    setring base_r;
    list L=list();
    L[1]=aff_r;
    L[2]=Proj_R;
    L[3]=degX;
    kill aff_r,Proj_R;
    return(L);
  }
  ERROR("basering must have 2 or 3 variables");
}
///////////////////////////////////////////////////////////////////////////////
static proc Aff_SL (ideal ISL)
{
  // computes : affine singular (closed) points as a list of lists of
  //     prime ideals and intvec (for storing the places over each point)
  // the ideal ISL=s_locus(CHI) is assumed to be computed in advance for
  //     a plane curve CHI, and it must be given by a standard basis
  // for our purpose the function must called with the "global" ideal
  //     "Aff_SLocus"
  list SL=list();
  ideal I=ISL;
  if ( I[1] != 1 )
  {
    list L=list();
    ideal aux;
    intvec iv;
    int i,s;
    L=closed_points(I);
    s=size(L);
    for (i=1;i<=s;i=i+1)
    {
      aux=std(L[i]);
      SL[i]=list(aux,iv);
    }
  }
  return(SL);
}
///////////////////////////////////////////////////////////////////////////////
static proc inf_P (poly f)
{
  // computes : all (closed) points at infinity as homogeneous polynomials
  // output : two lists with respectively singular and non-singular points
  intvec iv;
  def base_r=basering;
  ring r_auxz=char(basering),(x,y,z),lp;
  poly f=imap(base_r,f);
  poly F=homog(f,z);                    // equation of projective curve
  poly f_inf=subst(F,z,0);
  setring base_r;
  poly f_inf=imap(r_auxz,f_inf);
  ideal I=factorize(f_inf,1);           // points at infinity as homogeneous
                                        // polynomials
  int s=size(I);
  int i;
  list IP_S=list();                     // for singular points at infinity
  list IP_NS=list();                    // for non-singular points at infinity
  int counter_S;
  int counter_NS;
  poly aux;
  for (i=1;i<=s;i=i+1)
  {
    aux=subst(I[i],y,1);
    if (aux==1)
    {
      // the point is (1:0:0)
      setring r_auxz;
      poly f_yz=subst(F,x,1);
      if ( subst(subst(diff(f_yz,y),y,0),z,0)==0 &&
           subst(subst(diff(f_yz,z),y,0),z,0)==0 )
      {
        // the point is singular
        counter_S=counter_S+1;
        kill f_yz;
        setring base_r;
        IP_S[counter_S]=list(I[i],iv);
      }
      else
      {
        // the point is non-singular
        counter_NS=counter_NS+1;
        kill f_yz;
        setring base_r;
        IP_NS[counter_NS]=list(I[i],iv);
      }
    }
    else
    {
      // the point is (a:1:0) | a is root of aux
      if (deg(aux)==1)
      {
        // the point is rational and no field extension is needed
        setring r_auxz;
        poly f_xz=subst(F,y,1);
        poly aux=imap(base_r,aux);
        number A=-number(subst(aux,x,0));
        map phi=r_auxz,x+A,0,z;
        poly f_origin=phi(f_xz);
        if ( subst(subst(diff(f_origin,x),x,0),z,0)==0 &&
             subst(subst(diff(f_origin,z),x,0),z,0)==0 )
        {
          // the point is singular
          counter_S=counter_S+1;
          kill f_xz,aux,A,phi,f_origin;
          setring base_r;
          IP_S[counter_S]=list(I[i],iv);
        }
        else
        {
          // the point is non-singular
          counter_NS=counter_NS+1;
          kill f_xz,aux,A,phi,f_origin;
          setring base_r;
          IP_NS[counter_NS]=list(I[i],iv);
        }
      }
      else
      {
        // the point is non-rational and a field extension with minpoly=aux
        // is needed
        ring r_ext=(char(basering),@a),(x,y,z),lp;
        poly F=imap(r_auxz,F);
        poly f_xz=subst(F,y,1);
        poly aux=imap(base_r,aux);
        minpoly=number(subst(aux,x,@a));
        map phi=r_ext,x+@a,0,z;
        poly f_origin=phi(f_xz);
        if ( subst(subst(diff(f_origin,x),x,0),z,0)==0 &&
             subst(subst(diff(f_origin,z),x,0),z,0)==0 )
        {
          // the point is singular
          counter_S=counter_S+1;
          setring base_r;
          kill r_ext;
          IP_S[counter_S]=list(I[i],iv);
        }
        else
        {
          // the point is non-singular
          counter_NS=counter_NS+1;
          setring base_r;
          kill r_ext;
          IP_NS[counter_NS]=list(I[i],iv);
        }
      }
    }
  }
  kill r_auxz;
  return(list(IP_S,IP_NS));
}
///////////////////////////////////////////////////////////////////////////////
static proc closed_points_ext (poly f,int d,ideal SL)
{
  // computes : (closed affine non-singular) points over an extension of
  //            degree d
  // remark(1) : singular points are supposed to be listed appart
  // remark(2) : std SL=s_locus(f) is supposed to be computed in advance
  // remark(3) : ideal SL is used to remove those points which are singular
  // output : list of list of prime ideals with an intvec for storing the
  //          places
  int Q=char(basering)^d;             // cardinality of the extension field
  ideal I=f,x^Q-x,y^Q-y;              // ideal of the searched points
  I=std(I);
  if (I==1)
  {
    return(list());
  }
  list LP=list();
  int m=size(SL);
  list L=list();
  ideal aux;
  intvec iv;
  int i,s,j,counter;
  L=closed_points(I);
  s=size(L);
  for (i=1;i<=s;i=i+1)
  {
    aux=std(L[i]);
    for (j=1;j<=m;j=j+1)
    {
      // check if singular i.e. if SL is contained in aux
      if ( NF(SL[j],aux) != 0 )
      {
        counter=counter+1;
        LP[counter]=list(aux,iv);
        break;
      }
    }
  }
  return(LP);
}
///////////////////////////////////////////////////////////////////////////////
static proc degree_P (list P)
"USAGE:      degree_P(P), where P is either a polynomial or an ideal
RETURN:     integer with the degree of the closed point given by P
SEE ALSO:   closed_points
NOTE:       If P is a (homogeneous irreducible) polynomial the point is at
            infinity, and if P is a (prime) ideal the points is affine, and
            the ideal must be given by 2 generators: the first one irreducible
            and depending only on y, and the second one irreducible over the
            extension given by y with the first generator as minimal polynomial
"
{
  // computes : the degree of a given point
  // remark(1) : if the input is (irreducible homogeneous) polynomial => the point
  //             is at infinity
  // remark(2) : it the input is (std. resp. lp. prime) ideal => the point is
  //             affine
  if (typeof(P[1])=="ideal")
  {
    if (size(P[1])==2)
    {
      int d=deg(P[1][1]);
      poly aux=subst(P[1][2],y,1);
      d=d*deg(aux);
      return(d);
    }
    else
    {
      // this should not happen in principle
      ERROR("non-valid parameter");
    }
  }
  else
  {
    if (typeof(P[1])=="poly")
    {
      return(deg(P[1]));
    }
    else
    {
      ERROR("parameter must have a polynomial or ideal in the first component");
    }
  }
}
///////////////////////////////////////////////////////////////////////////////
static proc closed_points_deg (poly f,int d,ideal SL)
{
  // computes : (closed affine non-singular) points of degree d
  // remark(1) : singular points are supposed to be listed appart
  // remark(2) : std SL=s_locus(f) is supposed to be computed in advance
  list L=closed_points_ext(f,d,SL);
  int s=size(L);
  int i,counter;
  list LP=list();
  for (i=1;i<=s;i=i+1)
  {
    if (degree_P(L[i])==d)
    {
      counter=counter+1;
      LP[counter]=L[i];
    }
  }
  return(LP);
}
///////////////////////////////////////////////////////////////////////////////
static proc idealSubset (ideal I,ideal J)
{
  // checks wether I is contained in J and returns a boolean
  // remark : J is assumed to be given by a standard basis
  int s=ncols(I);
  int i;
  for (i=1;i<=s;i=i+1)
  {
    if ( NF(I[i],std(J)) != 0 )
    {
      return(0);
    }
  }
  return(1);
}
///////////////////////////////////////////////////////////////////////////////
static proc belongs (list P,ideal I)
{
  // checks if affine point P is contained in V(I) and returns a boolean
  // remark : P[1] is assumed to be an ideal given by a standard basis
  if (typeof(P[1])=="ideal")
  {
    return(idealSubset(I,P[1]));
  }
  else
  {
    ERROR("first argument must be an affine point");
  }
}
///////////////////////////////////////////////////////////////////////////////
static proc equals (ideal I,ideal J)
{
  // checks if I is equal to J and returns a boolean
  // remark : I and J are assumed to be given by a standard basis
  int answer=0;
  if (idealSubset(I,J)==1)
  {
    if (idealSubset(J,I)==1)
    {
      answer=1;
    }
  }
  return(answer);
}
///////////////////////////////////////////////////////////////////////////////
static proc isInLP (ideal P,list LP)
{
  // checks if affine point P is a list LP and returns either its position or
  //   zero
  // remark : all points in LP and P itself are assumed to be given by a
  //   standard basis
  // warning : the procedure does not check whether the points are affine or
  //   not
  int s=size(LP);
  if (s==0)
  {
    return(0);
  }
  int i;
  for (i=1;i<=s;i=i+1)
  {
    if (equals(P,LP[i][1])==1)
    {
      return(i);
    }
  }
  return(0);
}
///////////////////////////////////////////////////////////////////////////////
static proc res_deg ()
{
  // computes the residual degree of the basering with respect to its prime
  //   field
  // warning : minpoly must depend on a parameter called "a"
  int ext;
  string s_m=string(minpoly);
  if (s_m=="0")
  {
    ext=1;
  }
  else
  {
    ring auxr=char(basering),a,lp;
    execute("poly minp="+s_m+";");
    ext=deg(minp);
  }
  return(ext);
}
///////////////////////////////////////////////////////////////////////////////
static proc Frobenius (etwas,int r)
{
  // applies the Frobenius map over F_{p^r} to an object defined over an
  //   extension of such field
  // usually it is called with r=1, i.e. the Frobenius map over the prime
  //   field F_p
  // returns always an object of the same type, and works correctly on
  // numbers, polynomials, ideals, matrices or lists of the above types
  // maybe : types vector and module should be added in the future, but they
  //   are not needed now
  int q=char(basering)^r;
  if (typeof(etwas)=="number")
  {
    return(etwas^q);
  }
  if (typeof(etwas)=="poly")
  {
    int s=size(etwas);
    poly f;
    int i;
    for (i=1;i<=s;i=i+1)
    {
      f=f+(leadcoef(etwas[i])^q)*leadmonom(etwas[i]);
    }
    return(f);
  }
  if (typeof(etwas)=="ideal")
  {
    int s=ncols(etwas);
    ideal I;
    int i;
    for (i=1;i<=s;i=i+1)
    {
      I[i]=Frobenius(etwas[i],r);
    }
    return(I);
  }
  if (typeof(etwas)=="matrix")
  {
    int m=nrows(etwas);
    int n=ncols(etwas);
    matrix A[m][n];
    int i,j;
    for (i=1;i<=m;i=i+1)
    {
      for (j=1;j<=n;j=j+1)
      {
        A[i,j]=Frobenius(etwas[i,j],r);
      }
    }
    return(A);
  }
  if (typeof(etwas)=="list")
  {
    int s=size(etwas);
    list L=list();
    int i;
    for (i=1;i<=s;i=i+1)
    {
      if (typeof(etwas[i])<>"none")
      {
        L[i]=Frobenius(etwas[i],r);
      }
    }
    return(L);
  }
  return(etwas);
}
///////////////////////////////////////////////////////////////////////////////
static proc conj_b (list L,int r)
{
  // applies the Frobenius map over F_{p^r} to a list of type HNE defined over
  //   a larger extension
  // when r=1 it turns to be the Frobenius map over the prime field F_{p}
  // returns : a list of type HNE which is either conjugate of the input or
  //     the same list in case of L being actually defined over the base field
  //     F_{p^r}
  int p=char(basering);
  int Q=p^r;
  list LL=list();
  int m=nrows(L[1]);
  int n=ncols(L[1]);
  matrix A[m][n];
  poly f;
  poly aux;
  int i,j;
  for (i=1;i<=m;i=i+1)
  {
    for (j=1;j<=n;j=j+1)
    {
      aux=L[1][i,j];
      if (aux<>x)
      {
        A[i,j]=aux^Q;
      }
      else
      {
        A[i,j]=aux;
        break;
      }
    }
  }
  m=size(L[4]);
  for (i=1;i<=m;i=i+1)
  {
    f=f+(leadcoef(L[4][i])^Q)*leadmonom(L[4][i]);
  }
  LL[1]=A;
  LL[2]=L[2];
  LL[3]=L[3];
  LL[4]=f;
  return(LL);
}
///////////////////////////////////////////////////////////////////////////////
static proc grad_b (list L,int r)
{
  // computes the degree of a list of type HNE which is actually defined over
  //    F_{p^r} eventhough it is given in an extension of such field
  int gr=1;
  int rd=res_deg() div r;
  list LL=L;
  int i;
  for (i=1;i<=rd;i=i+1)
  {
    LL=conj_b(LL,r);
    if ( LL[1]==L[1] && LL[4]==L[4] )
    {
      break;
    }
    else
    {
      gr=gr+1;
    }
  }
  return(gr);
}
///////////////////////////////////////////////////////////////////////////////
static proc conj_bs (list L,int r)
{
  // computes all the conjugates over F_{p^r} of a list of type HNE defined
  //   over an extension
  // returns : a list of lists of type HNE, where the first one is the input
  //   list
  // remark : notice that the degree of the branch is then the size of the
  //   output
  list branches=list();
  int gr=1;
  branches[1]=L;
  int rd=res_deg() div r;
  list LL=L;
  int i;
  for (i=1;i<=rd;i=i+1)
  {
    LL=conj_b(LL,r);
    if ( LL[1]==L[1] && LL[4]==L[4] )
    {
      break;
    }
    else
    {
      gr=gr+1;
      branches[gr]=LL;
    }
  }
  return(branches);
}
///////////////////////////////////////////////////////////////////////////////
static proc subfield (sf)
{
  // writes the generator "a" of a subfield of the coefficients field of
  //   basering in terms of the current generator (also called "a") as a
  //   string sf is an existing ring whose coefficient field is such a subfield
  // warning : in basering there must be a variable called "x" and subfield
  //   must not be prime
  def base_r=basering;
  string new_m=string(minpoly);
  setring sf;
  string old_m=string(minpoly);
  if (old_m==new_m)
  {
    setring base_r;
    return("a");
  }
  else
  {
    if (old_m<>string(0))
    {
      ring auxring=char(basering),(a,x),lp;
      execute("poly mpol="+old_m+";");
      mpol=subst(mpol,a,x);
      setring base_r;
      poly mpol=imap(auxring,mpol);
      kill auxring;
      string answer="? error : non-primitive element";
      int r=res_deg();
      int q=char(basering)^r;
      int i;
      number b;
      for (i=1;i<=q-2;i=i+1)
      {
        b=a^i;
        if (subst(mpol,x,b)==0)
        {
          answer=string(b);
          break;
        }
      }
      if (answer<>"? error : non-primitive element")
      {
        return(answer);
      }
      else
      {
        // list all the elements of the finite field F_q
        int p=char(basering);
        list FF1,FF2;
        for (i=0;i<p;i=i+1)
        {
          FF1[i+1]=number(i);
        }
        int s,j,k;
        for (i=1;i<r;i=i+1)
        {
          s=size(FF1);
          for (j=1;j<=s;j=j+1)
          {
            FF1[j]=FF1[j]*a;
          }
          FF2=FF1;
          for (k=1;k<p;k=k+1)
          {
            for (j=1;j<=s;j=j+1)
            {
              FF1[j]=FF1[j]+number(1);
            }
            FF2=FF2+FF1;
          }
          FF1=FF2;
          FF2=list();
        }
        kill FF2;
        for (i=1;i<=q-1;i=i+1)
        {
          b=FF1[i];
          if (subst(mpol,x,b)==0)
          {
            answer=string(b);
            break;
          }
        }
      }
      kill FF1;
      return(answer);
    }
    else
    {
      dbprint(printlevel+1,"warning : minpoly=0 in the subfield;
                              you should check that nothing is wrong");
      return(string(1));
    }
  }
}
///////////////////////////////////////////////////////////////////////////////
static proc importdatum (sf,string datum,string rel)
{
  // fetchs a poly with name "datum" to the current basering from the ring sf
  //   such that the generator is given by string "rel"
  // warning : ring sf must have only variables (x,y) and basering must have
  //   at least (x,y)
  // warning : the case of minpoly=0 is not regarded; there you can use "imap"
  //   instead
  def base_r=basering;
  if (rel=="a")
  {
    setring sf;
    execute("poly pdatum="+datum+";");
    setring base_r;
    poly pdatum=imap(sf,pdatum);
    return(pdatum);
  }
  else
  {
    setring sf;
    execute("string sdatum=string("+datum+");");
    ring auxring=char(basering),(a,x,y),lp;
    execute("poly pdatum="+sdatum+";");
    execute("map phi=basering,"+rel+",x,y;");
    pdatum=phi(pdatum);
    string snewdatum=string(pdatum);
    setring base_r;
    execute("poly pdatum="+snewdatum+";");
    kill auxring;
    return(pdatum);
  }
}
///////////////////////////////////////////////////////////////////////////////
static proc rationalize (lf,string datum,string rel)
{
  // fetchs a poly with name "datum" to the current basering from the ring lf
  //   and larger coefficients field, where the generator of current ring is
  //   given by string "rel" and "datum" is actually defined over the small
  //   field
  // warning : ring lf must have only variables (x,y) and basering must have
  //           at least (x,y)
  // warning : the case of minpoly=0 is supposed unnecessary, since then
  //           "datum" should be already written in the right way, i.e. in
  //           terms of the prime field
  def base_r=basering;
  if (rel=="a")
  {
    setring lf;
    execute("poly pdatum="+datum+";");
    setring base_r;
    poly pdatum=imap(lf,pdatum);
    return(pdatum);
  }
  else
  {
    setring lf;
    execute("string sdatum=string("+datum+");");
    ring auxring=char(basering),(a,b,x,y,t),lp;
//  execute("poly pdatum="+datum+";");
    execute("poly pdatum="+sdatum+";");
    execute("poly prel=b-("+rel+");");
// was before:
//    ideal I=pdatum,prel;
//    I=eliminate(I,a);
//    if (size(I)<=1){ poly newdatum=I[1]; }
//    else { ERROR("Elimination error"); }
// CL 09/05 new:
    poly newdatum=reduce(pdatum,std(prel));
newdatum;
//
    newdatum=subst(newdatum,b,a);
    string snewdatum=string(newdatum);
    setring base_r;
    execute("poly newdatum="+snewdatum+";");
    kill auxring;
    return(newdatum);
  }
}
///////////////////////////////////////////////////////////////////////////////
static proc place (intvec Pp,int sing,list CURVE)
{
  // computes the "rational" places which are defined over a (closed) point
  // Pp points to an appropriate point of the given curve
  // creates : local rings (if they do not exist yet) and then add the
  //   places to a list
  //   each place is given basically by the coordinates of the point and a
  //   list HNdevelop
  // returns : list with all updated data of the curve
  // if the places already exist they are not computed again
  // if sing==1 the point is assumed singular and computes the local conductor
  //   for all places using the local invariants of the branches
  // if sing==2 the point is assumed singular and computes the local conductor
  //   for all places using the Dedekind formula and local parametrizations
  //   of the branches
  // if sing<>1&2 the point is assumed non-singular and the local conductor
  //   should be zero
  list PP=list();
  if (Pp[1]==0)
  {
    if (Pp[2]==0)
    {
      PP=Aff_SPoints[Pp[3]];
    }
    if (Pp[2]==1)
    {
      PP=Inf_Points[1][Pp[3]];
    }
    if (Pp[2]==2)
    {
      PP=Inf_Points[2][Pp[3]];
    }
  }
  else
  {
    PP=Aff_Points(Pp[2])[Pp[3]];
  }
  if (PP[2]<>0)
  {
    return(CURVE);
  }
  intvec PtoPl;
  def base_r=basering;
  int ext1;
  list Places=CURVE[3];
  intvec Conductor=CURVE[4];
  list update_CURVE=CURVE;
  if (typeof(PP[1])=="ideal")
  {
    ideal P=PP[1];
    if (size(P)==2)
    {
      int d=deg(P[1]);
      poly aux=subst(P[2],y,1);
      d=d*deg(aux);
      ext1=d;
      // the point is (A:B:1) but one must distinguish several cases
      // P is assumed to be a std. resp. "(x,y),lp" and thus P[1] depends
      // only on "y"
      if (d==1)
      {
        // the point is rational
        number B=-number(subst(P[1],y,0));
        poly aux2=subst(P[2],y,B);
        number A=-number(subst(aux2,x,0));
        // the point is (A:B:1)
        ring local_aux=char(basering),(x,y),ls;
        number coord@1=imap(base_r,A);
        number coord@2=imap(base_r,B);
        number coord@3=number(1);
        map phi=base_r,x+coord@1,y+coord@2;
        poly CHI=phi(CHI);
      }
      else
      {
        if (deg(P[1])==1)
        {
          // the point is non-rational but the second component needs no
          // field extension
          number B=-number(subst(P[1],y,0));
          poly aux2=subst(P[2],y,B);
          // the point has degree d>1
          // careful : the parameter will be called "a" anyway
          ring local_aux=(char(basering),a),(x,y),ls;
          map psi=base_r,a,0;
          minpoly=number(psi(aux2));
          number coord@1=a;
          number coord@2=imap(base_r,B);
          number coord@3=number(1);
          // the point is (a:B:1)
          map phi=base_r,x+a,y+coord@2;
          poly CHI=phi(CHI);
        }
        else
        {
          if (deg(subst(P[2],y,1))==1)
          {
            // the point is non-rational but the needed minpoly is just P[1]
            // careful : the parameter will be called "a" anyway
            poly P1=P[1];
            poly P2=P[2];
            ring local_aux=(char(basering),a),(x,y),ls;
            map psi=base_r,0,a;
            minpoly=number(psi(P1));
            // the point looks like (A:a:1)
            // A is computed by substituting y=a in P[2]
            poly aux1=imap(base_r,P2);
            poly aux2=subst(aux1,y,a);
            number coord@1=-number(subst(aux2,x,0));
            number coord@2=a;
            number coord@3=number(1);
            map phi=base_r,x+coord@1,y+a;
            poly CHI=phi(CHI);
          }
          else
          {
            // this is the most complicated case of non-rational point
            // firstly : construct an extension of degree d and guess the
            //           minpoly
            poly P1=P[1];
            poly P2=P[2];
            int p=char(basering);
            int Q=p^d;
            ring aux_r=(Q,a),(x,y,t),ls;
            string minpoly_string=string(minpoly);
            ring local_aux=(char(basering),a),(x,y),ls;
            execute("minpoly="+minpoly_string+";");
            // secondly : compute one root of P[1]
            poly P_1=imap(base_r,P1);
            poly P_2=imap(base_r,P2);
            ideal factors1=factorize(P_1,1);       // hopefully this works !!!!
            number coord@2=-number(subst(factors1[1],y,0));
            // thirdly : compute one of the first components for the above root
            poly P_0=subst(P_2,y,coord@2);
            ideal factors2=factorize(P_0,1);       // hopefully this works !!!!
            number coord@1=-number(subst(factors2[1],x,0));
            number coord@3=number(1);
            map phi=base_r,x+coord@1,y+coord@2;
            poly CHI=phi(CHI);
            kill aux_r;
          }
        }
      }
    }
    else
    {
      // this should not happen in principle
      ERROR("non-valid parameter");
    }
  }
  else
  {
    if (typeof(PP[1])=="poly")
    {
      poly P=PP[1];
      ring r_auxz=char(basering),(x,y,z),lp;
      poly CHI=imap(base_r,CHI);
      CHI=homog(CHI,z);
      setring base_r;
      poly aux=subst(P,y,1);
      if (aux==1)
      {
        // the point is (1:0:0)
        ring local_aux=char(basering),(x,y),ls;
        number coord@1=number(1);
        number coord@2=number(0);
        number coord@3=number(0);
        map Phi=r_auxz,1,x,y;
        poly CHI=Phi(CHI);
        ext1=1;
      }
      else
      {
        // the point is (A:1:0) where A is a root of aux
        int d=deg(aux);
        ext1=d;
        if (d==1)
        {
          // the point is rational
          number A=-number(subst(aux,x,0));
          ring local_aux=char(basering),(x,y),ls;
          number coord@1=imap(base_r,A);
          number coord@2=number(1);
          number coord@3=number(0);
          map Phi=r_auxz,x+coord@1,1,y;
          poly CHI=Phi(CHI);
        }
        else
        {
          // the point has degree d>1
          // careful : the parameter will be called "a" anyway
          ring local_aux=(char(basering),a),(x,y),ls;
          map psi=base_r,a,1;
          minpoly=number(psi(P));
          number coord@1=a;
          number coord@2=number(1);
          number coord@3=number(0);
          map Phi=r_auxz,x+a,1,y;
          poly CHI=Phi(CHI);
        }
      }
      kill r_auxz;
    }
    else
    {
      ERROR("a point must have a polynomial or ideal in the first component");
    }
  }
  export(coord@1);
  export(coord@2);
  export(coord@3);
  export(CHI);
  int i,j,k;
  int m,n;
  list LLL=ratdevelop(CHI);
//  list LLL=hnexpansion(CHI,"ess");
  if (typeof(LLL[1])=="ring") {
    def altring=basering;
    def HNEring = LLL[1]; setring HNEring;
    def L@HNE = HND;
    kill HND;
//    def L@HNE = hne;
//    kill hne;

  }
  else {
    def L@HNE=LLL;
//    def L@HNE=LLL[1];
    def HNEring=basering;
    setring HNEring;
  }
  kill LLL;
  export(L@HNE);
  int n_branches=size(L@HNE);
  list Li_aux=list();
  int N_branches;
  int N=size(Places);
  if (sing==1)
  {
    list delta2=list();
    for (i=1;i<=n_branches;i=i+1)
    {
      delta2[i]=invariants(L@HNE[i])[5];
    }
    int dq;
  }
  int ext2=res_deg();
  list dgs=list();
  int ext_0;
  int check;
  string sss,olda,newa;
  if (defined(Q)==0)
  {
    int Q;
  }
  if (ext1==1)
  {
    if (ext2==1)
    {
      if (sing==1)
      {
        intmat I_mult[n_branches][n_branches];
        if (n_branches>1)
        {
          for (i=1;i<=n_branches-1;i=i+1)
          {
            for (j=i+1;j<=n_branches;j=j+1)
            {
              I_mult[i,j]=intersection(L@HNE[i],L@HNE[j]);
              I_mult[j,i]=I_mult[i,j];
            }
          }
        }
      }
////////
check=0;
////////
      if (size(update_CURVE[5])>0)
      {
        if (typeof(update_CURVE[5][1])=="list")
        {
          check=1;
        }
      }
      if (check==0)
      {
        intvec dgs_points(1);
        ring S(1)=char(basering),(x,y,t),ls;
        list BRANCHES=list();
        list POINTS=list();
        list LOC_EQS=list();
        list PARAMETRIZATIONS=list();
        export(BRANCHES);
        export(POINTS);
        export(LOC_EQS);
        export(PARAMETRIZATIONS);
      }
      else
      {
        intvec dgs_points(1)=update_CURVE[5][1][2];
        def S1=update_CURVE[5][1][1];
        execute("ring S(1)="+string(update_CURVE[5][1][1])+";");
        fetchall(S1);
        kill S1;
      }
      N_branches=size(BRANCHES);
      for (i=1;i<=n_branches;i=i+1)
      {
        dgs_points(1)[N_branches+i]=1;
        POINTS[N_branches+i]=list();
        POINTS[N_branches+i][1]=imap(local_aux,coord@1);
        POINTS[N_branches+i][2]=imap(local_aux,coord@2);
        POINTS[N_branches+i][3]=imap(local_aux,coord@3);
        LOC_EQS[N_branches+i]=imap(local_aux,CHI);
        setring HNEring;
        Li_aux=L@HNE[i];
        setring S(1);
        BRANCHES=insert(BRANCHES,imap(HNEring,Li_aux),N_branches+i-1);
        PARAMETRIZATIONS[N_branches+i]=param(BRANCHES[N_branches+i],0);
        N=N+1;
        intvec iw=1,N_branches+i;
        Places[N]=iw;
        if (sing==1)
        {
          dq=delta2[i];
          for (j=1;j<=n_branches;j=j+1)
          {
            dq=dq+I_mult[i,j];
          }
          Conductor[N]=dq;
        }
        if (sing==2)
        {
          Conductor[N]=local_conductor(iw[2],S(1));
        }
        kill iw;
        PtoPl[i]=N;
      }
      setring base_r;
      update_CURVE[5][1]=list();
      update_CURVE[5][1][1]=S(1);
      update_CURVE[5][1][2]=dgs_points(1);
    }
    else
    {
      // we start with a rational point but we get non-rational branches
      // they may have different degrees and then we may need to reduce the
      // field extensions for each one, and finally check if the minpoly
      // fetches with S(i) or not
      // if one of the branches is rational, we may trust that it is written
      // correctly
      if (sing==1)
      {
        int n_geobrs;
        int counter_c;
        list auxgb=list();
        list geobrs=list();
        for (i=1;i<=n_branches;i=i+1)
        {
          auxgb=conj_bs(L@HNE[i],1);
          dgs[i]=size(auxgb);
          n_geobrs=n_geobrs+dgs[i];
          for (j=1;j<=dgs[i];j=j+1)
          {
            counter_c=counter_c+1;
            geobrs[counter_c]=auxgb[j];
          }
        }
        intmat I_mult[n_geobrs][n_geobrs];
        for (i=1;i<n_geobrs;i=i+1)
        {
          for (j=i+1;j<=n_geobrs;j=j+1)
          {
            I_mult[i,j]=intersection(geobrs[i],geobrs[j]);
            I_mult[j,i]=I_mult[i,j];
          }
        }
        kill auxgb,geobrs;
      }
      else
      {
        for (i=1;i<=n_branches;i=i+1)
        {
          dgs[i]=grad_b(L[i],1);
        }
      }
      // the actual degree of each branch is computed and now check if the
      // local ring exists
      for (i=1;i<=n_branches;i=i+1)
      {
        ext_0=dgs[i];
////////
check=0;
////////
        if (size(update_CURVE[5])>=ext_0)
        {
          if (typeof(update_CURVE[5][ext_0])=="list")
          {
            check=1;
          }
        }
        if (check==0)
        {
          if (ext_0>1)
          {
            if (ext_0==ext2)
            {
              sss=string(minpoly);
            }
            else
            {
              Q=char(basering)^ext_0;
              ring auxxx=(Q,a),z,lp;
              sss=string(minpoly);
              setring base_r;
              kill auxxx;
            }
            ring S(ext_0)=(char(basering),a),(x,y,t),ls;
            execute("minpoly="+sss+";");
          }
          else
          {
            ring S(ext_0)=char(basering),(x,y,t),ls;
          }
          intvec dgs_points(ext_0);
          list BRANCHES=list();
          list POINTS=list();
          list LOC_EQS=list();
          list PARAMETRIZATIONS=list();
          export(BRANCHES);
          export(POINTS);
          export(LOC_EQS);
          export(PARAMETRIZATIONS);
        }
        else
        {
          intvec dgs_points(ext_0)=update_CURVE[5][ext_0][2];
          def Sext_0=update_CURVE[5][ext_0][1];
          setring Sext_0;
          string SM=string(minpoly);
          string SR=string(update_CURVE[5][ext_0][1]);
          execute("ring S("+string(ext_0)+")="+SR+";");
          execute("minpoly="+SM+";");
          kill SM,SR;
          fetchall(Sext_0);
          kill Sext_0;
        }
        N_branches=size(BRANCHES);
        dgs_points(ext_0)[N_branches+1]=1;
        POINTS[N_branches+1]=list();
        POINTS[N_branches+1][1]=imap(local_aux,coord@1);
        POINTS[N_branches+1][2]=imap(local_aux,coord@2);
        POINTS[N_branches+1][3]=imap(local_aux,coord@3);
        LOC_EQS[N_branches+1]=imap(local_aux,CHI);
        // now fetch the branches into the new local ring
        if (ext_0==1)
        {
          setring HNEring;
          Li_aux=L@HNE[i];
          setring S(1);
          BRANCHES=insert(BRANCHES,imap(HNEring,Li_aux),N_branches);
        }
        else
        {
          // rationalize branch
          setring HNEring;
          newa=subfield(S(ext_0));
          m=nrows(L@HNE[i][1]);
          n=ncols(L@HNE[i][1]);
          setring S(ext_0);
          list Laux=list();
          poly paux=rationalize(HNEring,"L@HNE["+string(i)+"][4]",newa);
          matrix Maux[m][n];
          for (j=1;j<=m;j=j+1)
          {
            for (k=1;k<=n;k=k+1)
            {
              Maux[j,k]=rationalize(HNEring,"L@HNE["+string(i)+"][1]["+
                           string(j)+","+string(k)+"]",newa);
            }
          }
          setring HNEring;
          intvec Li2=L@HNE[i][2];
          int Li3=L@HNE[i][3];
          setring S(ext_0);
          Laux[1]=Maux;
          Laux[2]=Li2;
          Laux[3]=Li3;
          Laux[4]=paux;
          BRANCHES=insert(BRANCHES,Laux,N_branches);
          kill Laux,Maux,paux,Li2,Li3;
        }
        PARAMETRIZATIONS[N_branches+1]=param(BRANCHES[N_branches+1],0);
        N=N+1;
        intvec iw=ext_0,N_branches+1;
        Places[N]=iw;
        if (sing==2)
        {
          Conductor[N]=local_conductor(iw[2],S(ext_0));
        }
        kill iw;
        PtoPl[i]=N;
        setring HNEring;
        update_CURVE[5][ext_0]=list();
        update_CURVE[5][ext_0][1]=S(ext_0);
        update_CURVE[5][ext_0][2]=dgs_points(ext_0);
      }
      if (sing==1)
      {
        int N_ini=N-n_branches;
        counter_c=1;
        for (i=1;i<=n_branches;i=i+1)
        {
          dq=delta2[i];
          for (j=1;j<=n_geobrs;j=j+1)
          {
            dq=dq+I_mult[counter_c,j];
          }
          Conductor[N_ini+i]=dq;
          counter_c=counter_c+dgs[i];
        }
      }
      setring base_r;
    }
  }
  else
  {
    if (ext1==ext2)
    {
      // the degree of the point equals to the degree of all branches
      // one must just fetch the minpoly's of local_aux, HNEring and S(ext2)
      if (sing==1)
      {
        intmat I_mult[n_branches][n_branches];
        if (n_branches>1)
        {
          for (i=1;i<=n_branches-1;i=i+1)
          {
            for (j=i+1;j<=n_branches;j=j+1)
            {
              I_mult[i,j]=intersection(L@HNE[i],L@HNE[j]);
              I_mult[j,i]=I_mult[i,j];
            }
          }
        }
      }
////////
check=0;
////////
      if (size(update_CURVE[5])>=ext2)
      {
        if (typeof(update_CURVE[5][ext2])=="list")
        {
          check=1;
        }
      }
      if (check==0)
      {
        sss=string(minpoly);
        ring S(ext2)=(char(basering),a),(x,y,t),ls;
        execute("minpoly="+sss+";");
        intvec dgs_points(ext2);
        list BRANCHES=list();
        list POINTS=list();
        list LOC_EQS=list();
        list PARAMETRIZATIONS=list();
        export(BRANCHES);
        export(POINTS);
        export(LOC_EQS);
        export(PARAMETRIZATIONS);
      }
      else
      {
        intvec dgs_points(ext2)=update_CURVE[5][ext2][2];
        def Sext2=update_CURVE[5][ext2][1];
        setring Sext2;
        string SM=string(minpoly);
        string SR=string(update_CURVE[5][ext2][1]);
        execute("ring S("+string(ext2)+")="+SR+";");
        execute("minpoly="+SM+";");
        kill SM,SR;
        fetchall(Sext2);
        kill Sext2;
      }
      N_branches=size(BRANCHES);
      for (i=1;i<=n_branches;i=i+1)
      {
        // fetch all the data into the new local ring
        olda=subfield(local_aux);
        dgs_points(ext2)[N_branches+i]=ext1;
        POINTS[N_branches+i]=list();
        POINTS[N_branches+i][1]=number(importdatum(local_aux,"coord@1",olda));
        POINTS[N_branches+i][2]=number(importdatum(local_aux,"coord@2",olda));
        POINTS[N_branches+i][3]=number(importdatum(local_aux,"coord@3",olda));
        LOC_EQS[N_branches+i]=importdatum(local_aux,"CHI",olda);
        newa=subfield(HNEring);
        setring HNEring;
        m=nrows(L@HNE[i][1]);
        n=ncols(L@HNE[i][1]);
        setring S(ext2);
        list Laux=list();
        poly paux=importdatum(HNEring,"L@HNE["+string(i)+"][4]",newa);
        matrix Maux[m][n];
        for (j=1;j<=m;j=j+1)
        {
          for (k=1;k<=n;k=k+1)
          {
            Maux[j,k]=importdatum(HNEring,"L@HNE["+string(i)+"][1]["+
                                   string(j)+","+string(k)+"]",newa);
          }
        }
        setring HNEring;
        intvec Li2=L@HNE[i][2];
        int Li3=L@HNE[i][3];
        setring S(ext2);
        Laux[1]=Maux;
        Laux[2]=Li2;
        Laux[3]=Li3;
        Laux[4]=paux;
        BRANCHES=insert(BRANCHES,Laux,N_branches+i-1);
        kill Laux,Maux,paux,Li2,Li3;
        PARAMETRIZATIONS[N_branches+i]=param(BRANCHES[N_branches+i],0);
        N=N+1;
        intvec iw=ext2,N_branches+i;
        Places[N]=iw;
        if (sing==1)
        {
          dq=delta2[i];
          for (j=1;j<=n_branches;j=j+1)
          {
            dq=dq+I_mult[i,j];
          }
          Conductor[N]=dq;
        }
        if (sing==2)
        {
          Conductor[N]=local_conductor(iw[2],S(ext2));
        }
        kill iw;
        PtoPl[i]=N;
      }
      setring base_r;
      update_CURVE[5][ext2]=list();
      update_CURVE[5][ext2][1]=S(ext2);
      update_CURVE[5][ext2][2]=dgs_points(ext2);
    }
    else
    {
      // this is the most complicated case
      if (sing==1)
      {
        int n_geobrs;
        int counter_c;
        list auxgb=list();
        list geobrs=list();
        for (i=1;i<=n_branches;i=i+1)
        {
          auxgb=conj_bs(L@HNE[i],ext1);
          dgs[i]=size(auxgb);
          n_geobrs=n_geobrs+dgs[i];
          for (j=1;j<=dgs[i];j=j+1)
          {
            counter_c=counter_c+1;
            geobrs[counter_c]=auxgb[j];
          }
        }
        intmat I_mult[n_geobrs][n_geobrs];
        for (i=1;i<n_geobrs;i=i+1)
        {
          for (j=i+1;j<=n_geobrs;j=j+1)
          {
            I_mult[i,j]=intersection(geobrs[i],geobrs[j]);
            I_mult[j,i]=I_mult[i,j];
          }
        }
        kill auxgb,geobrs;
      }
      else
      {
        for (i=1;i<=n_branches;i=i+1)
        {
          dgs[i]=grad_b(L@HNE[i],ext1);
        }
      }
      for (i=1;i<=n_branches;i=i+1)
      {
        // first compute the actual degree of each branch and check if the
        // local ring exists
        ext_0=ext1*dgs[i];
////////
check=0;
////////
        if (size(update_CURVE[5])>=ext_0)
        {
          if (typeof(update_CURVE[5][ext_0])=="list")
          {
            check=1;
          }
        }
        if (check==0)
        {
          if (ext_0>ext1)
          {
            if (ext_0==ext2)
            {
              sss=string(minpoly);
            }
            else
            {
              Q=char(basering)^ext_0;
              ring auxxx=(Q,a),z,lp;
              sss=string(minpoly);
              setring base_r;
              kill auxxx;
            }
          }
          else
          {
            setring local_aux;
            sss=string(minpoly);
          }
          ring S(ext_0)=(char(basering),a),(x,y,t),ls;
          execute("minpoly="+sss+";");
          intvec dgs_points(ext_0);
          list BRANCHES=list();
          list POINTS=list();
          list LOC_EQS=list();
          list PARAMETRIZATIONS=list();
          export(BRANCHES);
          export(POINTS);
          export(LOC_EQS);
          export(PARAMETRIZATIONS);
        }
        else
        {
          intvec dgs_points(ext_0)=update_CURVE[5][ext_0][2];
          def Sext_0=update_CURVE[5][ext_0][1];
          setring Sext_0;
          string SM=string(minpoly);
          string SR=string(update_CURVE[5][ext_0][1]);
          execute("ring S("+string(ext_0)+")="+SR+";");
          execute("minpoly="+SM+";");
          kill SM,SR;
          fetchall(badring);
          kill badring;
        }
        N_branches=size(BRANCHES);
        // now fetch all the data into the new local ring
        olda=subfield(local_aux);
        dgs_points(ext_0)[N_branches+1]=ext1;
        POINTS[N_branches+1]=list();
        POINTS[N_branches+1][1]=number(importdatum(local_aux,"coord@1",olda));
        POINTS[N_branches+1][2]=number(importdatum(local_aux,"coord@2",olda));
        POINTS[N_branches+1][3]=number(importdatum(local_aux,"coord@3",olda));
        LOC_EQS[N_branches+1]=importdatum(local_aux,"CHI",olda);
        setring HNEring;
        newa=subfield(S(ext_0));
        m=nrows(L@HNE[i][1]);
        n=ncols(L@HNE[i][1]);
        setring S(ext_0);
        list Laux=list();
        poly paux=rationalize(HNEring,"L@HNE["+string(i)+"][4]",newa);
        matrix Maux[m][n];
        for (j=1;j<=m;j=j+1)
        {
          for (k=1;k<=n;k=k+1)
          {
            Maux[j,k]=rationalize(HNEring,"L@HNE["+string(i)+"][1]["+
               string(j)+","+string(k)+"]",newa);
          }
        }
        setring HNEring;
        intvec Li2=L@HNE[i][2];
        int Li3=L@HNE[i][3];
        setring S(ext_0);
        Laux[1]=Maux;
        Laux[2]=Li2;
        Laux[3]=Li3;
        Laux[4]=paux;
        BRANCHES=insert(BRANCHES,Laux,N_branches);
        kill Laux,Maux,paux,Li2,Li3;
        PARAMETRIZATIONS[N_branches+1]=param(BRANCHES[N_branches+1],0);
        N=N+1;
        intvec iw=ext_0,N_branches+1;
        Places[N]=iw;
        if (sing==2)
        {
          Conductor[N]=local_conductor(iw[2],S(ext_0));
        }
        kill iw;
        PtoPl[i]=N;
        setring HNEring;
        update_CURVE[5][ext_0]=list();
        update_CURVE[5][ext_0][1]=S(ext_0);
        update_CURVE[5][ext_0][2]=dgs_points(ext_0);
      }
      if (sing==1)
      {
        int N_ini=N-n_branches;
        counter_c=1;
        for (i=1;i<=n_branches;i=i+1)
        {
          dq=delta2[i];
          for (j=1;j<=n_geobrs;j=j+1)
          {
            dq=dq+I_mult[counter_c,j];
          }
          Conductor[N_ini+i]=dq;
          counter_c=counter_c+dgs[i];
        }
      }
      setring base_r;
    }
  }
  update_CURVE[3]=Places;
  update_CURVE[4]=Conductor;
  PP[2]=PtoPl;
  if (Pp[1]==0)
  {
    if (Pp[2]==0)
    {
      Aff_SPoints[Pp[3]]=PP;
    }
    if (Pp[2]==1)
    {
      Inf_Points[1][Pp[3]]=PP;
    }
    if (Pp[2]==2)
    {
      Inf_Points[2][Pp[3]]=PP;
    }
  }
  else
  {
    Aff_Points(Pp[2])[Pp[3]]=PP;
  }
  update_CURVE[1][1]=base_r;
  kill HNEring;
  kill local_aux;
  return(update_CURVE);
}
///////////////////////////////////////////////////////////////////////////////
static proc local_conductor (int k,SS)
{
  // computes the degree of the local conductor at a place of a plane curve
  // if the point is non-singular the result will be zero
  // the computation is carried out with the "Dedekind formula" via
  // parametrizations
  int a,b,Cq;
  def b_ring=basering;
  setring SS;
  poly fx=diff(LOC_EQS[k],x);
  poly fy=diff(LOC_EQS[k],y);
  int nr=ncols(BRANCHES[k][1]);
  poly xt=PARAMETRIZATIONS[k][1][1];
  poly yt=PARAMETRIZATIONS[k][1][2];
  int ordx=PARAMETRIZATIONS[k][2][1];
  int ordy=PARAMETRIZATIONS[k][2][2];
  map phi_t=basering,xt,yt,1;
  poly derf;
  if (fx<>0)
  {
    derf=fx;
    poly tt=diff(yt,t);
    b=mindeg(tt);
    if (ordy>-1)
    {
      while (b>=ordy)
      {
        BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
        nr=ncols(BRANCHES[k][1]);
        PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
        ordy=PARAMETRIZATIONS[k][2][2];
        yt=PARAMETRIZATIONS[k][1][2];
        tt=diff(yt,t);
        b=mindeg(tt);
      }
      xt=PARAMETRIZATIONS[k][1][1];
      ordx=PARAMETRIZATIONS[k][2][1];
    }
    poly ft=phi_t(derf);
  }
  else
  {
    derf=fy;
    poly tt=diff(xt,t);
    b=mindeg(tt);
    if (ordx>-1)
    {
      while (b>=ordx)
      {
        BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
        nr=ncols(BRANCHES[k][1]);
        PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
        ordx=PARAMETRIZATIONS[k][2][1];
        xt=PARAMETRIZATIONS[k][1][1];
        tt=diff(xt,t);
        b=mindeg(tt);
      }
      yt=PARAMETRIZATIONS[k][1][2];
      ordy=PARAMETRIZATIONS[k][2][2];
    }
    poly ft=phi_t(derf);
  }
  a=mindeg(ft);
  if ( ordx>-1 || ordy>-1 )
  {
    if (ordy==-1)
    {
      while (a>ordx)
      {
        BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
        nr=ncols(BRANCHES[k][1]);
        PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
        ordx=PARAMETRIZATIONS[k][2][1];
        xt=PARAMETRIZATIONS[k][1][1];
        ft=phi_t(derf);
        a=mindeg(ft);
      }
    }
    else
    {
      if (ordx==-1)
      {
        while (a>ordy)
        {
          BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
          nr=ncols(BRANCHES[k][1]);
          PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
          ordy=PARAMETRIZATIONS[k][2][2];
          yt=PARAMETRIZATIONS[k][1][2];
          ft=phi_t(derf);
          a=mindeg(ft);
        }
      }
      else
      {
        int ordf=ordx;
        if (ordx>ordy)
        {
          ordf=ordy;
        }
        while (a>ordf)
        {
          BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
          nr=ncols(BRANCHES[k][1]);
          PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
          ordx=PARAMETRIZATIONS[k][2][1];
          ordy=PARAMETRIZATIONS[k][2][2];
          ordf=ordx;
          if (ordx>ordy)
          {
            ordf=ordy;
          }
          xt=PARAMETRIZATIONS[k][1][1];
          yt=PARAMETRIZATIONS[k][1][2];
          ft=phi_t(derf);
          a=mindeg(ft);
        }
      }
    }
  }
  Cq=a-b;
  setring b_ring;
  return(Cq);
}
///////////////////////////////////////////////////////////////////////////////
static proc max_D (intvec D1,intvec D2)
{
  // computes the maximum of two divisors (intvec)
  int s1=size(D1);
  int s2=size(D2);
  int i;
  if (s1>s2)
  {
    for (i=1;i<=s2;i=i+1)
    {
      if (D2[i]>D1[i])
      {
        D1[i]=D2[i];
      }
    }
    for (i=s2+1;i<=s1;i=i+1)
    {
      if (D1[i]<0)
      {
        D1[i]=0;
      }
    }
    return(D1);
  }
  else
  {
    for (i=1;i<=s1;i=i+1)
    {
      if (D1[i]>D2[i])
      {
        D2[i]=D1[i];
      }
    }
    for (i=s1+1;i<=s2;i=i+1)
    {
      if (D2[i]<0)
      {
        D2[i]=0;
      }
    }
    return(D2);
  }
}
///////////////////////////////////////////////////////////////////////////////
static proc deg_D (intvec D,list PP)
{
  // computes the degree of a divisor (intvec or list of integers)
  int i;
  int d=0;
  int s=size(D);
  for (i=1;i<=s;i=i+1)
  {
    d=d+D[i]*PP[i][1];
  }
  return(d);
}
///////////////////////////////////////////////////////////////////////////////

// ============================================================================
// *******  MAIN PROCEDURES for the "preprocessing" of Brill-Noether   ********
// ============================================================================

proc Adj_div (poly f,list #)
"USAGE:    Adj_div( f [,l] ); f a poly, [l a list]
RETURN:   list L with the computed data:
  @format
  L[1] a list of rings: L[1][1]=aff_r (affine), L[1][2]=Proj_R (projective),
  L[2] an intvec with 2 entries (degree, genus),
  L[3] a list of intvec (closed places),
  L[4] an intvec (conductor),
  L[5] a list of lists:
     L[5][d][1] a (local) ring over an extension of degree d,
     L[5][d][2] an intvec (degrees of base points of places of degree d)
  @end format
NOTE:     @code{Adj_div(f);} computes and stores the fundamental data of the
          plane curve defined by f as needed for AG codes.
          In the affine ring you can find the following data:
   @format
   poly CHI:  affine equation of the curve,
   ideal Aff_SLocus:  affine singular locus (std),
   list Inf_Points:  points at infinity
            Inf_Points[1]:  singular points
            Inf_Points[2]:  non-singular points,
   list Aff_SPoints:  affine singular points (if not empty).
   @end format
          In the projective ring you can find the projective equation
          CHI of the curve (poly).
          In the local rings L[5][d][1] you find:
   @format
   list POINTS:  base points of the places of degree d,
   list LOC_EQS:  local equations of the curve at the base points,
   list BRANCHES:  Hamburger-Noether developments of the places,
   list PARAMETRIZATIONS:  local parametrizations of the places,
   @end format
          Each entry of the list L[3] corresponds to one closed place (i.e.,
          a place and all its conjugates) which is represented by an intvec
          of size two, the first entry is the degree of the place (in
          particular, it tells the local ring where to find the data
          describing one representative of the closed place), and the
          second one is the position of those data in the lists POINTS, etc.,
          inside this local ring.@*
          In the intvec L[4] (conductor) the i-th entry corresponds to the
          i-th entry in the list of places L[3].

          With no optional arguments, the conductor is computed by
          local invariants of the singularities; otherwise it is computed
          by the Dedekind formula. @*
          An affine point is represented by a list P where P[1] is std
          of a prime ideal and P[2] is an intvec containing the position
          of the places above P in the list of closed places L[3]. @*
          If the point is at infinity, P[1] is a homogeneous irreducible
          polynomial in two variables.

          If @code{printlevel>=0} additional comments are displayed (default:
          @code{printlevel=0}).
WARNING:  The parameter of the needed field extensions is called 'a'. Thus,
          there should be no global object named 'a' when executing Adj_div.
KEYWORDS: Hamburger-Noether expansions; adjunction divisor
SEE ALSO: closed_points, NSplaces
EXAMPLE:  example Adj_div; shows an example
"
{
  // computes the adjunction divisor and the genus of a (singular) plane curve
  // as a byproduct, it computes all the singular points with the corresponding
  //      places and the genus of the curve
  // the adjunction divisor is stored in an intvec
  // also the non-singular places at infinity are computed
  // returns a list with all the computed data
  if (char(basering)==0)
  {
    ERROR("Base field not implemented");
  }
  if (npars(basering)>0)
  {
    ERROR("Base field not implemented");
  }
  if (defined(a)==1)
  {
    if ((typeof(a)=="int") or (typeof(a)=="intvec") or (typeof(a)=="intmat") or
        (typeof(a)=="string") or (typeof(a)=="proc") or (typeof(a)=="list"))
    {
      print("WARNING: There is an object named 'a' of the global type "
             + typeof(a)+".");
      print("This will cause problems when "
             + "executing Adj_div, as 'a' is used for "
             + "the name of the parameter of the field extensions needed.");
      ERROR("Please rename or kill the object named 'a'");
    }
  }
  intvec opgt=option(get);
  option(redSB);
  def Base_R=basering;
  list CURVE=curve(f);
  def aff_r=CURVE[1];
  def Proj_R=CURVE[2];
  int degX=CURVE[3];
  int genusX=(degX-1)*(degX-2);
  genusX = genusX div 2;
  intvec iivv=degX,genusX;
  intvec Conductor;
  setring aff_r;
  dbprint(printlevel+1,"Computing affine singular points ... ");
  list Aff_SPoints=Aff_SL(Aff_SLocus);
  int s=size(Aff_SPoints);
  if (s>0)
  {
    export(Aff_SPoints);
  }
  dbprint(printlevel+1,"Computing all points at infinity ... ");
  list Inf_Points=inf_P(CHI);
  export(Inf_Points);
  list update_CURVE=list();
  update_CURVE[1]=list();
  update_CURVE[1][1]=aff_r;
  update_CURVE[1][2]=Proj_R;
  update_CURVE[2]=iivv;
  update_CURVE[3]=list();
  update_CURVE[4]=Conductor;
  update_CURVE[5]=list();
  int i;
  intvec pP=0,0,0;
  if (size(#)==0)
  {
    dbprint(printlevel+1,"Computing affine singular places ... ");
    if (s>0)
    {
      for (i=1;i<=s;i=i+1)
      {
        pP[3]=i;
        update_CURVE=place(pP,1,update_CURVE);
      }
    }
    dbprint(printlevel+1,"Computing singular places at infinity ... ");
    s=size(Inf_Points[1]);
    if (s>0)
    {
      pP[2]=1;
      for (i=1;i<=s;i=i+1)
      {
        pP[3]=i;
        update_CURVE=place(pP,1,update_CURVE);
      }
    }
  }
  else
  {
    dbprint(printlevel+1,"Computing affine singular places ... ");
    s=size(Aff_SPoints);
    if (s>0)
    {
      for (i=1;i<=s;i=i+1)
      {
        pP[3]=i;
        update_CURVE=place(pP,2,update_CURVE);
      }
    }
    dbprint(printlevel+1,"Computing singular places at infinity ... ");
    s=size(Inf_Points[1]);
    if (s>0)
    {
      pP[2]=1;
      for (i=1;i<=s;i=i+1)
      {
        pP[3]=i;
        update_CURVE=place(pP,2,update_CURVE);
      }
    }
  }
  dbprint(printlevel+1,"Computing non-singular places at infinity ... ");
  s=size(Inf_Points[2]);
  if (s>0)
  {
    pP[2]=2;
    for (i=1;i<=s;i=i+1)
    {
      pP[3]=i;
      update_CURVE=place(pP,0,update_CURVE);
    }
  }
  dbprint(printlevel+1,"Adjunction divisor computed successfully");
  list Places=update_CURVE[3];
  Conductor=update_CURVE[4];
  genusX = genusX - (deg_D(Conductor,Places) div 2);
  update_CURVE[2][2]=genusX;
  setring Base_R;
  dbprint(printlevel+1," ");
  dbprint(printlevel+2,"The genus of the curve is "+string(genusX));
  option(set,opgt);
  return(update_CURVE);
}
example
{
  "EXAMPLE:";  echo = 2;
  int plevel=printlevel;
  printlevel=-1;
  ring s=2,(x,y),lp;
  list C=Adj_div(y9+y8+xy6+x2y3+y2+x3);
  def aff_R=C[1][1];      // the affine ring
  setring aff_R;
  listvar(aff_R);         // data in the affine ring
  CHI;                    // affine equation of the curve
  Aff_SLocus;             // ideal of the affine singular locus
  Aff_SPoints[1];         // 1st affine singular point: (1:1:1), no.1
  Inf_Points[1];          // singular point(s) at infinity: (1:0:0), no.4
  Inf_Points[2];          // list of non-singular points at infinity
  //
  pause("press RETURN");
  def proj_R=C[1][2];     // the projective ring
  setring proj_R;
  CHI;                    // projective equation of the curve
  C[2][1];                // degree of the curve
  C[2][2];                // genus of the curve
  C[3];                   // list of computed places
  C[4];                   // adjunction divisor (all points are singular!)
  //
  pause("press RETURN");
  // we look at the place(s) of degree 2 by changing to the ring
  C[5][2][1];
  def S(2)=C[5][2][1];
  setring S(2);
  POINTS;                // base point(s) of place(s) of degree 2: (1:a:1)
  LOC_EQS;               // local equation(s)
  PARAMETRIZATIONS;      // parametrization(s) and exactness
  BRANCHES;              // Hamburger-Noether development
  printlevel=plevel;
}
///////////////////////////////////////////////////////////////////////////////

proc NSplaces (intvec h,list CURVE)
"USAGE:   NSplaces( h, CURVE ), where h is an intvec and CURVE is a list
RETURN:   list L with updated data of CURVE after computing all non-singular
          affine closed places whose degrees are in the intvec h: @*
   @format
   in L[1][1]: (affine ring) lists Aff_Points(d) with affine non-singular
               (closed) points of degree d (if non-empty),
   in L[3]:    the newly computed closed places are added,
   in L[5]:    local rings created/updated to store (repres. of) new places.
   @end format
          See @ref{Adj_div} for a description of the entries in L.
NOTE:     The list_expression should be the output of the procedure Adj_div.@*
          Raising @code{printlevel}, additional comments are displayed
          (default: @code{printlevel=0}).
WARNING:  The parameter of the needed field extensions is called 'a'. Thus,
          there should be no global object named 'a' when executing NSplaces.
SEE ALSO: closed_points, Adj_div
EXAMPLE:  example NSplaces; shows an example
"
{
  // computes all the non-singular affine (closed) places with fixed degrees
  // creates lists of points and the corresponding places
  // list CURVE must be the output of the procedure "Adj_div"
  // notice : if h==0 then nothing is done
  // warning : non-positive entries are forgotten
  if (defined(a)==1)
  {
    if ((typeof(a)=="int") or (typeof(a)=="intvec") or (typeof(a)=="intmat") or
        (typeof(a)=="string") or (typeof(a)=="proc") or (typeof(a)=="list"))
    {
      print("WARNING: There is an object named 'a' of the global type "
             + typeof(a)+".");
      print("This will cause problems when "
             + "executing Adj_div, as 'a' is used for "
             + "the name of the parameter of the field extensions needed.");
      ERROR("Please rename or kill the object named 'a'");
    }
  }
  if (h==0)
  {
    return(CURVE);
  }
  intvec opgt=option(get);
  option(redSB);
  def Base_R=basering;
  def aff_r=CURVE[1][1];
  int M=size(h);
  list update_CURVE=CURVE;
  int i,j,l,s;
  setring aff_r;
  intvec pP=1,0,0;
  for (i=1;i<=M;i=i+1)
  {
    l=h[i];
    if (l>0)
    {
      dbprint(printlevel+1,"Computing non-singular affine places of degree "
                             +string(l)+" ... ");
      if (defined(Aff_Points(l))==0)
      {
        list Aff_Points(l)=closed_points_deg(CHI,l,Aff_SLocus);
        s=size(Aff_Points(l));
        if (s>0)
        {
          export(Aff_Points(l));
          pP[2]=l;
          for (j=1;j<=s;j=j+1)
          {
            dbprint(printlevel-voice-2,"Place No.\ "+string(j));
            pP[3]=j;
            update_CURVE=place(pP,0,update_CURVE);
          }
        }
      }
    }
  }
  setring Base_R;
  option(set,opgt);
  return(update_CURVE);
}
example
{
  "EXAMPLE:";  echo = 2;
  int plevel=printlevel;
  printlevel=-1;
  ring s=2,(x,y),lp;
  list C=Adj_div(x3y+y3+x);
  // The list of computed places:
  C[3];
  // create places up to degree 4
  list L=NSplaces(1..4,C);
  // The list of computed places is now:
  L[3];
  // e.g., affine non-singular points of degree 4 :
  def aff_r=L[1][1];
  setring aff_r;
  Aff_Points(4);
  // e.g., base point of the 1st place of degree 4 :
  def S(4)=L[5][4][1];
  setring S(4);
  POINTS[1];
  printlevel=plevel;
}
///////////////////////////////////////////////////////////////////////////////

// ** SPECIAL PROCEDURES FOR LINEAR ALGEBRA **

static proc Ker (matrix A)
{
  // warning : "lp" ordering is necessary
  intvec opgt=option(get);
  option(redSB);
  matrix M=transpose(syz(A));
  option(set,opgt);
  return(M);
}
///////////////////////////////////////////////////////////////////////////////
static proc get_NZsol (matrix A)
{
  matrix sol=Ker(A);
  return(submat(sol,1..1,1..ncols(sol)));
}
///////////////////////////////////////////////////////////////////////////////
static proc supplement (matrix W,matrix V)
"USAGE:     supplement(W,V), where W,V are matrices of numbers such that the
            vector space generated by the rows of W is contained in that
            generated by the rows of V
RETURN:     matrix whose rows generate a supplementary vector space of W in V,
            or a zero row-matrix if <W>==<V>
NOTE:       W,V must be given with maximal rank
"
{
  // W and V represent independent sets of vectors and <W> is assumed to be
  // contained in <V>
  // computes matrix S whose rows are l.i. vectors s.t. <W> union <S> is a
  // basis of <V>
  // careful : the size of all vectors is assumed to be the same but it is
  //  not checked and neither the linear independence of the vectors is checked
  // the trivial case W=0 is not covered by this procedure (and thus V<>0)
  // if <W>=<V> then a zero row-matrix is returned
  // warning : option(redSB) must be set in advance
  int n1=nrows(W);
  int n2=nrows(V);
  int s=n2-n1;
  if (s==0)
  {
    int n=ncols(W);
    matrix HH[1][n];
    return(HH);
  }
  matrix H=transpose(lift(transpose(V),transpose(W)));
  H=supplem(H);
  return(H*V);
}
///////////////////////////////////////////////////////////////////////////////
static proc supplem (matrix M)
"USAGE:      supplem(M), where M is a matrix of numbers with maximal rank
RETURN:     matrix whose rows generate a supplementary vector space of <M> in
            k^n, where k is the base field and n is the number of columns
SEE ALSO:   supplement
NOTE:       The rank r is assumed to be 1<r<n.
"
{
  // warning : the linear independence of the rows is not checked
  int r=nrows(M);
  int n=ncols(M);
  int s=n-r;
  matrix A=M;
  matrix supl[s][n];
  int counter=0;
  int h=r+1;
  int i;
  for (i=1;i<=n;i=i+1)
  {
    matrix TT[1][n];
    TT[1,i]=1;
    A=transpose(concat(transpose(A),transpose(TT)));
    r=mat_rank(A);
    if (r==h)
    {
      h=h+1;
      counter=counter+1;
      supl=transpose(concat(transpose(supl),transpose(TT)));
      if (counter==s)
      {
        break;
      }
    }
    kill TT;
  }
  supl=transpose(compress(transpose(supl)));
  return(supl);
}
///////////////////////////////////////////////////////////////////////////////
static proc mat_rank (matrix A)
{
  // warning : "lp" ordering is necessary
  intvec opgt=option(get);
  option(redSB);
  int r=size(std(module(transpose(A))));
  option(set,opgt);
  return(r);
}
///////////////////////////////////////////////////////////////////////////////

// ***************************************************************
// * PROCEDURES FOR INTERPOLATION, INTERSECTION AND EXTRA PLACES *
// ***************************************************************

static proc estim_n (intvec Dplus,int dgX,list PL)
{
  // computes an estimate for the degree n in the Brill-Noether algorithm
  int estim=2*deg_D(Dplus,PL)+dgX*(dgX-3);
  estim=estim div (2*dgX);
  estim=estim+1;
  if (estim<dgX)
  {
    estim=dgX;
  }
  return(estim);
}
///////////////////////////////////////////////////////////////////////////////
static proc nforms (int n)
{
  // computes the list of all homogeneous monomials of degree n>=0
  // exports ideal nFORMS(n) whose generators are ranged with lp order
  //   in Proj_R and returns size(nFORMS(n))
  // warning : it is supposed to be called inside Proj_R
  // if n<=0 then nFORMS(0) is "computed fast"
  ideal nFORMS(n);
  int N;
  if (n>0)
  {
    N=(n+1)*(n+2);
    N=N div 2;
    N=N+1;
    int i,j,k;
    for (i=0;i<=n;i=i+1)
    {
      for (j=0;j<=n-i;j=j+1)
      {
        k=k+1;
        nFORMS(n)[N-k]=x^i*y^j*z^(n-i-j);
      }
    }
    export(nFORMS(n));
  }
  else
  {
    N=2;
    nFORMS(0)=1;
    export(nFORMS(0));
  }
  return(N-1);
}
///////////////////////////////////////////////////////////////////////////////
static proc nmultiples (int n,int dgX,poly f)
{
  // computes a basis of the space of forms of degree n which are multiple of
  //     CHI
  // returns a matrix whose rows are the coordinates (related to nFORMS(n))
  //     of such a basis
  // warning : it is supposed to be called inside Proj_R
  // warning : nFORMS(n) is created in the way, together with nFORMS(n-degX)
  // warning : n must be greater or equal than the degree of the curve
  if (defined(nFORMS(n))==0)
  {
    dbprint(printlevel+1,string(nforms(n)));
  }
  int m=n-dgX;
  if (defined(nFORMS(m))==0)
  {
    int k=nforms(m);
  }
  else
  {
    int k=size(nFORMS(m));
  }
  ideal nmults;
  int i;
  for (i=1;i<=k;i=i+1)
  {
    nmults[i]=f*nFORMS(m)[i];
  }
  return(transpose(lift(nFORMS(n),nmults)));
}
///////////////////////////////////////////////////////////////////////////////
static proc interpolating_forms (intvec D,int n,list CURVE)
{
  // computes a vector basis of the space of forms of degree n whose
  //     intersection divisor with the curve is greater or equal than D>=0
  // the procedure is supposed to be called inside the ring Proj_R and
  //     assumes that the forms nFORMS(n) are previously computed
  // the output is a matrix whose rows are the coordinates in nFORMS(n) of
  //     such a basis
  // remark : the support of D may contain "extra" places
  def BR=basering;
  def aff_r=CURVE[1][1];
  int N=size(nFORMS(n));
  matrix totalM[1][N];
  int s=size(D);
  list Places=CURVE[3];
  int NPls=size(Places);
  int i,j,k,kk,id,ip,RR,ordx,ordy,nr,NR;
  if (s<=NPls)
  {
    for (i=1;i<=s;i=i+1)
    {
      if (D[i]>0)
      {
        id=Places[i][1];
        ip=Places[i][2];
        RR=D[i];
        def SS=CURVE[5][id][1];
        setring SS;
        poly xt=PARAMETRIZATIONS[ip][1][1];
        poly yt=PARAMETRIZATIONS[ip][1][2];
        ordx=PARAMETRIZATIONS[ip][2][1];
        ordy=PARAMETRIZATIONS[ip][2][2];
        nr=ncols(BRANCHES[ip][1]);
        if ( ordx>-1 || ordy>-1 )
        {
          while ( ( RR>=ordx && ordx>-1 ) || ( RR>=ordy && ordy>-1 ) )
          {
            BRANCHES[ip]=extdevelop(BRANCHES[ip],2*nr);
            nr=ncols(BRANCHES[ip][1]);
            PARAMETRIZATIONS[ip]=param(BRANCHES[ip],0);
            xt=PARAMETRIZATIONS[ip][1][1];
            yt=PARAMETRIZATIONS[ip][1][2];
            ordx=PARAMETRIZATIONS[ip][2][1];
            ordy=PARAMETRIZATIONS[ip][2][2];
          }
        }
        if (POINTS[ip][3]==number(1))
        {
          number A=POINTS[ip][1];
          number B=POINTS[ip][2];
          map Mt=BR,A+xt,B+yt,1;
          kill A,B;
        }
        else
        {
          if (POINTS[ip][2]==number(1))
          {
            number A=POINTS[ip][1];
            map Mt=BR,A+xt,1,yt;
            kill A;
          }
          else
          {
            map Mt=BR,1,xt,yt;
          }
        }
        ideal nFORMS(n)=Mt(nFORMS(n));
        // rewrite properly the above conditions to obtain the local equations
        matrix partM[RR][N];
        matrix auxMC=coeffs(nFORMS(n),t);
        NR=nrows(auxMC);
        if (RR<=NR)
        {
          for (j=1;j<=RR;j=j+1)
          {
            for (k=1;k<=N;k=k+1)
            {
              partM[j,k]=number(auxMC[j,k]);
            }
          }
        }
        else
        {
          for (j=1;j<=NR;j=j+1)
          {
            for (k=1;k<=N;k=k+1)
            {
              partM[j,k]=number(auxMC[j,k]);
            }
          }
          for (j=NR+1;j<=RR;j=j+1)
          {
            for (k=1;k<=N;k=k+1)
            {
              partM[j,k]=number(0);
            }
          }
        }
        matrix localM=partM;
        matrix conjM=partM;
        for (j=2;j<=id;j=j+1)
        {
          conjM=Frobenius(conjM,1);
          localM=transpose(concat(transpose(localM),transpose(conjM)));
        }
        localM=gauss_row(localM);
        setring BR;
        totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM))));
        totalM=transpose(compress(transpose(totalM)));
        setring SS;
        kill xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM;
        setring BR;
        kill SS;
      }
    }
  }
  else
  {
    // distinguish between "standard" places and "extra" places
    for (i=1;i<=NPls;i=i+1)
    {
      if (D[i]>0)
      {
        id=Places[i][1];
        ip=Places[i][2];
        RR=D[i];
        def SS=CURVE[5][id][1];
        setring SS;
        poly xt=PARAMETRIZATIONS[ip][1][1];
        poly yt=PARAMETRIZATIONS[ip][1][2];
        ordx=PARAMETRIZATIONS[ip][2][1];
        ordy=PARAMETRIZATIONS[ip][2][2];
        nr=ncols(BRANCHES[ip][1]);
        if ( ordx>-1 || ordy>-1 )
        {
          while ( ( RR>=ordx && ordx>-1 ) || ( RR>=ordy && ordy>-1 ) )
          {
            BRANCHES[ip]=extdevelop(BRANCHES[ip],2*nr);
            nr=ncols(BRANCHES[ip][1]);
            PARAMETRIZATIONS[ip]=param(BRANCHES[ip],0);
            xt=PARAMETRIZATIONS[ip][1][1];
            yt=PARAMETRIZATIONS[ip][1][2];
            ordx=PARAMETRIZATIONS[ip][2][1];
            ordy=PARAMETRIZATIONS[ip][2][2];
          }
        }
        if (POINTS[ip][3]==number(1))
        {
          number A=POINTS[ip][1];
          number B=POINTS[ip][2];
          map Mt=BR,A+xt,B+yt,1;
          kill A,B;
        }
        else
        {
          if (POINTS[ip][2]==number(1))
          {
            number A=POINTS[ip][1];
            map Mt=BR,A+xt,1,yt;
            kill A;
          }
          else
          {
            map Mt=BR,1,xt,yt;
          }
        }
        ideal nFORMS(n)=Mt(nFORMS(n));
        // rewrite properly the above conditions to obtain the local equations
        matrix partM[RR][N];
        matrix auxMC=coeffs(nFORMS(n),t);
        NR=nrows(auxMC);
        if (RR<=NR)
        {
          for (j=1;j<=RR;j=j+1)
          {
            for (k=1;k<=N;k=k+1)
            {
              partM[j,k]=number(auxMC[j,k]);
            }
          }
        }
        else
        {
          for (j=1;j<=NR;j=j+1)
          {
            for (k=1;k<=N;k=k+1)
            {
              partM[j,k]=number(auxMC[j,k]);
            }
          }
          for (j=NR+1;j<=RR;j=j+1)
          {
            for (k=1;k<=N;k=k+1)
            {
              partM[j,k]=number(0);
            }
          }
        }
        matrix localM=partM;
        matrix conjM=partM;
        for (j=2;j<=id;j=j+1)
        {
          conjM=Frobenius(conjM,1);
          localM=transpose(concat(transpose(localM),transpose(conjM)));
        }
        localM=gauss_row(localM);
        setring BR;
        totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM))));
        totalM=transpose(compress(transpose(totalM)));
        setring SS;
        kill xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM;
        setring BR;
        kill SS;
      }
    }
    k=s-NPls;
    int l;
    for (i=1;i<=k;i=i+1)
    {
      // in this case D[NPls+i]>0 is assumed to be true during the
      // Brill-Noether algorithm
      RR=D[NPls+i];
      setring aff_r;
      def SS=@EXTRA@[2][i];
      def extra_dgs=@EXTRA@[3];
      setring SS;
      poly xt=PARAMETRIZATION[1][1];
      poly yt=PARAMETRIZATION[1][2];
      ordx=PARAMETRIZATION[2][1];
      ordy=PARAMETRIZATION[2][2];
      nr=ncols(BRANCH[1]);
      if ( ordx>-1 || ordy>-1 )
      {
        while ( ( RR>=ordx && ordx>-1 ) || ( RR>=ordy && ordy>-1 ) )
        {
          BRANCH=extdevelop(BRANCH,2*nr);
          nr=ncols(BRANCH[1]);
          PARAMETRIZATION=param(BRANCH,0);
          xt=PARAMETRIZATION[1][1];
          yt=PARAMETRIZATION[1][2];
          ordx=PARAMETRIZATION[2][1];
          ordy=PARAMETRIZATION[2][2];
        }
      }
      number A=POINT[1];
      number B=POINT[2];
      map Mt=BR,A+xt,B+yt,1;
      kill A,B;
      ideal nFORMS(n)=Mt(nFORMS(n));
      // rewrite properly the above conditions to obtain the local equations
      matrix partM[RR][N];
      matrix auxMC=coeffs(nFORMS(n),t);
      NR=nrows(auxMC);
      if (RR<=NR)
      {
        for (j=1;j<=RR;j=j+1)
        {
          for (kk=1;kk<=N;kk=kk+1)
          {
            partM[j,kk]=number(auxMC[j,kk]);
          }
        }
      }
      else
      {
        for (j=1;j<=NR;j=j+1)
        {
          for (kk=1;kk<=N;kk=kk+1)
          {
            partM[j,kk]=number(auxMC[j,kk]);
          }
        }
        for (j=NR+1;j<=RR;j=j+1)
        {
          for (kk=1;kk<=N;kk=kk+1)
          {
            partM[j,kk]=number(0);
          }
        }
      }
      matrix localM=partM;
      matrix conjM=partM;
      l=extra_dgs[i];
      for (j=2;j<=l;j=j+1)
      {
        conjM=Frobenius(conjM,1);
        localM=transpose(concat(transpose(localM),transpose(conjM)));
      }
      localM=gauss_row(localM);
      setring BR;
      totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM))));
      totalM=transpose(compress(transpose(totalM)));
      setring SS;
      kill xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM;
      setring BR;
      kill SS;
    }
  }
  return(Ker(totalM));
}
///////////////////////////////////////////////////////////////////////////////
static proc local_IN (poly h,int m)
{
  // computes the intersection number of h and the curve CHI at a certain place
  // returns a list with the intersection number and the "leading coefficient"
  // the procedure must be called inside a local ring, h must be a local
  //     equation with respect to the desired place, and m indicates the
  //     number of place inside that local ring, containing lists
  //     POINT(S)/BRANCH(ES)/PARAMETRIZATION(S) when m=0 an "extra place" is
  //     considered
  def BR=basering;
  if (m>0)
  {
    int nr=ncols(BRANCHES[m][1]);
    poly xt=PARAMETRIZATIONS[m][1][1];
    poly yt=PARAMETRIZATIONS[m][1][2];
    int ordx=PARAMETRIZATIONS[m][2][1];
    int ordy=PARAMETRIZATIONS[m][2][2];
    map phi=BR,xt,yt,1;
    poly ht=phi(h);
    int inum=mindeg(ht);
    if ( ordx>-1 || ordy>-1 )
    {
      while (((inum>ordx||inum==-1) && ordx>-1) || ((inum>ordy||inum==-1) && ordy>-1))
      {
        BRANCHES[m]=extdevelop(BRANCHES[m],2*nr);
        nr=ncols(BRANCHES[m][1]);
        PARAMETRIZATIONS[m]=param(BRANCHES[m],0);
        xt=PARAMETRIZATIONS[m][1][1];
        yt=PARAMETRIZATIONS[m][1][2];
        ordx=PARAMETRIZATIONS[m][2][1];
        ordy=PARAMETRIZATIONS[m][2][2];
        phi=BR,xt,yt,1;
        ht=phi(h);
        inum=mindeg(ht);
      }
    }
  }
  else
  {
    int nr=ncols(BRANCH[1]);
    poly xt=PARAMETRIZATION[1][1];
    poly yt=PARAMETRIZATION[1][2];
    int ordx=PARAMETRIZATION[2][1];
    int ordy=PARAMETRIZATION[2][2];
    map phi=basering,xt,yt,1;
    poly ht=phi(h);
    int inum=mindeg(ht);
    if ( ordx>-1 || ordy>-1 )
    {
      while (((inum>ordx||inum==-1) && ordx>-1) || ((inum>ordy||inum==-1) && ordy>-1))
      {
        BRANCH=extdevelop(BRANCH,2*nr);
        nr=ncols(BRANCH[1]);
        PARAMETRIZATION=param(BRANCH,0);
        xt=PARAMETRIZATION[1][1];
        yt=PARAMETRIZATION[1][2];
        ordx=PARAMETRIZATION[2][1];
        ordy=PARAMETRIZATION[2][2];
        phi=BR,xt,yt,1;
        ht=phi(h);
        inum=mindeg(ht);
      }
    }
  }
  list answer=list();
  answer[1]=inum;
  number AA=number(coeffs(ht,t)[inum+1,1]);
  answer[2]=AA;
  return(answer);
}
///////////////////////////////////////////////////////////////////////////////
static proc extra_place (ideal P)
{
  // computes the "rational" place which is defined over a (closed) "extra"
  //     point
  // an "extra" point will be necessarily affine, non-singular and non-rational
  // creates : a specific local ring to deal with the (unique) place above it
  // returns : list with the above local ring and the degree of the point/place
  // warning : the procedure must be called inside the affine ring aff_r
  def base_r=basering;
  int ext=deg(P[1]);
  poly aux=subst(P[2],y,1);
  ext=ext*deg(aux);
  // P is assumed to be a std. resp. "(x,y),lp" and thus P[1] depends only
  // on "y"
  if (deg(P[1])==1)
  {
    // the point is non-rational but the second component needs no field
    // extension
    number B=-number(subst(P[1],y,0));
    poly aux2=subst(P[2],y,B);
    // careful : the parameter will be called "a" anyway
    ring ES=(char(basering),a),(x,y,t),ls;
    map psi=base_r,a,0;
    minpoly=number(psi(aux2));
    kill psi;
    number A=a;
    number B=imap(base_r,B);
  }
  else
  {
    if (deg(subst(P[2],y,1))==1)
    {
      // the point is non-rational but the needed minpoly is just P[1]
      // careful : the parameter will be called "a" anyway
      poly P1=P[1];
      poly P2=P[2];
      ring ES=(char(basering),a),(x,y,t),ls;
      map psi=base_r,0,a;
      minpoly=number(psi(P1));
      kill psi;
      poly aux1=imap(base_r,P2);
      poly aux2=subst(aux1,y,a);
      number A=-number(subst(aux2,x,0));
      number B=a;
    }
    else
    {
      // this is the most complicated case of non-rational point
      poly P1=P[1];
      poly P2=P[2];
      int p=char(basering);
      int Q=p^ext;
      ring aux_r=(Q,a),(x,y,t),ls;
      string minpoly_string=string(minpoly);
      ring ES=(char(basering),a),(x,y,t),ls;
      execute("minpoly="+minpoly_string+";");
      poly P_1=imap(base_r,P1);
      poly P_2=imap(base_r,P2);
      ideal factors1=factorize(P_1,1);
      number B=-number(subst(factors1[1],y,0));
      poly P_0=subst(P_2,y,B);
      ideal factors2=factorize(P_0,1);
      number A=-number(subst(factors2[1],x,0));
      kill aux_r;
    }
  }
  list POINT=list();
  POINT[1]=A;
  POINT[2]=B;
  export(POINT);
  map phi=base_r,x+A,y+B;
  poly LOC_EQ=phi(CHI);
  kill A,B,phi;
  list LLL=hnexpansion(LOC_EQ,"ess");
  if (typeof(LLL[1])=="ring") {
    def altring=basering;
    def HNEring = LLL[1];
    setring HNEring;
    def L@HNE = hne[1];
    kill hne;
    export(L@HNE);
    int m=nrows(L@HNE[1]);
    int n=ncols(L@HNE[1]);
    intvec Li2=L@HNE[2];
    int Li3=L@HNE[3];
    setring ES;
    string newa=subfield(HNEring);
    poly paux=importdatum(HNEring,"L@HNE[4]",newa);
    matrix Maux[m][n];
    int i,j;
    for (i=1;i<=m;i=i+1)
    {
      for (j=1;j<=n;j=j+1)
      {
        Maux[i,j]=importdatum(HNEring,"L@HNE[1]["+string(i)+","+
                                                  string(j)+"]",newa);
      }
    }
  }
  else {
    def L@HNE=LLL[1][1];
    matrix Maux=L@HNE[1];
    intvec Li2=L@HNE[2];
    int Li3=L@HNE[3];
    poly paux=L@HNE[4];
    def HNEring=basering;
    setring HNEring;
  }
  kill LLL;

  list BRANCH=list();
  BRANCH[1]=Maux;
  BRANCH[2]=Li2;
  BRANCH[3]=Li3;
  BRANCH[4]=paux;
  export(BRANCH);
  list PARAMETRIZATION=param(BRANCH,0);
  export(PARAMETRIZATION);
  kill HNEring;
  setring base_r;
  list answer=list();
  answer[1]=ES;
  answer[2]=ext;
  kill ES;
  return(answer);
}
///////////////////////////////////////////////////////////////////////////////
static proc intersection_div (poly H,list CURVE)
"USAGE:     intersection_div(H,CURVE), where H is a homogeneous polynomial
            in ring Proj_R=p,(x,y,z),lp and CURVE is the list of data for
            the given curve
CREATE:     new places which had not been computed in advance if necessary
            those places are stored each one in a local ring where you find
            lists POINT,BRANCH,PARAMETRIZATION for the place living in that
            ring; the degree of the point/place in such a ring is stored in
            an intvec, and the base points in the remaining list
            Everything is exported in a list @EXTRA@ inside the ring
            aff_r=CURVE[1][1] and returned with the updated CURVE
RETURN:     list with the intersection divisor (intvec) between the underlying
            curve and H=0, and the list CURVE updated
SEE ALSO:   Adj_div, NSplaces, closed_points
NOTE:       The procedure must be called from the ring Proj_R=CURVE[1][2]
            (projective).
            If @EXTRA@ already exists, the new places are added to the
            previous data.
"
{
  // computes the intersection divisor of H and the curve CHI
  // returns a list with (possibly) "extra places" and it must be called
  //     inside Proj_R
  // in case of extra places, some local rings ES(1) ... ES(m) are created
  //     together with an integer list "extra_dgs" containing the degrees of
  //     those places
  intvec opgt=option(get);
  option(redSB);
  intvec interdiv;
  def BRing=basering;
  int Tz1=deg(H);
  list Places=CURVE[3];
  int N=size(Places);
  def aff_r=CURVE[1][1];
  setring aff_r;
  if (defined(@EXTRA@)==0)
  {
    list @EXTRA@=list();
    list EP=list();
    list ES=list();
    list extra_dgs=list();
  }
  else
  {
    list EP=@EXTRA@[1];
    list ES=@EXTRA@[2];
    list extra_dgs=@EXTRA@[3];
  }
  int NN=size(extra_dgs);
  int counterEPl=0;
  setring BRing;
  poly h=subst(H,z,1);
  int Tz2=deg(h);
  int Tz3=Tz1-Tz2;
  int i,j,k,l,m,n,s,np,NP,I_N;
  if (Tz3==0)
  {
    // if this still does not work -> try always with ALL points in
    // Inf_Points !!!!
    poly Hinf=subst(H,z,0);
    setring aff_r;
    // compute the points at infinity of H and see which of them are in
    // Inf_Points
    poly h=imap(BRing,h);
    poly hinf=imap(BRing,Hinf);
    ideal pinf=factorize(hinf,1);
    list TIP=Inf_Points[1]+Inf_Points[2];
    s=size(TIP);
    NP=size(pinf);
    for (i=1;i<=NP;i=i+1)
    {
      for (j=1;j<=s;j=j+1)
      {
        if (pinf[i]==TIP[j][1])
        {
          np=size(TIP[j][2]);
          for (k=1;k<=np;k=k+1)
          {
            n=TIP[j][2][k];
            l=Places[n][1];
            m=Places[n][2];
            def SS=CURVE[5][l][1];
            setring SS;
            // local equation h of H
            if (POINTS[m][2]==number(1))
            {
              number A=POINTS[m][1];
              map psi=BRing,x+A,1,y;
              kill A;
            }
            else
            {
              map psi=BRing,1,x,y;
            }
            poly h=psi(H);
            I_N=local_IN(h,m)[1];
            interdiv[n]=I_N;
            kill h,psi;
            setring aff_r;
            kill SS;
          }
          break;
        }
      }
    }
    kill hinf,pinf;
  }
  else
  {
    // H is a multiple of z and hence all the points in Inf_Points intersect
    // with H
    setring aff_r;
    poly h=imap(BRing,h);
    list TIP=Inf_Points[1]+Inf_Points[2];
    s=size(TIP);
    for (j=1;j<=s;j=j+1)
    {
      np=size(TIP[j][2]);
      for (k=1;k<=np;k=k+1)
      {
        n=TIP[j][2][k];
        l=Places[n][1];
        m=Places[n][2];
        def SS=CURVE[5][l][1];
        setring SS;
        // local equation h of H
        if (POINTS[m][2]==number(1))
        {
          number A=POINTS[m][1];
          map psi=BRing,x+A,1,y;
          kill A;
        }
        else
        {
          map psi=BRing,1,x,y;
        }
        poly h=psi(H);
        I_N=local_IN(h,m)[1];
        interdiv[n]=I_N;
        kill h,psi;
        setring aff_r;
        kill SS;
      }
    }
  }
  // compute common affine points of H and CHI
  ideal CAL=h,CHI;
  CAL=std(CAL);
  if (CAL<>1)
  {
    list TAP=list();
    TAP=closed_points(CAL);
    NP=size(TAP);
    list auxP=list();
    int dP;
    for (i=1;i<=NP;i=i+1)
    {
      if (belongs(TAP[i],Aff_SLocus)==1)
      {
        // search the point in the list Aff_SPoints
        j=isInLP(TAP[i],Aff_SPoints);
        np=size(Aff_SPoints[j][2]);
        for (k=1;k<=np;k=k+1)
        {
          n=Aff_SPoints[j][2][k];
          l=Places[n][1];
          m=Places[n][2];
          def SS=CURVE[5][l][1];
          setring SS;
          // local equation h of H
          number A=POINTS[m][1];
          number B=POINTS[m][2];
          map psi=BRing,x+A,y+B,1;
          poly h=psi(H);
          I_N=local_IN(h,m)[1];
          interdiv[n]=I_N;
          kill A,B,h,psi;
          setring aff_r;
          kill SS;
        }
      }
      else
      {
        auxP=list();
        auxP[1]=TAP[i];
        dP=degree_P(auxP);
        if (defined(Aff_Points(dP))<>0)
        {
          // search the point in the list Aff_Points(dP)
          j=isInLP(TAP[i],Aff_Points(dP));
          n=Aff_Points(dP)[j][2][1];
          l=Places[n][1];
          m=Places[n][2];
          def SS=CURVE[5][l][1];
          setring SS;
          // local equation h of H
          number A=POINTS[m][1];
          number B=POINTS[m][2];
          map psi=BRing,x+A,y+B,1;
          poly h=psi(H);
          I_N=local_IN(h,m)[1];
          interdiv[n]=I_N;
          kill A,B,h,psi;
          setring aff_r;
          kill SS;
        }
        else
        {
          // previously check if it is an existing "extra place"
          j=isInLP(TAP[i],EP);
          if (j>0)
          {
            def SS=ES[j];
            setring SS;
            // local equation h of H
            number A=POINT[1];
            number B=POINT[2];
            map psi=BRing,x+A,y+B,1;
            poly h=psi(H);
            I_N=local_IN(h,0)[1];
            interdiv[N+j]=I_N;
            setring aff_r;
            kill SS;
          }
          else
          {
            // then we must create a new "extra place"
            counterEPl=counterEPl+1;
            list EXTRA_PLACE=extra_place(TAP[i]);
            def SS=EXTRA_PLACE[1];
            ES[NN+counterEPl]=SS;
            extra_dgs[NN+counterEPl]=EXTRA_PLACE[2];
            EP[NN+counterEPl]=list();
            EP[NN+counterEPl][1]=TAP[i];
            EP[NN+counterEPl][2]=0;
            setring SS;
            // local equation h of H
            number A=POINT[1];
            number B=POINT[2];
            map psi=BRing,x+A,y+B,1;
            poly h=psi(H);
            I_N=local_IN(h,0)[1];
            kill A,B,h,psi;
            interdiv[N+NN+counterEPl]=I_N;
            setring aff_r;
            kill SS;
          }
        }
      }
    }
    kill TAP,auxP;
  }
  kill h,CAL,TIP;
  @EXTRA@[1]=EP;
  @EXTRA@[2]=ES;
  @EXTRA@[3]=extra_dgs;
  kill EP;
  list update_CURVE=CURVE;
  if (size(extra_dgs)>0)
  {
    export(@EXTRA@);
    update_CURVE[1][1]=basering;
  }
  else
  {
    kill @EXTRA@;
  }
  setring BRing;
  kill h;
  kill aff_r;
  list answer=list();
  answer[1]=interdiv;
  answer[2]=update_CURVE;
  option(set,opgt);
  return(answer);
}
///////////////////////////////////////////////////////////////////////////////
static proc local_eq (poly H,SS,int m)
{
  // computes a local equation of polynomial H in the ring SS related to the place
  //     "m"
  // list POINT/POINTS is searched depending on wether m=0 or m>0 respectively
  // warning : the procedure must be called from ring "Proj_R" and returns a
  //     string
  def BRing=basering;
  setring SS;
  if (m>0)
  {
    if (POINTS[m][3]==number(1))
    {
      number A=POINTS[m][1];
      number B=POINTS[m][2];
      map psi=BRing,x+A,y+B,1;
    }
    else
    {
      if (POINTS[m][2]==number(1))
      {
        number A=POINTS[m][1];
        map psi=BRing,x+A,1,y;
      }
      else
      {
        map psi=BRing,1,x,y;
      }
    }
  }
  else
  {
    number A=POINT[1];
    number B=POINT[2];
    map psi=BRing,x+A,y+B,1;
  }
  poly h=psi(H);
  string str_h=string(h);
  setring BRing;
  return(str_h);
}
///////////////////////////////////////////////////////////////////////////////
static proc min_wt_rmat (matrix M)
{
  // finds the row of M with minimum non-zero entries, i.e. minimum
  // "Hamming-weight"
  int m=nrows(M);
  int n=ncols(M);
  int i,j;
  int Hwt=0;
  for (j=1;j<=n;j=j+1)
  {
    if (M[1,j]<>0)
    {
      Hwt=Hwt+1;
    }
  }
  int minHwt=Hwt;
  int k=1;
  for (i=2;i<=m;i=i+1)
  {
    Hwt=0;
    for (j=1;j<=n;j=j+1)
    {
      if (M[i,j]<>0)
      {
        Hwt=Hwt+1;
      }
    }
    if (Hwt<minHwt)
    {
      minHwt=Hwt;
      k=i;
    }
  }
  return(k);
}
///////////////////////////////////////////////////////////////////////////////

// ============================================================================
// *******    MAIN PROCEDURE : the Brill-Noether algorithm             ********
// ============================================================================

proc BrillNoether (intvec G,list CURVE)
"USAGE:    BrillNoether(G,CURVE); G an intvec, CURVE a list
RETURN:   list of ideals (each of them with two homogeneous generators,
          which represent the numerator, resp. denominator, of a rational
          function).@*
          The corresponding rational functions form a vector space basis of the
          linear system L(G), G a rational divisor over a non-singular curve.
NOTE:     The procedure must be called from the ring CURVE[1][2], where
          CURVE is the output of the procedure @code{NSplaces}. @*
          The intvec G represents a rational divisor supported on the closed
          places of CURVE[3] (e.g. @code{G=2,0,-1;} means 2 times the closed
          place 1 minus 1 times the closed place 3).
SEE ALSO: Adj_div, NSplaces, Weierstrass
EXAMPLE:  example BrillNoether; shows an example
"
{
  // computes a vector basis for the space L(G),
  //     where G is a given rational divisor over the non-singular curve
  // returns : list of ideals in R each with 2 elements H,Ho such that
  //           the set of functions {H/Ho} is the searched basis
  // warning : the conductor and sufficiently many points of the plane
  //           curve should be computed in advance, in list CURVE
  // the algorithm of Brill-Noether is carried out in the procedure
  def BRing=basering;
  int degX=CURVE[2][1];
  list Places=CURVE[3];
  intvec Conductor=CURVE[4];
  if (deg_D(G,Places)<0)
  {
    return(list());
  }
  intvec nuldiv;
  if (G==nuldiv)
  {
    list quickL=list();
    ideal quickId;
    quickId[1]=1;
    quickId[2]=1;
    quickL[1]=quickId;
    return(quickL);
  }
  intvec J=max_D(G,nuldiv)+Conductor;
  int n=estim_n(J,degX,Places);
  dbprint(printlevel+1,"Forms of degree "+string(n)+" : ");
  matrix W=nmultiples(n,degX,CHI);
  kill nFORMS(n-degX);
  list update_CURVE=CURVE;
  matrix V=interpolating_forms(J,n,update_CURVE);
  matrix VmW=supplement(W,V);
  int k=min_wt_rmat(VmW);
  int N=size(nFORMS(n));
  matrix H0[1][N];
  int i,j;
  for (i=1;i<=N;i=i+1) { H0[1,i]=VmW[k,i]; }
  poly Ho;
  for (i=1;i<=N;i=i+1) { Ho=Ho+(H0[1,i]*nFORMS(n)[i]); }
  list INTERD=intersection_div(Ho,update_CURVE);
  intvec NHo=INTERD[1];
  update_CURVE=INTERD[2];
  intvec AR=NHo-G;
  matrix V2=interpolating_forms(AR,n,update_CURVE);
  def aux_RING=update_CURVE[1][1];
  setring aux_RING;
  if (defined(@EXTRA@)<>0)
  {
    kill @EXTRA@;
  }
  setring BRing;
  update_CURVE[1][1]=aux_RING;
  kill aux_RING;
  matrix B0=supplement(W,V2);
  if (Hamming_wt(B0)==0)
  {
    return(list());
  }
  int ld=nrows(B0);
  list Bres=list();
  ideal HH;
  poly H;
  for (j=1;j<=ld;j=j+1)
  {
    H=0;
    for (i=1;i<=N;i=i+1)
    {
      H=H+(B0[j,i]*nFORMS(n)[i]);
    }
    HH=H,Ho;
    Bres[j]=simplifyRF(HH);
  }
  kill nFORMS(n);
  dbprint(printlevel+1," ");
  dbprint(printlevel+2,"Vector basis successfully computed ");
  dbprint(printlevel+1," ");
  return(Bres);
}
example
{
  "EXAMPLE:";  echo = 2;
  int plevel=printlevel;
  printlevel=-1;
  ring s=2,(x,y),lp;
  list C=Adj_div(x3y+y3+x);
  C=NSplaces(1..4,C);
  // the first 3 Places in C[3] are of degree 1.
  // we define the rational divisor G = 4*C[3][1]+4*C[3][3] (of degree 8):
  intvec G=4,0,4;
  def R=C[1][2];
  setring R;
  list LG=BrillNoether(G,C);
  // here is the vector basis of L(G):
  LG;
  printlevel=plevel;
}
///////////////////////////////////////////////////////////////////////////////

// *** procedures for dealing with "RATIONAL FUNCTIONS" over a plane curve ***
// a rational function F may be given by (homogeneous) ideal or (affine) poly
// (or number)

static proc polytoRF (F)
{
  // converts a polynomial (or number) into a "rational function" of type "ideal"
  // warning : it must be called inside "R" and the polynomial is expected to be affine
  ideal RF;
  RF[1]=homog(F,z);
  RF[2]=z^(deg(F));
  return(RF);
}
///////////////////////////////////////////////////////////////////////////////
static proc simplifyRF (ideal F)
{
  // simplifies a rational function f/g extracting the gcd(f,g)
  // maybe add a "restriction" to the curve "CHI" but it is not easy to
  // programm
  poly auxp=gcd(F[1],F[2]);
  return(ideal(division(F,auxp)[1]));
}
///////////////////////////////////////////////////////////////////////////////
static proc sumRF (F,G)
{
  // sum of two "rational functions" F,G given by either a pair
  // numerator/denominator or a poly
  if ( typeof(F)=="ideal" && typeof(G)=="ideal" )
  {
    ideal FG;
    FG[1]=F[1]*G[2]+F[2]*G[1];
    FG[2]=F[2]*G[2];
    return(simplifyRF(FG));
  }
  else
  {
    if (typeof(F)=="ideal")
    {
      ideal GG=polytoRF(G);
      ideal FG;
      FG[1]=F[1]*GG[2]+F[2]*GG[1];
      FG[2]=F[2]*GG[2];
      return(simplifyRF(FG));
    }
    else
    {
      if (typeof(G)=="ideal")
      {
        ideal FF=polytoRF(F);
        ideal FG;
        FG[1]=FF[1]*G[2]+FF[2]*G[1];
        FG[2]=FF[2]*G[2];
        return(simplifyRF(FG));
      }
      else
      {
        return(F+G);
      }
    }
  }
}
///////////////////////////////////////////////////////////////////////////////
static proc negRF (F)
{
  // returns -F as rational function
  if (typeof(F)=="ideal")
  {
    ideal FF=F;
    FF[1]=-F[1];
    return(FF);
  }
  else
  {
    return(-F);
  }
}
///////////////////////////////////////////////////////////////////////////////
static proc escprodRF (l,F)
{
  // computes l*F as rational function
  // l should be either a number or a polynomial of degree zero
  if (typeof(F)=="ideal")
  {
    ideal lF=F;
    lF[1]=l*F[1];
    return(lF);
  }
  else
  {
    return(l*F);
  }
}
///////////////////////////////////////////////////////////////////////////////

// ******** procedures to compute Weierstrass semigroups ********

static proc orderRF (ideal F,SS,int m)
"USAGE:  orderRF(F,SS,m); F an ideal, SS a ring and m an integer
RETURN: list with the order (int) and the leading coefficient (number)
NOTE:   F represents a rational function, thus the procedure must be
        called from global ring R or R(d).
        SS contains the name of a local ring where rational places are
        stored, and then we take that which is in position m in the
        corresponding lists of data.
        The order of F at the place given by SS,m is returned together
        with the coefficient of minimum degree in the corresponding power
        series.
"
{
  // computes the order of a rational function F at a RATIONAL place given by
  //     a local ring SS and a position "m" inside SS
  // warning : the procedure must be called from global projective ring "R" or
  //     "R(i)"
  // returns a list with the order (int) and the "leading coefficient" (number)
  def BR=basering;
  poly f=F[1];
  string sf=local_eq(f,SS,m);
  poly g=F[2];
  string sg=local_eq(g,SS,m);
  setring SS;
  execute("poly ff="+sf+";");
  execute("poly gg="+sg+";");
  list o1=local_IN(ff,m);
  list o2=local_IN(gg,m);
  int oo=o1[1]-o2[1];
  number lc=o1[2]/o2[2];
  setring BR;
  number LC=number(imap(SS,lc));
  return(list(oo,LC));
}
///////////////////////////////////////////////////////////////////////////////

proc Weierstrass (int P,int m,list CURVE)
"USAGE:    Weierstrass( i, m, CURVE );  i,m integers and CURVE a list
RETURN:   list WS of two lists:
  @format
  WS[1] list of integers (Weierstr. semigroup of the curve at place i up to m)
  WS[2] list of ideals (the associated rational functions)
  @end format
NOTE:     The procedure must be called from the ring CURVE[1][2],
          where CURVE is the output of the procedure @code{NSplaces}.
@*        i represents the place CURVE[3][i].
@*        Rational functions are represented by numerator/denominator
          in form of ideals with two homogeneous generators.
WARNING:  The place must be rational, i.e., necessarily CURVE[3][i][1]=1. @*
SEE ALSO: Adj_div, NSplaces, BrillNoether
EXAMPLE:  example Weierstrass; shows an example
"
{
  // computes the Weierstrass semigroup at a RATIONAL place P up to a bound "m"
  //   together with the functions achieving each value up to m, via
  //   Brill-Noether
  // returns 2 lists : the first consists of all the poles up to m in
  //   increasing order and the second consists of the corresponging rational
  //   functions
  list Places=CURVE[3];
  intvec pl=Places[P];
  int dP=pl[1];
  int nP=pl[2];
  if (dP<>1)
  {
    ERROR("the given place is not defined over the prime field");
  }
  if (m<=0)
  {
    if (m==0)
    {
      list semig=list();
      int auxint=0;
      semig[1]=auxint;
      list funcs=list();
      ideal auxF;
      auxF[1]=1;
      auxF[2]=1;
      funcs[1]=auxF;
      return(list(semig,funcs));
    }
    else
    {
      ERROR("second argument must be non-negative");
    }
  }
  int auxint=0;
  ideal auxF;
  auxF[1]=1;
  auxF[2]=1;
  // Brill-Noether algorithm
  intvec mP;
  mP[P]=m;
  list LmP=BrillNoether(mP,CURVE);
  int lmP=size(LmP);
  if (lmP==1)
  {
    list semig=list();
    semig[1]=auxint;
    list funcs=list();
    funcs[1]=auxF;
    return(list(semig,funcs));
  }
  def SS=CURVE[5][dP][1];
  list ordLmP=list();
  list polLmP=list();
  list sortpol=list();
  int maxpol;
  int i,j,k;
  for (i=1;i<=lmP-1;i=i+1)
  {
    for (j=1;j<=lmP-i+1;j=j+1)
    {
      ordLmP[j]=orderRF(LmP[j],SS,nP);
      polLmP[j]=-ordLmP[j][1];
    }
    sortpol=sort(polLmP);
    polLmP=sortpol[1];
    maxpol=polLmP[lmP-i+1];
    LmP=permute_L(LmP,sortpol[2]);
    ordLmP=permute_L(ordLmP,sortpol[2]);
    // triangulate LmP
    for (k=1;k<=lmP-i;k=k+1)
    {
      if (polLmP[lmP-i+1-k]==maxpol)
      {
        LmP[lmP-i+1-k]=sumRF(LmP[lmP-i+1-k],negRF(escprodRF(
                       ordLmP[lmP-i+1-k][2]/ordLmP[lmP-i+1][2],LmP[lmP-i+1])));
      }
      else
      {
        break;
      }
    }
  }
  polLmP[1]=auxint;
  LmP[1]=auxF;
  return(list(polLmP,LmP));
}
example
{
  "EXAMPLE:";  echo = 2;
  int plevel=printlevel;
  printlevel=-1;
  ring s=2,(x,y),lp;
  list C=Adj_div(x3y+y3+x);
  C=NSplaces(1..4,C);
  def R=C[1][2];
  setring R;
  // Place C[3][1] has degree 1 (i.e it is rational);
  list WS=Weierstrass(1,7,C);
  // the first part of the list is the Weierstrass semigroup up to 7 :
  WS[1];
  // and the second part are the corresponding functions :
  WS[2];
  printlevel=plevel;
}
///////////////////////////////////////////////////////////////////////////////

// axiliary procedure for permuting a list or intvec

proc permute_L (L,P)
"USAGE:    permute_L( L, P ); L,P either intvecs or lists
RETURN:   list obtained from L by applying the permutation given by P.
NOTE:     If P is a list, all entries must be integers.
SEE ALSO: sys_code, AGcode_Omega, prepSV
EXAMPLE:  example permute_L; shows an example
"
{
  // permutes the list L according to the permutation P (both intvecs or
  //  lists of integers)
  int s=size(L);
  int n=size(P);
  int i;
  if (s<n)
  {
    for (i=s+1;i<=n;i=i+1)
    {
      L[i]=0;
    }
    s=size(L);
  }
  list auxL=L;
  for (i=1;i<=n;i=i+1)
  {
    auxL[i]=L[P[i]];
  }
  return(auxL);
}
example
{
  "EXAMPLE:";  echo = 2;
  list L=list();
  L[1]="a";
  L[2]="b";
  L[3]="c";
  L[4]="d";
  intvec P=1,3,4,2;
  // the list L is permuted according to P :
  permute_L(L,P);
}
///////////////////////////////////////////////////////////////////////////////
static proc evalRF (ideal F,SS,int m)
"USAGE:     evalRF(F,SS,m), where F is an ideal, SS is a ring and m is an
            integer
RETURN:     the evaluation (number) of F at the place given by SS,m if it is
            well-defined
NOTE:       F represents a rational function, thus the procedure must be
            called from R or R(d).
            SS contains the name of a local ring where rational places are
            stored, and then we take that which is in position m in the
            corresponding lists of data.
"
{
  // evaluates a rational function F at a RATIONAL place given by
  //     a local ring SS and a position "m" inside SS
  list olc=orderRF(F,SS,m);
  int oo=olc[1];
  if (oo==0)
  {
    return(olc[2]);
  }
  else
  {
    if (oo>0)
    {
      return(number(0));
    }
    else
    {
      ERROR("the function is not well-defined at the given place");
    }
  }
}
///////////////////////////////////////////////////////////////////////////////
//
// ******** procedures for constructing AG codes ********
//
///////////////////////////////////////////////////////////////////////////////

static proc gen_mat (list LF,intvec LP,RP)
"USAGE:      gen_mat(LF,LP,RP); LF list of rational functions,
            LP intvec of rational places and RP a local ring
RETURN:     a generator matrix of the evaluation code given by LF and LP
SEE ALSO:   extcurve
KEYWORDS:   evaluation codes
NOTE:       Rational places are searched in local ring RP
            The procedure must be called from R or R(d) fromlist CURVE
            after having executed extcurve(d,CURVE)
"
{
  // computes a generator matrix (with numbers) of the evaluation code given
  //    by a list of rational functions LF and a list of RATIONAL places LP
  int m=size(LF);
  int n=size(LP);
  matrix GM[m][n];
  int i,j;
  for (i=1;i<=m;i=i+1)
  {
    for (j=1;j<=n;j=j+1)
    {
      GM[i,j]=evalRF(LF[i],RP,LP[j]);
    }
  }
  return(GM);
}
///////////////////////////////////////////////////////////////////////////////

proc dual_code (matrix G)
"USAGE:   dual_code(G); G a matrix of numbers
RETURN:   a generator matrix of the dual code generated by G
NOTE:     The input should be a matrix G of numbers. @*
          The output is also a parity check matrix for the code defined by G
KEYWORDS: linear code, dual
EXAMPLE:  example dual_code; shows an example
"
{
  // computes the dual code of C given by a generator matrix G
  //     i.e. computes a parity-check matrix H of C
  // conversely : computes also G if H is given
  return(Ker(G));
}
example
{
  "EXAMPLE:";  echo = 2;
  ring s=2,T,lp;
  // here is the Hamming code of length 7 and dimension 3
  matrix G[3][7]=1,0,1,0,1,0,1,0,1,1,0,0,1,1,0,0,0,1,1,1,1;
  print(G);
  matrix H=dual_code(G);
  print(H);
}
///////////////////////////////////////////////////////////////////////////////

// ======================================================================
// *********** initial test for disjointness ***************
// ======================================================================

static proc disj_divs (intvec H,intvec P,list EC)
{
  int s1=size(H);
  int s2=size(P);
  list PLACES=EC[3];
  def BRing=basering;
  def auxR=EC[1][5];
  setring auxR;
  int s=res_deg();
  setring BRing;
  kill auxR;
  int i,j,k,d,l,N,M;
  intvec auxIV,iw;
  for (i=1;i<=s;i=i+1)
  {
    if ( (s mod i) == 0 )
    {
      if (typeof(EC[5][i])=="list")
      {
        def auxR=EC[5][i][1];
        setring auxR;
        auxIV[i]=size(POINTS);
        setring BRing;
        kill auxR;
      }
      else
      {
        auxIV[i]=0;
      }
    }
    else
    {
      auxIV[i]=0;
    }
  }
  for (i=1;i<=s1;i=i+1)
  {
    if (H[i]<>0)
    {
      iw=PLACES[i];
      d=iw[1];
      if ( (s mod d) == 0 )
      {
        l=iw[2];
        // check that this place(s) is/are not in sup(D)
        if (d==1)
        {
          for (j=1;j<=s2;j=j+1)
          {
            if (l==P[j])
            {
              return(0);
            }
          }
        }
        else
        {
          N=0;
          for (j=1;j<d;j=j+1)
          {
            N=N+j*auxIV[j];
          }
          N=N+d*(l-1);
          M=N+d;
          for (k=N+1;k<=M;k=k+1)
          {
            for (j=1;j<=s2;j=j+1)
            {
              if (k==P[j])
              {
                return(0);
              }
            }
          }
        }
      }
    }
  }
  kill auxIV,iw;
  return(1);
}
///////////////////////////////////////////////////////////////////////////////

proc AGcode_L (intvec G,intvec D,list EC)
"USAGE:    AGcode_L( G, D, EC );  G,D intvec, EC a list
RETURN:   a generator matrix for the evaluation AG code defined by the
          divisors G and D.
NOTE:     The procedure must be called within the ring EC[1][4],
          where EC is the output of @code{extcurve(d)} (or within
          the ring EC[1][2] if d=1). @*
          The entry i in the intvec D refers to the i-th rational
          place in EC[1][5] (i.e., to POINTS[i], etc., see @ref{extcurve}).@*
          The intvec G represents a rational divisor (see @ref{BrillNoether}
          for more details).@*
          The code evaluates the vector space basis of L(G) at the rational
          places given by D.
WARNINGS: G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
          not checked by the algorithm.
          G and D should have disjoint supports (checked by the algorithm).
SEE ALSO: Adj_div, BrillNoether, extcurve, AGcode_Omega
EXAMPLE:  example AGcode_L; shows an example
"
{
  // returns a generator matrix for the evaluation AG code given by G and D
  // G must be a divisor defined over the prime field and D an intvec of
  // "rational places"
  // it must be called inside R or R(d) and requires previously "extcurve(d)"
  def BR=basering;
  if (disj_divs(G,D,EC)==0)
  {
    dbprint(printlevel+3,"? the divisors do not have disjoint supports,
                                 0-matrix returned ?");
    matrix answer;
    return(answer);
  }
  if (res_deg()>1)
  {
    def R=EC[1][2];
    setring R;
    list LG=BrillNoether(G,EC);
    setring BR;
    list LG=imap(R,LG);
    setring R;
    kill LG;
    setring BR;
    kill R;
  }
  else
  {
    list LG=BrillNoether(G,EC);
  }
  def RP=EC[1][5];
  matrix M=gen_mat(LG,D,RP);
  kill LG;
  return(M);
}
example
{
  "EXAMPLE:";  echo = 2;
  int plevel=printlevel;
  printlevel=-1;
  ring s=2,(x,y),lp;
  list HC=Adj_div(x3+y2+y);
  HC=NSplaces(1..2,HC);
  HC=extcurve(2,HC);
  def ER=HC[1][4];
  setring ER;
  intvec G=5;      // the rational divisor G = 5*HC[3][1]
  intvec D=2..9;   // D = sum of the rational places no. 2..9 over F_4
  // let us construct the corresponding evaluation AG code :
  matrix C=AGcode_L(G,D,HC);
  // here is a linear code of type [8,5,>=3] over F_4
  print(C);
  printlevel=plevel;
}
///////////////////////////////////////////////////////////////////////////////

proc AGcode_Omega (intvec G,intvec D,list EC)
"USAGE:    AGcode_Omega( G, D, EC ); G,D intvec, EC a list
RETURN:   a generator matrix for the residual AG code defined by the
          divisors G and D.
NOTE:     The procedure must be called within the ring EC[1][4],
          where EC is the output of @code{extcurve(d)} (or within
          the ring EC[1][2] if d=1). @*
          The entry i in the intvec D refers to the i-th rational
          place in EC[1][5] (i.e., to POINTS[i], etc., see @ref{extcurve}).@*
          The intvec G represents a rational divisor (see @ref{BrillNoether}
          for more details).@*
          The code computes the residues of a vector space basis of
          @math{\Omega(G-D)} at the rational places given by D.
WARNINGS: G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
          not checked by the algorithm.
          G and D should have disjoint supports (checked by the algorithm).
SEE ALSO: Adj_div, BrillNoether, extcurve, AGcode_L
EXAMPLE:  example AGcode_Omega; shows an example
"
{
  // returns a generator matrix for the residual AG code given by G and D
  // G must be a divisor defined over the prime field and D an intvec or
  // "rational places"
  // it must be called inside R or R(d) and requires previously "extcurve(d)"
  return(dual_code(AGcode_L(G,D,EC)));
}
example
{
  "EXAMPLE:";  echo = 2;
  int plevel=printlevel;
  printlevel=-1;
  ring s=2,(x,y),lp;
  list HC=Adj_div(x3+y2+y);
  HC=NSplaces(1..2,HC);
  HC=extcurve(2,HC);
  def ER=HC[1][4];
  setring ER;
  intvec G=5;      // the rational divisor G = 5*HC[3][1]
  intvec D=2..9;   // D = sum of the rational places no. 2..9 over F_4
  // let us construct the corresponding residual AG code :
  matrix C=AGcode_Omega(G,D,HC);
  // here is a linear code of type [8,3,>=5] over F_4
  print(C);
  printlevel=plevel;
}
///////////////////////////////////////////////////////////////////////////////

// ============================================================================
// *******   auxiliary procedure to define AG codes over extensions    ********
// ============================================================================

proc extcurve (int d,list CURVE)
"USAGE:    extcurve( d, CURVE ); d an integer, CURVE a list
RETURN:   list L which is the update of the list CURVE with additional entries
   @format
   L[1][3]: ring (p,a),(x,y),lp (affine),
   L[1][4]: ring (p,a),(x,y,z),lp (projective),
   L[1][5]: ring (p,a),(x,y,t),ls (local),
   L[2][3]: int  (the number of rational places),
   @end format
          the rings being defined over a field extension of degree d. @*
          If d<2 then @code{extcurve(d,CURVE);} creates a list L which
          is the update of the list CURVE with additional entries
   @format
   L[1][5]: ring p,(x,y,t),ls,
   L[2][3]: int  (the number of computed places over the base field).
   @end format
          In both cases, in the ring L[1][5] lists with the data for all the
          computed rational places (after a field extension of degree d) are
          created (see @ref{Adj_div}):
   @format
   lists POINTS, LOC_EQS, BRANCHES, PARAMETRIZATIONS.
   @end format
NOTE:     The list CURVE should be the output of @code{NSplaces},
          and must contain (at least) one place of degree d. @*
          You actually need all the places with degree dividing d.
          Otherwise, not all the places are computed, but only part of them. @*
          This procedure must be executed before constructing AG codes,
          even if no extension is needed. The ring L[1][4] must be active
          when constructing codes over the field extension.@*
SEE ALSO: closed_points, Adj_div, NSplaces, AGcode_L, AGcode_Omega
EXAMPLE:  example extcurve; shows an example
"
{
  // extends the underlying curve and all its associated objects to a larger
  //     base field in order to evaluate points over such a extension
  // if d<2 then the only change is that a local ring "RatPl" (which is a
  //     copy of "S(1)") is created in order to store the rational places
  //     where we can do evaluations
  // otherwise, such a ring contains all places which are rational over the
  //     extension
  // warning : list Places does not change so that only divisors which are
  //     "rational over the prime field" are allowed; this probably will
  //     change in the future
  // warning : the places in RatPl are ranged in increasing degree, respecting
  //     the order from list Places and placing the conjugate branches all
  //     together
  def BR=basering;
  list ext_CURVE=CURVE;
  if (d<2)
  {
    def SS=CURVE[5][1][1];
    ring RatPl=char(basering),(x,y,t),ls;
    list POINTS=imap(SS,POINTS);
    list LOC_EQS=imap(SS,LOC_EQS);
    list BRANCHES=imap(SS,BRANCHES);
    list PARAMETRIZATIONS=imap(SS,PARAMETRIZATIONS);
    export(POINTS);
    export(LOC_EQS);
    export(BRANCHES);
    export(PARAMETRIZATIONS);
    int NrRatPl=size(POINTS);
    ext_CURVE[2][3]=NrRatPl;
    setring BR;
    ext_CURVE[1][5]=RatPl;
    dbprint(printlevel+1,"");
    dbprint(printlevel+2,"Total number of rational places : "+string(NrRatPl));
    dbprint(printlevel+1,"");
    kill RatPl;
    return(ext_CURVE);
  }
  else
  {
    // exclude the case when no place of degree d was previously computed/found
    int dd=size(CURVE[5]);
    if (dd<d)
    {
      ERROR("you did not compute/find any place of degree "+string(d));
    }
    if (typeof(CURVE[5][d])<>"list")
    {
      ERROR("you did not compute/find any place of degree "+string(d));
    }

    // Define the ring "RatPl" :
    def S(d)=CURVE[5][d][1];
    setring S(d);
    string smin=string(minpoly);
    setring BR;
    ring RatPl=(char(basering),a),(x,y,t),ls;
    execute("minpoly="+smin+";");

    // import data from local ring S(d)
    list POINTS=imap(S(d),POINTS);
    list LOC_EQS=imap(S(d),LOC_EQS);
    list BRANCHES=imap(S(d),BRANCHES);
    list PARAMETRIZATIONS=imap(S(d),PARAMETRIZATIONS);
    kill S(d);

    // conjugate data from S(d)
    int s=size(POINTS);
    int counter=0;
    int piv=0;
    int i,j,k;
    for (j=1;j<=s;j=j+1)
    {
      counter=counter+1;
      piv=counter;
      for (k=1;k<d;k=k+1)
      {
        POINTS=insert(POINTS,Frobenius(POINTS[piv],1),counter);
        LOC_EQS=insert(LOC_EQS,Frobenius(LOC_EQS[piv],1),counter);
        BRANCHES=insert(BRANCHES,conj_b(BRANCHES[piv],1),counter);
        PARAMETRIZATIONS=insert(PARAMETRIZATIONS,Frobenius(
                                  PARAMETRIZATIONS[piv],1),counter);
        counter=counter+1;
        piv=counter;
      }
    }

    // now the same thing for smaller base fields
    string olda;
    poly paux;
    intvec iv,iw;
    int ii,jj,kk,m,n;
    for (i=d-1;i>=2;i=i-1)
    {
      if ( (d mod i)==0 )
      {
        if (typeof(CURVE[5][i])=="list")
        {
          def S(i)=CURVE[5][i][1];
          // embedd S(i) inside basering == RatPl " ==S(d) "
          olda=subfield(S(i));
          setring S(i);
          s=size(POINTS);
          setring RatPl;
          // import data from S(i)
          for (j=s;j>=1;j=j-1)
          {
            counter=0;
            POINTS=insert(POINTS,list(),0);
            POINTS[1][1]=number(importdatum(S(i),"POINTS["+string(j)
                                                          +"][1]",olda));
            POINTS[1][2]=number(importdatum(S(i),"POINTS["+string(j)
                                                          +"][2]",olda));
            POINTS[1][3]=number(importdatum(S(i),"POINTS["+string(j)
                                                          +"][3]",olda));
            LOC_EQS=insert(LOC_EQS,importdatum(S(i),"LOC_EQS["+string(j)
                                                          +"]",olda),0);
            BRANCHES=insert(BRANCHES,list(),0);
            setring S(i);
            m=nrows(BRANCHES[j][1]);
            n=ncols(BRANCHES[j][1]);
            iv=BRANCHES[j][2];
            kk=BRANCHES[j][3];
            poly par@1=subst(PARAMETRIZATIONS[j][1][1],t,x);
            poly par@2=subst(PARAMETRIZATIONS[j][1][2],t,x);
            export(par@1);
            export(par@2);
            iw=PARAMETRIZATIONS[j][2];
            setring RatPl;
            paux=importdatum(S(i),"BRANCHES["+string(j)+"][4]",olda);
            matrix Maux[m][n];
            for (ii=1;ii<=m;ii=ii+1)
            {
              for (jj=1;jj<=n;jj=jj+1)
              {
                Maux[ii,jj]=importdatum(S(i),"BRANCHES["+string(j)
                               +"][1]["+string(ii)+","+string(jj)+"]",olda);
              }
            }
            BRANCHES[1][1]=Maux;
            BRANCHES[1][2]=iv;
            BRANCHES[1][3]=kk;
            BRANCHES[1][4]=paux;
            kill Maux;
            PARAMETRIZATIONS=insert(PARAMETRIZATIONS,list(),0);
            PARAMETRIZATIONS[1][1]=ideal(0);
            PARAMETRIZATIONS[1][1][1]=importdatum(S(i),"par@1",olda);
            PARAMETRIZATIONS[1][1][2]=importdatum(S(i),"par@2",olda);
            PARAMETRIZATIONS[1][1][1]=subst(PARAMETRIZATIONS[1][1][1],x,t);
            PARAMETRIZATIONS[1][1][2]=subst(PARAMETRIZATIONS[1][1][2],x,t);
            PARAMETRIZATIONS[1][2]=iw;
            for (k=1;k<i;k=k+1)
            {
              counter=counter+1;
              piv=counter;
              POINTS=insert(POINTS,Frobenius(POINTS[piv],1),counter);
              LOC_EQS=insert(LOC_EQS,Frobenius(LOC_EQS[piv],1),counter);
              BRANCHES=insert(BRANCHES,conj_b(BRANCHES[piv],1),counter);
              PARAMETRIZATIONS=insert(PARAMETRIZATIONS,Frobenius(
                                        PARAMETRIZATIONS[piv],1),counter);
            }
            setring S(i);
            kill par@1,par@2;
            setring RatPl;
          }
          kill S(i);
        }
      }
    }
    kill iw;
    kill paux;

    // finally add points from S(1)
    if (typeof(CURVE[5][1])=="list")
    {
      def S(1)=CURVE[5][1][1];
      POINTS=imap(S(1),POINTS)+POINTS;
      LOC_EQS=imap(S(1),LOC_EQS)+LOC_EQS;
      BRANCHES=imap(S(1),BRANCHES)+BRANCHES;
      PARAMETRIZATIONS=imap(S(1),PARAMETRIZATIONS)+PARAMETRIZATIONS;
    }

    // prepare data for output
    export(POINTS);
    export(LOC_EQS);
    export(BRANCHES);
    export(PARAMETRIZATIONS);
    int NrRatPl=size(POINTS);
    ext_CURVE[2][3]=NrRatPl;
    setring BR;
    ext_CURVE[1][5]=RatPl;
    ring r(d)=(char(basering),a),(x,y),lp;
    execute("minpoly="+smin+";");
    setring BR;
    ext_CURVE[1][3]=r(d);
    ring R(d)=(char(basering),a),(x,y,z),lp;
    execute("minpoly="+smin+";");
    setring BR;
    ext_CURVE[1][4]=R(d);
    dbprint(printlevel+1,"");
    dbprint(printlevel+2,"Total number of rational places : NrRatPl = "
                                   +string(NrRatPl));
    dbprint(printlevel+1,"");
    kill S(1);
    kill R(d);
    kill RatPl;
    return(ext_CURVE);
  }
}
example
{
  "EXAMPLE:";  echo = 2;
  int plevel=printlevel;
  printlevel=-1;
  ring s=2,(x,y),lp;
  list C=Adj_div(x5+y2+y);
  C=NSplaces(1..4,C);
  // since we have all points up to degree 4, we can extend the curve
  // to that extension, in order to get rational points over F_16;
  C=extcurve(4,C);
  // e.g., display the basepoint of place no. 32:
  def R=C[1][5];
  setring R;
  POINTS[32];
  printlevel=plevel;
}
///////////////////////////////////////////////////////////////////////////////

// specific procedures for linear/AG codes

static proc Hamming_wt (matrix A)
"USAGE:     Hamming_wt(A), where A is any matrix
RETURN:     the Hamming weight (number of non-zero entries) of the matrix A
"
{
  // computes the Hamming weight (number of non-zero entries) of any matrix
  // notice that "words" are represented by matrices of size 1xn
  // computing the Hamming distance between two matrices can be done by
  // Hamming_wt(A-B)
  int m=nrows(A);
  int n=ncols(A);
  int i,j;
  int w=0;
  for (i=1;i<=m;i=i+1)
  {
    for (j=1;j<=n;j=j+1)
    {
      if (A[i,j]<>0)
      {
        w=w+1;
      }
    }
  }
  return(w);
}
///////////////////////////////////////////////////////////////////////////////

// Basic Algorithm of Skorobogatov and Vladut for decoding AG codes
// warning : the user must choose carefully the parameters for the code and
//     the decoding since they will never be checked by the procedures

proc prepSV (intvec G,intvec D,intvec F,list EC)
"USAGE:    prepSV( G, D, F, EC ); G,D,F intvecs and EC a list
RETURN:   list E of size n+3, where n=size(D). All its entries but E[n+3]
          are matrices:
   @format
   E[1]:  parity check matrix for the current AG code
   E[2] ... E[n+2]:  matrices used in the procedure decodeSV
   E[n+3]:  intvec with
       E[n+3][1]: correction capacity @math{epsilon} of the algorithm
       E[n+3][2]: designed Goppa distance @math{delta} of the current AG code
   @end format
NOTE:     Computes the preprocessing for the basic (Skorobogatov-Vladut)
          decoding algorithm.@*
          The procedure must be called within the ring EC[1][4], where EC is
          the output of @code{extcurve(d)} (or in the ring EC[1][2] if d=1) @*
          The intvec G and F represent rational divisors (see
          @ref{BrillNoether} for more details).@*
          The intvec D refers to rational places (see @ref{AGcode_Omega}
          for more details.).
          The current AG code is @code{AGcode_Omega(G,D,EC)}.@*
          If you know the exact minimum distance d and you want to use it in
          @code{decodeSV} instead of @math{delta}, you can change the value
          of E[n+3][2] to d before applying decodeSV.
          If you have a systematic encoding for the current code and want to
          keep it during the decoding, you must previously permute D (using
          @code{permute_L(D,P);}), e.g., according to the permutation
          P=L[3], L being the output of @code{sys_code}.
WARNINGS: F must be a divisor with support disjoint from the support of D and
          with degree @math{epsilon + genus}, where
          @math{epsilon:=[(deg(G)-3*genus+1)/2]}.@*
          G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
          not checked by the algorithm.
          G and D should also have disjoint supports (checked by the
          algorithm).
KEYWORDS: SV-decoding algorithm, preprocessing
SEE ALSO: extcurve, AGcode_Omega, decodeSV, sys_code, permute_L
EXAMPLE:  example prepSV; shows an example
"
{
  if (disj_divs(F,D,EC)==0)
  {
    dbprint(printlevel+3,"? the divisors do not have disjoint supports,
                            empty list returned ?");
    return(list());
  }
  list E=list();
  list Places=EC[3];
  int m=deg_D(G,Places);
  int genusX=EC[2][2];
  int e=(m+1-3*genusX)/2;
  if (e<1)
  {
    dbprint(printlevel+3,"? the correction capacity of the basic algorithm
                            is zero, empty list returned ?");
    return(list());
  }
  // deg(F)==e+genusX should be satisfied, and sup(D),sup(F) should be
  // disjoint !!!!
  int n=size(D);
  // 2*genusX-2<m<n should also be satisfied !!!!
  matrix EE=AGcode_L(G,D,EC);
  int l=nrows(EE);
  E[1]=EE;
  matrix GP=AGcode_L(F,D,EC);
  int r=nrows(GP);
  intvec H=G-F;
  matrix HP=AGcode_L(H,D,EC);
  int s=nrows(HP);
  int i,j,k;
  kill EE;
  for (k=1;k<=n;k=k+1)
  {
    E[k+1]=GP[1,k]*submat(HP,1..s,k..k);
    for (i=2;i<=r;i=i+1)
    {
      E[k+1]=concat(E[k+1],GP[i,k]*submat(HP,1..s,k..k));
    }
  }
  E[n+2]=GP;
  intvec IW=e,m+2-2*genusX;
  E[n+3]=IW;
  kill IW;
  return(E);
}
example
{
  "EXAMPLE:";  echo = 2;
  int plevel=printlevel;
  printlevel=-1;
  ring s=2,(x,y),lp;
  list HC=Adj_div(x3+y2+y);
  HC=NSplaces(1..2,HC);
  HC=extcurve(2,HC);
  def ER=HC[1][4];
  setring ER;
  intvec G=5;      // the rational divisor G = 5*HC[3][1]
  intvec D=2..9;   // D = sum of the rational places no. 2..9 over F_4
  // construct the corresp. residual AG code of type [8,3,>=5] over F_4:
  matrix C=AGcode_Omega(G,D,HC);
  // we can correct 1 error and the genus is 1, thus F must have degree 2
  // and support disjoint from that of D;
  intvec F=2;
  list SV=prepSV(G,D,F,HC);
  // now everything is prepared to decode with the basic algorithm;
  // for example, here is a parity check matrix to compute the syndrome :
  print(SV[1]);
  // and here you have the correction capacity of the algorithm :
  int epsilon=SV[size(D)+3][1];
  epsilon;
  printlevel=plevel;
}
///////////////////////////////////////////////////////////////////////////////

proc decodeSV (matrix y,list K)
"USAGE:     decodeSV( y, K ); y a row-matrix and K a list
RETURN:     a codeword (row-matrix) if possible, resp. the 0-matrix (of size
            1) if decoding is impossible.
            For decoding the basic (Skorobogatov-Vladut) decoding algorithm
            is applied.
NOTE:       The list_expression should be the output K of the procedure
            @code{prepSV}.@*
            The matrix_expression should be a (1 x n)-matrix,  where
            n = ncols(K[1]).@*
            The decoding may fail if the number of errors is greater than
            the correction capacity of the algorithm.
KEYWORDS:   SV-decoding algorithm
SEE ALSO:   extcurve, AGcode_Omega, prepSV
EXAMPLE:    example decodeSV; shows an example
"
{
  // decodes y with the "basic decoding algorithm", if possible
  // requires the preprocessing done by the procedure "prepSV"
  // the procedure must be called from ring R or R(d)
  // returns either a codeword (matrix) of none (in case of too many errors)
  matrix syndr=K[1]*transpose(y);
  if (Hamming_wt(syndr)==0)
  {
    return(y);
  }
  matrix Ey=y[1,1]*K[2];
  int n=ncols(y);
  int i;
  for (i=2;i<=n;i=i+1)
  {
    Ey=Ey+y[1,i]*K[i+1];
  }
  matrix Ky=get_NZsol(Ey);
  if (Hamming_wt(Ky)==0)
  {
    dbprint(printlevel+3,"? no error-locator found ?");
    dbprint(printlevel+3,"? too many errors occur, 0-matrix returned ?");
    matrix answer;
    return(answer);
  }
  int r=nrows(K[n+2]);
  matrix ErrLoc[1][n];
  list Z=list();
  list NZ=list();
  int j;
  for (j=1;j<=n;j=j+1)
  {
    for (i=1;i<=r;i=i+1)
    {
      ErrLoc[1,j]=ErrLoc[1,j]+K[n+2][i,j]*Ky[1,i];
    }
    if (ErrLoc[1,j]==0)
    {
      Z=insert(Z,j,size(Z));
    }
    else
    {
      NZ=insert(NZ,j,size(NZ));
    }
  }
  int k=size(NZ);
  int l=nrows(K[1]);
  int s=l+k;
  matrix A[s][n];
  matrix b[s][1];
  for (i=1;i<=l;i=i+1)
  {
    for (j=1;j<=n;j=j+1)
    {
      A[i,j]=K[1][i,j];
    }
    b[i,1]=syndr[i,1];
  }
  for (i=1;i<=k;i=i+1)
  {
    A[l+i,NZ[i]]=number(1);
  }
  intvec opgt=option(get);
  option(redSB);
  matrix L=transpose(syz(concat(A,-b)));
  if (nrows(L)==1)
  {
    if (L[1,n+1]<>0)
    {
      poly pivote=L[1,n+1];
      matrix sol=submat(L,1..1,1..n);
      if (pivote<>1)
      {
        sol=(number(1)/number(pivote))*sol;
      }
      // check at least that the number of comitted errors is less than half
      //     the Goppa distance
      // imposing Hamming_wt(sol)<=K[n+3][1] would be more correct, but maybe
      //     is too strong
      // on the other hand, if Hamming_wt(sol) is too large the decoding may
      //     not be acceptable
      if ( Hamming_wt(sol) <= (K[n+3][2]-1)/2 )
      {
        option(set,opgt);
        return(y-sol);
      }
      else
      {
        dbprint(printlevel+3,"? non-acceptable decoding ?");
      }
    }
    else
    {
      dbprint(printlevel+3,"? no solution found ?");
    }
  }
  else
  {
    dbprint(printlevel+3,"? non-unique solution ?");
  }
  option(set,opgt);
  dbprint(printlevel+3,"? too many errors occur, 0-matrix returned ?");
  matrix answer;
  return(answer);
}
example
{
  "EXAMPLE:";  echo = 2;
  int plevel=printlevel;
  printlevel=-1;
  ring s=2,(x,y),lp;
  list HC=Adj_div(x3+y2+y);
  HC=NSplaces(1..2,HC);
  HC=extcurve(2,HC);
  def ER=HC[1][4];
  setring ER;
  intvec G=5;      // the rational divisor G = 5*HC[3][1]
  intvec D=2..9;   // D = sum of the rational places no. 2..9 over F_4
  // construct the corresp. residual AG code of type [8,3,>=5] over F_4:
  matrix C=AGcode_Omega(G,D,HC);
  // we can correct 1 error and the genus is 1, thus F must have degree 2
  // and support disjoint from that of D
  intvec F=2;
  list SV=prepSV(G,D,F,HC);
  // now we produce 1 error on the zero-codeword :
  matrix y[1][8];
  y[1,3]=a;
  // and then we decode :
  print(decodeSV(y,SV));
  printlevel=plevel;
}
///////////////////////////////////////////////////////////////////////////////

proc sys_code (matrix C)
"USAGE:    sys_code(C); C is a matrix of constants
RETURN:   list L with:
   @format
   L[1] is the generator matrix in standard form of an equivalent code,
   L[2] is the parity check matrix in standard form of such code,
   L[3] is an intvec which represents the needed permutation.
   @end format
NOTE:     Computes a systematic code which is equivalent to the given one.@*
          The input should be a matrix of numbers.@*
          The output has to be interpreted as follows: if the input was
          the generator matrix of an AG code then one should apply the
          permutation L[3] to the divisor D of rational points by means
          of @code{permute_L(D,L[3]);} before continuing to work with the
          code (for instance, if you want to use the systematic encoding
          together with a decoding algorithm).
KEYWORDS: linear code, systematic
SEE ALSO: permute_L, AGcode_Omega, prepSV
EXAMPLE:  example sys_code; shows an example
"
{
  // computes a systematic code which is equivalent to that given by C
  int i,j,k,l,h,r;
  int m=nrows(C);
  int n=ncols(C);
  int mr=m;
  matrix A[m][n]=C;
  poly c,p;
  list corners=list();
  if(m>n)
  {
    mr=n;
  }
  // first the matrix A will be reduced with elementary operations by rows
  for(i=1;i<=mr;i=i+1)
  {
    if((i+l)>n)
    {
      // the process is finished
      break;
    }
    // look for a non-zero element
    if(A[i,i+l]==0)
    {
      h=i;
      p=0;
      for (j=i+1;j<=m;j=j+1)
      {
        c=A[j,i+l];
        if (c!=0)
        {
          p=c;
          h=j;
          break;
        }
      }
      if (h!=i)
      {
        // permutation of rows i and h
        for (j=1;j<=n;j=j+1)
        {
          c=A[i,j];
          A[i,j]=A[h,j];
          A[h,j]=c;
        }
      }
      if(p==0)
      {
        // non-zero element not found in the current column
        l=l+1;
        continue;
      }
    }
    // non-zero element was found in "strategic" position
    corners[i]=i+l;
    // make zeros below that position
    for(j=i+1;j<=m;j=j+1)
    {
      c=A[j,i+l]/A[i,i+l];
      for(k=i+l+1;k<=n;k=k+1)
      {
        A[j,k]=A[j,k]-A[i,k]*c;
      }
      A[j,i+l]=0;
      A[j,i]=0;
    }
    // the rank is at least r
    // when the process stops the last r is actually the true rank of A=a
    r=i;
  }
  if (r<m)
  {
    ERROR("the given matrix does not have maximum rank");
  }
  // set the corners to the beginning and construct the permutation
  intvec PCols=1..n;
  for (j=1;j<=m;j=j+1)
  {
    if (corners[j]>j)
    {
      // interchange columns j and corners[j]
      for (i=1;i<=m;i=i+1)
      {
        c=A[i,j];
        A[i,j]=A[i,corners[j]];
        A[i,corners[j]]=c;
      }
      k=PCols[j];
      PCols[j]=PCols[corners[j]];
      PCols[corners[j]]=k;
    }
  }
  // convert the diagonal into ones
  for (i=1;i<=m;i=i+1)
  {
    for (j=i+1;j<=n;j=j+1)
    {
      A[i,j]=A[i,j]/A[i,i];
    }
    A[i,i]=1;
  }
  // complete a block with the identity matrix
  for (k=1;k<m;k=k+1)
  {
    for (i=k+1;i<=m;i=i+1)
    {
      for (j=i+1;j<=n;j=j+1)
      {
        A[k,j]=A[k,j]-A[i,j]*A[k,i];
      }
      A[k,i]=0;
    }
  }
  // compute a parity-check matrix in standard form
  matrix B=concat(-transpose(submat(A,1..m,m+1..n)),diag(1,n-m));
  list L=list();
  L[1]=A;
  L[2]=B;
  L[3]=PCols;
  return(L);
}
example
{
  "EXAMPLE:";  echo = 2;
  ring s=3,T,lp;
  matrix C[2][5]=0,1,0,1,1,0,1,0,0,1;
  print(C);
  list L=sys_code(C);
  L[3];
  // here is the generator matrix in standard form
  print(L[1]);
  // here is the control matrix in standard form
  print(L[2]);
  // we can check that both codes are dual to each other
  print(L[1]*transpose(L[2]));
}

////////////// this procedure should be moved to hnoether.lib /////////
static proc ratdevelop (poly chi)
{
  int ring_is_changed;
  int k1=res_deg();
  list L=hnexpansion(chi,"ess");
  if (typeof(L[1])=="ring") {
    def altring=basering;
    def HNring = L[1]; setring HNring;
    def HND = hne;
    kill hne;
    ring_is_changed=1;
  }
  else {
    def HND=L[1];
  }
  int k2=res_deg();
  if (k1<k2)
  {
    int N=size(HND);
    list CB=list();
    int i,j,k,d;
    intvec dgs;
    list auxL=list();
    intvec auxiv;
    for (i=1;i<=N;i++)
    {
      CB[i]=conj_bs(HND[i],k1);
      dgs[i]=size(CB[i]);
    }
    for (i=1;i<N;i++)
    {
      for (j=i+1;j<=N;j++)
      {
        if (dgs[i]==dgs[j])
        {
          d=dgs[i];
          for (k=1;k<=d;k++)
          {
            if (intersection(CB[i][1],CB[j][k])==-1)
            {
              HND=delete(HND,j);
              CB=delete(CB,j);
//              dgs=delete(dgs,j);
              execute("auxL="+string(dgs)+";");
              auxL=delete(auxL,j);
              execute("dgs="+string(auxL)+";");
              N=N-1;
              j=j-1;
              break;
            }
          }
        }
      }
    }
  }
  if (ring_is_changed==0) { return(HND); }
  else { export HND; setring altring; return(HNring); }
}



///////////////////////////////////////////////////////////////////////////////

/*
// ============================================================================
// *******       ADDITIONAL INFORMATION ABOUT THE LIBRARY              ********
// ============================================================================

A SINGULAR library for plane curves, Weierstrass semigroups and AG codes
Also available via http://wmatem.eis.uva.es/~ignfar/singular/

PREVIOUS WARNINGS :

   (1) The procedures will work only in positive characteristic
   (2) The base field must be prime (this may change in the future)
           This limitation is not too serious, since in coding theory
           the used curves are usually (if not always) defined over a
           prime field, and extensions are only considered for
           evaluating functions in a field with many points;
           by the same reason, auxiliary divisors are usually defined
           over the prime field,
           with the exception of that consisting of "rational points"
   (3) The curve must be absolutely irreducible (but it is not checked)
   (4) Only (algebraic projective) plane curves are considered

GENERAL CONCEPTS :

   (1) An affine point P is represented by a std of a prime ideal,
       and an intvec containing the position of the places above P
       in the list of Places; if the point is at infinity, the ideal is
       changed by a homogeneous irreducible polynomial in two variables
   (2) A place is represented by :
       a base point (list of homogeneous coordinates),
       a local equation for the curve at the base point,
       a Hamburger-Noether expansion of the corresponding branch,
       and a local parametrization (in "t") of such branch; everything is
       stored in a local ring "_[5][d][1]", d being the degree of the place,
       by means of lists "POINTS,LOC_EQS,BRANCHES,PARAMETRIZATIONS", and
       the degrees of the base points corresponding to the places in the
       ring "_[5][d][1]" are stored in an intvec "_[5][d][2]"
   (3) A divisor is represented by an intvec, where the integer at the
       position i means the coefficient of the i-th place in the divisor
   (4) Rational functions are represented by numerator/denominator
       in form of ideals with two homogeneous generators

OUTLINE/EXAMPLE OF THE USE OF THE LIBRARY :

Plane curves :

      (1.0) ring s=p,(x,y[,z]),lp;

            Notice that if you use 3 variables, then the equation
            of the curve is assumed to be a homogeneous polynomial.

      (1.1) list CURVE=Adj_div(equation);

            In CURVE[3] are listed all the (singular closed) places
            with their corresponding degrees; thus, you can now decide
            how many other points you want to compute with NSplaces.

      (1.2) CURVE=NSplaces(degrees,CURVE);

            "degrees" is an intvec with the degrees we want to compute

      (1.3) CURVE=extcurve(extension,CURVE);

            The rational places over the extension are ranged in
            the ring CURVE[1][5] with the following rules:

                (i)    all the representatives of the same closed point
                       are listed in consecutive positions;
                (ii)   if deg(P)<deg(Q), then the representatives of P
                       are listed before those of Q;
                (iii)  if two closed points P,Q have the same degree,
                       then the representatives of P are listed before
                       if P appears before in the list CURVE[3].

Rational functions :

      (2.0) def R=CURVE[1][2];
            setring R;
      (2.1) list LG=BrillNoether(intvec divisor,CURVE);
      (2.2) list WS=Weierstrass(int place,int bound,CURVE);

Algebraic Geometry codes :

      (3.0) def ER=CURVE[1][4];    // if extension>1; else use R instead
            setring ER;

            Now you have to decide the divisor G and the sequence of
            rational points D to use for constructing the codes;
            first notice that the syntax for G and D is different:

                (a) G is a divisor supported on the closed places of
                    CURVE[3], and you must say just the coefficient
                    of each such a point; for example, G=2,0,-1 means
                    2 times the place 1 minus 1 times the place 3.

                (b) D is a sequence of rational points (all different
                    and counted 1 time each one), whose data are read
                    from the lists inside CURVE[1][5] and now you must
                    say just the order how you range the chosen point;
                    for example, D=2,4,1 means that you choose the
                    rational places 1,2,4 and you range them as 2,4,1.

      (3.1) matrix C=AGcode_L(divisor,places,CURVE);

      (3.2) AGcode_Omega(divisor,places,CURVE);

            In the same way as for defining the codes, the auxiliary
            divisor F must have disjoint support to that of D, and
            its degree has to be given by a formula (see help prepSV).

      (3.3) list SV=prepSV(divisor,places,auxiliary_divisor,CURVE);

      (3.4) decodeSV(word,SV);

Special Issues :

      (I)  AG codes with systematic encoding :

              matrix C=AGcode_Omega(G,D,CURVE);
              list CODE=sys_code(G);
              C=CODE[1];               // generator matrix in standard form
              D=permute_L(D,CODE[3]);  // suitable permutation of coordinates
              list SV=prepSV(G,D,F,CURVE);
              SV[1]=CODE[2];           // parity-check matrix in standard form

      (II) Using the true minimum distance d for decoding :

              matrix C=AGcode_Omega(G,D,CURVE);
              int n=size(D);
              list SV=prepSV(G,D,F,CURVE);
              SV[n+3][2]=d;            // then use decodeSV(y,SV);


// ============================================================================
// ***   Some "macros" with typical examples of curves in Coding Theory    ****
// ============================================================================


proc Klein ()
{
  list KLEIN=Adj_div(x3y+y3+x);
  KLEIN=NSplaces(1..3,KLEIN);
  KLEIN=extcurve(3,KLEIN);
  dbprint(printlevel+1,"Klein quartic over F_8 successfully constructed");
  return(KLEIN);
}

proc Hermite (int m)
{
  int p=char(basering);
  int r=p^m;
  list HERMITE=Adj_div(y^r+y-x^(r+1));
  HERMITE=NSplaces(1..2*m,HERMITE);
  HERMITE=extcurve(2*m,HERMITE);
  dbprint(printlevel+1,"Hermitian curve over F_("+string(r)+"^2)
                          successfully constructed");
  return(HERMITE);
}

proc Suzuki ()
{
  list SUZUKI=Adj_div(x10+x3+y8+y);
  SUZUKI=NSplaces(1..3,SUZUKI);
  SUZUKI=extcurve(3,SUZUKI);
  dbprint(printlevel+1,"Suzuki curve over F_8 successfully constructed");
  return(SUZUKI);
}

// **** Other interesting examples :

// A hyperelliptic curve with 33 rational points over F_16

list CURVE=Adj_div(x5+y2+y);
CURVE=NSplaces(1..4,CURVE);
CURVE=extcurve(4,CURVE);

// A nice curve with 113 rational points over F_64

list CURVE=Adj_div(y9+y8+xy6+x2y3+y2+x3);
intvec ww=1,2,3,6;
CURVE=NSplaces(ww,CURVE);
CURVE=extcurve(6,CURVE);

*/