source: git/Singular/LIB/brnoeth.lib @ 489a49

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