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

version="$Id: brnoeth.lib,v 1.4 2000-12-15 12:22:56 obachman Exp $";
info="
LIBRARY:   brnoeth.lib  PROCEDURES FOR THE BRILL-NOETHER ALGORITHM
           
AUTHORS:   Jose Ignacio Farran Martin,    ignfar@eis.uva.es
           Christoph Lossen,              lossen@mathematik.uni-kl.de

SEE ALSO:  hnoether_lib, triang_lib

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

OVERVIEW:  Implementation of the Brill-Noether algorithm for solving the
           Riemann-Roch problem and applications in 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 charateristic.@* 
           You can get more information about the library just
           by reading the end of the file brnoeth.lib.

MAIN PROCEDURES:
 Adj_div(f);            computes the conductor of a curve 
 NSplaces(h,A);         computes non-singular places up to given degree
 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 given divisors 
                        G and D
 AGcode_Omega(G,D,E);   computes the residual AG code with given 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

AUXILIARY PROCEDURES:
 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
";


// ===========================================================================


LIB "matrix.lib";
LIB "triang.lib";
LIB "hnoether.lib";
// -> LIB "general.lib","ring.lib";
// maybe useful : LIB "linalg.lib","primdec.lib","normal.lib";


// ===========================================================================


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


proc closed_points (ideal I)
"USAGE:     closed_points(I), where I is an ideal

RETURN:     list of prime ideals (std), 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.@*
            The option(redSB) is convenient to be set 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;
  }
  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 of projective)
CREATE:     poly CHI in both ring aff_r=p,(x,y),lp and ring 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 aff_r,Proj_R and 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 )
  {
    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);
        }
      }
    }
  }
  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) poly => 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 poly 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 subset (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=size(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(subset(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 (subset(I,J)==1)
  {
    if (subset(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 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);
      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 Fq;
        ideal facs=factorize(x^(q-1)-1,1);
        for (i=1;i<=q-1;i=i+1)
        {
          Fq[i]=number(subst(facs[i],x,0));
        }
        for (i=1;i<=q-1;i=i+1)
        {
          b=Fq[i];
          if (subst(mpol,x,b)==0)
          {
            answer=string(b);
            break;
          }
        }
      }
      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+";");
    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 only int 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="+sdatum+";");
    execute("poly prel=b-"+rel+";");
    ideal I=pdatum,prel;
    I=eliminate(I,a);
    poly newdatum=I[1];   // hopefully it was done correctly and size(I)=1 !!!!
    newdatum=subst(newdatum,b,a);
    string snewdatum=string(newdatum);
    setring base_r;
    execute("poly newdatum="+snewdatum+";");
    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 poly 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 L@HNE=essdevelop(CHI);
  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];
            }
          }
        }
      }
      if (size(update_CURVE[5])>0)
      {
        if (typeof(update_CURVE[5][1])=="list")
        {
          check=1;
        }
      }      
      if (check==0)
      {
        ring S(1)=char(basering),(x,y,t),ls;
        intvec dgs_points(1);
        list BRANCHES=list();
        list POINTS=list();
        list LOC_EQS=list();
        list PARAMETRIZATIONS=list();
        export(BRANCHES);
        export(POINTS);
        export(LOC_EQS);
        export(PARAMETRIZATIONS);
      }
      else
      {
        def S(1)=update_CURVE[5][1][1];
        def dgs_points(1)=update_CURVE[5][1][2];
      }
      setring S(1);
      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));
        }
        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 reduce the 
      // field extensions for each one, and finally check if the minpoly 
      // fetchs with S(i) or not
      // if one of the branches is rational, we may trust that is 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];
        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
        {
          def S(ext_0)=update_CURVE[5][ext_0][1];
          def dgs_points(ext_0)=update_CURVE[5][ext_0][2];
        }
        setring S(ext_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 branche
          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;
        kill(iw);
        if (sing==2)
        {
          Conductor[N]=local_conductor(iw[2],S(ext_0));
        }
        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];
            }
          }
        }
      }
      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
      {
        def S(ext2)=update_CURVE[5][ext2][1];
        def dgs_points(ext2)=update_CURVE[5][ext2][2];
      }
      setring S(ext2);
      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;
        kill(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));
        }
        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];
        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
        {
          def S(ext_0)=update_CURVE[5][ext_0][1];
          def dgs_points(ext_0)=update_CURVE[5][ext_0][2];
        }
        setring S(ext_0);
        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));
        }
        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);
  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 [,#] ), where f is a poly, [ # a list ]

RETURN:      list L with the computed data:
  @format 
  L[1] is a list of rings: L[1][1]=aff_r (affine), L[1][2]=Proj_R (projective),
  L[2] is an intvec with 2 entries (degree, genus),
  L[3] is a list of intvec (closed places),
  L[4] is an intvec (conductor),
  L[5] is a list of lists: 
         L[5][d][1] is a (local) ring over an extension of degree d, 
         L[5][d][2] is 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.

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");
  }
  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);
  C;
  // affine ring
  def aff_R=C[1][1];
  setring aff_R;
  // the affine equation of the curve
  CHI;
  // the ideal of affine singular locus
  Aff_SLocus;
  // the list of affine singular points
  Aff_SPoints;
  // the list of singular/non-singular points at infinity
  Inf_Points;
  // the projective ring
  def proj_R=C[1][2];
  setring proj_R;
  // the projective equation of the curve
  CHI;
  // the degree of the curve :
  C[2][1];
  // the genus of the curve :
  C[2][2];
  // the adjunction divisor :
  C[4];
  // the list of computed places
  C[3];
  // the list of local rings and degrees of base points
  C[5];
  // we look at some places
  def S(1)=C[5][1][1];
  setring S(1);
  POINTS;
  LOC_EQS;
  PARAMETRIZATIONS;
  BRANCHES;
  printlevel=plevel;
}


// ============================================================================


proc NSplaces (int h,list CURVE)
"USAGE:    NSplaces( h, CURVE ), where h is an integer and CURVE is a list

RETURN:    list L with updated data of CURVE after computing all 
           points up to degree H+h (H the maximum degree of the previously 
           computed places: @*
   @format 
   in affine ring L[1][1]: 
       lists Aff_Points(d) with affine non-singular 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 (represent. 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.@*

SEE ALSO:  closed_points, Adj_div
EXAMPLE:   example NSplaces; shows an example
"
{
  // computes all the non-singular closed places with degree up to a certain 
  // bound; this bound is the maximum degree of an existing singular place or 
  // non-singular place at infinity plus an increment h>=0 which is given as 
  // input
  // creates lists of points and the corresponding places
  // list CURVE must be the output of the procedure "Adj_div"
  // warning : if h<0 then it will be replaced by h=0
  intvec opgt=option(get);
  option(redSB);
  def Base_R=basering;
  def aff_r=CURVE[1][1];
  int M=size(CURVE[5]);
  if (h>0)
  {
    M=M+h;
  }
  list update_CURVE=CURVE;
  int i,j,s;
  setring aff_r;
  intvec pP=1,0,0;
  for (i=1;i<=M;i=i+1)
  {
    dbprint(printlevel+1,"Computing non-singular affine places of degree "
                           +string(i)+" ... ");
    list Aff_Points(i)=closed_points_deg(CHI,i,Aff_SLocus);
    s=size(Aff_Points(i));
    if (s>0)
    {
      export(Aff_Points(i));
      pP[2]=i;
      for (j=1;j<=s;j=j+1)
      {
        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);
  // create places up to degree 1+3
  list L=NSplaces(3,C);
  L;
  // for example, here is the list with the affine non-singular 
  // points of degree 4 :
  def aff_r=L[1][1];
  setring aff_r;
  Aff_Points(4);
  // for example, we check the places of degree 4 :
  def S(4)=L[5][4][1];
  setring S(4);
  // and for example, the base points of such places :
  POINTS;
  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:      suplement(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=rowred(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=rowred(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[1]=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=rowred(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 && ordx>-1 ) || ( inum>ordy && 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];
        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 && ordx>-1 ) || ( inum>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];
        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 L@HNE=essdevelop(LOC_EQ)[1];
  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);
    }
  }
  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;
  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);
  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 poly 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), where G is an intvec and CURVE is a list

RETURN:     list of ideals (each of them with two homogeneous generators, 
            which represent the nominator, resp. denominator, of a rational 
            function).@*
            The corresponding rational functions form a vector 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(3,C);
  // Places C[3][1] and C[3][2] are rational,
  // so that we define a divisor of degree 8
  intvec G=4,4;
  def R=C[1][2];
  setring R;
  list LG=BrillNoether(G,C);
  // here is the vector basis :
  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 poly (or number) into a "rational function" of type "ideal"
  // warning : it must be called inside "R" and poly should 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(auxp,F)[1]));
}


static proc sumRF (F,G)
{
  // sum of two "rational functions" F,G given by either a pair 
  // nominator/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 poly 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), where F is an ideal, SS is a ring and m is 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 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 ), where i,m are integers and CURVE a list

RETURN:     list WS of two lists:
  @format 
  WS[1] is a list of integers (the Weierstrass semigroup of the curve 
                                                      at the place i up to m)
  WS[2] is a 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]. 

            WARNING: the place must be rational, i.e., necessarily 
            CURVE[3][P][1]=1.
            
            Rational functions are represented by nominator/denominator
            in form of ideals with two homogeneous generators.

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(3,C);
  def R=C[1][2];
  setring R;
  // Place C[3][1] has degree 1 (i.e it is rational);
  list WS=Weierstrass(1,10,C);
  // the first part of the list is the Weierstrass semigroup up to 10 :
  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 ), where L,P are 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), where G is 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 )
    {
      def auxR=EC[5][i][1];
      setring auxR;
      auxIV[i]=size(POINTS);
      setring BRing;
      kill(auxR);
    }
    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) 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);
              }
            }
          }
        }
      }
    }
  }
  return(1);
}




// ============================================================================


proc AGcode_L (intvec G,intvec D,list EC)
"USAGE:      AGcode_L( G, D, EC ), where G,D are intvec and EC is 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 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,HC);
  HC=extcurve(2,HC);
  def ER=HC[1][4];
  setring ER;
  intvec G=5;
  intvec D=2..9;
  // we already have a rational divisor G and 8 more points 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 ), where G,D are intvec and EC is 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,HC);
  HC=extcurve(2,HC);
  def ER=HC[1][4];
  setring ER;
  intvec G=5;
  intvec D=2..9;
  // we already have a rational divisor G and 8 more points 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 ), where d is an integer and CURVE is 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 places over the base field).
   @end format
            In both cases, in the ring L[1][5] lists with the data for all the 
            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 has 
            to contain (at least) all places up to degree d. @*
            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,"");
    return(ext_CURVE);
  }
  else
  {
    if (size(CURVE[5])>=d)
    {
      int i,j,k;
      for (i=1;i<=d;i=i+1)
      {
        if (typeof(CURVE[5][i])=="list")
        {
          def S(i)=CURVE[5][i][1];
        }
        else
        {
          ring S(i)=char(basering),t,ls;
          list POINTS=list();
          setring BR;
        }
      }
      setring S(d);
      string smin=string(minpoly);
      setring BR;
      ring RatPl=(char(basering),a),(x,y,t),ls;
      execute("minpoly="+smin+";");
      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);
      int s=size(POINTS);
      int counter=0;
      int piv=0;
      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;
        }
      }
      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 )
        {
          olda=subfield(S(i));
          setring S(i);
          s=size(POINTS);
          setring RatPl;
          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(paux);
      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;
      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,"");
      return(ext_CURVE);
    }
    else
    {
      ERROR("you must compute first all the places up to degree "+string(d));
      return();
    }
  }
}
example
{
  "EXAMPLE:";  echo = 2;
  int plevel=printlevel;
  printlevel=-1;
  ring s=2,(x,y),lp;
  list C=Adj_div(x5+y2+y);
  C=NSplaces(3,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);
  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 ), where G,D,F are intvec and EC is 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 within 
            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 to 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;
  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,HC);
  HC=extcurve(2,HC);
  def ER=HC[1][4];
  setring ER;
  intvec G=5;
  intvec D=2..9;
  // we already have a rational divisor G and 8 more points over F_4;
  // let us construct the corresponding 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 to 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[11][1];
  epsilon;
  printlevel=plevel; 
}


// ============================================================================


proc decodeSV (matrix y,list K)
"USAGE:     decodeSV( y, K ), where y is a row-matrix and K is 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,HC);
  HC=extcurve(2,HC);
  def ER=HC[1][4];
  setring ER;
  intvec G=5;
  intvec D=2..9;
  // we already have a rational divisor G and 8 more points over F_4;
  // let us construct the corresponding 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 to 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) where 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;j<=n;j=j+1)
    {
      A[i,j]=A[i,j]/A[i,i];
    }
  }
  // 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;j<=n;j=j+1)
      {
        A[k,j]=A[k,j]-A[i,j]*A[k,i];
      }
    }
  }
  // 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 each other
  print(L[1]*transpose(L[2]));
}


/* 


// ============================================================================
// *******       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(range,CURVE);

            See help NSplaces tp know the meaning of "range"; 
            for instance, if you have that the singular places 
            have degree 2 at most and you want all the places 
            up to degree 5, you must write range=3. 

      (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(2,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(2*m-1,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(2,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(3,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);
CURVE=NSplaces(4,CURVE);
CURVE=extcurve(6,CURVE);


*/


;