source: git/Singular/LIB/brnoeth.lib @ 50cbdc

spielwiese
Last change on this file since 50cbdc was 50cbdc, checked in by Hans Schönemann <hannes@…>, 23 years ago
*hannes: merge-2-0-2 git-svn-id: file:///usr/local/Singular/svn/trunk@5619 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 136.5 KB
Line 
1version="$Id: brnoeth.lib,v 1.12 2001-08-27 14:47:46 Singular Exp $";
2category="Miscellaneous";
3info="
4LIBRARY:   brnoeth.lib  Brill-Noether Algorithm, Weierstrass-SG and AG-codes
5AUTHORS:   Jose Ignacio Farran Martin, ignfar@eis.uva.es
6           Christoph Lossen,           lossen@mathematik.uni-kl.de
7
8OVERVIEW:
9 Implementation of the Brill-Noether algorithm for solving the
10 Riemann-Roch problem and applications in Algebraic Geometry codes.
11 The computation of Weierstrass semigroups is also implemented.@*
12 The procedures are intended only for plane (singular) curves defined over
13 a prime field of positive charateristic.@*
14 For more information about the library see the end of the file brnoeth.lib.
15
16MAIN PROCEDURES:
17 Adj_div(f);            computes the conductor of a curve
18 NSplaces(h,A);         computes non-singular places up to given degree
19 BrillNoether(D,C);     computes a vector space basis of the linear system L(D)
20 Weierstrass(P,m,C);    computes the Weierstrass semigroup of C at P up to m
21 extcurve(d,C);         extends the curve C to an extension of degree d
22 AGcode_L(G,D,E);       computes the evaluation AG code with divisors G and D
23 AGcode_Omega(G,D,E);   computes the residual AG code with divisors G and D
24 prepSV(G,D,F,E);       preprocessing for the basic decoding algorithm
25 decodeSV(y,K);         decoding of a word with the basic decoding algorithm
26
27AUXILIARY PROCEDURES:
28 closed_points(I);      computes the zero-set of a zero-dim. ideal in 2 vars
29 dual_code(C);          computes the dual code
30 sys_code(C);           computes an equivalent systematic code
31 permute_L(L,P);        applies a permutation to a list
32
33SEE ALSO:  hnoether_lib, triang_lib
34
35KEYWORDS:  Weierstrass semigroup; Algebraic Geometry codes;
36           Brill-Noether algorithm
37";
38
39LIB "matrix.lib";
40LIB "triang.lib";
41LIB "hnoether.lib";
42LIB "inout.lib";
43
44///////////////////////////////////////////////////////////////////////////////
45
46// **********************************************************
47// * POINTS, PLACES AND DIVISORS OF (SINGULAR) PLANE CURVES *
48// **********************************************************
49
50proc closed_points (ideal I)
51"USAGE:     closed_points(I); I an ideal
52RETURN:     list of prime ideals (each a Groebner basis), corresponding to
53            the (distinct affine closed) points of V(I)
54NOTE:       The ideal must have dimension 0, the basering must have 2
55            variables, the ordering must be lp, and the base field must
56            be finite and prime.@*
57            It might be convenient to set the option(redSB) in advance.
58SEE ALSO:   triang_lib
59EXAMPLE:    example closed_points; shows an example
60"
61{
62  ideal II=std(I);
63  if (II==1)
64  {
65    return(list());
66  }
67  list TL=triangMH(II);
68  int s=size(TL);
69  list L=list();
70  int i,j,k;
71  ideal Facts;
72  poly G2;
73  for (i=1;i<=s;i=i+1)
74  {
75    Facts=factorize(TL[i][1],1);
76    k=size(Facts);
77    G2=TL[i][2];
78    for (j=1;j<=k;j=j+1)
79    {
80      L=L+pd2(Facts[j],G2);
81    }
82  }
83  // eliminate possible repetitions
84  s=size(L);
85  list LP=list();
86  LP[1]=std(L[1]);
87  int counter=1;
88  for (i=2;i<=s;i=i+1)
89  {
90    if (isPinL(L[i],LP)==0)
91    {
92      counter=counter+1;
93      LP[counter]=std(L[i]);
94    }
95  }
96  return(LP);
97}
98example
99{
100  "EXAMPLE:";  echo = 2;
101  ring s=2,(x,y),lp;
102  // this is just the affine plane over F_4 :
103  ideal I=x4+x,y4+y;
104  list L=closed_points(I);
105  // and here you have all the points :
106  L;
107}
108///////////////////////////////////////////////////////////////////////////////
109
110static proc pd2 (poly g1,poly g2)
111{
112  // If g1,g2 is a std. resp. lex. in (x,y) then the procedure
113  //    factorizes g2 in the "extension given by g1"
114  //    (then g1 must be irreducible) and returns a list of
115  //    ideals with always g1 as first component and the
116  //    distinct factors of g2 as second components
117  list L=list();
118  ideal J=g1;
119  int i,s;
120  if (deg(g1)==1)
121  {
122    poly A=-subst(g1,var(2),0);
123    poly B=subst(g2,var(2),A);
124    ideal facts=factorize(B,1);
125    s=size(facts);
126    for (i=1;i<=s;i=i+1)
127    {
128      J[2]=facts[i];
129      L[i]=J;
130    }
131  }
132  else
133  {
134    def BR=basering;
135    poly A=g1;
136    poly B=g2;
137    ring raux1=char(basering),(x,y,a),lp;
138    poly G;
139    ring raux2=(char(basering),a),(x,y),lp;
140    map psi=BR,x,a;
141    minpoly=number(psi(A));
142    poly f=psi(B);
143    ideal facts=factorize(f,1);
144    s=size(facts);
145    poly g;
146    string sg;
147    for (i=1;i<=s;i=i+1)
148    {
149      g=facts[i];
150      sg=string(g);
151      setring raux1;
152      execute("G="+sg+";");
153      G=subst(G,a,y);
154      setring BR;
155      map ppssii=raux1,var(1),var(2),0;
156      J[2]=ppssii(G);
157      L[i]=J;
158      kill(ppssii);
159      setring raux2;
160    }
161    setring BR;
162    kill(raux1,raux2);
163  }
164  return(L);
165}
166///////////////////////////////////////////////////////////////////////////////
167static proc isPinL (ideal P,list L)
168{
169  // checks if a (plane) point P is in a list of (plane) points L
170  //    by just comparing generators
171  // it works only if all (prime) ideals are given in a "canonical way",
172  // namely:
173  //    the first generator is monic and irreducible,
174  //        and depends only on the second variable,
175  //    and the second one is monic in the first variable
176  //        and irreducible over the field extension determined by
177  //        the second variable and the first generator as minpoly
178  int s=size(L);
179  int i;
180  for (i=1;i<=s;i=i+1)
181  {
182    if ( P[1]==L[i][1] && P[2]==L[i][2] )
183    {
184      return(1);
185    }
186  }
187  return(0);
188}
189///////////////////////////////////////////////////////////////////////////////
190static proc s_locus (poly f)
191{
192  // computes : ideal of affine singular locus
193  // the equation f must be affine
194  // warning : if there is an error message then the output is "none"
195  // option(redSB) is convenient to be set in advance
196  ideal I=f,jacob(f);
197  I=std(I);
198  if (dim(I)>0)
199  {
200    // dimension check (has to be 0)
201    ERROR("something was wrong; possibly non-reduced curve");
202  }
203  else
204  {
205    return(I);
206  }
207}
208///////////////////////////////////////////////////////////////////////////////
209static proc curve (poly f)
210"USAGE:      curve(f), where f is a polynomial (affine or projective)
211CREATE:     poly CHI in both rings aff_r=p,(x,y),lp and Proj_R=p,(x,y,z),lp
212            also ideal (std) Aff_SLocus of affine singular locus in the ring
213            aff_r
214RETURN:     list (size 3) with two rings aff_r,Proj_R and an integer deg(f)
215NOTE:       f must be absolutely irreducible, but this is not checked
216            it is not implemented yet for extensions of prime fields
217"
218{
219  def base_r=basering;
220  ring aff_r=char(basering),(x,y),lp;
221  ring Proj_R=char(basering),(x,y,z),lp;
222  setring base_r;
223  int degX=deg(f);
224  if (nvars(basering)==2)
225  {
226    setring aff_r;
227    map embpol=base_r,x,y;
228    poly CHI=embpol(f);
229    export(CHI);
230    kill(embpol);
231    ideal Aff_SLocus=s_locus(CHI);
232    export(Aff_SLocus);
233    setring Proj_R;
234    poly CHI=homog(imap(aff_r,CHI),z);
235    export(CHI);
236    setring base_r;
237    list L=list();
238    L[1]=aff_r;
239    L[2]=Proj_R;
240    L[3]=degX;
241    kill(aff_r,Proj_R);
242    return(L);
243  }
244  if (nvars(basering)==3)
245  {
246    setring Proj_R;
247    map embpol=base_r,x,y,z;
248    poly CHI=embpol(f);
249    export(CHI);
250    kill(embpol);
251    string s=string(subst(CHI,z,1));
252    setring aff_r;
253    execute("poly CHI="+s+";");
254    export(CHI);
255    ideal Aff_SLocus=s_locus(CHI);
256    export(Aff_SLocus);
257    setring base_r;
258    list L=list();
259    L[1]=aff_r;
260    L[2]=Proj_R;
261    L[3]=degX;
262    kill(aff_r,Proj_R);
263    return(L);
264  }
265  ERROR("basering must have 2 or 3 variables");
266}
267///////////////////////////////////////////////////////////////////////////////
268static proc Aff_SL (ideal ISL)
269{
270  // computes : affine singular (closed) points as a list of lists of
271  //     prime ideals and intvec (for storing the places over each point)
272  // the ideal ISL=s_locus(CHI) is assumed to be computed in advance for
273  //     a plane curve CHI, and it must be given by a standard basis
274  // for our purpose the function must called with the "global" ideal
275  //     "Aff_SLocus"
276  list SL=list();
277  ideal I=ISL;
278  if ( I != 1 )
279  {
280    list L=list();
281    ideal aux;
282    intvec iv;
283    int i,s;
284    L=closed_points(I);
285    s=size(L);
286    for (i=1;i<=s;i=i+1)
287    {
288      aux=std(L[i]);
289      SL[i]=list(aux,iv);
290    }
291  }
292  return(SL);
293}
294///////////////////////////////////////////////////////////////////////////////
295static proc inf_P (poly f)
296{
297  // computes : all (closed) points at infinity as homogeneous polynomials
298  // output : two lists with respectively singular and non-singular points
299  intvec iv;
300  def base_r=basering;
301  ring r_auxz=char(basering),(x,y,z),lp;
302  poly f=imap(base_r,f);
303  poly F=homog(f,z);                    // equation of projective curve
304  poly f_inf=subst(F,z,0);
305  setring base_r;
306  poly f_inf=imap(r_auxz,f_inf);
307  ideal I=factorize(f_inf,1);           // points at infinity as homogeneous
308                                        // polynomials
309  int s=size(I);
310  int i;
311  list IP_S=list();                     // for singular points at infinity
312  list IP_NS=list();                    // for non-singular points at infinity
313  int counter_S;
314  int counter_NS;
315  poly aux;
316  for (i=1;i<=s;i=i+1)
317  {
318    aux=subst(I[i],y,1);
319    if (aux==1)
320    {
321      // the point is (1:0:0)
322      setring r_auxz;
323      poly f_yz=subst(F,x,1);
324      if ( subst(subst(diff(f_yz,y),y,0),z,0)==0 &&
325           subst(subst(diff(f_yz,z),y,0),z,0)==0 )
326      {
327        // the point is singular
328        counter_S=counter_S+1;
329        kill(f_yz);
330        setring base_r;
331        IP_S[counter_S]=list(I[i],iv);
332      }
333      else
334      {
335        // the point is non-singular
336        counter_NS=counter_NS+1;
337        kill(f_yz);
338        setring base_r;
339        IP_NS[counter_NS]=list(I[i],iv);
340      }
341    }
342    else
343    {
344      // the point is (a:1:0) | a is root of aux
345      if (deg(aux)==1)
346      {
347        // the point is rational and no field extension is needed
348        setring r_auxz;
349        poly f_xz=subst(F,y,1);
350        poly aux=imap(base_r,aux);
351        number A=-number(subst(aux,x,0));
352        map phi=r_auxz,x+A,0,z;
353        poly f_origin=phi(f_xz);
354        if ( subst(subst(diff(f_origin,x),x,0),z,0)==0 &&
355             subst(subst(diff(f_origin,z),x,0),z,0)==0 )
356        {
357          // the point is singular
358          counter_S=counter_S+1;
359          kill(f_xz,aux,A,phi,f_origin);
360          setring base_r;
361          IP_S[counter_S]=list(I[i],iv);
362        }
363        else
364        {
365          // the point is non-singular
366          counter_NS=counter_NS+1;
367          kill(f_xz,aux,A,phi,f_origin);
368          setring base_r;
369          IP_NS[counter_NS]=list(I[i],iv);
370        }
371      }
372      else
373      {
374        // the point is non-rational and a field extension with minpoly=aux
375        // is needed
376        ring r_ext=(char(basering),a),(x,y,z),lp;
377        poly F=imap(r_auxz,F);
378        poly f_xz=subst(F,y,1);
379        poly aux=imap(base_r,aux);
380        minpoly=number(subst(aux,x,a));
381        map phi=r_ext,x+a,0,z;
382        poly f_origin=phi(f_xz);
383        if ( subst(subst(diff(f_origin,x),x,0),z,0)==0 &&
384             subst(subst(diff(f_origin,z),x,0),z,0)==0 )
385        {
386          // the point is singular
387          counter_S=counter_S+1;
388          setring base_r;
389          kill(r_ext);
390          IP_S[counter_S]=list(I[i],iv);
391        }
392        else
393        {
394          // the point is non-singular
395          counter_NS=counter_NS+1;
396          setring base_r;
397          kill(r_ext);
398          IP_NS[counter_NS]=list(I[i],iv);
399        }
400      }
401    }
402  }
403  kill(r_auxz);
404  return(list(IP_S,IP_NS));
405}
406///////////////////////////////////////////////////////////////////////////////
407static proc closed_points_ext (poly f,int d,ideal SL)
408{
409  // computes : (closed affine non-singular) points over an extension of
410  //            degree d
411  // remark(1) : singular points are supposed to be listed appart
412  // remark(2) : std SL=s_locus(f) is supposed to be computed in advance
413  // remark(3) : ideal SL is used to remove those points which are singular
414  // output : list of list of prime ideals with an intvec for storing the
415  //          places
416  int Q=char(basering)^d;             // cardinality of the extension field
417  ideal I=f,x^Q-x,y^Q-y;              // ideal of the searched points
418  I=std(I);
419  if (I==1)
420  {
421    return(list());
422  }
423  list LP=list();
424  int m=size(SL);
425  list L=list();
426  ideal aux;
427  intvec iv;
428  int i,s,j,counter;
429  L=closed_points(I);
430  s=size(L);
431  for (i=1;i<=s;i=i+1)
432  {
433    aux=std(L[i]);
434    for (j=1;j<=m;j=j+1)
435    {
436      // check if singular i.e. if SL is contained in aux
437      if ( NF(SL[j],aux) != 0 )
438      {
439        counter=counter+1;
440        LP[counter]=list(aux,iv);
441        break;
442      }
443    }
444  }
445  return(LP);
446}
447///////////////////////////////////////////////////////////////////////////////
448static proc degree_P (list P)
449"USAGE:      degree_P(P), where P is either a polynomial or an ideal
450RETURN:     integer with the degree of the closed point given by P
451SEE ALSO:   closed_points
452NOTE:       If P is a (homogeneous irreducible) polynomial the point is at
453            infinity, and if P is a (prime) ideal the points is affine, and
454            the ideal must be given by 2 generators: the first one irreducible
455            and depending only on y, and the second one irreducible over the
456            extension given by y with the first generator as minimal polynomial
457"
458{
459  // computes : the degree of a given point
460  // remark(1) : if the input is (irreducible homogeneous) poly => the point
461  //             is at infinity
462  // remark(2) : it the input is (std. resp. lp. prime) ideal => the point is
463  //             affine
464  if (typeof(P[1])=="ideal")
465  {
466    if (size(P[1])==2)
467    {
468      int d=deg(P[1][1]);
469      poly aux=subst(P[1][2],y,1);
470      d=d*deg(aux);
471      return(d);
472    }
473    else
474    {
475      // this should not happen in principle
476      ERROR("non-valid parameter");
477    }
478  }
479  else
480  {
481    if (typeof(P[1])=="poly")
482    {
483      return(deg(P[1]));
484    }
485    else
486    {
487      ERROR("parameter must have a poly or ideal in the first component");
488    }
489  }
490}
491///////////////////////////////////////////////////////////////////////////////
492static proc closed_points_deg (poly f,int d,ideal SL)
493{
494  // computes : (closed affine non-singular) points of degree d
495  // remark(1) : singular points are supposed to be listed appart
496  // remark(2) : std SL=s_locus(f) is supposed to be computed in advance
497  list L=closed_points_ext(f,d,SL);
498  int s=size(L);
499  int i,counter;
500  list LP=list();
501  for (i=1;i<=s;i=i+1)
502  {
503    if (degree_P(L[i])==d)
504    {
505      counter=counter+1;
506      LP[counter]=L[i];
507    }
508  }
509  return(LP);
510}
511///////////////////////////////////////////////////////////////////////////////
512static proc subset (ideal I,ideal J)
513{
514  // checks wether I is contained in J and returns a boolean
515  // remark : J is assumed to be given by a standard basis
516  int s=size(I);
517  int i;
518  for (i=1;i<=s;i=i+1)
519  {
520    if ( NF(I[i],std(J)) != 0 )
521    {
522      return(0);
523    }
524  }
525  return(1);
526}
527///////////////////////////////////////////////////////////////////////////////
528static proc belongs (list P,ideal I)
529{
530  // checks if affine point P is contained in V(I) and returns a boolean
531  // remark : P[1] is assumed to be an ideal given by a standard basis
532  if (typeof(P[1])=="ideal")
533  {
534    return(subset(I,P[1]));
535  }
536  else
537  {
538    ERROR("first argument must be an affine point");
539  }
540}
541///////////////////////////////////////////////////////////////////////////////
542static proc equals (ideal I,ideal J)
543{
544  // checks if I is equal to J and returns a boolean
545  // remark : I and J are assumed to be given by a standard basis
546  int answer=0;
547  if (subset(I,J)==1)
548  {
549    if (subset(J,I)==1)
550    {
551      answer=1;
552    }
553  }
554  return(answer);
555}
556///////////////////////////////////////////////////////////////////////////////
557static proc isInLP (ideal P,list LP)
558{
559  // checks if affine point P is a list LP and returns either its position or
560  //   zero
561  // remark : all points in LP and P itself are assumed to be given by a
562  //   standard basis
563  // warning : the procedure does not check whether the points are affine or
564  //   not
565  int s=size(LP);
566  if (s==0)
567  {
568    return(0);
569  }
570  int i;
571  for (i=1;i<=s;i=i+1)
572  {
573    if (equals(P,LP[i][1])==1)
574    {
575      return(i);
576    }
577  }
578  return(0);
579}
580///////////////////////////////////////////////////////////////////////////////
581static proc res_deg ()
582{
583  // computes the residual degree of the basering with respect to its prime
584  //   field
585  // warning : minpoly must depend on a parameter called "a"
586  int ext;
587  string s_m=string(minpoly);
588  if (s_m=="0")
589  {
590    ext=1;
591  }
592  else
593  {
594    ring auxr=char(basering),a,lp;
595    execute("poly minp="+s_m+";");
596    ext=deg(minp);
597  }
598  return(ext);
599}
600///////////////////////////////////////////////////////////////////////////////
601static proc Frobenius (etwas,int r)
602{
603  // applies the Frobenius map over F_{p^r} to an object defined over an
604  //   extension of such field
605  // usually it is called with r=1, i.e. the Frobenius map over the prime
606  //   field F_p
607  // returns always an object of the same type, and works correctly on
608  // numbers, polynomials, ideals, matrices or lists of the above types
609  // maybe : types vector and module should be added in the future, but they
610  //   are not needed now
611  int q=char(basering)^r;
612  if (typeof(etwas)=="number")
613  {
614    return(etwas^q);
615  }
616  if (typeof(etwas)=="poly")
617  {
618    int s=size(etwas);
619    poly f;
620    int i;
621    for (i=1;i<=s;i=i+1)
622    {
623      f=f+(leadcoef(etwas[i])^q)*leadmonom(etwas[i]);
624    }
625    return(f);
626  }
627  if (typeof(etwas)=="ideal")
628  {
629    int s=ncols(etwas);
630    ideal I;
631    int i;
632    for (i=1;i<=s;i=i+1)
633    {
634      I[i]=Frobenius(etwas[i],r);
635    }
636    return(I);
637  }
638  if (typeof(etwas)=="matrix")
639  {
640    int m=nrows(etwas);
641    int n=ncols(etwas);
642    matrix A[m][n];
643    int i,j;
644    for (i=1;i<=m;i=i+1)
645    {
646      for (j=1;j<=n;j=j+1)
647      {
648        A[i,j]=Frobenius(etwas[i,j],r);
649      }
650    }
651    return(A);
652  }
653  if (typeof(etwas)=="list")
654  {
655    int s=size(etwas);
656    list L=list();
657    int i;
658    for (i=1;i<=s;i=i+1)
659    {
660      if (typeof(etwas[i])<>"none")
661      {
662        L[i]=Frobenius(etwas[i],r);
663      }
664    }
665    return(L);
666  }
667  return(etwas);
668}
669///////////////////////////////////////////////////////////////////////////////
670static proc conj_b (list L,int r)
671{
672  // applies the Frobenius map over F_{p^r} to a list of type HNE defined over
673  //   a larger extension
674  // when r=1 it turns to be the Frobenius map over the prime field F_{p}
675  // returns : a list of type HNE which is either conjugate of the input or
676  //     the same list in case of L being actually defined over the base field
677  //     F_{p^r}
678  int p=char(basering);
679  int Q=p^r;
680  list LL=list();
681  int m=nrows(L[1]);
682  int n=ncols(L[1]);
683  matrix A[m][n];
684  poly f;
685  poly aux;
686  int i,j;
687  for (i=1;i<=m;i=i+1)
688  {
689    for (j=1;j<=n;j=j+1)
690    {
691      aux=L[1][i,j];
692      if (aux<>x)
693      {
694        A[i,j]=aux^Q;
695      }
696      else
697      {
698        A[i,j]=aux;
699        break;
700      }
701    }
702  }
703  m=size(L[4]);
704  for (i=1;i<=m;i=i+1)
705  {
706    f=f+(leadcoef(L[4][i])^Q)*leadmonom(L[4][i]);
707  }
708  LL[1]=A;
709  LL[2]=L[2];
710  LL[3]=L[3];
711  LL[4]=f;
712  return(LL);
713}
714///////////////////////////////////////////////////////////////////////////////
715static proc grad_b (list L,int r)
716{
717  // computes the degree of a list of type HNE which is actually defined over
718  //    F_{p^r} eventhough it is given in an extension of such field
719  int gr=1;
720  int rd=res_deg() div r;
721  list LL=L;
722  int i;
723  for (i=1;i<=rd;i=i+1)
724  {
725    LL=conj_b(LL,r);
726    if ( LL[1]==L[1] && LL[4]==L[4] )
727    {
728      break;
729    }
730    else
731    {
732      gr=gr+1;
733    }
734  }
735  return(gr);
736}
737///////////////////////////////////////////////////////////////////////////////
738static proc conj_bs (list L,int r)
739{
740  // computes all the conjugates over F_{p^r} of a list of type HNE defined
741  //   over an extension
742  // returns : a list of lists of type HNE, where the first one is the input
743  //   list
744  // remark : notice that the degree of the branch is then the size of the
745  //   output
746  list branches=list();
747  int gr=1;
748  branches[1]=L;
749  int rd=res_deg() div r;
750  list LL=L;
751  int i;
752  for (i=1;i<=rd;i=i+1)
753  {
754    LL=conj_b(LL,r);
755    if ( LL[1]==L[1] && LL[4]==L[4] )
756    {
757      break;
758    }
759    else
760    {
761      gr=gr+1;
762      branches[gr]=LL;
763    }
764  }
765  return(branches);
766}
767///////////////////////////////////////////////////////////////////////////////
768static proc subfield (sf)
769{
770  // writes the generator "a" of a subfield of the coefficients field of
771  //   basering in terms of the the current generator (also called "a") as a
772  //   string sf is an existing ring whose coefficient field is such a subfield
773  // warning : in basering there must be a variable called "x" and subfield
774  //   must not be prime
775  def base_r=basering;
776  string new_m=string(minpoly);
777  setring sf;
778  string old_m=string(minpoly);
779  if (old_m==new_m)
780  {
781    setring base_r;
782    return("a");
783  }
784  else
785  {
786    if (old_m<>string(0))
787    {
788      ring auxring=char(basering),(a,x),lp;
789      execute("poly mpol="+old_m+";");
790      mpol=subst(mpol,a,x);
791      setring base_r;
792      poly mpol=imap(auxring,mpol);
793      kill(auxring);
794      string answer="? error : non-primitive element";
795      int r=res_deg();
796      int q=char(basering)^r;
797      int i;
798      number b;
799      for (i=1;i<=q-2;i=i+1)
800      {
801        b=a^i;
802        if (subst(mpol,x,b)==0)
803        {
804          answer=string(b);
805          break;
806        }
807      }
808      if (answer<>"? error : non-primitive element")
809      {
810        return(answer);
811      }
812      else
813      {
814        // list all the elements of the finite field F_q
815        int p=char(basering);
816        list FF1,FF2;
817        for (i=0;i<p;i=i+1)
818        {
819          FF1[i+1]=number(i);
820        }
821        int s,j,k;
822        for (i=1;i<r;i=i+1)
823        {
824          s=size(FF1);
825          for (j=1;j<=s;j=j+1)
826          {
827            FF1[j]=FF1[j]*a;
828          }
829          FF2=FF1;
830          for (k=1;k<p;k=k+1)
831          {
832            for (j=1;j<=s;j=j+1)
833            {
834              FF1[j]=FF1[j]+number(1);
835            }
836            FF2=FF2+FF1;
837          }
838          FF1=FF2;
839          FF2=list();
840        }
841        kill(FF2);
842//        list Fp;
843//        ideal facs=factorize(x^(q-1)-1,1);
844//        for (i=1;i<=q-1;i=i+1)
845//        {
846//          Fq[i]=number(subst(facs[i],x,0));
847//        }
848        for (i=1;i<=q-1;i=i+1)
849        {
850//          b=Fq[i];
851          b=FF1[i];
852          if (subst(mpol,x,b)==0)
853          {
854            answer=string(b);
855            break;
856          }
857        }
858      }
859      kill(FF1);
860      return(answer);
861    }
862    else
863    {
864      dbprint(printlevel+1,"warning : minpoly=0 in the subfield;
865                              you should check that nothing is wrong");
866      return(string(1));
867    }
868  }
869}
870///////////////////////////////////////////////////////////////////////////////
871static proc importdatum (sf,string datum,string rel)
872{
873  // fetchs a poly with name "datum" to the current basering from the ring sf
874  //   such that the generator is given by string "rel"
875  // warning : ring sf must have only variables (x,y) and basering must have
876  //   at least (x,y)
877  // warning : the case of minpoly=0 is not regarded; there you can use "imap"
878  //   instead
879  def base_r=basering;
880  if (rel=="a")
881  {
882    setring sf;
883    execute("poly pdatum="+datum+";");
884    setring base_r;
885    poly pdatum=imap(sf,pdatum);
886    return(pdatum);
887  }
888  else
889  {
890    setring sf;
891    execute("string sdatum=string("+datum+");");
892    ring auxring=char(basering),(a,x,y),lp;
893    execute("poly pdatum="+sdatum+";");
894    execute("map phi=basering,"+rel+",x,y;");
895    pdatum=phi(pdatum);
896    string snewdatum=string(pdatum);
897    setring base_r;
898    execute("poly pdatum="+snewdatum+";");
899    kill(auxring);
900    return(pdatum);
901  }
902}
903///////////////////////////////////////////////////////////////////////////////
904static proc rationalize (lf,string datum,string rel)
905{
906  // fetchs a poly with name "datum" to the current basering from the ring lf
907  //   and larger coefficients field, where the generator of current ring is
908  //   given by string "rel" and "datum" is actually defined over the small
909  //   field
910  // warning : "ring lf must have only variables (x,y) and basering must have
911  //           at least (x,y)
912  // warning : the case of minpoly=0 is supposed unnecessary, since then
913  //           "datum" should be
914  //   already written only int the right way, i.e. in terms of the prime field
915  def base_r=basering;
916  if (rel=="a")
917  {
918    setring lf;
919    execute("poly pdatum="+datum+";");
920    setring base_r;
921    poly pdatum=imap(lf,pdatum);
922    return(pdatum);
923  }
924  else
925  {
926    setring lf;
927    execute("string sdatum=string("+datum+");");
928    ring auxring=char(basering),(a,b,x,y,t),lp;
929    execute("poly pdatum="+sdatum+";");
930    execute("poly prel=b-"+rel+";");
931    ideal I=pdatum,prel;
932    I=eliminate(I,a);
933    poly newdatum=I[1];   // hopefully it was done correctly and size(I)=1 !!!!
934    newdatum=subst(newdatum,b,a);
935    string snewdatum=string(newdatum);
936    setring base_r;
937    execute("poly newdatum="+snewdatum+";");
938    kill(auxring);
939    return(newdatum);
940  }
941}
942///////////////////////////////////////////////////////////////////////////////
943static proc place (intvec Pp,int sing,list CURVE)
944{
945  // computes the "rational" places which are defined over a (closed) point
946  // Pp points to an appropriate point of the given curve
947  // creates : local rings (if they do not exist yet) and then add the
948  //   places to a list
949  //   each place is given basically by the coordinates of the point and a
950  //   list HNdevelop
951  // returns : list with all updated data of the curve
952  // if the places already exist they are not computed again
953  // if sing==1 the point is assumed singular and computes the local conductor
954  //   for all places using the local invariants of the branches
955  // if sing==2 the point is assumed singular and computes the local conductor
956  //   for all places using the Dedekind formula and local parametrizations
957  //   of the branches
958  // if sing<>1&2 the point is assumed non-singular and the local conductor
959  //   should be zero
960  list PP=list();
961  if (Pp[1]==0)
962  {
963    if (Pp[2]==0)
964    {
965      PP=Aff_SPoints[Pp[3]];
966    }
967    if (Pp[2]==1)
968    {
969      PP=Inf_Points[1][Pp[3]];
970    }
971    if (Pp[2]==2)
972    {
973      PP=Inf_Points[2][Pp[3]];
974    }
975  }
976  else
977  {
978    PP=Aff_Points(Pp[2])[Pp[3]];
979  }
980  if (PP[2]<>0)
981  {
982    return(CURVE);
983  }
984  intvec PtoPl;
985  def base_r=basering;
986  int ext1;
987  list Places=CURVE[3];
988  intvec Conductor=CURVE[4];
989  list update_CURVE=CURVE;
990  if (typeof(PP[1])=="ideal")
991  {
992    ideal P=PP[1];
993    if (size(P)==2)
994    {
995      int d=deg(P[1]);
996      poly aux=subst(P[2],y,1);
997      d=d*deg(aux);
998      ext1=d;
999      // the point is (A:B:1) but one must distinguish several cases
1000      // P is assumed to be a std. resp. "(x,y),lp" and thus P[1] depends
1001      // only on "y"
1002      if (d==1)
1003      {
1004        // the point is rational
1005        number B=-number(subst(P[1],y,0));
1006        poly aux2=subst(P[2],y,B);
1007        number A=-number(subst(aux2,x,0));
1008        // the point is (A:B:1)
1009        ring local_aux=char(basering),(x,y),ls;
1010        number coord@1=imap(base_r,A);
1011        number coord@2=imap(base_r,B);
1012        number coord@3=number(1);
1013        map phi=base_r,x+coord@1,y+coord@2;
1014        poly CHI=phi(CHI);
1015      }
1016      else
1017      {
1018        if (deg(P[1])==1)
1019        {
1020          // the point is non-rational but the second component needs no
1021          // field extension
1022          number B=-number(subst(P[1],y,0));
1023          poly aux2=subst(P[2],y,B);
1024          // the point has degree d>1
1025          // careful : the parameter will be called "a" anyway
1026          ring local_aux=(char(basering),a),(x,y),ls;
1027          map psi=base_r,a,0;
1028          minpoly=number(psi(aux2));
1029          number coord@1=a;
1030          number coord@2=imap(base_r,B);
1031          number coord@3=number(1);
1032          // the point is (a:B:1)
1033          map phi=base_r,x+a,y+coord@2;
1034          poly CHI=phi(CHI);
1035        }
1036        else
1037        {
1038          if (deg(subst(P[2],y,1))==1)
1039          {
1040            // the point is non-rational but the needed minpoly is just P[1]
1041            // careful : the parameter will be called "a" anyway
1042            poly P1=P[1];
1043            poly P2=P[2];
1044            ring local_aux=(char(basering),a),(x,y),ls;
1045            map psi=base_r,0,a;
1046            minpoly=number(psi(P1));
1047            // the point looks like (A:a:1)
1048            // A is computed by substituting y=a in P[2]
1049            poly aux1=imap(base_r,P2);
1050            poly aux2=subst(aux1,y,a);
1051            number coord@1=-number(subst(aux2,x,0));
1052            number coord@2=a;
1053            number coord@3=number(1);
1054            map phi=base_r,x+coord@1,y+a;
1055            poly CHI=phi(CHI);
1056          }
1057          else
1058          {
1059            // this is the most complicated case of non-rational point
1060            // firstly : construct an extension of degree d and guess the
1061            //           minpoly
1062            poly P1=P[1];
1063            poly P2=P[2];
1064            int p=char(basering);
1065            int Q=p^d;
1066            ring aux_r=(Q,a),(x,y,t),ls;
1067            string minpoly_string=string(minpoly);
1068            ring local_aux=(char(basering),a),(x,y),ls;
1069            execute("minpoly="+minpoly_string+";");
1070            // secondly : compute one root of P[1]
1071            poly P_1=imap(base_r,P1);
1072            poly P_2=imap(base_r,P2);
1073            ideal factors1=factorize(P_1,1);       // hopefully this works !!!!
1074            number coord@2=-number(subst(factors1[1],y,0));
1075            // thirdly : compute one of the first components for the above root
1076            poly P_0=subst(P_2,y,coord@2);
1077            ideal factors2=factorize(P_0,1);       // hopefully this works !!!!
1078            number coord@1=-number(subst(factors2[1],x,0));
1079            number coord@3=number(1);
1080            map phi=base_r,x+coord@1,y+coord@2;
1081            poly CHI=phi(CHI);
1082            kill(aux_r);
1083          }
1084        }
1085      }
1086    }
1087    else
1088    {
1089      // this should not happen in principle
1090      ERROR("non-valid parameter");
1091    }
1092  }
1093  else
1094  {
1095    if (typeof(PP[1])=="poly")
1096    {
1097      poly P=PP[1];
1098      ring r_auxz=char(basering),(x,y,z),lp;
1099      poly CHI=imap(base_r,CHI);
1100      CHI=homog(CHI,z);
1101      setring base_r;
1102      poly aux=subst(P,y,1);
1103      if (aux==1)
1104      {
1105        // the point is (1:0:0)
1106        ring local_aux=char(basering),(x,y),ls;
1107        number coord@1=number(1);
1108        number coord@2=number(0);
1109        number coord@3=number(0);
1110        map Phi=r_auxz,1,x,y;
1111        poly CHI=Phi(CHI);
1112        ext1=1;
1113      }
1114      else
1115      {
1116        // the point is (A:1:0) where A is a root of aux
1117        int d=deg(aux);
1118        ext1=d;
1119        if (d==1)
1120        {
1121          // the point is rational
1122          number A=-number(subst(aux,x,0));
1123          ring local_aux=char(basering),(x,y),ls;
1124          number coord@1=imap(base_r,A);
1125          number coord@2=number(1);
1126          number coord@3=number(0);
1127          map Phi=r_auxz,x+coord@1,1,y;
1128          poly CHI=Phi(CHI);
1129        }
1130        else
1131        {
1132          // the point has degree d>1
1133          // careful : the parameter will be called "a" anyway
1134          ring local_aux=(char(basering),a),(x,y),ls;
1135          map psi=base_r,a,1;
1136          minpoly=number(psi(P));
1137          number coord@1=a;
1138          number coord@2=number(1);
1139          number coord@3=number(0);
1140          map Phi=r_auxz,x+a,1,y;
1141          poly CHI=Phi(CHI);
1142        }
1143      }
1144      kill(r_auxz);
1145    }
1146    else
1147    {
1148      ERROR("a point must have a poly or ideal in the first component");
1149    }
1150  }
1151  export(coord@1);
1152  export(coord@2);
1153  export(coord@3);
1154  export(CHI);
1155  int i,j,k;
1156  int m,n;
1157  list L@HNE=essdevelop(CHI);
1158  export(L@HNE);
1159  int n_branches=size(L@HNE);
1160  list Li_aux=list();
1161  int N_branches;
1162  int N=size(Places);
1163  if (sing==1)
1164  {
1165    list delta2=list();
1166    for (i=1;i<=n_branches;i=i+1)
1167    {
1168      delta2[i]=invariants(L@HNE[i])[5];
1169    }
1170    int dq;
1171  }
1172  int ext2=res_deg();
1173  list dgs=list();
1174  int ext_0;
1175  int check;
1176  string sss,olda,newa;
1177  if (defined(Q)==0)
1178  {
1179    int Q;
1180  }
1181  if (ext1==1)
1182  {
1183    if (ext2==1)
1184    {
1185      if (sing==1)
1186      {
1187        intmat I_mult[n_branches][n_branches];
1188        if (n_branches>1)
1189        {
1190          for (i=1;i<=n_branches-1;i=i+1)
1191          {
1192            for (j=i+1;j<=n_branches;j=j+1)
1193            {
1194              I_mult[i,j]=intersection(L@HNE[i],L@HNE[j]);
1195              I_mult[j,i]=I_mult[i,j];
1196            }
1197          }
1198        }
1199      }
1200      if (size(update_CURVE[5])>0)
1201      {
1202        if (typeof(update_CURVE[5][1])=="list")
1203        {
1204          check=1;
1205        }
1206      }
1207      if (check==0)
1208      {
1209        intvec dgs_points(1);
1210        ring S(1)=char(basering),(x,y,t),ls;
1211        list BRANCHES=list();
1212        list POINTS=list();
1213        list LOC_EQS=list();
1214        list PARAMETRIZATIONS=list();
1215        export(BRANCHES);
1216        export(POINTS);
1217        export(LOC_EQS);
1218        export(PARAMETRIZATIONS);
1219      }
1220      else
1221      {
1222        intvec dgs_points(1)=update_CURVE[5][1][2];
1223        def S1=update_CURVE[5][1][1];
1224        execute("ring S(1)="+string(update_CURVE[5][1][1])+";");
1225        fetchall(S1);
1226        kill(S1);
1227      }
1228      N_branches=size(BRANCHES);
1229      for (i=1;i<=n_branches;i=i+1)
1230      {
1231        dgs_points(1)[N_branches+i]=1;
1232        POINTS[N_branches+i]=list();
1233        POINTS[N_branches+i][1]=imap(local_aux,coord@1);
1234        POINTS[N_branches+i][2]=imap(local_aux,coord@2);
1235        POINTS[N_branches+i][3]=imap(local_aux,coord@3);
1236        LOC_EQS[N_branches+i]=imap(local_aux,CHI);
1237        setring HNEring;
1238        Li_aux=L@HNE[i];
1239        setring S(1);
1240        BRANCHES=insert(BRANCHES,imap(HNEring,Li_aux),N_branches+i-1);
1241        PARAMETRIZATIONS[N_branches+i]=param(BRANCHES[N_branches+i],0);
1242        N=N+1;
1243        intvec iw=1,N_branches+i;
1244        Places[N]=iw;
1245        if (sing==1)
1246        {
1247          dq=delta2[i];
1248          for (j=1;j<=n_branches;j=j+1)
1249          {
1250            dq=dq+I_mult[i,j];
1251          }
1252          Conductor[N]=dq;
1253        }
1254        if (sing==2)
1255        {
1256          Conductor[N]=local_conductor(iw[2],S(1));
1257        }
1258        kill(iw);
1259        PtoPl[i]=N;
1260      }
1261      setring base_r;
1262      update_CURVE[5][1]=list();
1263      update_CURVE[5][1][1]=S(1);
1264      update_CURVE[5][1][2]=dgs_points(1);
1265    }
1266    else
1267    {
1268      // we start with a rational point but we get non-rational branches
1269      // they may have different degrees and then we may need reduce the
1270      // field extensions for each one, and finally check if the minpoly
1271      // fetchs with S(i) or not
1272      // if one of the branches is rational, we may trust that is is written
1273      // correctly
1274      if (sing==1)
1275      {
1276        int n_geobrs;
1277        int counter_c;
1278        list auxgb=list();
1279        list geobrs=list();
1280        for (i=1;i<=n_branches;i=i+1)
1281        {
1282          auxgb=conj_bs(L@HNE[i],1);
1283          dgs[i]=size(auxgb);
1284          n_geobrs=n_geobrs+dgs[i];
1285          for (j=1;j<=dgs[i];j=j+1)
1286          {
1287            counter_c=counter_c+1;
1288            geobrs[counter_c]=auxgb[j];
1289          }
1290        }
1291        intmat I_mult[n_geobrs][n_geobrs];
1292        for (i=1;i<n_geobrs;i=i+1)
1293        {
1294          for (j=i+1;j<=n_geobrs;j=j+1)
1295          {
1296            I_mult[i,j]=intersection(geobrs[i],geobrs[j]);
1297            I_mult[j,i]=I_mult[i,j];
1298          }
1299        }
1300        kill(auxgb,geobrs);
1301      }
1302      else
1303      {
1304        for (i=1;i<=n_branches;i=i+1)
1305        {
1306          dgs[i]=grad_b(L[i],1);
1307        }
1308      }
1309      // the actual degree of each branch is computed and now check if the
1310      // local ring exists
1311      for (i=1;i<=n_branches;i=i+1)
1312      {
1313        ext_0=dgs[i];
1314        if (size(update_CURVE[5])>=ext_0)
1315        {
1316          if (typeof(update_CURVE[5][ext_0])=="list")
1317          {
1318            check=1;
1319          }
1320        }
1321        if (check==0)
1322        {
1323          if (ext_0>1)
1324          {
1325            if (ext_0==ext2)
1326            {
1327              sss=string(minpoly);
1328            }
1329            else
1330            {
1331              Q=char(basering)^ext_0;
1332              ring auxxx=(Q,a),z,lp;
1333              sss=string(minpoly);
1334              setring base_r;
1335              kill(auxxx);
1336            }
1337            ring S(ext_0)=(char(basering),a),(x,y,t),ls;
1338            execute("minpoly="+sss+";");
1339          }
1340          else
1341          {
1342            ring S(ext_0)=char(basering),(x,y,t),ls;
1343          }
1344          intvec dgs_points(ext_0);
1345          list BRANCHES=list();
1346          list POINTS=list();
1347          list LOC_EQS=list();
1348          list PARAMETRIZATIONS=list();
1349          export(BRANCHES);
1350          export(POINTS);
1351          export(LOC_EQS);
1352          export(PARAMETRIZATIONS);
1353        }
1354        else
1355        {
1356          intvec dgs_points(ext_0)=update_CURVE[5][ext_0][2];
1357          def Sext_0=update_CURVE[5][ext_0][1];
1358          setring Sext_0;
1359          string SM=string(minpoly);
1360          string SR=string(update_CURVE[5][ext_0][1]);
1361          execute("ring S("+string(ext_0)+")="+SR+";");
1362          execute("minpoly="+SM+";");
1363          kill(SM,SR);
1364          fetchall(Sext_0);
1365          kill(Sext_0);
1366        }
1367        N_branches=size(BRANCHES);
1368        dgs_points(ext_0)[N_branches+1]=1;
1369        POINTS[N_branches+1]=list();
1370        POINTS[N_branches+1][1]=imap(local_aux,coord@1);
1371        POINTS[N_branches+1][2]=imap(local_aux,coord@2);
1372        POINTS[N_branches+1][3]=imap(local_aux,coord@3);
1373        LOC_EQS[N_branches+1]=imap(local_aux,CHI);
1374        // now fetch the branches into the new local ring
1375        if (ext_0==1)
1376        {
1377          setring HNEring;
1378          Li_aux=L@HNE[i];
1379          setring S(1);
1380          BRANCHES=insert(BRANCHES,imap(HNEring,Li_aux),N_branches);
1381        }
1382        else
1383        {
1384          // rationalize branche
1385          setring HNEring;
1386          newa=subfield(S(ext_0));
1387          m=nrows(L@HNE[i][1]);
1388          n=ncols(L@HNE[i][1]);
1389          setring S(ext_0);
1390          list Laux=list();
1391          poly paux=rationalize(HNEring,"L@HNE["+string(i)+"][4]",newa);
1392          matrix Maux[m][n];
1393          for (j=1;j<=m;j=j+1)
1394          {
1395            for (k=1;k<=n;k=k+1)
1396            {
1397              Maux[j,k]=rationalize(HNEring,"L@HNE["+string(i)+"][1]["+
1398                           string(j)+","+string(k)+"]",newa);
1399            }
1400          }
1401          setring HNEring;
1402          intvec Li2=L@HNE[i][2];
1403          int Li3=L@HNE[i][3];
1404          setring S(ext_0);
1405          Laux[1]=Maux;
1406          Laux[2]=Li2;
1407          Laux[3]=Li3;
1408          Laux[4]=paux;
1409          BRANCHES=insert(BRANCHES,Laux,N_branches);
1410          kill(Laux,Maux,paux,Li2,Li3);
1411        }
1412        PARAMETRIZATIONS[N_branches+1]=param(BRANCHES[N_branches+1],0);
1413        N=N+1;
1414        intvec iw=ext_0,N_branches+1;
1415        Places[N]=iw;
1416        if (sing==2)
1417        {
1418          Conductor[N]=local_conductor(iw[2],S(ext_0));
1419        }
1420        kill(iw);
1421        PtoPl[i]=N;
1422        setring HNEring;
1423        update_CURVE[5][ext_0]=list();
1424        update_CURVE[5][ext_0][1]=S(ext_0);
1425        update_CURVE[5][ext_0][2]=dgs_points(ext_0);
1426      }
1427      if (sing==1)
1428      {
1429        int N_ini=N-n_branches;
1430        counter_c=1;
1431        for (i=1;i<=n_branches;i=i+1)
1432        {
1433          dq=delta2[i];
1434          for (j=1;j<=n_geobrs;j=j+1)
1435          {
1436            dq=dq+I_mult[counter_c,j];
1437          }
1438          Conductor[N_ini+i]=dq;
1439          counter_c=counter_c+dgs[i];
1440        }
1441      }
1442      setring base_r;
1443    }
1444  }
1445  else
1446  {
1447    if (ext1==ext2)
1448    {
1449      // the degree of the point equals to the degree of all branches
1450      // one must just fetch the minpoly's of local_aux, HNEring and S(ext2)
1451      if (sing==1)
1452      {
1453        intmat I_mult[n_branches][n_branches];
1454        if (n_branches>1)
1455        {
1456          for (i=1;i<=n_branches-1;i=i+1)
1457          {
1458            for (j=i+1;j<=n_branches;j=j+1)
1459            {
1460              I_mult[i,j]=intersection(L@HNE[i],L@HNE[j]);
1461              I_mult[j,i]=I_mult[i,j];
1462            }
1463          }
1464        }
1465      }
1466      if (size(update_CURVE[5])>=ext2)
1467      {
1468        if (typeof(update_CURVE[5][ext2])=="list")
1469        {
1470          check=1;
1471        }
1472      }
1473      if (check==0)
1474      {
1475        sss=string(minpoly);
1476        ring S(ext2)=(char(basering),a),(x,y,t),ls;
1477        execute("minpoly="+sss+";");
1478        intvec dgs_points(ext2);
1479        list BRANCHES=list();
1480        list POINTS=list();
1481        list LOC_EQS=list();
1482        list PARAMETRIZATIONS=list();
1483        export(BRANCHES);
1484        export(POINTS);
1485        export(LOC_EQS);
1486        export(PARAMETRIZATIONS);
1487      }
1488      else
1489      {
1490        intvec dgs_points(ext2)=update_CURVE[5][ext2][2];
1491        def Sext2=update_CURVE[5][ext2][1];
1492        setring Sext2;
1493        string SM=string(minpoly);
1494        string SR=string(update_CURVE[5][ext2][1]);
1495        execute("ring S("+string(ext2)+")="+SR+";");
1496        execute("minpoly="+SM+";");
1497        kill(SM,SR);
1498        fetchall(Sext2);
1499        kill(Sext2);
1500      }
1501      N_branches=size(BRANCHES);
1502      for (i=1;i<=n_branches;i=i+1)
1503      {
1504        // fetch all the data into the new local ring
1505        olda=subfield(local_aux);
1506        dgs_points(ext2)[N_branches+i]=ext1;
1507        POINTS[N_branches+i]=list();
1508        POINTS[N_branches+i][1]=number(importdatum(local_aux,"coord@1",olda));
1509        POINTS[N_branches+i][2]=number(importdatum(local_aux,"coord@2",olda));
1510        POINTS[N_branches+i][3]=number(importdatum(local_aux,"coord@3",olda));
1511        LOC_EQS[N_branches+i]=importdatum(local_aux,"CHI",olda);
1512        newa=subfield(HNEring);
1513        setring HNEring;
1514        m=nrows(L@HNE[i][1]);
1515        n=ncols(L@HNE[i][1]);
1516        setring S(ext2);
1517        list Laux=list();
1518        poly paux=importdatum(HNEring,"L@HNE["+string(i)+"][4]",newa);
1519        matrix Maux[m][n];
1520        for (j=1;j<=m;j=j+1)
1521        {
1522          for (k=1;k<=n;k=k+1)
1523          {
1524            Maux[j,k]=importdatum(HNEring,"L@HNE["+string(i)+"][1]["+
1525                                   string(j)+","+string(k)+"]",newa);
1526          }
1527        }
1528        setring HNEring;
1529        intvec Li2=L@HNE[i][2];
1530        int Li3=L@HNE[i][3];
1531        setring S(ext2);
1532        Laux[1]=Maux;
1533        Laux[2]=Li2;
1534        Laux[3]=Li3;
1535        Laux[4]=paux;
1536        BRANCHES=insert(BRANCHES,Laux,N_branches+i-1);
1537        kill(Laux,Maux,paux,Li2,Li3);
1538        PARAMETRIZATIONS[N_branches+i]=param(BRANCHES[N_branches+i],0);
1539        N=N+1;
1540        intvec iw=ext2,N_branches+i;
1541        Places[N]=iw;
1542        if (sing==1)
1543        {
1544          dq=delta2[i];
1545          for (j=1;j<=n_branches;j=j+1)
1546          {
1547            dq=dq+I_mult[i,j];
1548          }
1549          Conductor[N]=dq;
1550        }
1551        if (sing==2)
1552        {
1553          Conductor[N]=local_conductor(iw[2],S(ext2));
1554        }
1555        kill(iw);
1556        PtoPl[i]=N;
1557      }
1558      setring base_r;
1559      update_CURVE[5][ext2]=list();
1560      update_CURVE[5][ext2][1]=S(ext2);
1561      update_CURVE[5][ext2][2]=dgs_points(ext2);
1562    }
1563    else
1564    {
1565      // this is the most complicated case
1566      if (sing==1)
1567      {
1568        int n_geobrs;
1569        int counter_c;
1570        list auxgb=list();
1571        list geobrs=list();
1572        for (i=1;i<=n_branches;i=i+1)
1573        {
1574          auxgb=conj_bs(L@HNE[i],ext1);
1575          dgs[i]=size(auxgb);
1576          n_geobrs=n_geobrs+dgs[i];
1577          for (j=1;j<=dgs[i];j=j+1)
1578          {
1579            counter_c=counter_c+1;
1580            geobrs[counter_c]=auxgb[j];
1581          }
1582        }
1583        intmat I_mult[n_geobrs][n_geobrs];
1584        for (i=1;i<n_geobrs;i=i+1)
1585        {
1586          for (j=i+1;j<=n_geobrs;j=j+1)
1587          {
1588            I_mult[i,j]=intersection(geobrs[i],geobrs[j]);
1589            I_mult[j,i]=I_mult[i,j];
1590          }
1591        }
1592        kill(auxgb,geobrs);
1593      }
1594      else
1595      {
1596        for (i=1;i<=n_branches;i=i+1)
1597        {
1598          dgs[i]=grad_b(L@HNE[i],ext1);
1599        }
1600      }
1601      for (i=1;i<=n_branches;i=i+1)
1602      {
1603        // first compute the actual degree of each branch and check if the
1604        // local ring exists
1605        ext_0=ext1*dgs[i];
1606        if (size(update_CURVE[5])>=ext_0)
1607        {
1608          if (typeof(update_CURVE[5][ext_0])=="list")
1609          {
1610            check=1;
1611          }
1612        }
1613        if (check==0)
1614        {
1615          if (ext_0>ext1)
1616          {
1617            if (ext_0==ext2)
1618            {
1619              sss=string(minpoly);
1620            }
1621            else
1622            {
1623              Q=char(basering)^ext_0;
1624              ring auxxx=(Q,a),z,lp;
1625              sss=string(minpoly);
1626              setring base_r;
1627              kill(auxxx);
1628            }
1629          }
1630          else
1631          {
1632            setring local_aux;
1633            sss=string(minpoly);
1634          }
1635          ring S(ext_0)=(char(basering),a),(x,y,t),ls;
1636          execute("minpoly="+sss+";");
1637          intvec dgs_points(ext_0);
1638          list BRANCHES=list();
1639          list POINTS=list();
1640          list LOC_EQS=list();
1641          list PARAMETRIZATIONS=list();
1642          export(BRANCHES);
1643          export(POINTS);
1644          export(LOC_EQS);
1645          export(PARAMETRIZATIONS);
1646        }
1647        else
1648        {
1649          intvec dgs_points(ext_0)=update_CURVE[5][ext_0][2];
1650          def Sext_0=update_CURVE[5][ext_0][1];
1651          setring Sext_0;
1652          string SM=string(minpoly);
1653          string SR=string(update_CURVE[5][ext_0][1]);
1654          execute("ring S("+string(ext_0)+")="+SR+";");
1655          execute("minpoly="+SM+";");
1656          kill(SM,SR);
1657          fetchall(badring);
1658          kill(badring);
1659        }
1660        N_branches=size(BRANCHES);
1661        // now fetch all the data into the new local ring
1662        olda=subfield(local_aux);
1663        dgs_points(ext_0)[N_branches+1]=ext1;
1664        POINTS[N_branches+1]=list();
1665        POINTS[N_branches+1][1]=number(importdatum(local_aux,"coord@1",olda));
1666        POINTS[N_branches+1][2]=number(importdatum(local_aux,"coord@2",olda));
1667        POINTS[N_branches+1][3]=number(importdatum(local_aux,"coord@3",olda));
1668        LOC_EQS[N_branches+1]=importdatum(local_aux,"CHI",olda);
1669        setring HNEring;
1670        newa=subfield(S(ext_0));
1671        m=nrows(L@HNE[i][1]);
1672        n=ncols(L@HNE[i][1]);
1673        setring S(ext_0);
1674        list Laux=list();
1675        poly paux=rationalize(HNEring,"L@HNE["+string(i)+"][4]",newa);
1676        matrix Maux[m][n];
1677        for (j=1;j<=m;j=j+1)
1678        {
1679          for (k=1;k<=n;k=k+1)
1680          {
1681            Maux[j,k]=rationalize(HNEring,"L@HNE["+string(i)+"][1]["+
1682               string(j)+","+string(k)+"]",newa);
1683          }
1684        }
1685        setring HNEring;
1686        intvec Li2=L@HNE[i][2];
1687        int Li3=L@HNE[i][3];
1688        setring S(ext_0);
1689        Laux[1]=Maux;
1690        Laux[2]=Li2;
1691        Laux[3]=Li3;
1692        Laux[4]=paux;
1693        BRANCHES=insert(BRANCHES,Laux,N_branches);
1694        kill(Laux,Maux,paux,Li2,Li3);
1695        PARAMETRIZATIONS[N_branches+1]=param(BRANCHES[N_branches+1],0);
1696        N=N+1;
1697        intvec iw=ext_0,N_branches+1;
1698        Places[N]=iw;
1699        if (sing==2)
1700        {
1701          Conductor[N]=local_conductor(iw[2],S(ext_0));
1702        }
1703        kill(iw);
1704        PtoPl[i]=N;
1705        setring HNEring;
1706        update_CURVE[5][ext_0]=list();
1707        update_CURVE[5][ext_0][1]=S(ext_0);
1708        update_CURVE[5][ext_0][2]=dgs_points(ext_0);
1709      }
1710      if (sing==1)
1711      {
1712        int N_ini=N-n_branches;
1713        counter_c=1;
1714        for (i=1;i<=n_branches;i=i+1)
1715        {
1716          dq=delta2[i];
1717          for (j=1;j<=n_geobrs;j=j+1)
1718          {
1719            dq=dq+I_mult[counter_c,j];
1720          }
1721          Conductor[N_ini+i]=dq;
1722          counter_c=counter_c+dgs[i];
1723        }
1724      }
1725      setring base_r;
1726    }
1727  }
1728  update_CURVE[3]=Places;
1729  update_CURVE[4]=Conductor;
1730  PP[2]=PtoPl;
1731  if (Pp[1]==0)
1732  {
1733    if (Pp[2]==0)
1734    {
1735      Aff_SPoints[Pp[3]]=PP;
1736    }
1737    if (Pp[2]==1)
1738    {
1739      Inf_Points[1][Pp[3]]=PP;
1740    }
1741    if (Pp[2]==2)
1742    {
1743      Inf_Points[2][Pp[3]]=PP;
1744    }
1745  }
1746  else
1747  {
1748    Aff_Points(Pp[2])[Pp[3]]=PP;
1749  }
1750  update_CURVE[1][1]=base_r;
1751  kill(HNEring);
1752  kill(local_aux);
1753  return(update_CURVE);
1754}
1755///////////////////////////////////////////////////////////////////////////////
1756static proc local_conductor (int k,SS)
1757{
1758  // computes the degree of the local conductor at a place of a plane curve
1759  // if the point is non-singular the result will be zero
1760  // the computation is carried out with the "Dedekind formula" via
1761  // parametrizations
1762  int a,b,Cq;
1763  def b_ring=basering;
1764  setring SS;
1765  poly fx=diff(LOC_EQS[k],x);
1766  poly fy=diff(LOC_EQS[k],y);
1767  int nr=ncols(BRANCHES[k][1]);
1768  poly xt=PARAMETRIZATIONS[k][1][1];
1769  poly yt=PARAMETRIZATIONS[k][1][2];
1770  int ordx=PARAMETRIZATIONS[k][2][1];
1771  int ordy=PARAMETRIZATIONS[k][2][2];
1772  map phi_t=basering,xt,yt,1;
1773  poly derf;
1774  if (fx<>0)
1775  {
1776    derf=fx;
1777    poly tt=diff(yt,t);
1778    b=mindeg(tt);
1779    if (ordy>-1)
1780    {
1781      while (b>=ordy)
1782      {
1783        BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
1784        nr=ncols(BRANCHES[k][1]);
1785        PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
1786        ordy=PARAMETRIZATIONS[k][2][2];
1787        yt=PARAMETRIZATIONS[k][1][2];
1788        tt=diff(yt,t);
1789        b=mindeg(tt);
1790      }
1791      xt=PARAMETRIZATIONS[k][1][1];
1792      ordx=PARAMETRIZATIONS[k][2][1];
1793    }
1794    poly ft=phi_t(derf);
1795  }
1796  else
1797  {
1798    derf=fy;
1799    poly tt=diff(xt,t);
1800    b=mindeg(tt);
1801    if (ordx>-1)
1802    {
1803      while (b>=ordx)
1804      {
1805        BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
1806        nr=ncols(BRANCHES[k][1]);
1807        PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
1808        ordx=PARAMETRIZATIONS[k][2][1];
1809        xt=PARAMETRIZATIONS[k][1][1];
1810        tt=diff(xt,t);
1811        b=mindeg(tt);
1812      }
1813      yt=PARAMETRIZATIONS[k][1][2];
1814      ordy=PARAMETRIZATIONS[k][2][2];
1815    }
1816    poly ft=phi_t(derf);
1817  }
1818  a=mindeg(ft);
1819  if ( ordx>-1 || ordy>-1 )
1820  {
1821    if (ordy==-1)
1822    {
1823      while (a>ordx)
1824      {
1825        BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
1826        nr=ncols(BRANCHES[k][1]);
1827        PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
1828        ordx=PARAMETRIZATIONS[k][2][1];
1829        xt=PARAMETRIZATIONS[k][1][1];
1830        ft=phi_t(derf);
1831        a=mindeg(ft);
1832      }
1833    }
1834    else
1835    {
1836      if (ordx==-1)
1837      {
1838        while (a>ordy)
1839        {
1840          BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
1841          nr=ncols(BRANCHES[k][1]);
1842          PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
1843          ordy=PARAMETRIZATIONS[k][2][2];
1844          yt=PARAMETRIZATIONS[k][1][2];
1845          ft=phi_t(derf);
1846          a=mindeg(ft);
1847        }
1848      }
1849      else
1850      {
1851        int ordf=ordx;
1852        if (ordx>ordy)
1853        {
1854          ordf=ordy;
1855        }
1856        while (a>ordf)
1857        {
1858          BRANCHES[k]=extdevelop(BRANCHES[k],2*nr);
1859          nr=ncols(BRANCHES[k][1]);
1860          PARAMETRIZATIONS[k]=param(BRANCHES[k],0);
1861          ordx=PARAMETRIZATIONS[k][2][1];
1862          ordy=PARAMETRIZATIONS[k][2][2];
1863          ordf=ordx;
1864          if (ordx>ordy)
1865          {
1866            ordf=ordy;
1867          }
1868          xt=PARAMETRIZATIONS[k][1][1];
1869          yt=PARAMETRIZATIONS[k][1][2];
1870          ft=phi_t(derf);
1871          a=mindeg(ft);
1872        }
1873      }
1874    }
1875  }
1876  Cq=a-b;
1877  setring b_ring;
1878  return(Cq);
1879}
1880///////////////////////////////////////////////////////////////////////////////
1881static proc max_D (intvec D1,intvec D2)
1882{
1883  // computes the maximum of two divisors (intvec)
1884  int s1=size(D1);
1885  int s2=size(D2);
1886  int i;
1887  if (s1>s2)
1888  {
1889    for (i=1;i<=s2;i=i+1)
1890    {
1891      if (D2[i]>D1[i])
1892      {
1893        D1[i]=D2[i];
1894      }
1895    }
1896    for (i=s2+1;i<=s1;i=i+1)
1897    {
1898      if (D1[i]<0)
1899      {
1900        D1[i]=0;
1901      }
1902    }
1903    return(D1);
1904  }
1905  else
1906  {
1907    for (i=1;i<=s1;i=i+1)
1908    {
1909      if (D1[i]>D2[i])
1910      {
1911        D2[i]=D1[i];
1912      }
1913    }
1914    for (i=s1+1;i<=s2;i=i+1)
1915    {
1916      if (D2[i]<0)
1917      {
1918        D2[i]=0;
1919      }
1920    }
1921    return(D2);
1922  }
1923}
1924///////////////////////////////////////////////////////////////////////////////
1925static proc deg_D (intvec D,list PP)
1926{
1927  // computes the degree of a divisor (intvec or list of integers)
1928  int i;
1929  int d=0;
1930  int s=size(D);
1931  for (i=1;i<=s;i=i+1)
1932  {
1933    d=d+D[i]*PP[i][1];
1934  }
1935  return(d);
1936}
1937///////////////////////////////////////////////////////////////////////////////
1938
1939// ============================================================================
1940// *******  MAIN PROCEDURES for the "preprocessing" of Brill-Noether   ********
1941// ============================================================================
1942
1943proc Adj_div (poly f,list #)
1944"USAGE:    Adj_div( f [,l] ); f a poly, [l a list]
1945RETURN:   list L with the computed data:
1946  @format
1947  L[1] a list of rings: L[1][1]=aff_r (affine), L[1][2]=Proj_R (projective),
1948  L[2] an intvec with 2 entries (degree, genus),
1949  L[3] a list of intvec (closed places),
1950  L[4] an intvec (conductor),
1951  L[5] a list of lists:
1952     L[5][d][1] a (local) ring over an extension of degree d,
1953     L[5][d][2] an intvec (degrees of base points of places of degree d)
1954  @end format
1955NOTE:     @code{Adj_div(f);} computes and stores the fundamental data of the
1956          plane curve defined by f as needed for AG codes.
1957          In the affine ring you can find the following data:
1958   @format
1959   poly CHI:  affine equation of the curve,
1960   ideal Aff_SLocus:  affine singular locus (std),
1961   list Inf_Points:  points at infinity
1962            Inf_Points[1]:  singular points
1963            Inf_Points[2]:  non-singular points,
1964   list Aff_SPoints:  affine singular points (if not empty).
1965   @end format
1966          In the projective ring you can find the projective equation
1967          CHI of the curve (poly).
1968          In the local rings L[5][d][1] you find:
1969   @format
1970   list POINTS:  base points of the places of degree d,
1971   list LOC_EQS:  local equations of the curve at the base points,
1972   list BRANCHES:  Hamburger-Noether developments of the places,
1973   list PARAMETRIZATIONS:  local parametrizations of the places,
1974   @end format
1975          Each entry of the list L[3] corresponds to one closed place (i.e.,
1976          a place and all its conjugates) which is represented by an intvec
1977          of size two, the first entry is the degree of the place (in
1978          particular, it tells the local ring where to find the data
1979          describing one representative of the closed place), and the
1980          second one is the position of those data in the lists POINTS, etc.,
1981          inside this local ring.@*
1982          In the intvec L[4] (conductor) the i-th entry corresponds to the
1983          i-th entry in the list of places L[3].
1984
1985          With no optional arguments, the conductor is computed by
1986          local invariants of the singularities; otherwise it is computed
1987          by the Dedekind formula. @*
1988          An affine point is represented by a list P where P[1] is std
1989          of a prime ideal and P[2] is an intvec containing the position
1990          of the places above P in the list of closed places L[3]. @*
1991          If the point is at infinity, P[1] is a homogeneous irreducible
1992          polynomial in two variables.
1993
1994          If @code{printlevel>=0} additional comments are displayed (default:
1995          @code{printlevel=0}).
1996KEYWORDS: Hamburger-Noether expansions; adjunction divisor
1997SEE ALSO: closed_points, NSplaces
1998EXAMPLE:  example Adj_div; shows an example
1999"
2000{
2001  // computes the adjunction divisor and the genus of a (singular) plane curve
2002  // as a byproduct, it computes all the singular points with the corresponding
2003  //      places and the genus of the curve
2004  // the adjunction divisor is stored in an intvec
2005  // also the non-singular places at infinity are computed
2006  // returns a list with all the computed data
2007  if (char(basering)==0)
2008  {
2009    ERROR("Base field not implemented");
2010  }
2011  if (npars(basering)>0)
2012  {
2013    ERROR("Base field not implemented");
2014  }
2015  intvec opgt=option(get);
2016  option(redSB);
2017  def Base_R=basering;
2018  list CURVE=curve(f);
2019  def aff_r=CURVE[1];
2020  def Proj_R=CURVE[2];
2021  int degX=CURVE[3];
2022  int genusX=(degX-1)*(degX-2);
2023  genusX = genusX div 2;
2024  intvec iivv=degX,genusX;
2025  intvec Conductor;
2026  setring aff_r;
2027  dbprint(printlevel+1,"Computing affine singular points ... ");
2028  list Aff_SPoints=Aff_SL(Aff_SLocus);
2029  int s=size(Aff_SPoints);
2030  if (s>0)
2031  {
2032    export(Aff_SPoints);
2033  }
2034  dbprint(printlevel+1,"Computing all points at infinity ... ");
2035  list Inf_Points=inf_P(CHI);
2036  export(Inf_Points);
2037  list update_CURVE=list();
2038  update_CURVE[1]=list();
2039  update_CURVE[1][1]=aff_r;
2040  update_CURVE[1][2]=Proj_R;
2041  update_CURVE[2]=iivv;
2042  update_CURVE[3]=list();
2043  update_CURVE[4]=Conductor;
2044  update_CURVE[5]=list();
2045  int i;
2046  intvec pP=0,0,0;
2047  if (size(#)==0)
2048  {
2049    dbprint(printlevel+1,"Computing affine singular places ... ");
2050    if (s>0)
2051    {
2052      for (i=1;i<=s;i=i+1)
2053      {
2054        pP[3]=i;
2055        update_CURVE=place(pP,1,update_CURVE);
2056      }
2057    }
2058    dbprint(printlevel+1,"Computing singular places at infinity ... ");
2059    s=size(Inf_Points[1]);
2060    if (s>0)
2061    {
2062      pP[2]=1;
2063      for (i=1;i<=s;i=i+1)
2064      {
2065        pP[3]=i;
2066        update_CURVE=place(pP,1,update_CURVE);
2067      }
2068    }
2069  }
2070  else
2071  {
2072    dbprint(printlevel+1,"Computing affine singular places ... ");
2073    s=size(Aff_SPoints);
2074    if (s>0)
2075    {
2076      for (i=1;i<=s;i=i+1)
2077      {
2078        pP[3]=i;
2079        update_CURVE=place(pP,2,update_CURVE);
2080      }
2081    }
2082    dbprint(printlevel+1,"Computing singular places at infinity ... ");
2083    s=size(Inf_Points[1]);
2084    if (s>0)
2085    {
2086      pP[2]=1;
2087      for (i=1;i<=s;i=i+1)
2088      {
2089        pP[3]=i;
2090        update_CURVE=place(pP,2,update_CURVE);
2091      }
2092    }
2093  }
2094  dbprint(printlevel+1,"Computing non-singular places at infinity ... ");
2095  s=size(Inf_Points[2]);
2096  if (s>0)
2097  {
2098    pP[2]=2;
2099    for (i=1;i<=s;i=i+1)
2100    {
2101      pP[3]=i;
2102      update_CURVE=place(pP,0,update_CURVE);
2103    }
2104  }
2105  dbprint(printlevel+1,"Adjunction divisor computed successfully");
2106  list Places=update_CURVE[3];
2107  Conductor=update_CURVE[4];
2108  genusX = genusX - (deg_D(Conductor,Places) div 2);
2109  update_CURVE[2][2]=genusX;
2110  setring Base_R;
2111  dbprint(printlevel+1," ");
2112  dbprint(printlevel+2,"The genus of the curve is "+string(genusX));
2113  option(set,opgt);
2114  return(update_CURVE);
2115}
2116example
2117{
2118  "EXAMPLE:";  echo = 2;
2119  int plevel=printlevel;
2120  printlevel=-1;
2121  ring s=2,(x,y),lp;
2122  list C=Adj_div(y9+y8+xy6+x2y3+y2+x3);
2123  def aff_R=C[1][1];      // the affine ring
2124  setring aff_R;
2125  listvar(aff_R);         // data in the affine ring
2126  CHI;                    // affine equation of the curve
2127  Aff_SLocus;             // ideal of the affine singular locus
2128  Aff_SPoints[1];         // 1st affine singular point: (1:1:1), no.1
2129  Inf_Points[1];          // singular point(s) at infinity: (1:0:0), no.4
2130  Inf_Points[2];          // list of non-singular points at infinity
2131  //
2132  pause("press RETURN");
2133  def proj_R=C[1][2];     // the projective ring
2134  setring proj_R;
2135  CHI;                    // projective equation of the curve
2136  C[2][1];                // degree of the curve
2137  C[2][2];                // genus of the curve
2138  C[3];                   // list of computed places
2139  C[4];                   // adjunction divisor (all points are singular!)
2140  //
2141  pause("press RETURN");
2142  // we look at the place(s) of degree 2 by changing to the ring
2143  C[5][2][1];
2144  def S(2)=C[5][2][1];
2145  setring S(2);
2146  POINTS;                // base point(s) of place(s) of degree 2: (1:a:1)
2147  LOC_EQS;               // local equation(s)
2148  PARAMETRIZATIONS;      // parametrization(s) and exactness
2149  BRANCHES;              // Hamburger-Noether development
2150  printlevel=plevel;
2151}
2152///////////////////////////////////////////////////////////////////////////////
2153
2154proc NSplaces (int h,list CURVE)
2155"USAGE:   NSplaces( h, CURVE ), where h is an integer and CURVE is a list
2156RETURN:   list L with updated data of CURVE after computing all
2157          points up to degree H+h (H the maximum degree of the previously
2158          computed places): @*
2159   @format
2160   in L[1][1]: (affine ring) lists Aff_Points(d) with affine non-singular
2161               points of degree d (if non-empty)
2162   in L[3]:    the newly computed closed places are added,
2163   in L[5]:    local rings created/updated to store (repres. of) new places.
2164   @end format
2165          See @ref{Adj_div} for a description of the entries in L.
2166NOTE:     The list_expression should be the output of the procedure Adj_div.@*
2167          If @code{printlevel>=0} additional comments are displayed (default:
2168          @code{printlevel=0}).
2169SEE ALSO: closed_points, Adj_div
2170EXAMPLE:  example NSplaces; shows an example
2171"
2172{
2173  // computes all the non-singular closed places with degree up to a certain
2174  // bound; this bound is the maximum degree of an existing singular place or
2175  // non-singular place at infinity plus an increment h>=0 which is given as
2176  // input
2177  // creates lists of points and the corresponding places
2178  // list CURVE must be the output of the procedure "Adj_div"
2179  // warning : if h<0 then it will be replaced by h=0
2180  intvec opgt=option(get);
2181  option(redSB);
2182  def Base_R=basering;
2183  def aff_r=CURVE[1][1];
2184  int M=size(CURVE[5]);
2185  if (h>0)
2186  {
2187    M=M+h;
2188  }
2189  list update_CURVE=CURVE;
2190  int i,j,s;
2191  setring aff_r;
2192  intvec pP=1,0,0;
2193  for (i=1;i<=M;i=i+1)
2194  {
2195    dbprint(printlevel+1,"Computing non-singular affine places of degree "
2196                           +string(i)+" ... ");
2197    if (defined(Aff_Points(i))==0)
2198    {
2199      list Aff_Points(i)=closed_points_deg(CHI,i,Aff_SLocus);
2200      s=size(Aff_Points(i));
2201      if (s>0)
2202      {
2203        export(Aff_Points(i));
2204        pP[2]=i;
2205        for (j=1;j<=s;j=j+1)
2206        {
2207          pP[3]=j;
2208          update_CURVE=place(pP,0,update_CURVE);
2209        }
2210      }
2211    }
2212  }
2213  setring Base_R;
2214  option(set,opgt);
2215  return(update_CURVE);
2216}
2217example
2218{
2219  "EXAMPLE:";  echo = 2;
2220  int plevel=printlevel;
2221  printlevel=-1;
2222  ring s=2,(x,y),lp;
2223  list C=Adj_div(x3y+y3+x);
2224  // The list of computed places:
2225  C[3];
2226  // create places up to degree 1+3
2227  list L=NSplaces(3,C);
2228  // The list of computed places is now:
2229  L[3];
2230  // e.g., affine non-singular points of degree 4 :
2231  def aff_r=L[1][1];
2232  setring aff_r;
2233  Aff_Points(4);
2234  // e.g., base point of the 1st place of degree 4 :
2235  def S(4)=L[5][4][1];
2236  setring S(4);
2237  POINTS[1];
2238  printlevel=plevel;
2239}
2240///////////////////////////////////////////////////////////////////////////////
2241
2242// ** SPECIAL PROCEDURES FOR LINEAR ALGEBRA **
2243
2244static proc Ker (matrix A)
2245{
2246  // warning : "lp" ordering is necessary
2247  intvec opgt=option(get);
2248  option(redSB);
2249  matrix M=transpose(syz(A));
2250  option(set,opgt);
2251  return(M);
2252}
2253///////////////////////////////////////////////////////////////////////////////
2254static proc get_NZsol (matrix A)
2255{
2256  matrix sol=Ker(A);
2257  return(submat(sol,1..1,1..ncols(sol)));
2258}
2259///////////////////////////////////////////////////////////////////////////////
2260static proc supplement (matrix W,matrix V)
2261"USAGE:     supplement(W,V), where W,V are matrices of numbers such that the
2262            vector space generated by the rows of W is contained in that
2263            generated by the rows of V
2264RETURN:     matrix whose rows generate a supplementary vector space of W in V,
2265            or a zero row-matrix if <W>==<V>
2266NOTE:       W,V must be given with maximal rank
2267"
2268{
2269  // W and V represent independent sets of vectors and <W> is assumed to be
2270  // contained in <V>
2271  // computes matrix S whose rows are l.i. vectors s.t. <W> union <S> is a
2272  // basis of <V>
2273  // careful : the size of all vectors is assumed to be the same but it is
2274  //  not checked and neither the linear independence of the vectors is checked
2275  // the trivial case W=0 is not covered by this procedure (and thus V<>0)
2276  // if <W>=<V> then a zero row-matrix is returned
2277  // warning : option(redSB) must be set in advance
2278  int n1=nrows(W);
2279  int n2=nrows(V);
2280  int s=n2-n1;
2281  if (s==0)
2282  {
2283    int n=ncols(W);
2284    matrix HH[1][n];
2285    return(HH);
2286  }
2287  matrix H=transpose(lift(transpose(V),transpose(W)));
2288  H=supplem(H);
2289  return(H*V);
2290}
2291///////////////////////////////////////////////////////////////////////////////
2292static proc supplem (matrix M)
2293"USAGE:      suplement(M), where M is a matrix of numbers with maximal rank
2294RETURN:     matrix whose rows generate a supplementary vector space of <M> in
2295            k^n, where k is the base field and n is the number of columns
2296SEE ALSO:   supplement
2297NOTE:       The rank r is assumed to be 1<r<n.
2298"
2299{
2300  // warning : the linear independence of the rows is not checked
2301  int r=nrows(M);
2302  int n=ncols(M);
2303  int s=n-r;
2304  matrix A=M;
2305  matrix supl[s][n];
2306  int counter=0;
2307  int h=r+1;
2308  int i;
2309  for (i=1;i<=n;i=i+1)
2310  {
2311    matrix TT[1][n];
2312    TT[1,i]=1;
2313    A=transpose(concat(transpose(A),transpose(TT)));
2314    r=mat_rank(A);
2315    if (r==h)
2316    {
2317      h=h+1;
2318      counter=counter+1;
2319      supl=transpose(concat(transpose(supl),transpose(TT)));
2320      if (counter==s)
2321      {
2322        break;
2323      }
2324    }
2325    kill(TT);
2326  }
2327  supl=transpose(compress(transpose(supl)));
2328  return(supl);
2329}
2330///////////////////////////////////////////////////////////////////////////////
2331static proc mat_rank (matrix A)
2332{
2333  // warning : "lp" ordering is necessary
2334  intvec opgt=option(get);
2335  option(redSB);
2336  int r=size(std(module(transpose(A))));
2337  option(set,opgt);
2338  return(r);
2339}
2340///////////////////////////////////////////////////////////////////////////////
2341
2342// ***************************************************************
2343// * PROCEDURES FOR INTERPOLATION, INTERSECTION AND EXTRA PLACES *
2344// ***************************************************************
2345
2346static proc estim_n (intvec Dplus,int dgX,list PL)
2347{
2348  // computes an estimate for the degree n in the Brill-Noether algorithm
2349  int estim=2*deg_D(Dplus,PL)+dgX*(dgX-3);
2350  estim=estim div (2*dgX);
2351  estim=estim+1;
2352  if (estim<dgX)
2353  {
2354    estim=dgX;
2355  }
2356  return(estim);
2357}
2358///////////////////////////////////////////////////////////////////////////////
2359static proc nforms (int n)
2360{
2361  // computes the list of all homogeneous monomials of degree n>=0
2362  // exports ideal nFORMS(n) whose generators are ranged with lp order
2363  //   in Proj_R and returns size(nFORMS(n))
2364  // warning : it is supposed to be called inside Proj_R
2365  // if n<=0 then nFORMS(0) is "computed fast"
2366  ideal nFORMS(n);
2367  int N;
2368  if (n>0)
2369  {
2370    N=(n+1)*(n+2);
2371    N=N div 2;
2372    N=N+1;
2373    int i,j,k;
2374    for (i=0;i<=n;i=i+1)
2375    {
2376      for (j=0;j<=n-i;j=j+1)
2377      {
2378        k=k+1;
2379        nFORMS(n)[N-k]=x^i*y^j*z^(n-i-j);
2380      }
2381    }
2382    export(nFORMS(n));
2383  }
2384  else
2385  {
2386    N=2;
2387    nFORMS(0)=1;
2388    export(nFORMS(0));
2389  }
2390  return(N-1);
2391}
2392///////////////////////////////////////////////////////////////////////////////
2393static proc nmultiples (int n,int dgX,poly f)
2394{
2395  // computes a basis of the space of forms of degree n which are multiple of
2396  //     CHI
2397  // returns a matrix whose rows are the coordinates (related to nFORMS(n))
2398  //     of such a basis
2399  // warning : it is supposed to be called inside Proj_R
2400  // warning : nFORMS(n) is created in the way, together with nFORMS(n-degX)
2401  // warning : n must be greater or equal than the degree of the curve
2402  if (defined(nFORMS(n))==0)
2403  {
2404    dbprint(printlevel+1,string(nforms(n)));
2405  }
2406  int m=n-dgX;
2407  if (defined(nFORMS(m))==0)
2408  {
2409    int k=nforms(m);
2410  }
2411  else
2412  {
2413    int k=size(nFORMS(m));
2414  }
2415  ideal nmults;
2416  int i;
2417  for (i=1;i<=k;i=i+1)
2418  {
2419    nmults[i]=f*nFORMS(m)[i];
2420  }
2421  return(transpose(lift(nFORMS(n),nmults)));
2422}
2423///////////////////////////////////////////////////////////////////////////////
2424static proc interpolating_forms (intvec D,int n,list CURVE)
2425{
2426  // computes a vector basis of the space of forms of degree n whose
2427  //     intersection divisor with the curve is greater or equal than D>=0
2428  // the procedure is supposed to be called inside the ring Proj_R and
2429  //     assumes that the forms nFORMS(n) are previously computed
2430  // the output is a matrix whose rows are the coordinates in nFORMS(n) of
2431  //     such a basis
2432  // remark : the support of D may contain "extra" places
2433  def BR=basering;
2434  def aff_r=CURVE[1][1];
2435  int N=size(nFORMS(n));
2436  matrix totalM[1][N];
2437  int s=size(D);
2438  list Places=CURVE[3];
2439  int NPls=size(Places);
2440  int i,j,k,kk,id,ip,RR,ordx,ordy,nr,NR;
2441  if (s<=NPls)
2442  {
2443    for (i=1;i<=s;i=i+1)
2444    {
2445      if (D[i]>0)
2446      {
2447        id=Places[i][1];
2448        ip=Places[i][2];
2449        RR=D[i];
2450        def SS=CURVE[5][id][1];
2451        setring SS;
2452        poly xt=PARAMETRIZATIONS[ip][1][1];
2453        poly yt=PARAMETRIZATIONS[ip][1][2];
2454        ordx=PARAMETRIZATIONS[ip][2][1];
2455        ordy=PARAMETRIZATIONS[ip][2][2];
2456        nr=ncols(BRANCHES[ip][1]);
2457        if ( ordx>-1 || ordy>-1 )
2458        {
2459          while ( ( RR>ordx && ordx>-1 ) || ( RR>ordy && ordy>-1 ) )
2460          {
2461            BRANCHES[ip]=extdevelop(BRANCHES[ip],2*nr);
2462            nr=ncols(BRANCHES[ip][1]);
2463            PARAMETRIZATIONS[ip]=param(BRANCHES[ip],0);
2464            xt=PARAMETRIZATIONS[ip][1][1];
2465            yt=PARAMETRIZATIONS[ip][1][2];
2466            ordx=PARAMETRIZATIONS[ip][2][1];
2467            ordy=PARAMETRIZATIONS[ip][2][2];
2468          }
2469        }
2470        if (POINTS[ip][3]==number(1))
2471        {
2472          number A=POINTS[ip][1];
2473          number B=POINTS[ip][2];
2474          map Mt=BR,A+xt,B+yt,1;
2475          kill(A,B);
2476        }
2477        else
2478        {
2479          if (POINTS[ip][2]==number(1))
2480          {
2481            number A=POINTS[ip][1];
2482            map Mt=BR,A+xt,1,yt;
2483            kill(A);
2484          }
2485          else
2486          {
2487            map Mt=BR,1,xt,yt;
2488          }
2489        }
2490        ideal nFORMS(n)=Mt(nFORMS(n));
2491        // rewrite properly the above conditions to obtain the local equations
2492        matrix partM[RR][N];
2493        matrix auxMC=coeffs(nFORMS(n),t);
2494        NR=nrows(auxMC);
2495        if (RR<=NR)
2496        {
2497          for (j=1;j<=RR;j=j+1)
2498          {
2499            for (k=1;k<=N;k=k+1)
2500            {
2501              partM[j,k]=number(auxMC[j,k]);
2502            }
2503          }
2504        }
2505        else
2506        {
2507          for (j=1;j<=NR;j=j+1)
2508          {
2509            for (k=1;k<=N;k=k+1)
2510            {
2511              partM[j,k]=number(auxMC[j,k]);
2512            }
2513          }
2514          for (j=NR+1;j<=RR;j=j+1)
2515          {
2516            for (k=1;k<=N;k=k+1)
2517            {
2518              partM[j,k]=number(0);
2519            }
2520          }
2521        }
2522        matrix localM=partM;
2523        matrix conjM=partM;
2524        for (j=2;j<=id;j=j+1)
2525        {
2526          conjM=Frobenius(conjM,1);
2527          localM=transpose(concat(transpose(localM),transpose(conjM)));
2528        }
2529        localM=gauss_row(localM);
2530        setring BR;
2531        totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM))));
2532        totalM=transpose(compress(transpose(totalM)));
2533        setring SS;
2534        kill(xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM);
2535        setring BR;
2536        kill(SS);
2537      }
2538    }
2539  }
2540  else
2541  {
2542    // distinguish between "standard" places and "extra" places
2543    for (i=1;i<=NPls;i=i+1)
2544    {
2545      if (D[i]>0)
2546      {
2547        id=Places[i][1];
2548        ip=Places[i][2];
2549        RR=D[i];
2550        def SS=CURVE[5][id][1];
2551        setring SS;
2552        poly xt=PARAMETRIZATIONS[ip][1][1];
2553        poly yt=PARAMETRIZATIONS[ip][1][2];
2554        ordx=PARAMETRIZATIONS[ip][2][1];
2555        ordy=PARAMETRIZATIONS[ip][2][2];
2556        nr=ncols(BRANCHES[ip][1]);
2557        if ( ordx>-1 || ordy>-1 )
2558        {
2559          while ( ( RR>ordx && ordx>-1 ) || ( RR>ordy && ordy>-1 ) )
2560          {
2561            BRANCHES[ip]=extdevelop(BRANCHES[ip],2*nr);
2562            nr=ncols(BRANCHES[ip][1]);
2563            PARAMETRIZATIONS[ip]=param(BRANCHES[ip],0);
2564            xt=PARAMETRIZATIONS[ip][1][1];
2565            yt=PARAMETRIZATIONS[ip][1][2];
2566            ordx=PARAMETRIZATIONS[ip][2][1];
2567            ordy=PARAMETRIZATIONS[ip][2][2];
2568          }
2569        }
2570        if (POINTS[ip][3]==number(1))
2571        {
2572          number A=POINTS[ip][1];
2573          number B=POINTS[ip][2];
2574          map Mt=BR,A+xt,B+yt,1;
2575          kill(A,B);
2576        }
2577        else
2578        {
2579          if (POINTS[ip][2]==number(1))
2580          {
2581            number A=POINTS[ip][1];
2582            map Mt=BR,A+xt,1,yt;
2583            kill(A);
2584          }
2585          else
2586          {
2587            map Mt=BR,1,xt,yt;
2588          }
2589        }
2590        ideal nFORMS(n)=Mt(nFORMS(n));
2591        // rewrite properly the above conditions to obtain the local equations
2592        matrix partM[RR][N];
2593        matrix auxMC=coeffs(nFORMS(n),t);
2594        NR=nrows(auxMC);
2595        if (RR<=NR)
2596        {
2597          for (j=1;j<=RR;j=j+1)
2598          {
2599            for (k=1;k<=N;k=k+1)
2600            {
2601              partM[j,k]=number(auxMC[j,k]);
2602            }
2603          }
2604        }
2605        else
2606        {
2607          for (j=1;j<=NR;j=j+1)
2608          {
2609            for (k=1;k<=N;k=k+1)
2610            {
2611              partM[j,k]=number(auxMC[j,k]);
2612            }
2613          }
2614          for (j=NR+1;j<=RR;j=j+1)
2615          {
2616            for (k=1;k<=N;k=k+1)
2617            {
2618              partM[j,k]=number(0);
2619            }
2620          }
2621        }
2622        matrix localM=partM;
2623        matrix conjM=partM;
2624        for (j=2;j<=id;j=j+1)
2625        {
2626          conjM=Frobenius(conjM,1);
2627          localM=transpose(concat(transpose(localM),transpose(conjM)));
2628        }
2629        localM=gauss_row(localM);
2630        setring BR;
2631        totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM))));
2632        totalM=transpose(compress(transpose(totalM)));
2633        setring SS;
2634        kill(xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM);
2635        setring BR;
2636        kill(SS);
2637      }
2638    }
2639    k=s-NPls;
2640    int l;
2641    for (i=1;i<=k;i=i+1)
2642    {
2643      // in this case D[NPls+i]>0 is assumed to be true during the
2644      // Brill-Noether algorithm
2645      RR=D[NPls+i];
2646      setring aff_r;
2647      def SS=@EXTRA@[2][i];
2648      def extra_dgs=@EXTRA@[3];
2649      setring SS;
2650      poly xt=PARAMETRIZATION[1][1];
2651      poly yt=PARAMETRIZATION[1][2];
2652      ordx=PARAMETRIZATION[2][1];
2653      ordy=PARAMETRIZATION[2][2];
2654      nr=ncols(BRANCH[1]);
2655      if ( ordx>-1 || ordy>-1 )
2656      {
2657        while ( ( RR>ordx && ordx>-1 ) || ( RR>ordy && ordy>-1 ) )
2658        {
2659          BRANCH[1]=extdevelop(BRANCH,2*nr);
2660          nr=ncols(BRANCH[1]);
2661          PARAMETRIZATION=param(BRANCH,0);
2662          xt=PARAMETRIZATION[1][1];
2663          yt=PARAMETRIZATION[1][2];
2664          ordx=PARAMETRIZATION[2][1];
2665          ordy=PARAMETRIZATION[2][2];
2666        }
2667      }
2668      number A=POINT[1];
2669      number B=POINT[2];
2670      map Mt=BR,A+xt,B+yt,1;
2671      kill(A,B);
2672      ideal nFORMS(n)=Mt(nFORMS(n));
2673      // rewrite properly the above conditions to obtain the local equations
2674      matrix partM[RR][N];
2675      matrix auxMC=coeffs(nFORMS(n),t);
2676      NR=nrows(auxMC);
2677      if (RR<=NR)
2678      {
2679        for (j=1;j<=RR;j=j+1)
2680        {
2681          for (kk=1;kk<=N;kk=kk+1)
2682          {
2683            partM[j,kk]=number(auxMC[j,kk]);
2684          }
2685        }
2686      }
2687      else
2688      {
2689        for (j=1;j<=NR;j=j+1)
2690        {
2691          for (kk=1;kk<=N;kk=kk+1)
2692          {
2693            partM[j,kk]=number(auxMC[j,kk]);
2694          }
2695        }
2696        for (j=NR+1;j<=RR;j=j+1)
2697        {
2698          for (kk=1;kk<=N;kk=kk+1)
2699          {
2700            partM[j,kk]=number(0);
2701          }
2702        }
2703      }
2704      matrix localM=partM;
2705      matrix conjM=partM;
2706      l=extra_dgs[i];
2707      for (j=2;j<=l;j=j+1)
2708      {
2709        conjM=Frobenius(conjM,1);
2710        localM=transpose(concat(transpose(localM),transpose(conjM)));
2711      }
2712      localM=gauss_row(localM);
2713      setring BR;
2714      totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM))));
2715      totalM=transpose(compress(transpose(totalM)));
2716      setring SS;
2717      kill(xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM);
2718      setring BR;
2719      kill(SS);
2720    }
2721  }
2722  return(Ker(totalM));
2723}
2724///////////////////////////////////////////////////////////////////////////////
2725static proc local_IN (poly h,int m)
2726{
2727  // computes the intersection number of h and the curve CHI at a certain place
2728  // returns a list with the intersection number and the "leading coefficient"
2729  // the procedure must be called inside a local ring, h must be a local
2730  //     equation with respect to the desired place, and m indicates the
2731  //     number of place inside that local ring, containing lists
2732  //     POINT(S)/BRANCH(ES)/PARAMETRIZATION(S) when m=0 an "extra place" is
2733  //     considered
2734  def BR=basering;
2735  if (m>0)
2736  {
2737    int nr=ncols(BRANCHES[m][1]);
2738    poly xt=PARAMETRIZATIONS[m][1][1];
2739    poly yt=PARAMETRIZATIONS[m][1][2];
2740    int ordx=PARAMETRIZATIONS[m][2][1];
2741    int ordy=PARAMETRIZATIONS[m][2][2];
2742    map phi=BR,xt,yt,1;
2743    poly ht=phi(h);
2744    int inum=mindeg(ht);
2745    if ( ordx>-1 || ordy>-1 )
2746    {
2747      while ( ( inum>ordx && ordx>-1 ) || ( inum>ordy && ordy>-1 ) )
2748      {
2749        BRANCHES[m]=extdevelop(BRANCHES[m],2*nr);
2750        nr=ncols(BRANCHES[m][1]);
2751        PARAMETRIZATIONS[m]=param(BRANCHES[m],0);
2752        xt=PARAMETRIZATIONS[m][1][1];
2753        yt=PARAMETRIZATIONS[m][1][2];
2754        ordx=PARAMETRIZATIONS[m][2][1];
2755        ordy=PARAMETRIZATIONS[m][2][2];
2756        ht=phi(h);
2757        inum=mindeg(ht);
2758      }
2759    }
2760  }
2761  else
2762  {
2763    int nr=ncols(BRANCH[1]);
2764    poly xt=PARAMETRIZATION[1][1];
2765    poly yt=PARAMETRIZATION[1][2];
2766    int ordx=PARAMETRIZATION[2][1];
2767    int ordy=PARAMETRIZATION[2][2];
2768    map phi=basering,xt,yt,1;
2769    poly ht=phi(h);
2770    int inum=mindeg(ht);
2771    if ( ordx>-1 || ordy>-1 )
2772    {
2773      while ( ( inum>ordx && ordx>-1 ) || ( inum>ordy && ordy>-1 ) )
2774      {
2775        BRANCH=extdevelop(BRANCH,2*nr);
2776        nr=ncols(BRANCH[1]);
2777        PARAMETRIZATION=param(BRANCH,0);
2778        xt=PARAMETRIZATION[1][1];
2779        yt=PARAMETRIZATION[1][2];
2780        ordx=PARAMETRIZATION[2][1];
2781        ordy=PARAMETRIZATION[2][2];
2782        ht=phi(h);
2783        inum=mindeg(ht);
2784      }
2785    }
2786  }
2787  list answer=list();
2788  answer[1]=inum;
2789  number AA=number(coeffs(ht,t)[inum+1,1]);
2790  answer[2]=AA;
2791  return(answer);
2792}
2793///////////////////////////////////////////////////////////////////////////////
2794static proc extra_place (ideal P)
2795{
2796  // computes the "rational" place which is defined over a (closed) "extra"
2797  //     point
2798  // an "extra" point will be necessarily affine, non-singular and non-rational
2799  // creates : a specific local ring to deal with the (unique) place above it
2800  // returns : list with the above local ring and the degree of the point/place
2801  // warning : the procedure must be called inside the affine ring aff_r
2802  def base_r=basering;
2803  int ext=deg(P[1]);
2804  poly aux=subst(P[2],y,1);
2805  ext=ext*deg(aux);
2806  // P is assumed to be a std. resp. "(x,y),lp" and thus P[1] depends only
2807  // on "y"
2808  if (deg(P[1])==1)
2809  {
2810    // the point is non-rational but the second component needs no field
2811    // extension
2812    number B=-number(subst(P[1],y,0));
2813    poly aux2=subst(P[2],y,B);
2814    // careful : the parameter will be called "a" anyway
2815    ring ES=(char(basering),a),(x,y,t),ls;
2816    map psi=base_r,a,0;
2817    minpoly=number(psi(aux2));
2818    kill(psi);
2819    number A=a;
2820    number B=imap(base_r,B);
2821  }
2822  else
2823  {
2824    if (deg(subst(P[2],y,1))==1)
2825    {
2826      // the point is non-rational but the needed minpoly is just P[1]
2827      // careful : the parameter will be called "a" anyway
2828      poly P1=P[1];
2829      poly P2=P[2];
2830      ring ES=(char(basering),a),(x,y,t),ls;
2831      map psi=base_r,0,a;
2832      minpoly=number(psi(P1));
2833      kill(psi);
2834      poly aux1=imap(base_r,P2);
2835      poly aux2=subst(aux1,y,a);
2836      number A=-number(subst(aux2,x,0));
2837      number B=a;
2838    }
2839    else
2840    {
2841      // this is the most complicated case of non-rational point
2842      poly P1=P[1];
2843      poly P2=P[2];
2844      int p=char(basering);
2845      int Q=p^ext;
2846      ring aux_r=(Q,a),(x,y,t),ls;
2847      string minpoly_string=string(minpoly);
2848      ring ES=(char(basering),a),(x,y,t),ls;
2849      execute("minpoly="+minpoly_string+";");
2850      poly P_1=imap(base_r,P1);
2851      poly P_2=imap(base_r,P2);
2852      ideal factors1=factorize(P_1,1);
2853      number B=-number(subst(factors1[1],y,0));
2854      poly P_0=subst(P_2,y,B);
2855      ideal factors2=factorize(P_0,1);
2856      number A=-number(subst(factors2[1],x,0));
2857      kill(aux_r);
2858    }
2859  }
2860  list POINT=list();
2861  POINT[1]=A;
2862  POINT[2]=B;
2863  export(POINT);
2864  map phi=base_r,x+A,y+B;
2865  poly LOC_EQ=phi(CHI);
2866  kill(A,B,phi);
2867  list L@HNE=essdevelop(LOC_EQ)[1];
2868  export(L@HNE);
2869  int m=nrows(L@HNE[1]);
2870  int n=ncols(L@HNE[1]);
2871  intvec Li2=L@HNE[2];
2872  int Li3=L@HNE[3];
2873  setring ES;
2874  string newa=subfield(HNEring);
2875  poly paux=importdatum(HNEring,"L@HNE[4]",newa);
2876  matrix Maux[m][n];
2877  int i,j;
2878  for (i=1;i<=m;i=i+1)
2879  {
2880    for (j=1;j<=n;j=j+1)
2881    {
2882      Maux[i,j]=importdatum(HNEring,"L@HNE[1]["+string(i)+","+
2883                                                string(j)+"]",newa);
2884    }
2885  }
2886  list BRANCH=list();
2887  BRANCH[1]=Maux;
2888  BRANCH[2]=Li2;
2889  BRANCH[3]=Li3;
2890  BRANCH[4]=paux;
2891  export(BRANCH);
2892  list PARAMETRIZATION=param(BRANCH,0);
2893  export(PARAMETRIZATION);
2894  kill(HNEring);
2895  setring base_r;
2896  list answer=list();
2897  answer[1]=ES;
2898  answer[2]=ext;
2899  kill(ES);
2900  return(answer);
2901}
2902///////////////////////////////////////////////////////////////////////////////
2903static proc intersection_div (poly H,list CURVE)
2904"USAGE:     intersection_div(H,CURVE), where H is a homogeneous polynomial
2905            in ring Proj_R=p,(x,y,z),lp and CURVE is the list of data for
2906            the given curve
2907CREATE:     new places which had not been computed in advance if necessary
2908            those places are stored each one in a local ring where you find
2909            lists POINT,BRANCH,PARAMETRIZATION for the place living in that
2910            ring; the degree of the point/place in such a ring is stored in
2911            an intvec, and the base points in the remaining list
2912            Everything is exported in a list @EXTRA@ inside the ring
2913            aff_r=CURVE[1][1] and returned with the updated CURVE
2914RETURN:     list with the intersection divisor (intvec) between the underlying
2915            curve and H=0, and the list CURVE updated
2916SEE ALSO:   Adj_div, NSplaces, closed_points
2917NOTE:       The procedure must be called from the ring Proj_R=CURVE[1][2]
2918            (projective).
2919            If @EXTRA@ already exists, the new places are added to the
2920            previous data.
2921"
2922{
2923  // computes the intersection divisor of H and the curve CHI
2924  // returns a list with (possibly) "extra places" and it must be called
2925  //     inside Proj_R
2926  // in case of extra places, some local rings ES(1) ... ES(m) are created
2927  //     together with an integer list "extra_dgs" containing the degrees of
2928  //     those places
2929  intvec opgt=option(get);
2930  option(redSB);
2931  intvec interdiv;
2932  def BRing=basering;
2933  int Tz1=deg(H);
2934  list Places=CURVE[3];
2935  int N=size(Places);
2936  def aff_r=CURVE[1][1];
2937  setring aff_r;
2938  if (defined(@EXTRA@)==0)
2939  {
2940    list @EXTRA@=list();
2941    list EP=list();
2942    list ES=list();
2943    list extra_dgs=list();
2944  }
2945  else
2946  {
2947    list EP=@EXTRA@[1];
2948    list ES=@EXTRA@[2];
2949    list extra_dgs=@EXTRA@[3];
2950  }
2951  int NN=size(extra_dgs);
2952  int counterEPl=0;
2953  setring BRing;
2954  poly h=subst(H,z,1);
2955  int Tz2=deg(h);
2956  int Tz3=Tz1-Tz2;
2957  int i,j,k,l,m,n,s,np,NP,I_N;
2958  if (Tz3==0)
2959  {
2960    // if this still does not work -> try always with ALL points in
2961    // Inf_Points !!!!
2962    poly Hinf=subst(H,z,0);
2963    setring aff_r;
2964    // compute the points at infinity of H and see which of them are in
2965    // Inf_Points
2966    poly h=imap(BRing,h);
2967    poly hinf=imap(BRing,Hinf);
2968    ideal pinf=factorize(hinf,1);
2969    list TIP=Inf_Points[1]+Inf_Points[2];
2970    s=size(TIP);
2971    NP=size(pinf);
2972    for (i=1;i<=NP;i=i+1)
2973    {
2974      for (j=1;j<=s;j=j+1)
2975      {
2976        if (pinf[i]==TIP[j][1])
2977        {
2978          np=size(TIP[j][2]);
2979          for (k=1;k<=np;k=k+1)
2980          {
2981            n=TIP[j][2][k];
2982            l=Places[n][1];
2983            m=Places[n][2];
2984            def SS=CURVE[5][l][1];
2985            setring SS;
2986            // local equation h of H
2987            if (POINTS[m][2]==number(1))
2988            {
2989              number A=POINTS[m][1];
2990              map psi=BRing,x+A,1,y;
2991              kill(A);
2992            }
2993            else
2994            {
2995              map psi=BRing,1,x,y;
2996            }
2997            poly h=psi(H);
2998            I_N=local_IN(h,m)[1];
2999            interdiv[n]=I_N;
3000            kill(h,psi);
3001            setring aff_r;
3002            kill(SS);
3003          }
3004          break;
3005        }
3006      }
3007    }
3008    kill(hinf,pinf);
3009  }
3010  else
3011  {
3012    // H is a multiple of z and hence all the points in Inf_Points intersect
3013    // with H
3014    setring aff_r;
3015    poly h=imap(BRing,h);
3016    list TIP=Inf_Points[1]+Inf_Points[2];
3017    s=size(TIP);
3018    for (j=1;j<=s;j=j+1)
3019    {
3020      np=size(TIP[j][2]);
3021      for (k=1;k<=np;k=k+1)
3022      {
3023        n=TIP[j][2][k];
3024        l=Places[n][1];
3025        m=Places[n][2];
3026        def SS=CURVE[5][l][1];
3027        setring SS;
3028        // local equation h of H
3029        if (POINTS[m][2]==number(1))
3030        {
3031          number A=POINTS[m][1];
3032          map psi=BRing,x+A,1,y;
3033          kill(A);
3034        }
3035        else
3036        {
3037          map psi=BRing,1,x,y;
3038        }
3039        poly h=psi(H);
3040        I_N=local_IN(h,m)[1];
3041        interdiv[n]=I_N;
3042        kill(h,psi);
3043        setring aff_r;
3044        kill(SS);
3045      }
3046    }
3047  }
3048  // compute common affine points of H and CHI
3049  ideal CAL=h,CHI;
3050  CAL=std(CAL);
3051  if (CAL<>1)
3052  {
3053    list TAP=list();
3054    TAP=closed_points(CAL);
3055    NP=size(TAP);
3056    list auxP=list();
3057    int dP;
3058    for (i=1;i<=NP;i=i+1)
3059    {
3060      if (belongs(TAP[i],Aff_SLocus)==1)
3061      {
3062        // search the point in the list Aff_SPoints
3063        j=isInLP(TAP[i],Aff_SPoints);
3064        np=size(Aff_SPoints[j][2]);
3065        for (k=1;k<=np;k=k+1)
3066        {
3067          n=Aff_SPoints[j][2][k];
3068          l=Places[n][1];
3069          m=Places[n][2];
3070          def SS=CURVE[5][l][1];
3071          setring SS;
3072          // local equation h of H
3073          number A=POINTS[m][1];
3074          number B=POINTS[m][2];
3075          map psi=BRing,x+A,y+B,1;
3076          poly h=psi(H);
3077          I_N=local_IN(h,m)[1];
3078          interdiv[n]=I_N;
3079          kill(A,B,h,psi);
3080          setring aff_r;
3081          kill(SS);
3082        }
3083      }
3084      else
3085      {
3086        auxP=list();
3087        auxP[1]=TAP[i];
3088        dP=degree_P(auxP);
3089        if (defined(Aff_Points(dP))<>0)
3090        {
3091          // search the point in the list Aff_Points(dP)
3092          j=isInLP(TAP[i],Aff_Points(dP));
3093          n=Aff_Points(dP)[j][2][1];
3094          l=Places[n][1];
3095          m=Places[n][2];
3096          def SS=CURVE[5][l][1];
3097          setring SS;
3098          // local equation h of H
3099          number A=POINTS[m][1];
3100          number B=POINTS[m][2];
3101          map psi=BRing,x+A,y+B,1;
3102          poly h=psi(H);
3103          I_N=local_IN(h,m)[1];
3104          interdiv[n]=I_N;
3105          kill(A,B,h,psi);
3106          setring aff_r;
3107          kill(SS);
3108        }
3109        else
3110        {
3111          // previously check if it is an existing "extra place"
3112          j=isInLP(TAP[i],EP);
3113          if (j>0)
3114          {
3115            def SS=ES[j];
3116            setring SS;
3117            // local equation h of H
3118            number A=POINT[1];
3119            number B=POINT[2];
3120            map psi=BRing,x+A,y+B,1;
3121            poly h=psi(H);
3122            I_N=local_IN(h,0)[1];
3123            interdiv[N+j]=I_N;
3124            setring aff_r;
3125            kill(SS);
3126          }
3127          else
3128          {
3129            // then we must create a new "extra place"
3130            counterEPl=counterEPl+1;
3131            list EXTRA_PLACE=extra_place(TAP[i]);
3132            def SS=EXTRA_PLACE[1];
3133            ES[NN+counterEPl]=SS;
3134            extra_dgs[NN+counterEPl]=EXTRA_PLACE[2];
3135            EP[NN+counterEPl]=list();
3136            EP[NN+counterEPl][1]=TAP[i];
3137            EP[NN+counterEPl][2]=0;
3138            setring SS;
3139            // local equation h of H
3140            number A=POINT[1];
3141            number B=POINT[2];
3142            map psi=BRing,x+A,y+B,1;
3143            poly h=psi(H);
3144            I_N=local_IN(h,0)[1];
3145            kill(A,B,h,psi);
3146            interdiv[N+NN+counterEPl]=I_N;
3147            setring aff_r;
3148            kill(SS);
3149          }
3150        }
3151      }
3152    }
3153    kill(TAP,auxP);
3154  }
3155  kill(h,CAL,TIP);
3156  @EXTRA@[1]=EP;
3157  @EXTRA@[2]=ES;
3158  @EXTRA@[3]=extra_dgs;
3159  kill(EP);
3160  list update_CURVE=CURVE;
3161  if (size(extra_dgs)>0)
3162  {
3163    export(@EXTRA@);
3164    update_CURVE[1][1]=basering;
3165  }
3166  else
3167  {
3168    kill(@EXTRA@);
3169  }
3170  setring BRing;
3171  kill(h);
3172  kill(aff_r);
3173  list answer=list();
3174  answer[1]=interdiv;
3175  answer[2]=update_CURVE;
3176  option(set,opgt);
3177  return(answer);
3178}
3179///////////////////////////////////////////////////////////////////////////////
3180static proc local_eq (poly H,SS,int m)
3181{
3182  // computes a local equation of poly H in the ring SS related to the place
3183  //     "m"
3184  // list POINT/POINTS is searched depending on wether m=0 or m>0 respectively
3185  // warning : the procedure must be called from ring "Proj_R" and returns a
3186  //     string
3187  def BRing=basering;
3188  setring SS;
3189  if (m>0)
3190  {
3191    if (POINTS[m][3]==number(1))
3192    {
3193      number A=POINTS[m][1];
3194      number B=POINTS[m][2];
3195      map psi=BRing,x+A,y+B,1;
3196    }
3197    else
3198    {
3199      if (POINTS[m][2]==number(1))
3200      {
3201        number A=POINTS[m][1];
3202        map psi=BRing,x+A,1,y;
3203      }
3204      else
3205      {
3206        map psi=BRing,1,x,y;
3207      }
3208    }
3209  }
3210  else
3211  {
3212    number A=POINT[1];
3213    number B=POINT[2];
3214    map psi=BRing,x+A,y+B,1;
3215  }
3216  poly h=psi(H);
3217  string str_h=string(h);
3218  setring BRing;
3219  return(str_h);
3220}
3221///////////////////////////////////////////////////////////////////////////////
3222static proc min_wt_rmat (matrix M)
3223{
3224  // finds the row of M with minimum non-zero entries, i.e. minimum
3225  // "Hamming-weight"
3226  int m=nrows(M);
3227  int n=ncols(M);
3228  int i,j;
3229  int Hwt=0;
3230  for (j=1;j<=n;j=j+1)
3231  {
3232    if (M[1,j]<>0)
3233    {
3234      Hwt=Hwt+1;
3235    }
3236  }
3237  int minHwt=Hwt;
3238  int k=1;
3239  for (i=2;i<=m;i=i+1)
3240  {
3241    Hwt=0;
3242    for (j=1;j<=n;j=j+1)
3243    {
3244      if (M[i,j]<>0)
3245      {
3246        Hwt=Hwt+1;
3247      }
3248    }
3249    if (Hwt<minHwt)
3250    {
3251      minHwt=Hwt;
3252      k=i;
3253    }
3254  }
3255  return(k);
3256}
3257///////////////////////////////////////////////////////////////////////////////
3258
3259// ============================================================================
3260// *******    MAIN PROCEDURE : the Brill-Noether algorithm             ********
3261// ============================================================================
3262
3263proc BrillNoether (intvec G,list CURVE)
3264"USAGE:    BrillNoether(G,CURVE); G an intvec, CURVE a list
3265RETURN:   list of ideals (each of them with two homogeneous generators,
3266          which represent the numerator, resp. denominator, of a rational
3267          function).@*
3268          The corresponding rational functions form a vector basis of the
3269          linear system L(G), G a rational divisor over a non-singular curve.
3270NOTE:     The procedure must be called from the ring CURVE[1][2], where
3271          CURVE is the output of the procedure @code{NSplaces}. @*
3272          The intvec G represents a rational divisor supported on the closed
3273          places of CURVE[3] (e.g. @code{G=2,0,-1;} means 2 times the closed
3274          place 1 minus 1 times the closed place 3).
3275SEE ALSO: Adj_div, NSplaces, Weierstrass
3276EXAMPLE:  example BrillNoether; shows an example
3277"
3278{
3279  // computes a vector basis for the space L(G),
3280  //     where G is a given rational divisor over the non-singular curve
3281  // returns : list of ideals in R each with 2 elements H,Ho such that
3282  //           the set of functions {H/Ho} is the searched basis
3283  // warning : the conductor and sufficiently many points of the plane
3284  //           curve should be computed in advance, in list CURVE
3285  // the algorithm of Brill-Noether is carried out in the procedure
3286  def BRing=basering;
3287  int degX=CURVE[2][1];
3288  list Places=CURVE[3];
3289  intvec Conductor=CURVE[4];
3290  if (deg_D(G,Places)<0)
3291  {
3292    return(list());
3293  }
3294  intvec nuldiv;
3295  if (G==nuldiv)
3296  {
3297    list quickL=list();
3298    ideal quickId;
3299    quickId[1]=1;
3300    quickId[2]=1;
3301    quickL[1]=quickId;
3302    return(quickL);
3303  }
3304  intvec J=max_D(G,nuldiv)+Conductor;
3305  int n=estim_n(J,degX,Places);
3306  dbprint(printlevel+1,"Forms of degree "+string(n)+" : ");
3307  matrix W=nmultiples(n,degX,CHI);
3308  kill(nFORMS(n-degX));
3309  list update_CURVE=CURVE;
3310  matrix V=interpolating_forms(J,n,update_CURVE);
3311  matrix VmW=supplement(W,V);
3312  int k=min_wt_rmat(VmW);
3313  int N=size(nFORMS(n));
3314  matrix H0[1][N];
3315  int i,j;
3316  for (i=1;i<=N;i=i+1)
3317  {
3318    H0[1,i]=VmW[k,i];
3319  }
3320  poly Ho;
3321  for (i=1;i<=N;i=i+1)
3322  {
3323    Ho=Ho+(H0[1,i]*nFORMS(n)[i]);
3324  }
3325  list INTERD=intersection_div(Ho,update_CURVE);
3326  intvec NHo=INTERD[1];
3327  update_CURVE=INTERD[2];
3328  intvec AR=NHo-G;
3329  matrix V2=interpolating_forms(AR,n,update_CURVE);
3330  def aux_RING=update_CURVE[1][1];
3331  setring aux_RING;
3332  if (defined(@EXTRA@)<>0)
3333  {
3334    kill(@EXTRA@);
3335  }
3336  setring BRing;
3337  update_CURVE[1][1]=aux_RING;
3338  kill(aux_RING);
3339  matrix B0=supplement(W,V2);
3340  if (Hamming_wt(B0)==0)
3341  {
3342    return(list());
3343  }
3344  int ld=nrows(B0);
3345  list Bres=list();
3346  ideal HH;
3347  poly H;
3348  for (j=1;j<=ld;j=j+1)
3349  {
3350    H=0;
3351    for (i=1;i<=N;i=i+1)
3352    {
3353      H=H+(B0[j,i]*nFORMS(n)[i]);
3354    }
3355    HH=H,Ho;
3356    Bres[j]=simplifyRF(HH);
3357  }
3358  kill(nFORMS(n));
3359  dbprint(printlevel+1," ");
3360  dbprint(printlevel+2,"Vector basis successfully computed ");
3361  dbprint(printlevel+1," ");
3362  return(Bres);
3363}
3364example
3365{
3366  "EXAMPLE:";  echo = 2;
3367  int plevel=printlevel;
3368  printlevel=-1;
3369  ring s=2,(x,y),lp;
3370  list C=Adj_div(x3y+y3+x);
3371  C=NSplaces(3,C);
3372  // the first 3 Places in C[3] are of degree 1.
3373  // we define the rational divisor G = 4*C[3][1]+4*C[3][3] (of degree 8):
3374  intvec G=4,0,4;
3375  def R=C[1][2];
3376  setring R;
3377  list LG=BrillNoether(G,C);
3378  // here is the vector basis of L(G):
3379  LG;
3380  printlevel=plevel;
3381}
3382///////////////////////////////////////////////////////////////////////////////
3383
3384// *** procedures for dealing with "RATIONAL FUNCTIONS" over a plane curve ***
3385// a rational function F may be given by (homogeneous) ideal or (affine) poly
3386// (or number)
3387
3388static proc polytoRF (F)
3389{
3390  // converts a poly (or number) into a "rational function" of type "ideal"
3391  // warning : it must be called inside "R" and poly should be affine
3392  ideal RF;
3393  RF[1]=homog(F,z);
3394  RF[2]=z^(deg(F));
3395  return(RF);
3396}
3397///////////////////////////////////////////////////////////////////////////////
3398static proc simplifyRF (ideal F)
3399{
3400  // simplifies a rational function f/g extracting the gcd(f,g)
3401  // maybe add a "restriction" to the curve "CHI" but it is not easy to
3402  // programm
3403  poly auxp=gcd(F[1],F[2]);
3404  return(ideal(division(F,auxp)[1]));
3405}
3406///////////////////////////////////////////////////////////////////////////////
3407static proc sumRF (F,G)
3408{
3409  // sum of two "rational functions" F,G given by either a pair
3410  // numerator/denominator or a poly
3411  if ( typeof(F)=="ideal" && typeof(G)=="ideal" )
3412  {
3413    ideal FG;
3414    FG[1]=F[1]*G[2]+F[2]*G[1];
3415    FG[2]=F[2]*G[2];
3416    return(simplifyRF(FG));
3417  }
3418  else
3419  {
3420    if (typeof(F)=="ideal")
3421    {
3422      ideal GG=polytoRF(G);
3423      ideal FG;
3424      FG[1]=F[1]*GG[2]+F[2]*GG[1];
3425      FG[2]=F[2]*GG[2];
3426      return(simplifyRF(FG));
3427    }
3428    else
3429    {
3430      if (typeof(G)=="ideal")
3431      {
3432        ideal FF=polytoRF(F);
3433        ideal FG;
3434        FG[1]=FF[1]*G[2]+FF[2]*G[1];
3435        FG[2]=FF[2]*G[2];
3436        return(simplifyRF(FG));
3437      }
3438      else
3439      {
3440        return(F+G);
3441      }
3442    }
3443  }
3444}
3445///////////////////////////////////////////////////////////////////////////////
3446static proc negRF (F)
3447{
3448  // returns -F as rational function
3449  if (typeof(F)=="ideal")
3450  {
3451    ideal FF=F;
3452    FF[1]=-F[1];
3453    return(FF);
3454  }
3455  else
3456  {
3457    return(-F);
3458  }
3459}
3460///////////////////////////////////////////////////////////////////////////////
3461static proc escprodRF (l,F)
3462{
3463  // computes l*F as rational function
3464  // l should be either a number or a poly of degree zero
3465  if (typeof(F)=="ideal")
3466  {
3467    ideal lF=F;
3468    lF[1]=l*F[1];
3469    return(lF);
3470  }
3471  else
3472  {
3473    return(l*F);
3474  }
3475}
3476///////////////////////////////////////////////////////////////////////////////
3477
3478// ******** procedures to compute Weierstrass semigroups ********
3479
3480static proc orderRF (ideal F,SS,int m)
3481"USAGE:  orderRF(F,SS,m); F an ideal, SS a ring and m an integer
3482RETURN: list with the order (int) and the leading coefficient (number)
3483NOTE:   F represents a rational function, thus the procedure must be
3484        called from global ring R or R(d).
3485        SS contains the name of a local ring where rational places are
3486        stored, and then we take that which is in position m in the
3487        corresponding lists of data.
3488        The order of F at the place given by SS,m is returned together
3489        with the coefficient of minimum degree in the corresponding power
3490        series.
3491"
3492{
3493  // computes the order of a rational function F at a RATIONAL place given by
3494  //     a local ring SS and a position "m" inside SS
3495  // warning : the procedure must be called from global projective ring "R" or
3496  //     "R(i)"
3497  // returns a list with the order (int) and the "leading coefficient" (number)
3498  def BR=basering;
3499  poly f=F[1];
3500  string sf=local_eq(f,SS,m);
3501  poly g=F[2];
3502  string sg=local_eq(g,SS,m);
3503  setring SS;
3504  execute("poly ff="+sf+";");
3505  execute("poly gg="+sg+";");
3506  list o1=local_IN(ff,m);
3507  list o2=local_IN(gg,m);
3508  int oo=o1[1]-o2[1];
3509  number lc=o1[2]/o2[2];
3510  setring BR;
3511  number LC=number(imap(SS,lc));
3512  return(list(oo,LC));
3513}
3514///////////////////////////////////////////////////////////////////////////////
3515
3516proc Weierstrass (int P,int m,list CURVE)
3517"USAGE:    Weierstrass( i, m, CURVE );  i,m integers and CURVE a list
3518RETURN:   list WS of two lists:
3519  @format
3520  WS[1] list of integers (Weierstr. semigroup of the curve at place i up to m)
3521  WS[2] list of ideals (the associated rational functions)
3522  @end format
3523NOTE:     The procedure must be called from the ring CURVE[1][2],
3524          where CURVE is the output of the procedure @code{NSplaces}.
3525@*        i represents the place CURVE[3][i].
3526@*        Rational functions are represented by numerator/denominator
3527          in form of ideals with two homogeneous generators.
3528WARNING:  The place must be rational, i.e., necessarily CURVE[3][i][1]=1. @*
3529SEE ALSO: Adj_div, NSplaces, BrillNoether
3530EXAMPLE:  example Weierstrass; shows an example
3531"
3532{
3533  // computes the Weierstrass semigroup at a RATIONAL place P up to a bound "m"
3534  //   together with the functions achieving each value up to m, via
3535  //   Brill-Noether
3536  // returns 2 lists : the first consists of all the poles up to m in
3537  //   increasing order and the second consists of the corresponging rational
3538  //   functions
3539  list Places=CURVE[3];
3540  intvec pl=Places[P];
3541  int dP=pl[1];
3542  int nP=pl[2];
3543  if (dP<>1)
3544  {
3545    ERROR("the given place is not defined over the prime field");
3546  }
3547  if (m<=0)
3548  {
3549    if (m==0)
3550    {
3551      list semig=list();
3552      int auxint=0;
3553      semig[1]=auxint;
3554      list funcs=list();
3555      ideal auxF;
3556      auxF[1]=1;
3557      auxF[2]=1;
3558      funcs[1]=auxF;
3559      return(list(semig,funcs));
3560    }
3561    else
3562    {
3563      ERROR("second argument must be non-negative");
3564    }
3565  }
3566  int auxint=0;
3567  ideal auxF;
3568  auxF[1]=1;
3569  auxF[2]=1;
3570  // Brill-Noether algorithm
3571  intvec mP;
3572  mP[P]=m;
3573  list LmP=BrillNoether(mP,CURVE);
3574  int lmP=size(LmP);
3575  if (lmP==1)
3576  {
3577    list semig=list();
3578    semig[1]=auxint;
3579    list funcs=list();
3580    funcs[1]=auxF;
3581    return(list(semig,funcs));
3582  }
3583  def SS=CURVE[5][dP][1];
3584  list ordLmP=list();
3585  list polLmP=list();
3586  list sortpol=list();
3587  int maxpol;
3588  int i,j,k;
3589  for (i=1;i<=lmP-1;i=i+1)
3590  {
3591    for (j=1;j<=lmP-i+1;j=j+1)
3592    {
3593      ordLmP[j]=orderRF(LmP[j],SS,nP);
3594      polLmP[j]=-ordLmP[j][1];
3595    }
3596    sortpol=sort(polLmP);
3597    polLmP=sortpol[1];
3598    maxpol=polLmP[lmP-i+1];
3599    LmP=permute_L(LmP,sortpol[2]);
3600    ordLmP=permute_L(ordLmP,sortpol[2]);
3601    // triangulate LmP
3602    for (k=1;k<=lmP-i;k=k+1)
3603    {
3604      if (polLmP[lmP-i+1-k]==maxpol)
3605      {
3606        LmP[lmP-i+1-k]=sumRF(LmP[lmP-i+1-k],negRF(escprodRF(
3607                       ordLmP[lmP-i+1-k][2]/ordLmP[lmP-i+1][2],LmP[lmP-i+1])));
3608      }
3609      else
3610      {
3611        break;
3612      }
3613    }
3614  }
3615  polLmP[1]=auxint;
3616  LmP[1]=auxF;
3617  return(list(polLmP,LmP));
3618}
3619example
3620{
3621  "EXAMPLE:";  echo = 2;
3622  int plevel=printlevel;
3623  printlevel=-1;
3624  ring s=2,(x,y),lp;
3625  list C=Adj_div(x3y+y3+x);
3626  C=NSplaces(3,C);
3627  def R=C[1][2];
3628  setring R;
3629  // Place C[3][1] has degree 1 (i.e it is rational);
3630  list WS=Weierstrass(1,7,C);
3631  // the first part of the list is the Weierstrass semigroup up to 7 :
3632  WS[1];
3633  // and the second part are the corresponding functions :
3634  WS[2];
3635  printlevel=plevel;
3636}
3637///////////////////////////////////////////////////////////////////////////////
3638
3639// axiliary procedure for permuting a list or intvec
3640
3641proc permute_L (L,P)
3642"USAGE:    permute_L( L, P ); L,P either intvecs or lists
3643RETURN:   list obtained from L by applying the permutation given by P.
3644NOTE:     If P is a list, all entries must be integers.
3645SEE ALSO: sys_code, AGcode_Omega, prepSV
3646EXAMPLE:  example permute_L; shows an example
3647"
3648{
3649  // permutes the list L according to the permutation P (both intvecs or
3650  //  lists of integers)
3651  int s=size(L);
3652  int n=size(P);
3653  int i;
3654  if (s<n)
3655  {
3656    for (i=s+1;i<=n;i=i+1)
3657    {
3658      L[i]=0;
3659    }
3660    s=size(L);
3661  }
3662  list auxL=L;
3663  for (i=1;i<=n;i=i+1)
3664  {
3665    auxL[i]=L[P[i]];
3666  }
3667  return(auxL);
3668}
3669example
3670{
3671  "EXAMPLE:";  echo = 2;
3672  list L=list();
3673  L[1]="a";
3674  L[2]="b";
3675  L[3]="c";
3676  L[4]="d";
3677  intvec P=1,3,4,2;
3678  // the list L is permuted according to P :
3679  permute_L(L,P);
3680}
3681///////////////////////////////////////////////////////////////////////////////
3682static proc evalRF (ideal F,SS,int m)
3683"USAGE:     evalRF(F,SS,m), where F is an ideal, SS is a ring and m is an
3684            integer
3685RETURN:     the evaluation (number) of F at the place given by SS,m if it is
3686            well-defined
3687NOTE:       F represents a rational function, thus the procedure must be
3688            called from R or R(d).
3689            SS contains the name of a local ring where rational places are
3690            stored, and then we take that which is in position m in the
3691            corresponding lists of data.
3692"
3693{
3694  // evaluates a rational function F at a RATIONAL place given by
3695  //     a local ring SS and a position "m" inside SS
3696  list olc=orderRF(F,SS,m);
3697  int oo=olc[1];
3698  if (oo==0)
3699  {
3700    return(olc[2]);
3701  }
3702  else
3703  {
3704    if (oo>0)
3705    {
3706      return(number(0));
3707    }
3708    else
3709    {
3710      ERROR("the function is not well-defined at the given place");
3711    }
3712  }
3713}
3714///////////////////////////////////////////////////////////////////////////////
3715//
3716// ******** procedures for constructing AG codes ********
3717//
3718///////////////////////////////////////////////////////////////////////////////
3719
3720static proc gen_mat (list LF,intvec LP,RP)
3721"USAGE:      gen_mat(LF,LP,RP); LF list of rational functions,
3722            LP intvec of rational places and RP a local ring
3723RETURN:     a generator matrix of the evaluation code given by LF and LP
3724SEE ALSO:   extcurve
3725KEYWORDS:   evaluation codes
3726NOTE:       Rational places are searched in local ring RP
3727            The procedure must be called from R or R(d) fromlist CURVE
3728            after having executed extcurve(d,CURVE)
3729"
3730{
3731  // computes a generator matrix (with numbers) of the evaluation code given
3732  //    by a list of rational functions LF and a list of RATIONAL places LP
3733  int m=size(LF);
3734  int n=size(LP);
3735  matrix GM[m][n];
3736  int i,j;
3737  for (i=1;i<=m;i=i+1)
3738  {
3739    for (j=1;j<=n;j=j+1)
3740    {
3741      GM[i,j]=evalRF(LF[i],RP,LP[j]);
3742    }
3743  }
3744  return(GM);
3745}
3746///////////////////////////////////////////////////////////////////////////////
3747
3748proc dual_code (matrix G)
3749"USAGE:   dual_code(G); G a matrix of numbers
3750RETURN:   a generator matrix of the dual code generated by G
3751NOTE:     The input should be a matrix G of numbers. @*
3752          The output is also a parity check matrix for the code defined by G
3753KEYWORDS: linear code, dual
3754EXAMPLE:  example dual_code; shows an example
3755"
3756{
3757  // computes the dual code of C given by a generator matrix G
3758  //     i.e. computes a parity-check matrix H of C
3759  // conversely : computes also G if H is given
3760  return(Ker(G));
3761}
3762example
3763{
3764  "EXAMPLE:";  echo = 2;
3765  ring s=2,T,lp;
3766  // here is the Hamming code of length 7 and dimension 3
3767  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;
3768  print(G);
3769  matrix H=dual_code(G);
3770  print(H);
3771}
3772///////////////////////////////////////////////////////////////////////////////
3773
3774// ======================================================================
3775// *********** initial test for disjointness ***************
3776// ======================================================================
3777
3778static proc disj_divs (intvec H,intvec P,list EC)
3779{
3780  int s1=size(H);
3781  int s2=size(P);
3782  list PLACES=EC[3];
3783  def BRing=basering;
3784  def auxR=EC[1][5];
3785  setring auxR;
3786  int s=res_deg();
3787  setring BRing;
3788  kill(auxR);
3789  int i,j,k,d,l,N,M;
3790  intvec auxIV,iw;
3791  for (i=1;i<=s;i=i+1)
3792  {
3793    if ( (s mod i) == 0 )
3794    {
3795      def auxR=EC[5][i][1];
3796      setring auxR;
3797      auxIV[i]=size(POINTS);
3798      setring BRing;
3799      kill(auxR);
3800    }
3801    else
3802    {
3803      auxIV[i]=0;
3804    }
3805  }
3806  for (i=1;i<=s1;i=i+1)
3807  {
3808    if (H[i]<>0)
3809    {
3810      iw=PLACES[i];
3811      d=iw[1];
3812      if ( (s mod d) == 0 )
3813      {
3814        l=iw[2];
3815        // check that this place(s) are not in sup(D)
3816        if (d==1)
3817        {
3818          for (j=1;j<=s2;j=j+1)
3819          {
3820            if (l==P[j])
3821            {
3822              return(0);
3823            }
3824          }
3825        }
3826        else
3827        {
3828          N=0;
3829          for (j=1;j<d;j=j+1)
3830          {
3831            N=N+j*auxIV[j];
3832          }
3833          N=N+d*(l-1);
3834          M=N+d;
3835          for (k=N+1;k<=M;k=k+1)
3836          {
3837            for (j=1;j<=s2;j=j+1)
3838            {
3839              if (k==P[j])
3840              {
3841                return(0);
3842              }
3843            }
3844          }
3845        }
3846      }
3847    }
3848  }
3849  kill(auxIV,iw);
3850  return(1);
3851}
3852///////////////////////////////////////////////////////////////////////////////
3853
3854proc AGcode_L (intvec G,intvec D,list EC)
3855"USAGE:    AGcode_L( G, D, EC );  G,D intvec, EC a list
3856RETURN:   a generator matrix for the evaluation AG code defined by the
3857          divisors G and D.
3858NOTE:     The procedure must be called within the ring EC[1][4],
3859          where EC is the output of @code{extcurve(d)} (or within
3860          the ring EC[1][2] if d=1). @*
3861          The entry i in the intvec D refers to the i-th rational
3862          place in EC[1][5] (i.e., to POINTS[i], etc., see @ref{extcurve}).@*
3863          The intvec G represents a rational divisor (see @ref{BrillNoether}
3864          for more details).@*
3865          The code evaluates the vector basis of L(G) at the rational
3866          places given by D.
3867WARNINGS: G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
3868          not checked by the algorithm.
3869          G and D should have disjoint supports (checked by the algorithm).
3870SEE ALSO: Adj_div, BrillNoether, extcurve, AGcode_Omega
3871EXAMPLE:  example AGcode_L; shows an example
3872"
3873{
3874  // returns a generator matrix for the evaluation AG code given by G and D
3875  // G must be a divisor defined over the prime field and D an intvec of
3876  // "rational places"
3877  // it must be called inside R or R(d) and requires previously "extcurve(d)"
3878  def BR=basering;
3879  if (disj_divs(G,D,EC)==0)
3880  {
3881    dbprint(printlevel+3,"? the divisors do not have disjoint supports,
3882                                 0-matrix returned ?");
3883    matrix answer;
3884    return(answer);
3885  }
3886  if (res_deg()>1)
3887  {
3888    def R=EC[1][2];
3889    setring R;
3890    list LG=BrillNoether(G,EC);
3891    setring BR;
3892    list LG=imap(R,LG);
3893    setring R;
3894    kill(LG);
3895    setring BR;
3896    kill(R);
3897  }
3898  else
3899  {
3900    list LG=BrillNoether(G,EC);
3901  }
3902  def RP=EC[1][5];
3903  matrix M=gen_mat(LG,D,RP);
3904  kill(LG);
3905  return(M);
3906}
3907example
3908{
3909  "EXAMPLE:";  echo = 2;
3910  int plevel=printlevel;
3911  printlevel=-1;
3912  ring s=2,(x,y),lp;
3913  list HC=Adj_div(x3+y2+y);
3914  HC=NSplaces(1,HC);
3915  HC=extcurve(2,HC);
3916  def ER=HC[1][4];
3917  setring ER;
3918  intvec G=5;      // the rational divisor G = 5*HC[3][1]
3919  intvec D=2..9;   // D = sum of the rational places no. 2..9 over F_4
3920  // let us construct the corresponding evaluation AG code :
3921  matrix C=AGcode_L(G,D,HC);
3922  // here is a linear code of type [8,5,>=3] over F_4
3923  print(C);
3924  printlevel=plevel;
3925}
3926///////////////////////////////////////////////////////////////////////////////
3927
3928proc AGcode_Omega (intvec G,intvec D,list EC)
3929"USAGE:    AGcode_Omega( G, D, EC ); G,D intvec, EC a list
3930RETURN:   a generator matrix for the residual AG code defined by the
3931          divisors G and D.
3932NOTE:     The procedure must be called within the ring EC[1][4],
3933          where EC is the output of @code{extcurve(d)} (or within
3934          the ring EC[1][2] if d=1). @*
3935          The entry i in the intvec D refers to the i-th rational
3936          place in EC[1][5] (i.e., to POINTS[i], etc., see @ref{extcurve}).@*
3937          The intvec G represents a rational divisor (see @ref{BrillNoether}
3938          for more details).@*
3939          The code computes the residues of a vector space basis of
3940          @math{\Omega(G-D)} at the rational places given by D.
3941WARNINGS: G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
3942          not checked by the algorithm.
3943          G and D should have disjoint supports (checked by the algorithm).
3944SEE ALSO: Adj_div, BrillNoether, extcurve, AGcode_L
3945EXAMPLE:  example AGcode_Omega; shows an example
3946"
3947{
3948  // returns a generator matrix for the residual AG code given by G and D
3949  // G must be a divisor defined over the prime field and D an intvec or
3950  // "rational places"
3951  // it must be called inside R or R(d) and requires previously "extcurve(d)"
3952  return(dual_code(AGcode_L(G,D,EC)));
3953}
3954example
3955{
3956  "EXAMPLE:";  echo = 2;
3957  int plevel=printlevel;
3958  printlevel=-1;
3959  ring s=2,(x,y),lp;
3960  list HC=Adj_div(x3+y2+y);
3961  HC=NSplaces(1,HC);
3962  HC=extcurve(2,HC);
3963  def ER=HC[1][4];
3964  setring ER;
3965  intvec G=5;      // the rational divisor G = 5*HC[3][1]
3966  intvec D=2..9;   // D = sum of the rational places no. 2..9 over F_4
3967  // let us construct the corresponding residual AG code :
3968  matrix C=AGcode_Omega(G,D,HC);
3969  // here is a linear code of type [8,3,>=5] over F_4
3970  print(C);
3971  printlevel=plevel;
3972}
3973///////////////////////////////////////////////////////////////////////////////
3974
3975// ============================================================================
3976// *******   auxiliary procedure to define AG codes over extensions    ********
3977// ============================================================================
3978
3979proc extcurve (int d,list CURVE)
3980"USAGE:    extcurve( d, CURVE ); d an integer, CURVE a list
3981RETURN:   list L which is the update of the list CURVE with additional entries
3982   @format
3983   L[1][3]: ring (p,a),(x,y),lp (affine),
3984   L[1][4]: ring (p,a),(x,y,z),lp (projective),
3985   L[1][5]: ring (p,a),(x,y,t),ls (local),
3986   L[2][3]: int  (the number of rational places),
3987   @end format
3988          the rings being defined over a field extension of degree d. @*
3989          If d<2 then @code{extcurve(d,CURVE);} creates a list L which
3990          is the update of the list CURVE with additional entries
3991   @format
3992   L[1][5]: ring p,(x,y,t),ls,
3993   L[2][3]: int  (the number of places over the base field).
3994   @end format
3995          In both cases, in the ring L[1][5] lists with the data for all the
3996          rational places (after a field extension of degree d) are
3997          created (see @ref{Adj_div}):
3998   @format
3999   lists POINTS, LOC_EQS, BRANCHES, PARAMETRIZATIONS.
4000   @end format
4001NOTE:     The list CURVE should be the output of @code{NSplaces} and has
4002          to contain (at least) all places up to degree d. @*
4003          This procedure must be executed before constructing AG codes,
4004          even if no extension is needed. The ring L[1][4] must be active
4005          when constructing codes over the field extension.@*
4006SEE ALSO: closed_points, Adj_div, NSplaces, AGcode_L, AGcode_Omega
4007EXAMPLE:  example extcurve; shows an example
4008"
4009{
4010  // extends the underlying curve and all its associated objects to a larger
4011  //     base field in order to evaluate points over such a extension
4012  // if d<2 then the only change is that a local ring "RatPl" (which is a
4013  //     copy of "S(1)") is created in order to store the rational places
4014  //     where we can do evaluations
4015  // otherwise, such a ring contains all places which are rational over the
4016  //     extension
4017  // warning : list Places does not change so that only divisors which are
4018  //     "rational over the prime field" are allowed; this probably will
4019  //     change in the future
4020  // warning : the places in RatPl are ranged in increasing degree, respecting
4021  //     the order from list Places and placing the conjugate branches all
4022  //     together
4023  def BR=basering;
4024  list ext_CURVE=CURVE;
4025  if (d<2)
4026  {
4027    def SS=CURVE[5][1][1];
4028    ring RatPl=char(basering),(x,y,t),ls;
4029    list POINTS=imap(SS,POINTS);
4030    list LOC_EQS=imap(SS,LOC_EQS);
4031    list BRANCHES=imap(SS,BRANCHES);
4032    list PARAMETRIZATIONS=imap(SS,PARAMETRIZATIONS);
4033    export(POINTS);
4034    export(LOC_EQS);
4035    export(BRANCHES);
4036    export(PARAMETRIZATIONS);
4037    int NrRatPl=size(POINTS);
4038    ext_CURVE[2][3]=NrRatPl;
4039    setring BR;
4040    ext_CURVE[1][5]=RatPl;
4041    dbprint(printlevel+1,"");
4042    dbprint(printlevel+2,"Total number of rational places : "+string(NrRatPl));
4043    dbprint(printlevel+1,"");
4044    kill(RatPl);
4045    return(ext_CURVE);
4046  }
4047  else
4048  {
4049    if (size(CURVE[5])>=d)
4050    {
4051      int i,j,k;
4052      if (typeof(CURVE[5][d])=="list")
4053      {
4054        def S(d)=CURVE[5][d][1];
4055        setring S(d);
4056      }
4057      else
4058      {
4059        ring S(d)=char(basering),t,ls;
4060        list POINTS=list();
4061      }
4062      string smin=string(minpoly);
4063      setring BR;
4064      ring RatPl=(char(basering),a),(x,y,t),ls;
4065      execute("minpoly="+smin+";");
4066      list POINTS=imap(S(d),POINTS);
4067      list LOC_EQS=imap(S(d),LOC_EQS);
4068      list BRANCHES=imap(S(d),BRANCHES);
4069      list PARAMETRIZATIONS=imap(S(d),PARAMETRIZATIONS);
4070      kill(S(d));
4071      int s=size(POINTS);
4072      int counter=0;
4073      int piv=0;
4074      for (j=1;j<=s;j=j+1)
4075      {
4076        counter=counter+1;
4077        piv=counter;
4078        for (k=1;k<d;k=k+1)
4079        {
4080          POINTS=insert(POINTS,Frobenius(POINTS[piv],1),counter);
4081          LOC_EQS=insert(LOC_EQS,Frobenius(LOC_EQS[piv],1),counter);
4082          BRANCHES=insert(BRANCHES,conj_b(BRANCHES[piv],1),counter);
4083          PARAMETRIZATIONS=insert(PARAMETRIZATIONS,Frobenius(
4084                                             PARAMETRIZATIONS[piv],1),counter);
4085          counter=counter+1;
4086          piv=counter;
4087        }
4088      }
4089      string olda;
4090      poly paux;
4091      intvec iv,iw;
4092      int ii,jj,kk,m,n;
4093      for (i=d-1;i>=2;i=i-1)
4094      {
4095        if ( (d mod i)==0 )
4096        {
4097          if (typeof(CURVE[5][i])=="list")
4098          {
4099            def S(i)=CURVE[5][i][1];
4100          }
4101          else
4102          {
4103            ring S(i)=char(basering),t,ls;
4104            list POINTS=list();
4105            setring BR;
4106           }
4107          olda=subfield(S(i));
4108          setring S(i);
4109          s=size(POINTS);
4110          setring RatPl;
4111          for (j=s;j>=1;j=j-1)
4112          {
4113            counter=0;
4114            POINTS=insert(POINTS,list(),0);
4115            POINTS[1][1]=number(importdatum(S(i),"POINTS["+string(j)
4116                                                            +"][1]",olda));
4117            POINTS[1][2]=number(importdatum(S(i),"POINTS["+string(j)
4118                                                            +"][2]",olda));
4119            POINTS[1][3]=number(importdatum(S(i),"POINTS["+string(j)
4120                                                            +"][3]",olda));
4121            LOC_EQS=insert(LOC_EQS,importdatum(S(i),"LOC_EQS["+string(j)
4122                                                            +"]",olda),0);
4123            BRANCHES=insert(BRANCHES,list(),0);
4124            setring S(i);
4125            m=nrows(BRANCHES[j][1]);
4126            n=ncols(BRANCHES[j][1]);
4127            iv=BRANCHES[j][2];
4128            kk=BRANCHES[j][3];
4129            poly par@1=subst(PARAMETRIZATIONS[j][1][1],t,x);
4130            poly par@2=subst(PARAMETRIZATIONS[j][1][2],t,x);
4131            export(par@1);
4132            export(par@2);
4133            iw=PARAMETRIZATIONS[j][2];
4134            setring RatPl;
4135            paux=importdatum(S(i),"BRANCHES["+string(j)+"][4]",olda);
4136            matrix Maux[m][n];
4137            for (ii=1;ii<=m;ii=ii+1)
4138            {
4139              for (jj=1;jj<=n;jj=jj+1)
4140              {
4141                Maux[ii,jj]=importdatum(S(i),"BRANCHES["+string(j)
4142                                 +"][1]["+string(ii)+","+string(jj)+"]",olda);
4143              }
4144            }
4145            BRANCHES[1][1]=Maux;
4146            BRANCHES[1][2]=iv;
4147            BRANCHES[1][3]=kk;
4148            BRANCHES[1][4]=paux;
4149            kill(Maux);
4150            PARAMETRIZATIONS=insert(PARAMETRIZATIONS,list(),0);
4151            PARAMETRIZATIONS[1][1]=ideal(0);
4152            PARAMETRIZATIONS[1][1][1]=importdatum(S(i),"par@1",olda);
4153            PARAMETRIZATIONS[1][1][2]=importdatum(S(i),"par@2",olda);
4154            PARAMETRIZATIONS[1][1][1]=subst(PARAMETRIZATIONS[1][1][1],x,t);
4155            PARAMETRIZATIONS[1][1][2]=subst(PARAMETRIZATIONS[1][1][2],x,t);
4156            PARAMETRIZATIONS[1][2]=iw;
4157            for (k=1;k<i;k=k+1)
4158            {
4159              counter=counter+1;
4160              piv=counter;
4161              POINTS=insert(POINTS,Frobenius(POINTS[piv],1),counter);
4162              LOC_EQS=insert(LOC_EQS,Frobenius(LOC_EQS[piv],1),counter);
4163              BRANCHES=insert(BRANCHES,conj_b(BRANCHES[piv],1),counter);
4164              PARAMETRIZATIONS=insert(PARAMETRIZATIONS,Frobenius(
4165                                             PARAMETRIZATIONS[piv],1),counter);
4166            }
4167            setring S(i);
4168            kill(par@1,par@2);
4169            setring RatPl;
4170          }
4171          kill(S(i));
4172        }
4173      }
4174      kill(iw);
4175      kill(paux);
4176      if (typeof(CURVE[5][1])=="list")
4177      {
4178        def S(1)=CURVE[5][1][1];
4179      }
4180      else
4181      {
4182        ring S(1)=char(basering),t,ls;
4183        list POINTS=list();
4184        setring RatPl;
4185      }
4186      POINTS=imap(S(1),POINTS)+POINTS;
4187      LOC_EQS=imap(S(1),LOC_EQS)+LOC_EQS;
4188      BRANCHES=imap(S(1),BRANCHES)+BRANCHES;
4189      PARAMETRIZATIONS=imap(S(1),PARAMETRIZATIONS)+PARAMETRIZATIONS;
4190      export(POINTS);
4191      export(LOC_EQS);
4192      export(BRANCHES);
4193      export(PARAMETRIZATIONS);
4194      int NrRatPl=size(POINTS);
4195      ext_CURVE[2][3]=NrRatPl;
4196      setring BR;
4197      ext_CURVE[1][5]=RatPl;
4198      ring r(d)=(char(basering),a),(x,y),lp;
4199      execute("minpoly="+smin+";");
4200      setring BR;
4201      ext_CURVE[1][3]=r(d);
4202      ring R(d)=(char(basering),a),(x,y,z),lp;
4203      execute("minpoly="+smin+";");
4204      setring BR;
4205      ext_CURVE[1][4]=R(d);
4206      dbprint(printlevel+1,"");
4207      dbprint(printlevel+2,"Total number of rational places : NrRatPl = "
4208                                   +string(NrRatPl));
4209      dbprint(printlevel+1,"");
4210      kill(S(1));
4211      kill(R(d));
4212      kill(RatPl);
4213      return(ext_CURVE);
4214    }
4215    else
4216    {
4217      ERROR("you must compute first all the places up to degree "+string(d));
4218      return();
4219    }
4220  }
4221}
4222example
4223{
4224  "EXAMPLE:";  echo = 2;
4225  int plevel=printlevel;
4226  printlevel=-1;
4227  ring s=2,(x,y),lp;
4228  list C=Adj_div(x5+y2+y);
4229  C=NSplaces(3,C);
4230  // since we have all points up to degree 4, we can extend the curve
4231  // to that extension, in order to get rational points over F_16;
4232  C=extcurve(4,C);
4233  // e.g., display the basepoint of place no. 32:
4234  def R=C[1][5];
4235  setring R;
4236  POINTS[32];
4237  printlevel=plevel;
4238}
4239///////////////////////////////////////////////////////////////////////////////
4240
4241// specific procedures for linear/AG codes
4242
4243static proc Hamming_wt (matrix A)
4244"USAGE:     Hamming_wt(A), where A is any matrix
4245RETURN:     the Hamming weight (number of non-zero entries) of the matrix A
4246"
4247{
4248  // computes the Hamming weight (number of non-zero entries) of any matrix
4249  // notice that "words" are represented by matrices of size 1xn
4250  // computing the Hamming distance between two matrices can be done by
4251  // Hamming_wt(A-B)
4252  int m=nrows(A);
4253  int n=ncols(A);
4254  int i,j;
4255  int w=0;
4256  for (i=1;i<=m;i=i+1)
4257  {
4258    for (j=1;j<=n;j=j+1)
4259    {
4260      if (A[i,j]<>0)
4261      {
4262        w=w+1;
4263      }
4264    }
4265  }
4266  return(w);
4267}
4268///////////////////////////////////////////////////////////////////////////////
4269
4270// Basic Algorithm of Skorobogatov and Vladut for decoding AG codes
4271// warning : the user must choose carefully the parameters for the code and
4272//     the decoding since they will never be checked by the procedures
4273
4274proc prepSV (intvec G,intvec D,intvec F,list EC)
4275"USAGE:    prepSV( G, D, F, EC ); G,D,F intvecs and EC a list
4276RETURN:   list E of size n+3, where n=size(D). All its entries but E[n+3]
4277          are matrices:
4278   @format
4279   E[1]:  parity check matrix for the current AG code
4280   E[2] ... E[n+2]:  matrices used in the procedure decodeSV
4281   E[n+3]:  intvec with
4282       E[n+3][1]: correction capacity @math{epsilon} of the algorithm
4283       E[n+3][2]: designed Goppa distance @math{delta} of the current AG code
4284   @end format
4285NOTE:     Computes the preprocessing for the basic (Skorobogatov-Vladut)
4286          decoding algorithm.@*
4287          The procedure must be called within the ring EC[1][4], where EC is
4288          the output of @code{extcurve(d)} (or in the ring EC[1][2] if d=1) @*
4289          The intvec G and F represent rational divisors (see
4290          @ref{BrillNoether} for more details).@*
4291          The intvec D refers to rational places (see @ref{AGcode_Omega}
4292          for more details.).
4293          The current AG code is @code{AGcode_Omega(G,D,EC)}.@*
4294          If you know the exact minimum distance d and you want to use it in
4295          @code{decodeSV} instead of @math{delta}, you can change the value
4296          of E[n+3][2] to d before applying decodeSV.
4297          If you have a systematic encoding for the current code and want to
4298          keep it during the decoding, you must previously permute D (using
4299          @code{permute_L(D,P);}), e.g., according to the permutation
4300          P=L[3], L being the output of @code{sys_code}.
4301WARNINGS: F must be a divisor with support disjoint from the support of D and
4302          with degree @math{epsilon + genus}, where
4303          @math{epsilon:=[(deg(G)-3*genus+1)/2]}.@*
4304          G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is
4305          not checked by the algorithm.
4306          G and D should also have disjoint supports (checked by the
4307          algorithm).
4308KEYWORDS: SV-decoding algorithm, preprocessing
4309SEE ALSO: extcurve, AGcode_Omega, decodeSV, sys_code, permute_L
4310EXAMPLE:  example prepSV; shows an example
4311"
4312{
4313  if (disj_divs(F,D,EC)==0)
4314  {
4315    dbprint(printlevel+3,"? the divisors do not have disjoint supports,
4316                            empty list returned ?");
4317    return(list());
4318  }
4319  list E=list();
4320  list Places=EC[3];
4321  int m=deg_D(G,Places);
4322  int genusX=EC[2][2];
4323  int e=(m+1-3*genusX)/2;
4324  if (e<1)
4325  {
4326    dbprint(printlevel+3,"? the correction capacity of the basic algorithm
4327                            is zero, empty list returned ?");
4328    return(list());
4329  }
4330  // deg(F)==e+genusX should be satisfied, and sup(D),sup(F) should be
4331  // disjoint !!!!
4332  int n=size(D);
4333  // 2*genusX-2<m<n should also be satisfied !!!!
4334  matrix EE=AGcode_L(G,D,EC);
4335  int l=nrows(EE);
4336  E[1]=EE;
4337  matrix GP=AGcode_L(F,D,EC);
4338  int r=nrows(GP);
4339  intvec H=G-F;
4340  matrix HP=AGcode_L(H,D,EC);
4341  int s=nrows(HP);
4342  int i,j,k;
4343  kill(EE);
4344  for (k=1;k<=n;k=k+1)
4345  {
4346    E[k+1]=GP[1,k]*submat(HP,1..s,k..k);
4347    for (i=2;i<=r;i=i+1)
4348    {
4349      E[k+1]=concat(E[k+1],GP[i,k]*submat(HP,1..s,k..k));
4350    }
4351  }
4352  E[n+2]=GP;
4353  intvec IW=e,m+2-2*genusX;
4354  E[n+3]=IW;
4355  kill(IW);
4356  return(E);
4357}
4358example
4359{
4360  "EXAMPLE:";  echo = 2;
4361  int plevel=printlevel;
4362  printlevel=-1;
4363  ring s=2,(x,y),lp;
4364  list HC=Adj_div(x3+y2+y);
4365  HC=NSplaces(1,HC);
4366  HC=extcurve(2,HC);
4367  def ER=HC[1][4];
4368  setring ER;
4369  intvec G=5;      // the rational divisor G = 5*HC[3][1]
4370  intvec D=2..9;   // D = sum of the rational places no. 2..9 over F_4
4371  // construct the corresp. residual AG code of type [8,3,>=5] over F_4:
4372  matrix C=AGcode_Omega(G,D,HC);
4373  // we can correct 1 error and the genus is 1, thus F must have degree 2
4374  // and support disjoint from that of D;
4375  intvec F=2;
4376  list SV=prepSV(G,D,F,HC);
4377  // now everything is prepared to decode with the basic algorithm;
4378  // for example, here is a parity check matrix to compute the syndrome :
4379  print(SV[1]);
4380  // and here you have the correction capacity of the algorithm :
4381  int epsilon=SV[size(D)+3][1];
4382  epsilon;
4383  printlevel=plevel;
4384}
4385///////////////////////////////////////////////////////////////////////////////
4386
4387proc decodeSV (matrix y,list K)
4388"USAGE:     decodeSV( y, K ); y a row-matrix and K a list
4389RETURN:     a codeword (row-matrix) if possible, resp. the 0-matrix (of size
4390            1) if decoding is impossible.
4391            For decoding the basic (Skorobogatov-Vladut) decoding algorithm
4392            is applied.
4393NOTE:       The list_expression should be the output K of the procedure
4394            @code{prepSV}.@*
4395            The matrix_expression should be a (1 x n)-matrix,  where
4396            n = ncols(K[1]).@*
4397            The decoding may fail if the number of errors is greater than
4398            the correction capacity of the algorithm.
4399KEYWORDS:   SV-decoding algorithm
4400SEE ALSO:   extcurve, AGcode_Omega, prepSV
4401EXAMPLE:    example decodeSV; shows an example
4402"
4403{
4404  // decodes y with the "basic decoding algorithm", if possible
4405  // requires the preprocessing done by the procedure "prepSV"
4406  // the procedure must be called from ring R or R(d)
4407  // returns either a codeword (matrix) of none (in case of too many errors)
4408  matrix syndr=K[1]*transpose(y);
4409  if (Hamming_wt(syndr)==0)
4410  {
4411    return(y);
4412  }
4413  matrix Ey=y[1,1]*K[2];
4414  int n=ncols(y);
4415  int i;
4416  for (i=2;i<=n;i=i+1)
4417  {
4418    Ey=Ey+y[1,i]*K[i+1];
4419  }
4420  matrix Ky=get_NZsol(Ey);
4421  if (Hamming_wt(Ky)==0)
4422  {
4423    dbprint(printlevel+3,"? no error-locator found ?");
4424    dbprint(printlevel+3,"? too many errors occur, 0-matrix returned ?");
4425    matrix answer;
4426    return(answer);
4427  }
4428  int r=nrows(K[n+2]);
4429  matrix ErrLoc[1][n];
4430  list Z=list();
4431  list NZ=list();
4432  int j;
4433  for (j=1;j<=n;j=j+1)
4434  {
4435    for (i=1;i<=r;i=i+1)
4436    {
4437      ErrLoc[1,j]=ErrLoc[1,j]+K[n+2][i,j]*Ky[1,i];
4438    }
4439    if (ErrLoc[1,j]==0)
4440    {
4441      Z=insert(Z,j,size(Z));
4442    }
4443    else
4444    {
4445      NZ=insert(NZ,j,size(NZ));
4446    }
4447  }
4448  int k=size(NZ);
4449  int l=nrows(K[1]);
4450  int s=l+k;
4451  matrix A[s][n];
4452  matrix b[s][1];
4453  for (i=1;i<=l;i=i+1)
4454  {
4455    for (j=1;j<=n;j=j+1)
4456    {
4457      A[i,j]=K[1][i,j];
4458    }
4459    b[i,1]=syndr[i,1];
4460  }
4461  for (i=1;i<=k;i=i+1)
4462  {
4463    A[l+i,NZ[i]]=number(1);
4464  }
4465  intvec opgt=option(get);
4466  option(redSB);
4467  matrix L=transpose(syz(concat(A,-b)));
4468  if (nrows(L)==1)
4469  {
4470    if (L[1,n+1]<>0)
4471    {
4472      poly pivote=L[1,n+1];
4473      matrix sol=submat(L,1..1,1..n);
4474      if (pivote<>1)
4475      {
4476        sol=(number(1)/number(pivote))*sol;
4477      }
4478      // check at least that the number of comitted errors is less than half
4479      //     the Goppa distance
4480      // imposing Hamming_wt(sol)<=K[n+3][1] would be more correct, but maybe
4481      //     is too strong
4482      // on the other hand, if Hamming_wt(sol) is too large the decoding may
4483      //     not be acceptable
4484      if ( Hamming_wt(sol) <= (K[n+3][2]-1)/2 )
4485      {
4486        option(set,opgt);
4487        return(y-sol);
4488      }
4489      else
4490      {
4491        dbprint(printlevel+3,"? non-acceptable decoding ?");
4492      }
4493    }
4494    else
4495    {
4496      dbprint(printlevel+3,"? no solution found ?");
4497    }
4498  }
4499  else
4500  {
4501    dbprint(printlevel+3,"? non-unique solution ?");
4502  }
4503  option(set,opgt);
4504  dbprint(printlevel+3,"? too many errors occur, 0-matrix returned ?");
4505  matrix answer;
4506  return(answer);
4507}
4508example
4509{
4510  "EXAMPLE:";  echo = 2;
4511  int plevel=printlevel;
4512  printlevel=-1;
4513  ring s=2,(x,y),lp;
4514  list HC=Adj_div(x3+y2+y);
4515  HC=NSplaces(1,HC);
4516  HC=extcurve(2,HC);
4517  def ER=HC[1][4];
4518  setring ER;
4519  intvec G=5;      // the rational divisor G = 5*HC[3][1]
4520  intvec D=2..9;   // D = sum of the rational places no. 2..9 over F_4
4521  // construct the corresp. residual AG code of type [8,3,>=5] over F_4:
4522  matrix C=AGcode_Omega(G,D,HC);
4523  // we can correct 1 error and the genus is 1, thus F must have degree 2
4524  // and support disjoint from that of D
4525  intvec F=2;
4526  list SV=prepSV(G,D,F,HC);
4527  // now we produce 1 error on the zero-codeword :
4528  matrix y[1][8];
4529  y[1,3]=a;
4530  // and then we decode :
4531  print(decodeSV(y,SV));
4532  printlevel=plevel;
4533}
4534///////////////////////////////////////////////////////////////////////////////
4535
4536proc sys_code (matrix C)
4537"USAGE:    sys_code(C); C is a matrix of constants
4538RETURN:   list L with:
4539   @format
4540   L[1] is the generator matrix in standard form of an equivalent code,
4541   L[2] is the parity check matrix in standard form of such code,
4542   L[3] is an intvec which represents the needed permutation.
4543   @end format
4544NOTE:     Computes a systematic code which is equivalent to the given one.@*
4545          The input should be a matrix of numbers.@*
4546          The output has to be interpreted as follows: if the input was
4547          the generator matrix of an AG code then one should apply the
4548          permutation L[3] to the divisor D of rational points by means
4549          of @code{permute_L(D,L[3]);} before continuing to work with the
4550          code (for instance, if you want to use the systematic encoding
4551          together with a decoding algorithm).
4552KEYWORDS: linear code, systematic
4553SEE ALSO: permute_L, AGcode_Omega, prepSV
4554EXAMPLE:  example sys_code; shows an example
4555"
4556{
4557  // computes a systematic code which is equivalent to that given by C
4558  int i,j,k,l,h,r;
4559  int m=nrows(C);
4560  int n=ncols(C);
4561  int mr=m;
4562  matrix A[m][n]=C;
4563  poly c,p;
4564  list corners=list();
4565  if(m>n)
4566  {
4567    mr=n;
4568  }
4569  // first the matrix A will be reduced with elementary operations by rows
4570  for(i=1;i<=mr;i=i+1)
4571  {
4572    if((i+l)>n)
4573    {
4574      // the process is finished
4575      break;
4576    }
4577    // look for a non-zero element
4578    if(A[i,i+l]==0)
4579    {
4580      h=i;
4581      p=0;
4582      for (j=i+1;j<=m;j=j+1)
4583      {
4584        c=A[j,i+l];
4585        if (c!=0)
4586        {
4587          p=c;
4588          h=j;
4589          break;
4590        }
4591      }
4592      if (h!=i)
4593      {
4594        // permutation of rows i and h
4595        for (j=1;j<=n;j=j+1)
4596        {
4597          c=A[i,j];
4598          A[i,j]=A[h,j];
4599          A[h,j]=c;
4600        }
4601      }
4602      if(p==0)
4603      {
4604        // non-zero element not found in the current column
4605        l=l+1;
4606        continue;
4607      }
4608    }
4609    // non-zero element was found in "strategic" position
4610    corners[i]=i+l;
4611    // make zeros below that position
4612    for(j=i+1;j<=m;j=j+1)
4613    {
4614      c=A[j,i+l]/A[i,i+l];
4615      for(k=i+l+1;k<=n;k=k+1)
4616      {
4617        A[j,k]=A[j,k]-A[i,k]*c;
4618      }
4619      A[j,i+l]=0;
4620      A[j,i]=0;
4621    }
4622    // the rank is at least r
4623    // when the process stops the last r is actually the true rank of A=a
4624    r=i;
4625  }
4626  if (r<m)
4627  {
4628    ERROR("the given matrix does not have maximum rank");
4629  }
4630  // set the corners to the beginning and construct the permutation
4631  intvec PCols=1..n;
4632  for (j=1;j<=m;j=j+1)
4633  {
4634    if (corners[j]>j)
4635    {
4636      // interchange columns j and corners[j]
4637      for (i=1;i<=m;i=i+1)
4638      {
4639        c=A[i,j];
4640        A[i,j]=A[i,corners[j]];
4641        A[i,corners[j]]=c;
4642      }
4643      k=PCols[j];
4644      PCols[j]=PCols[corners[j]];
4645      PCols[corners[j]]=k;
4646    }
4647  }
4648  // convert the diagonal into ones
4649  for (i=1;i<=m;i=i+1)
4650  {
4651    for (j=i;j<=n;j=j+1)
4652    {
4653      A[i,j]=A[i,j]/A[i,i];
4654    }
4655  }
4656  // complete a block with the identity matrix
4657  for (k=1;k<m;k=k+1)
4658  {
4659    for (i=k+1;i<=m;i=i+1)
4660    {
4661      for (j=i;j<=n;j=j+1)
4662      {
4663        A[k,j]=A[k,j]-A[i,j]*A[k,i];
4664      }
4665    }
4666  }
4667  // compute a parity-check matrix in standard form
4668  matrix B=concat(-transpose(submat(A,1..m,m+1..n)),diag(1,n-m));
4669  list L=list();
4670  L[1]=A;
4671  L[2]=B;
4672  L[3]=PCols;
4673  return(L);
4674}
4675example
4676{
4677  "EXAMPLE:";  echo = 2;
4678  ring s=3,T,lp;
4679  matrix C[2][5]=0,1,0,1,1,0,1,0,0,1;
4680  print(C);
4681  list L=sys_code(C);
4682  L[3];
4683  // here is the generator matrix in standard form
4684  print(L[1]);
4685  // here is the control matrix in standard form
4686  print(L[2]);
4687  // we can check that both codes are dual to each other
4688  print(L[1]*transpose(L[2]));
4689}
4690///////////////////////////////////////////////////////////////////////////////
4691
4692/*
4693// ============================================================================
4694// *******       ADDITIONAL INFORMATION ABOUT THE LIBRARY              ********
4695// ============================================================================
4696
4697A SINGULAR library for plane curves, Weierstrass semigroups and AG codes
4698Also available via http://wmatem.eis.uva.es/~ignfar/singular/
4699
4700PREVIOUS WARNINGS :
4701
4702   (1) The procedures will work only in positive characteristic
4703   (2) The base field must be prime (this may change in the future)
4704           This limitation is not too serious, since in coding theory
4705           the used curves are usually (if not always) defined over a
4706           prime field, and extensions are only considered for
4707           evaluating functions in a field with many points;
4708           by the same reason, auxiliary divisors are usually defined
4709           over the prime field,
4710           with the exception of that consisting of "rational points"
4711   (3) The curve must be absolutely irreducible (but it is not checked)
4712   (4) Only (algebraic projective) plane curves are considered
4713
4714GENERAL CONCEPTS :
4715
4716   (1) An affine point P is represented by a std of a prime ideal,
4717       and an intvec containing the position of the places above P
4718       in the list of Places; if the point is at infinity, the ideal is
4719       changed by a homogeneous irreducible polynomial in two variables
4720   (2) A place is represented by :
4721       a base point (list of homogeneous coordinates),
4722       a local equation for the curve at the base point,
4723       a Hamburger-Noether expansion of the corresponding branch,
4724       and a local parametrization (in "t") of such branch; everything is
4725       stored in a local ring "_[5][d][1]", d being the degree of the place,
4726       by means of lists "POINTS,LOC_EQS,BRANCHES,PARAMETRIZATIONS", and
4727       the degrees of the base points corresponding to the places in the
4728       ring "_[5][d][1]" are stored in an intvec "_[5][d][2]"
4729   (3) A divisor is represented by an intvec, where the integer at the
4730       position i means the coefficient of the i-th place in the divisor
4731   (4) Rational functions are represented by numerator/denominator
4732       in form of ideals with two homogeneous generators
4733
4734OUTLINE/EXAMPLE OF THE USE OF THE LIBRARY :
4735
4736Plane curves :
4737
4738      (1.0) ring s=p,(x,y[,z]),lp;
4739
4740            Notice that if you use 3 variables, then the equation
4741            of the curve is assumed to be a homogeneous polynomial.
4742
4743      (1.1) list CURVE=Adj_div(equation);
4744
4745            In CURVE[3] are listed all the (singular closed) places
4746            with their corresponding degrees; thus, you can now decide
4747            how many other points you want to compute with NSplaces.
4748
4749      (1.2) CURVE=NSplaces(range,CURVE);
4750
4751            See help NSplaces tp know the meaning of "range";
4752            for instance, if you have that the singular places
4753            have degree 2 at most and you want all the places
4754            up to degree 5, you must write range=3.
4755
4756      (1.3) CURVE=extcurve(extension,CURVE);
4757
4758            The rational places over the extension are ranged in
4759            the ring CURVE[1][5] with the following rules:
4760
4761                (i)    all the representatives of the same closed point
4762                       are listed in consecutive positions;
4763                (ii)   if deg(P)<deg(Q), then the representatives of P
4764                       are listed before those of Q;
4765                (iii)  if two closed points P,Q have the same degree,
4766                       then the representatives of P are listed before
4767                       if P appears before in the list CURVE[3].
4768
4769Rational functions :
4770
4771      (2.0) def R=CURVE[1][2];
4772            setring R;
4773      (2.1) list LG=BrillNoether(intvec divisor,CURVE);
4774      (2.2) list WS=Weierstrass(int place,int bound,CURVE);
4775
4776Algebraic Geometry codes :
4777
4778      (3.0) def ER=CURVE[1][4];    // if extension>1; else use R instead
4779            setring ER;
4780
4781            Now you have to decide the divisor G and the sequence of
4782            rational points D to use for constructing the codes;
4783            first notice that the syntax for G and D is different:
4784
4785                (a) G is a divisor supported on the closed places of
4786                    CURVE[3], and you must say just the coefficient
4787                    of each such a point; for example, G=2,0,-1 means
4788                    2 times the place 1 minus 1 times the place 3.
4789
4790                (b) D is a sequence of rational points (all different
4791                    and counted 1 time each one), whose data are read
4792                    from the lists inside CURVE[1][5] and now you must
4793                    say just the order how you range the chosen point;
4794                    for example, D=2,4,1 means that you choose the
4795                    rational places 1,2,4 and you range them as 2,4,1.
4796
4797      (3.1) matrix C=AGcode_L(divisor,places,CURVE);
4798
4799      (3.2) AGcode_Omega(divisor,places,CURVE);
4800
4801            In the same way as for defining the codes, the auxiliary
4802            divisor F must have disjoint support to that of D, and
4803            its degree has to be given by a formula (see help prepSV).
4804
4805      (3.3) list SV=prepSV(divisor,places,auxiliary_divisor,CURVE);
4806
4807      (3.4) decodeSV(word,SV);
4808
4809Special Issues :
4810
4811      (I)  AG codes with systematic encoding :
4812
4813              matrix C=AGcode_Omega(G,D,CURVE);
4814              list CODE=sys_code(G);
4815              C=CODE[1];               // generator matrix in standard form
4816              D=permute_L(D,CODE[3]);  // suitable permutation of coordinates
4817              list SV=prepSV(G,D,F,CURVE);
4818              SV[1]=CODE[2];           // parity-check matrix in standard form
4819
4820      (II) Using the true minimum distance d for decoding :
4821
4822              matrix C=AGcode_Omega(G,D,CURVE);
4823              int n=size(D);
4824              list SV=prepSV(G,D,F,CURVE);
4825              SV[n+3][2]=d;            // then use decodeSV(y,SV);
4826
4827
4828// ============================================================================
4829// ***   Some "macros" with typical examples of curves in Coding Theory    ****
4830// ============================================================================
4831
4832
4833proc Klein ()
4834{
4835  list KLEIN=Adj_div(x3y+y3+x);
4836  KLEIN=NSplaces(2,KLEIN);
4837  KLEIN=extcurve(3,KLEIN);
4838  dbprint(printlevel+1,"Klein quartic over F_8 successfully constructed");
4839  return(KLEIN);
4840}
4841
4842proc Hermite (int m)
4843{
4844  int p=char(basering);
4845  int r=p^m;
4846  list HERMITE=Adj_div(y^r+y-x^(r+1));
4847  HERMITE=NSplaces(2*m-1,HERMITE);
4848  HERMITE=extcurve(2*m,HERMITE);
4849  dbprint(printlevel+1,"Hermitian curve over F_("+string(r)+"^2)
4850                          successfully constructed");
4851  return(HERMITE);
4852}
4853
4854proc Suzuki ()
4855{
4856  list SUZUKI=Adj_div(x10+x3+y8+y);
4857  SUZUKI=NSplaces(2,SUZUKI);
4858  SUZUKI=extcurve(3,SUZUKI);
4859  dbprint(printlevel+1,"Suzuki curve over F_8 successfully constructed");
4860  return(SUZUKI);
4861}
4862
4863// **** Other interesting examples :
4864
4865// A hyperelliptic curve with 33 rational points over F_16
4866
4867list CURVE=Adj_div(x5+y2+y);
4868CURVE=NSplaces(3,CURVE);
4869CURVE=extcurve(4,CURVE);
4870
4871// A nice curve with 113 rational points over F_64
4872
4873list CURVE=Adj_div(y9+y8+xy6+x2y3+y2+x3);
4874CURVE=NSplaces(4,CURVE);
4875CURVE=extcurve(6,CURVE);
4876
4877*/
4878;
Note: See TracBrowser for help on using the repository browser.