Index: Singular/LIB/brnoeth.lib
===================================================================
--- Singular/LIB/brnoeth.lib	(revision 489a498405fb2726f955024c28acd54e5b19b8b3)
+++ Singular/LIB/brnoeth.lib	(revision 489a498405fb2726f955024c28acd54e5b19b8b3)
@@ -0,0 +1,4903 @@
+version="$Id: brnoeth.lib,v 1.1 2000-12-11 13:50:19 Singular Exp $";
+info="
+LIBRARY:   brnoeth.lib
+
+AUTHORS:   Jose Ignacio Farran Martin,    ignfar@eis.uva.es
+           Christoph Lossen,              lossen@mathematik.uni-kl.de
+
+PURPOSE:   Brill-Noether algorithm, Weierstrass semigroups and AG codes
+
+SEE ALSO:  hnoether_lib, triang_lib
+
+KEYWORDS:  Weierstrass semigroup; Algebraic Geometry codes; Brill-Noether algorithm
+
+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
+
+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.
+";
+
+
+// =============================================================================
+
+
+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)
+"SYNTAX:    closed_points( ideal_expression )
+
+TYPE:       list
+
+PURPOSE:    compute the affine closed points in the zero-set of an ideal
+
+CREATE:     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)
+    print("? something was wrong; possibly non-reduced curve");
+    return();
+  }
+  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);
+  }
+  print("? error : basering must have 2 or 3 variables");
+  return();
+}
+
+
+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
+      print("? error : non-valid parameter");
+      return();
+    }
+  }
+  else
+  {
+    if (typeof(P[1])=="poly")
+    {
+      return(deg(P[1]));
+    }
+    else
+    {
+      print("? error : parameter must have a poly or ideal in the first component");
+      return();
+    }
+  }
+}
+
+
+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
+  {
+    print("? error : first argument must be an affine point");
+    return();
+  }
+}
+
+
+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
+    {
+      print("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
+      print("? error : non-valid parameter");
+      return();
+    }
+  }
+  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
+    {
+      print("? error : a point must have a poly or ideal in the first component");
+      return();
+    }
+  }
+  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 #)
+"SYNTAX:      Adj_div( poly_expression [, list_expression] )
+
+TYPE:       list
+
+CREATE:     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
+
+PURPOSE:    computes and stores the fundamental data of a plane curve 
+            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].
+
+NOTE:       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)
+  {
+    print("? error : base field not implemented");
+    return();
+  }
+  if (npars(basering)>0)
+  {
+    print("? error : base field not implemented");
+    return();
+  }
+  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)
+"SYNTAX:      NSplaces( int_expression, list_expression )
+
+TYPE:       list 
+
+PURPOSE:    @code{NSplaces(h,CURVE);} returns a list L with updated data of 
+            @code{CURVE} after computing all points up to degree H+@code{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 (representatives 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)
+"SYNTAX:    BrillNoether( intvec_expression, list_expression )
+
+TYPE:       list
+
+PURPOSE:    @code{BrillNoether(G,CURVE);} computes a vector basis of the linear
+            system L(G), where G is a given rational divisor over a non-singular 
+            curve.
+
+CREATE:     list of ideals (each of them with two homogeneous generators, which 
+            represent the nominator, resp. denominator, of a rational function).
+
+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 numerator/denominator or a poly
+  if ( typeof(F)=="ideal" && typeof(G)=="ideal" )
+  {
+    ideal FG;
+    FG[1]=F[1]*G[2]+F[2]*G[1];
+    FG[2]=F[2]*G[2];
+    return(simplifyRF(FG));
+  }
+  else
+  {
+    if (typeof(F)=="ideal")
+    {
+      ideal GG=polytoRF(G);
+      ideal FG;
+      FG[1]=F[1]*GG[2]+F[2]*GG[1];
+      FG[2]=F[2]*GG[2];
+      return(simplifyRF(FG));
+    }
+    else
+    {
+      if (typeof(G)=="ideal")
+      {
+        ideal FF=polytoRF(F);
+        ideal FG;
+        FG[1]=FF[1]*G[2]+FF[2]*G[1];
+        FG[2]=FF[2]*G[2];
+        return(simplifyRF(FG));
+      }
+      else
+      {
+        return(F+G);
+      }
+    }
+  }
+}
+
+
+static proc negRF (F)
+{
+  // returns -F as rational function
+  if (typeof(F)=="ideal")
+  {
+    ideal FF=F;
+    FF[1]=-F[1];
+    return(FF);
+  }
+  else
+  {
+    return(-F);
+  }
+}
+
+
+static proc escprodRF (l,F)
+{
+  // computes l*F as rational function
+  // l should be either a number or a 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)
+"SYNTAX:    Weierstrass( int_expression, int_expression, list_expression )
+
+TYPE:       list
+
+PURPOSE:    compute the Weierstrass semigroup up to a given bound and the 
+            associated rational functions
+
+CREATE:    @code{Weierstrass(i,m,CURVE);} returns a 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 @code{CURVE[1][2]}, 
+            where @code{CURVE} is the output of the procedure @code{NSplaces}. 
+            @code{i} represents the place @code{CURVE[3][i]}. 
+
+            WARNING: the place must be rational, i.e., necessarily 
+            @code{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)
+  {
+    print("? error : the given place is not defined over the prime field");
+    return();
+  }
+  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
+    {
+      print("? error : second argument must be non-negative");
+      return();
+    }
+  }
+  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)
+"SYNTAX:    permute_L( intvec_expression, intvec_expression )
+            permute_L( intvec_expression, list_expression )
+            permute_L( list_expression, intvec_expression )
+            permute_L( list_expression, list_expression )
+
+TYPE:       list
+
+PURPOSE:    apply a permutation to an intvec/list
+
+CREATE:     @code{permute_L(L,P);} returns a list obtained from @code{L} by 
+            applying the permutation given by @code{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
+    {
+      print("? error : the function is not well-defined at the given place");
+      return();
+    }
+  }
+}
+
+
+// ******** 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)
+"SYNTAX:    dual_code( matrix_expression )
+
+TYPE:       matrix
+
+PURPOSE:    compute a generator matrix of the dual code
+ 
+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)
+"SYNTAX:    AGcode_L( intvec_expression, intvec_expression, list_expression )
+
+TYPE:       matrix
+
+PURPOSE:    @code{AGcode_L(G,D,EC)} computes a generator matrix for the 
+            evaluation AG code defined by the divisors G and D. 
+
+NOTE:       The procedure must be called within the ring @code{EC[1][4]}, 
+            where @code{EC} is the output of @code{extcurve(d)} (or within 
+            the ring @code{EC[1][2]} if d=1). @*
+            The entry i in the intvec @code{D} refers to the i-th rational 
+            place in @code{EC[1][5]} (i.e., to @code{POINTS[i]}, etc., see 
+            @ref{extcurve}).@*
+            The intvec @code{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)
+  {     
+    print("? 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_8
+  print(C);
+  printlevel=plevel;
+}
+
+
+// =============================================================================
+
+
+proc AGcode_Omega (intvec G,intvec D,list EC)
+"SYNTAX:    AGcode_Omega( intvec_expression, intvec_expression, list_expression )
+
+TYPE:       matrix
+
+PURPOSE:    @code{AGcode_Omega(G,D,EC)} computes a generator matrix for the 
+            residual AG code defined by the divisors G and D. 
+
+NOTE:       The procedure must be called within the ring @code{EC[1][4]}, 
+            where @code{EC} is the output of @code{extcurve(d)} (or within 
+            the ring @code{EC[1][2]} if d=1). @*
+            The entry i in the intvec @code{D} refers to the i-th rational place 
+            in @code{EC[1][5]} (i.e., to @code{POINTS[i]}, etc., see 
+            @ref{extcurve}).@*
+            The intvec @code{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_8
+  print(C);
+  printlevel=plevel;
+}
+
+
+// =============================================================================
+// ********   auxiliary procedure to define AG codes over extensions    ********
+// =============================================================================
+
+
+proc extcurve (int d,list CURVE)
+"SYNTAX:    extcurve( int_expression, list_expression )
+
+TYPE:       list
+
+PURPOSE:    create local rings over field extensions in order to have
+            sufficiently many rational places (as needed for the construction
+            of AG codes).
+
+CREATE:     If @code{d}>1 then @code{extcurve(d,CURVE);} creates a list L which
+            is the update of the list @code{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 @code{d}. 
+
+            If @code{d}<2 then @code{extcurve(d,CURVE);} creates a list L which
+            is the update of the list @code{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 @code{d}) are 
+            created (see @ref{Adj_div}):
+   @format 
+   lists POINTS, LOC_EQS, BRANCHES, PARAMETRIZATIONS.
+   @end format
+
+NOTE:       The list_expression @code{CURVE} should be the output of 
+            @code{NSplaces} and has to contain (at least) all places up to 
+            degree @code{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
+    {
+      print("? 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)
+"SYNTAX:    prepSV( intvec_expression, intvec_expression, intvec_expression, list_expression ) 
+
+TYPE:       list 
+
+PURPOSE:    preprocessing for the basic (Skorobogatov-Vladut) decoding 
+            algorithm 
+
+CREATE:     @code{prepSV(G,D,F,EC);} returns a 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:       The procedure must be called within the ring @code{EC[1][4]}, 
+            where @code{EC} is the output of @code{extcurve(d)} (or within 
+            the ring @code{EC[1][2]} if d=1). @*
+            The intvec @code{G} and @code{F} represent rational divisors (see 
+            @ref{BrillNoether} for more details).@*
+            The intvec @code{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 
+            @code{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)
+  {     
+    print("? 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)
+  {
+    print("? error : the correction capacity of the basic algorithm is zero");
+    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_8
+  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)
+"SYNTAX:    decodeSV( matrix_expression, list_expression )
+
+TYPE:       matrix
+
+PURPOSE:    basic (Skorobogatov-Vladut) decoding algorithm
+
+CREATE:     a codeword (row-matrix) if possible, resp. the 0-matrix (of size 1) if
+            decoding is impossible
+
+NOTE:       The list_expression should be the output K of the procedure 
+            @code{prepSV}.@* 
+            The matrix_expression should be a (1xn)-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)
+  {
+    print("? : no error-locator found");
+    print("? : too many errors occur");
+    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 that 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
+      {
+        print("? : non-acceptable decoding");
+      }
+    }
+    else
+    {
+      print("? : no solution found");
+    }
+  }
+  else
+  {
+    print("? : non-unique solution");
+  }
+  option(set,opgt);
+  print("? : too many errors occur");
+  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_8
+  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)
+"SYNTAX:    sys_code( matrix_expression )
+
+TYPE:       list
+
+PURPOSE:    computes a systematic code which is equivalent to the given one
+
+CREATE:     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:       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)
+  {
+    print("? error : the given matrix does not have maximum rank");
+    return(list());
+  }
+  // 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);
+
+
+*/
+
+
+;
