version="$Id: brnoeth.lib,v 1.11.2.6 2003/10/08 13:40:52 lossen Exp $"; category="Coding theory"; info=" LIBRARY: brnoeth.lib Brill-Noether Algorithm, Weierstrass-SG and AG-codes AUTHORS: Jose Ignacio Farran Martin, ignfar@eis.uva.es Christoph Lossen, lossen@mathematik.uni-kl.de OVERVIEW: Implementation of the Brill-Noether algorithm for solving the Riemann-Roch problem and applications in Algebraic Geometry codes. The computation of Weierstrass semigroups is also implemented.@* The procedures are intended only for plane (singular) curves defined over a prime field of positive characteristic.@* For more information about the library see the end of the file brnoeth.lib. MAIN PROCEDURES: Adj_div(f); computes the conductor of a curve NSplaces(h,A); computes non-singular places with given degrees BrillNoether(D,C); computes a vector space basis of the linear system L(D) Weierstrass(P,m,C); computes the Weierstrass semigroup of C at P up to m extcurve(d,C); extends the curve C to an extension of degree d AGcode_L(G,D,E); computes the evaluation AG code with divisors G and D AGcode_Omega(G,D,E); computes the residual AG code with divisors G and D prepSV(G,D,F,E); preprocessing for the basic decoding algorithm decodeSV(y,K); decoding of a word with the basic decoding algorithm AUXILIARY PROCEDURES: closed_points(I); computes the zero-set of a zero-dim. ideal in 2 vars dual_code(C); computes the dual code sys_code(C); computes an equivalent systematic code permute_L(L,P); applies a permutation to a list SEE ALSO: hnoether_lib, triang_lib KEYWORDS: Weierstrass semigroup; Algebraic Geometry codes; Brill-Noether algorithm "; LIB "matrix.lib"; LIB "triang.lib"; LIB "hnoether.lib"; LIB "inout.lib"; /////////////////////////////////////////////////////////////////////////////// // ********************************************************** // * POINTS, PLACES AND DIVISORS OF (SINGULAR) PLANE CURVES * // ********************************************************** proc closed_points (ideal I) "USAGE: closed_points(I); I an ideal RETURN: list of prime ideals (each a Groebner basis), corresponding to the (distinct affine closed) points of V(I) NOTE: The ideal must have dimension 0, the basering must have 2 variables, the ordering must be lp, and the base field must be finite and prime.@* It might be convenient to set the option(redSB) in advance. SEE ALSO: triang_lib EXAMPLE: example closed_points; shows an example " { ideal II=std(I); if (II==1) { return(list()); } list TL=triangMH(II); int s=size(TL); list L=list(); int i,j,k; ideal Facts; poly G2; for (i=1;i<=s;i=i+1) { Facts=factorize(TL[i][1],1); k=size(Facts); G2=TL[i][2]; for (j=1;j<=k;j=j+1) { L=L+pd2(Facts[j],G2); } } // eliminate possible repetitions s=size(L); list LP=list(); LP[1]=std(L[1]); int counter=1; for (i=2;i<=s;i=i+1) { if (isPinL(L[i],LP)==0) { counter=counter+1; LP[counter]=std(L[i]); } } return(LP); } example { "EXAMPLE:"; echo = 2; ring s=2,(x,y),lp; // this is just the affine plane over F_4 : ideal I=x4+x,y4+y; list L=closed_points(I); // and here you have all the points : L; } /////////////////////////////////////////////////////////////////////////////// static proc pd2 (poly g1,poly g2) { // If g1,g2 is a std. resp. lex. in (x,y) then the procedure // factorizes g2 in the "extension given by g1" // (then g1 must be irreducible) and returns a list of // ideals with always g1 as first component and the // distinct factors of g2 as second components list L=list(); ideal J=g1; int i,s; if (deg(g1)==1) { poly A=-subst(g1,var(2),0); poly B=subst(g2,var(2),A); ideal facts=factorize(B,1); s=size(facts); for (i=1;i<=s;i=i+1) { J[2]=facts[i]; L[i]=J; } } else { def BR=basering; poly A=g1; poly B=g2; ring raux1=char(basering),(x,y,a),lp; poly G; ring raux2=(char(basering),a),(x,y),lp; map psi=BR,x,a; minpoly=number(psi(A)); poly f=psi(B); ideal facts=factorize(f,1); s=size(facts); poly g; string sg; for (i=1;i<=s;i=i+1) { g=facts[i]; sg=string(g); setring raux1; execute("G="+sg+";"); G=subst(G,a,y); setring BR; map ppssii=raux1,var(1),var(2),0; J[2]=ppssii(G); L[i]=J; kill(ppssii); setring raux2; } setring BR; kill(raux1,raux2); } return(L); } /////////////////////////////////////////////////////////////////////////////// static proc isPinL (ideal P,list L) { // checks if a (plane) point P is in a list of (plane) points L // by just comparing generators // it works only if all (prime) ideals are given in a "canonical way", // namely: // the first generator is monic and irreducible, // and depends only on the second variable, // and the second one is monic in the first variable // and irreducible over the field extension determined by // the second variable and the first generator as minpoly int s=size(L); int i; for (i=1;i<=s;i=i+1) { if ( P[1]==L[i][1] && P[2]==L[i][2] ) { return(1); } } return(0); } /////////////////////////////////////////////////////////////////////////////// static proc s_locus (poly f) { // computes : ideal of affine singular locus // the equation f must be affine // warning : if there is an error message then the output is "none" // option(redSB) is convenient to be set in advance ideal I=f,jacob(f); I=std(I); if (dim(I)>0) { // dimension check (has to be 0) ERROR("something was wrong; possibly non-reduced curve"); } else { return(I); } } /////////////////////////////////////////////////////////////////////////////// static proc curve (poly f) "USAGE: curve(f), where f is a polynomial (affine or projective) CREATE: poly CHI in both rings aff_r=p,(x,y),lp and Proj_R=p,(x,y,z),lp also ideal (std) Aff_SLocus of affine singular locus in the ring aff_r RETURN: list (size 3) with two rings aff_r,Proj_R and an integer deg(f) NOTE: f must be absolutely irreducible, but this is not checked it is not implemented yet for extensions of prime fields " { def base_r=basering; ring aff_r=char(basering),(x,y),lp; ring Proj_R=char(basering),(x,y,z),lp; setring base_r; int degX=deg(f); if (nvars(basering)==2) { setring aff_r; map embpol=base_r,x,y; poly CHI=embpol(f); export(CHI); kill(embpol); ideal Aff_SLocus=s_locus(CHI); export(Aff_SLocus); setring Proj_R; poly CHI=homog(imap(aff_r,CHI),z); export(CHI); setring base_r; list L=list(); L[1]=aff_r; L[2]=Proj_R; L[3]=degX; kill(aff_r,Proj_R); return(L); } if (nvars(basering)==3) { setring Proj_R; map embpol=base_r,x,y,z; poly CHI=embpol(f); export(CHI); kill(embpol); string s=string(subst(CHI,z,1)); setring aff_r; execute("poly CHI="+s+";"); export(CHI); ideal Aff_SLocus=s_locus(CHI); export(Aff_SLocus); setring base_r; list L=list(); L[1]=aff_r; L[2]=Proj_R; L[3]=degX; kill(aff_r,Proj_R); return(L); } ERROR("basering must have 2 or 3 variables"); } /////////////////////////////////////////////////////////////////////////////// static proc Aff_SL (ideal ISL) { // computes : affine singular (closed) points as a list of lists of // prime ideals and intvec (for storing the places over each point) // the ideal ISL=s_locus(CHI) is assumed to be computed in advance for // a plane curve CHI, and it must be given by a standard basis // for our purpose the function must called with the "global" ideal // "Aff_SLocus" list SL=list(); ideal I=ISL; if ( I != 1 ) { list L=list(); ideal aux; intvec iv; int i,s; L=closed_points(I); s=size(L); for (i=1;i<=s;i=i+1) { aux=std(L[i]); SL[i]=list(aux,iv); } } return(SL); } /////////////////////////////////////////////////////////////////////////////// static proc inf_P (poly f) { // computes : all (closed) points at infinity as homogeneous polynomials // output : two lists with respectively singular and non-singular points intvec iv; def base_r=basering; ring r_auxz=char(basering),(x,y,z),lp; poly f=imap(base_r,f); poly F=homog(f,z); // equation of projective curve poly f_inf=subst(F,z,0); setring base_r; poly f_inf=imap(r_auxz,f_inf); ideal I=factorize(f_inf,1); // points at infinity as homogeneous // polynomials int s=size(I); int i; list IP_S=list(); // for singular points at infinity list IP_NS=list(); // for non-singular points at infinity int counter_S; int counter_NS; poly aux; for (i=1;i<=s;i=i+1) { aux=subst(I[i],y,1); if (aux==1) { // the point is (1:0:0) setring r_auxz; poly f_yz=subst(F,x,1); if ( subst(subst(diff(f_yz,y),y,0),z,0)==0 && subst(subst(diff(f_yz,z),y,0),z,0)==0 ) { // the point is singular counter_S=counter_S+1; kill(f_yz); setring base_r; IP_S[counter_S]=list(I[i],iv); } else { // the point is non-singular counter_NS=counter_NS+1; kill(f_yz); setring base_r; IP_NS[counter_NS]=list(I[i],iv); } } else { // the point is (a:1:0) | a is root of aux if (deg(aux)==1) { // the point is rational and no field extension is needed setring r_auxz; poly f_xz=subst(F,y,1); poly aux=imap(base_r,aux); number A=-number(subst(aux,x,0)); map phi=r_auxz,x+A,0,z; poly f_origin=phi(f_xz); if ( subst(subst(diff(f_origin,x),x,0),z,0)==0 && subst(subst(diff(f_origin,z),x,0),z,0)==0 ) { // the point is singular counter_S=counter_S+1; kill(f_xz,aux,A,phi,f_origin); setring base_r; IP_S[counter_S]=list(I[i],iv); } else { // the point is non-singular counter_NS=counter_NS+1; kill(f_xz,aux,A,phi,f_origin); setring base_r; IP_NS[counter_NS]=list(I[i],iv); } } else { // the point is non-rational and a field extension with minpoly=aux // is needed ring r_ext=(char(basering),a),(x,y,z),lp; poly F=imap(r_auxz,F); poly f_xz=subst(F,y,1); poly aux=imap(base_r,aux); minpoly=number(subst(aux,x,a)); map phi=r_ext,x+a,0,z; poly f_origin=phi(f_xz); if ( subst(subst(diff(f_origin,x),x,0),z,0)==0 && subst(subst(diff(f_origin,z),x,0),z,0)==0 ) { // the point is singular counter_S=counter_S+1; setring base_r; kill(r_ext); IP_S[counter_S]=list(I[i],iv); } else { // the point is non-singular counter_NS=counter_NS+1; setring base_r; kill(r_ext); IP_NS[counter_NS]=list(I[i],iv); } } } } kill(r_auxz); return(list(IP_S,IP_NS)); } /////////////////////////////////////////////////////////////////////////////// static proc closed_points_ext (poly f,int d,ideal SL) { // computes : (closed affine non-singular) points over an extension of // degree d // remark(1) : singular points are supposed to be listed appart // remark(2) : std SL=s_locus(f) is supposed to be computed in advance // remark(3) : ideal SL is used to remove those points which are singular // output : list of list of prime ideals with an intvec for storing the // places int Q=char(basering)^d; // cardinality of the extension field ideal I=f,x^Q-x,y^Q-y; // ideal of the searched points I=std(I); if (I==1) { return(list()); } list LP=list(); int m=size(SL); list L=list(); ideal aux; intvec iv; int i,s,j,counter; L=closed_points(I); s=size(L); for (i=1;i<=s;i=i+1) { aux=std(L[i]); for (j=1;j<=m;j=j+1) { // check if singular i.e. if SL is contained in aux if ( NF(SL[j],aux) != 0 ) { counter=counter+1; LP[counter]=list(aux,iv); break; } } } return(LP); } /////////////////////////////////////////////////////////////////////////////// static proc degree_P (list P) "USAGE: degree_P(P), where P is either a polynomial or an ideal RETURN: integer with the degree of the closed point given by P SEE ALSO: closed_points NOTE: If P is a (homogeneous irreducible) polynomial the point is at infinity, and if P is a (prime) ideal the points is affine, and the ideal must be given by 2 generators: the first one irreducible and depending only on y, and the second one irreducible over the extension given by y with the first generator as minimal polynomial " { // computes : the degree of a given point // remark(1) : if the input is (irreducible homogeneous) poly => the point // is at infinity // remark(2) : it the input is (std. resp. lp. prime) ideal => the point is // affine if (typeof(P[1])=="ideal") { if (size(P[1])==2) { int d=deg(P[1][1]); poly aux=subst(P[1][2],y,1); d=d*deg(aux); return(d); } else { // this should not happen in principle ERROR("non-valid parameter"); } } else { if (typeof(P[1])=="poly") { return(deg(P[1])); } else { ERROR("parameter must have a poly or ideal in the first component"); } } } /////////////////////////////////////////////////////////////////////////////// static proc closed_points_deg (poly f,int d,ideal SL) { // computes : (closed affine non-singular) points of degree d // remark(1) : singular points are supposed to be listed appart // remark(2) : std SL=s_locus(f) is supposed to be computed in advance list L=closed_points_ext(f,d,SL); int s=size(L); int i,counter; list LP=list(); for (i=1;i<=s;i=i+1) { if (degree_P(L[i])==d) { counter=counter+1; LP[counter]=L[i]; } } return(LP); } /////////////////////////////////////////////////////////////////////////////// static proc subset (ideal I,ideal J) { // checks wether I is contained in J and returns a boolean // remark : J is assumed to be given by a standard basis int s=ncols(I); int i; for (i=1;i<=s;i=i+1) { if ( NF(I[i],std(J)) != 0 ) { return(0); } } return(1); } /////////////////////////////////////////////////////////////////////////////// static proc belongs (list P,ideal I) { // checks if affine point P is contained in V(I) and returns a boolean // remark : P[1] is assumed to be an ideal given by a standard basis if (typeof(P[1])=="ideal") { return(subset(I,P[1])); } else { ERROR("first argument must be an affine point"); } } /////////////////////////////////////////////////////////////////////////////// static proc equals (ideal I,ideal J) { // checks if I is equal to J and returns a boolean // remark : I and J are assumed to be given by a standard basis int answer=0; if (subset(I,J)==1) { if (subset(J,I)==1) { answer=1; } } return(answer); } /////////////////////////////////////////////////////////////////////////////// static proc isInLP (ideal P,list LP) { // checks if affine point P is a list LP and returns either its position or // zero // remark : all points in LP and P itself are assumed to be given by a // standard basis // warning : the procedure does not check whether the points are affine or // not int s=size(LP); if (s==0) { return(0); } int i; for (i=1;i<=s;i=i+1) { if (equals(P,LP[i][1])==1) { return(i); } } return(0); } /////////////////////////////////////////////////////////////////////////////// static proc res_deg () { // computes the residual degree of the basering with respect to its prime // field // warning : minpoly must depend on a parameter called "a" int ext; string s_m=string(minpoly); if (s_m=="0") { ext=1; } else { ring auxr=char(basering),a,lp; execute("poly minp="+s_m+";"); ext=deg(minp); } return(ext); } /////////////////////////////////////////////////////////////////////////////// static proc Frobenius (etwas,int r) { // applies the Frobenius map over F_{p^r} to an object defined over an // extension of such field // usually it is called with r=1, i.e. the Frobenius map over the prime // field F_p // returns always an object of the same type, and works correctly on // numbers, polynomials, ideals, matrices or lists of the above types // maybe : types vector and module should be added in the future, but they // are not needed now int q=char(basering)^r; if (typeof(etwas)=="number") { return(etwas^q); } if (typeof(etwas)=="poly") { int s=size(etwas); poly f; int i; for (i=1;i<=s;i=i+1) { f=f+(leadcoef(etwas[i])^q)*leadmonom(etwas[i]); } return(f); } if (typeof(etwas)=="ideal") { int s=ncols(etwas); ideal I; int i; for (i=1;i<=s;i=i+1) { I[i]=Frobenius(etwas[i],r); } return(I); } if (typeof(etwas)=="matrix") { int m=nrows(etwas); int n=ncols(etwas); matrix A[m][n]; int i,j; for (i=1;i<=m;i=i+1) { for (j=1;j<=n;j=j+1) { A[i,j]=Frobenius(etwas[i,j],r); } } return(A); } if (typeof(etwas)=="list") { int s=size(etwas); list L=list(); int i; for (i=1;i<=s;i=i+1) { if (typeof(etwas[i])<>"none") { L[i]=Frobenius(etwas[i],r); } } return(L); } return(etwas); } /////////////////////////////////////////////////////////////////////////////// static proc conj_b (list L,int r) { // applies the Frobenius map over F_{p^r} to a list of type HNE defined over // a larger extension // when r=1 it turns to be the Frobenius map over the prime field F_{p} // returns : a list of type HNE which is either conjugate of the input or // the same list in case of L being actually defined over the base field // F_{p^r} int p=char(basering); int Q=p^r; list LL=list(); int m=nrows(L[1]); int n=ncols(L[1]); matrix A[m][n]; poly f; poly aux; int i,j; for (i=1;i<=m;i=i+1) { for (j=1;j<=n;j=j+1) { aux=L[1][i,j]; if (aux<>x) { A[i,j]=aux^Q; } else { A[i,j]=aux; break; } } } m=size(L[4]); for (i=1;i<=m;i=i+1) { f=f+(leadcoef(L[4][i])^Q)*leadmonom(L[4][i]); } LL[1]=A; LL[2]=L[2]; LL[3]=L[3]; LL[4]=f; return(LL); } /////////////////////////////////////////////////////////////////////////////// static proc grad_b (list L,int r) { // computes the degree of a list of type HNE which is actually defined over // F_{p^r} eventhough it is given in an extension of such field int gr=1; int rd=res_deg() div r; list LL=L; int i; for (i=1;i<=rd;i=i+1) { LL=conj_b(LL,r); if ( LL[1]==L[1] && LL[4]==L[4] ) { break; } else { gr=gr+1; } } return(gr); } /////////////////////////////////////////////////////////////////////////////// static proc conj_bs (list L,int r) { // computes all the conjugates over F_{p^r} of a list of type HNE defined // over an extension // returns : a list of lists of type HNE, where the first one is the input // list // remark : notice that the degree of the branch is then the size of the // output list branches=list(); int gr=1; branches[1]=L; int rd=res_deg() div r; list LL=L; int i; for (i=1;i<=rd;i=i+1) { LL=conj_b(LL,r); if ( LL[1]==L[1] && LL[4]==L[4] ) { break; } else { gr=gr+1; branches[gr]=LL; } } return(branches); } /////////////////////////////////////////////////////////////////////////////// static proc subfield (sf) { // writes the generator "a" of a subfield of the coefficients field of // basering in terms of the the current generator (also called "a") as a // string sf is an existing ring whose coefficient field is such a subfield // warning : in basering there must be a variable called "x" and subfield // must not be prime def base_r=basering; string new_m=string(minpoly); setring sf; string old_m=string(minpoly); if (old_m==new_m) { setring base_r; return("a"); } else { if (old_m<>string(0)) { ring auxring=char(basering),(a,x),lp; execute("poly mpol="+old_m+";"); mpol=subst(mpol,a,x); setring base_r; poly mpol=imap(auxring,mpol); kill(auxring); string answer="? error : non-primitive element"; int r=res_deg(); int q=char(basering)^r; int i; number b; for (i=1;i<=q-2;i=i+1) { b=a^i; if (subst(mpol,x,b)==0) { answer=string(b); break; } } if (answer<>"? error : non-primitive element") { return(answer); } else { // list all the elements of the finite field F_q int p=char(basering); list FF1,FF2; for (i=0;i1&2 the point is assumed non-singular and the local conductor // should be zero list PP=list(); if (Pp[1]==0) { if (Pp[2]==0) { PP=Aff_SPoints[Pp[3]]; } if (Pp[2]==1) { PP=Inf_Points[1][Pp[3]]; } if (Pp[2]==2) { PP=Inf_Points[2][Pp[3]]; } } else { PP=Aff_Points(Pp[2])[Pp[3]]; } if (PP[2]<>0) { return(CURVE); } intvec PtoPl; def base_r=basering; int ext1; list Places=CURVE[3]; intvec Conductor=CURVE[4]; list update_CURVE=CURVE; if (typeof(PP[1])=="ideal") { ideal P=PP[1]; if (size(P)==2) { int d=deg(P[1]); poly aux=subst(P[2],y,1); d=d*deg(aux); ext1=d; // the point is (A:B:1) but one must distinguish several cases // P is assumed to be a std. resp. "(x,y),lp" and thus P[1] depends // only on "y" if (d==1) { // the point is rational number B=-number(subst(P[1],y,0)); poly aux2=subst(P[2],y,B); number A=-number(subst(aux2,x,0)); // the point is (A:B:1) ring local_aux=char(basering),(x,y),ls; number coord@1=imap(base_r,A); number coord@2=imap(base_r,B); number coord@3=number(1); map phi=base_r,x+coord@1,y+coord@2; poly CHI=phi(CHI); } else { if (deg(P[1])==1) { // the point is non-rational but the second component needs no // field extension number B=-number(subst(P[1],y,0)); poly aux2=subst(P[2],y,B); // the point has degree d>1 // careful : the parameter will be called "a" anyway ring local_aux=(char(basering),a),(x,y),ls; map psi=base_r,a,0; minpoly=number(psi(aux2)); number coord@1=a; number coord@2=imap(base_r,B); number coord@3=number(1); // the point is (a:B:1) map phi=base_r,x+a,y+coord@2; poly CHI=phi(CHI); } else { if (deg(subst(P[2],y,1))==1) { // the point is non-rational but the needed minpoly is just P[1] // careful : the parameter will be called "a" anyway poly P1=P[1]; poly P2=P[2]; ring local_aux=(char(basering),a),(x,y),ls; map psi=base_r,0,a; minpoly=number(psi(P1)); // the point looks like (A:a:1) // A is computed by substituting y=a in P[2] poly aux1=imap(base_r,P2); poly aux2=subst(aux1,y,a); number coord@1=-number(subst(aux2,x,0)); number coord@2=a; number coord@3=number(1); map phi=base_r,x+coord@1,y+a; poly CHI=phi(CHI); } else { // this is the most complicated case of non-rational point // firstly : construct an extension of degree d and guess the // minpoly poly P1=P[1]; poly P2=P[2]; int p=char(basering); int Q=p^d; ring aux_r=(Q,a),(x,y,t),ls; string minpoly_string=string(minpoly); ring local_aux=(char(basering),a),(x,y),ls; execute("minpoly="+minpoly_string+";"); // secondly : compute one root of P[1] poly P_1=imap(base_r,P1); poly P_2=imap(base_r,P2); ideal factors1=factorize(P_1,1); // hopefully this works !!!! number coord@2=-number(subst(factors1[1],y,0)); // thirdly : compute one of the first components for the above root poly P_0=subst(P_2,y,coord@2); ideal factors2=factorize(P_0,1); // hopefully this works !!!! number coord@1=-number(subst(factors2[1],x,0)); number coord@3=number(1); map phi=base_r,x+coord@1,y+coord@2; poly CHI=phi(CHI); kill(aux_r); } } } } else { // this should not happen in principle ERROR("non-valid parameter"); } } else { if (typeof(PP[1])=="poly") { poly P=PP[1]; ring r_auxz=char(basering),(x,y,z),lp; poly CHI=imap(base_r,CHI); CHI=homog(CHI,z); setring base_r; poly aux=subst(P,y,1); if (aux==1) { // the point is (1:0:0) ring local_aux=char(basering),(x,y),ls; number coord@1=number(1); number coord@2=number(0); number coord@3=number(0); map Phi=r_auxz,1,x,y; poly CHI=Phi(CHI); ext1=1; } else { // the point is (A:1:0) where A is a root of aux int d=deg(aux); ext1=d; if (d==1) { // the point is rational number A=-number(subst(aux,x,0)); ring local_aux=char(basering),(x,y),ls; number coord@1=imap(base_r,A); number coord@2=number(1); number coord@3=number(0); map Phi=r_auxz,x+coord@1,1,y; poly CHI=Phi(CHI); } else { // the point has degree d>1 // careful : the parameter will be called "a" anyway ring local_aux=(char(basering),a),(x,y),ls; map psi=base_r,a,1; minpoly=number(psi(P)); number coord@1=a; number coord@2=number(1); number coord@3=number(0); map Phi=r_auxz,x+a,1,y; poly CHI=Phi(CHI); } } kill(r_auxz); } else { ERROR("a point must have a poly or ideal in the first component"); } } export(coord@1); export(coord@2); export(coord@3); export(CHI); int i,j,k; int m,n; // list L@HNE=essdevelop(CHI); list L@HNE=ratdevelop(CHI); export(L@HNE); int n_branches=size(L@HNE); list Li_aux=list(); int N_branches; int N=size(Places); if (sing==1) { list delta2=list(); for (i=1;i<=n_branches;i=i+1) { delta2[i]=invariants(L@HNE[i])[5]; } int dq; } int ext2=res_deg(); list dgs=list(); int ext_0; int check; string sss,olda,newa; if (defined(Q)==0) { int Q; } if (ext1==1) { if (ext2==1) { if (sing==1) { intmat I_mult[n_branches][n_branches]; if (n_branches>1) { for (i=1;i<=n_branches-1;i=i+1) { for (j=i+1;j<=n_branches;j=j+1) { I_mult[i,j]=intersection(L@HNE[i],L@HNE[j]); I_mult[j,i]=I_mult[i,j]; } } } } //////// check=0; //////// if (size(update_CURVE[5])>0) { if (typeof(update_CURVE[5][1])=="list") { check=1; } } if (check==0) { intvec dgs_points(1); ring S(1)=char(basering),(x,y,t),ls; list BRANCHES=list(); list POINTS=list(); list LOC_EQS=list(); list PARAMETRIZATIONS=list(); export(BRANCHES); export(POINTS); export(LOC_EQS); export(PARAMETRIZATIONS); } else { intvec dgs_points(1)=update_CURVE[5][1][2]; def S1=update_CURVE[5][1][1]; execute("ring S(1)="+string(update_CURVE[5][1][1])+";"); fetchall(S1); kill(S1); } N_branches=size(BRANCHES); for (i=1;i<=n_branches;i=i+1) { dgs_points(1)[N_branches+i]=1; POINTS[N_branches+i]=list(); POINTS[N_branches+i][1]=imap(local_aux,coord@1); POINTS[N_branches+i][2]=imap(local_aux,coord@2); POINTS[N_branches+i][3]=imap(local_aux,coord@3); LOC_EQS[N_branches+i]=imap(local_aux,CHI); setring HNEring; Li_aux=L@HNE[i]; setring S(1); BRANCHES=insert(BRANCHES,imap(HNEring,Li_aux),N_branches+i-1); PARAMETRIZATIONS[N_branches+i]=param(BRANCHES[N_branches+i],0); N=N+1; intvec iw=1,N_branches+i; Places[N]=iw; if (sing==1) { dq=delta2[i]; for (j=1;j<=n_branches;j=j+1) { dq=dq+I_mult[i,j]; } Conductor[N]=dq; } if (sing==2) { Conductor[N]=local_conductor(iw[2],S(1)); } kill(iw); PtoPl[i]=N; } setring base_r; update_CURVE[5][1]=list(); update_CURVE[5][1][1]=S(1); update_CURVE[5][1][2]=dgs_points(1); } else { // we start with a rational point but we get non-rational branches // they may have different degrees and then we may need reduce the // field extensions for each one, and finally check if the minpoly // fetchs with S(i) or not // if one of the branches is rational, we may trust that is is written // correctly if (sing==1) { int n_geobrs; int counter_c; list auxgb=list(); list geobrs=list(); for (i=1;i<=n_branches;i=i+1) { auxgb=conj_bs(L@HNE[i],1); dgs[i]=size(auxgb); n_geobrs=n_geobrs+dgs[i]; for (j=1;j<=dgs[i];j=j+1) { counter_c=counter_c+1; geobrs[counter_c]=auxgb[j]; } } intmat I_mult[n_geobrs][n_geobrs]; for (i=1;i=ext_0) { if (typeof(update_CURVE[5][ext_0])=="list") { check=1; } } if (check==0) { if (ext_0>1) { if (ext_0==ext2) { sss=string(minpoly); } else { Q=char(basering)^ext_0; ring auxxx=(Q,a),z,lp; sss=string(minpoly); setring base_r; kill(auxxx); } ring S(ext_0)=(char(basering),a),(x,y,t),ls; execute("minpoly="+sss+";"); } else { ring S(ext_0)=char(basering),(x,y,t),ls; } intvec dgs_points(ext_0); list BRANCHES=list(); list POINTS=list(); list LOC_EQS=list(); list PARAMETRIZATIONS=list(); export(BRANCHES); export(POINTS); export(LOC_EQS); export(PARAMETRIZATIONS); } else { intvec dgs_points(ext_0)=update_CURVE[5][ext_0][2]; def Sext_0=update_CURVE[5][ext_0][1]; setring Sext_0; string SM=string(minpoly); string SR=string(update_CURVE[5][ext_0][1]); execute("ring S("+string(ext_0)+")="+SR+";"); execute("minpoly="+SM+";"); kill(SM,SR); fetchall(Sext_0); kill(Sext_0); } N_branches=size(BRANCHES); dgs_points(ext_0)[N_branches+1]=1; POINTS[N_branches+1]=list(); POINTS[N_branches+1][1]=imap(local_aux,coord@1); POINTS[N_branches+1][2]=imap(local_aux,coord@2); POINTS[N_branches+1][3]=imap(local_aux,coord@3); LOC_EQS[N_branches+1]=imap(local_aux,CHI); // now fetch the branches into the new local ring if (ext_0==1) { setring HNEring; Li_aux=L@HNE[i]; setring S(1); BRANCHES=insert(BRANCHES,imap(HNEring,Li_aux),N_branches); } else { // rationalize branche setring HNEring; newa=subfield(S(ext_0)); m=nrows(L@HNE[i][1]); n=ncols(L@HNE[i][1]); setring S(ext_0); list Laux=list(); poly paux=rationalize(HNEring,"L@HNE["+string(i)+"][4]",newa); matrix Maux[m][n]; for (j=1;j<=m;j=j+1) { for (k=1;k<=n;k=k+1) { Maux[j,k]=rationalize(HNEring,"L@HNE["+string(i)+"][1]["+ string(j)+","+string(k)+"]",newa); } } setring HNEring; intvec Li2=L@HNE[i][2]; int Li3=L@HNE[i][3]; setring S(ext_0); Laux[1]=Maux; Laux[2]=Li2; Laux[3]=Li3; Laux[4]=paux; BRANCHES=insert(BRANCHES,Laux,N_branches); kill(Laux,Maux,paux,Li2,Li3); } PARAMETRIZATIONS[N_branches+1]=param(BRANCHES[N_branches+1],0); N=N+1; intvec iw=ext_0,N_branches+1; Places[N]=iw; if (sing==2) { Conductor[N]=local_conductor(iw[2],S(ext_0)); } kill(iw); PtoPl[i]=N; setring HNEring; update_CURVE[5][ext_0]=list(); update_CURVE[5][ext_0][1]=S(ext_0); update_CURVE[5][ext_0][2]=dgs_points(ext_0); } if (sing==1) { int N_ini=N-n_branches; counter_c=1; for (i=1;i<=n_branches;i=i+1) { dq=delta2[i]; for (j=1;j<=n_geobrs;j=j+1) { dq=dq+I_mult[counter_c,j]; } Conductor[N_ini+i]=dq; counter_c=counter_c+dgs[i]; } } setring base_r; } } else { if (ext1==ext2) { // the degree of the point equals to the degree of all branches // one must just fetch the minpoly's of local_aux, HNEring and S(ext2) if (sing==1) { intmat I_mult[n_branches][n_branches]; if (n_branches>1) { for (i=1;i<=n_branches-1;i=i+1) { for (j=i+1;j<=n_branches;j=j+1) { I_mult[i,j]=intersection(L@HNE[i],L@HNE[j]); I_mult[j,i]=I_mult[i,j]; } } } } //////// check=0; //////// if (size(update_CURVE[5])>=ext2) { if (typeof(update_CURVE[5][ext2])=="list") { check=1; } } if (check==0) { sss=string(minpoly); ring S(ext2)=(char(basering),a),(x,y,t),ls; execute("minpoly="+sss+";"); intvec dgs_points(ext2); list BRANCHES=list(); list POINTS=list(); list LOC_EQS=list(); list PARAMETRIZATIONS=list(); export(BRANCHES); export(POINTS); export(LOC_EQS); export(PARAMETRIZATIONS); } else { intvec dgs_points(ext2)=update_CURVE[5][ext2][2]; def Sext2=update_CURVE[5][ext2][1]; setring Sext2; string SM=string(minpoly); string SR=string(update_CURVE[5][ext2][1]); execute("ring S("+string(ext2)+")="+SR+";"); execute("minpoly="+SM+";"); kill(SM,SR); fetchall(Sext2); kill(Sext2); } N_branches=size(BRANCHES); for (i=1;i<=n_branches;i=i+1) { // fetch all the data into the new local ring olda=subfield(local_aux); dgs_points(ext2)[N_branches+i]=ext1; POINTS[N_branches+i]=list(); POINTS[N_branches+i][1]=number(importdatum(local_aux,"coord@1",olda)); POINTS[N_branches+i][2]=number(importdatum(local_aux,"coord@2",olda)); POINTS[N_branches+i][3]=number(importdatum(local_aux,"coord@3",olda)); LOC_EQS[N_branches+i]=importdatum(local_aux,"CHI",olda); newa=subfield(HNEring); setring HNEring; m=nrows(L@HNE[i][1]); n=ncols(L@HNE[i][1]); setring S(ext2); list Laux=list(); poly paux=importdatum(HNEring,"L@HNE["+string(i)+"][4]",newa); matrix Maux[m][n]; for (j=1;j<=m;j=j+1) { for (k=1;k<=n;k=k+1) { Maux[j,k]=importdatum(HNEring,"L@HNE["+string(i)+"][1]["+ string(j)+","+string(k)+"]",newa); } } setring HNEring; intvec Li2=L@HNE[i][2]; int Li3=L@HNE[i][3]; setring S(ext2); Laux[1]=Maux; Laux[2]=Li2; Laux[3]=Li3; Laux[4]=paux; BRANCHES=insert(BRANCHES,Laux,N_branches+i-1); kill(Laux,Maux,paux,Li2,Li3); PARAMETRIZATIONS[N_branches+i]=param(BRANCHES[N_branches+i],0); N=N+1; intvec iw=ext2,N_branches+i; Places[N]=iw; if (sing==1) { dq=delta2[i]; for (j=1;j<=n_branches;j=j+1) { dq=dq+I_mult[i,j]; } Conductor[N]=dq; } if (sing==2) { Conductor[N]=local_conductor(iw[2],S(ext2)); } kill(iw); PtoPl[i]=N; } setring base_r; update_CURVE[5][ext2]=list(); update_CURVE[5][ext2][1]=S(ext2); update_CURVE[5][ext2][2]=dgs_points(ext2); } else { // this is the most complicated case if (sing==1) { int n_geobrs; int counter_c; list auxgb=list(); list geobrs=list(); for (i=1;i<=n_branches;i=i+1) { auxgb=conj_bs(L@HNE[i],ext1); dgs[i]=size(auxgb); n_geobrs=n_geobrs+dgs[i]; for (j=1;j<=dgs[i];j=j+1) { counter_c=counter_c+1; geobrs[counter_c]=auxgb[j]; } } intmat I_mult[n_geobrs][n_geobrs]; for (i=1;i=ext_0) { if (typeof(update_CURVE[5][ext_0])=="list") { check=1; } } if (check==0) { if (ext_0>ext1) { if (ext_0==ext2) { sss=string(minpoly); } else { Q=char(basering)^ext_0; ring auxxx=(Q,a),z,lp; sss=string(minpoly); setring base_r; kill(auxxx); } } else { setring local_aux; sss=string(minpoly); } ring S(ext_0)=(char(basering),a),(x,y,t),ls; execute("minpoly="+sss+";"); intvec dgs_points(ext_0); list BRANCHES=list(); list POINTS=list(); list LOC_EQS=list(); list PARAMETRIZATIONS=list(); export(BRANCHES); export(POINTS); export(LOC_EQS); export(PARAMETRIZATIONS); } else { intvec dgs_points(ext_0)=update_CURVE[5][ext_0][2]; def Sext_0=update_CURVE[5][ext_0][1]; setring Sext_0; string SM=string(minpoly); string SR=string(update_CURVE[5][ext_0][1]); execute("ring S("+string(ext_0)+")="+SR+";"); execute("minpoly="+SM+";"); kill(SM,SR); fetchall(badring); kill(badring); } N_branches=size(BRANCHES); // now fetch all the data into the new local ring olda=subfield(local_aux); dgs_points(ext_0)[N_branches+1]=ext1; POINTS[N_branches+1]=list(); POINTS[N_branches+1][1]=number(importdatum(local_aux,"coord@1",olda)); POINTS[N_branches+1][2]=number(importdatum(local_aux,"coord@2",olda)); POINTS[N_branches+1][3]=number(importdatum(local_aux,"coord@3",olda)); LOC_EQS[N_branches+1]=importdatum(local_aux,"CHI",olda); setring HNEring; newa=subfield(S(ext_0)); m=nrows(L@HNE[i][1]); n=ncols(L@HNE[i][1]); setring S(ext_0); list Laux=list(); poly paux=rationalize(HNEring,"L@HNE["+string(i)+"][4]",newa); matrix Maux[m][n]; for (j=1;j<=m;j=j+1) { for (k=1;k<=n;k=k+1) { Maux[j,k]=rationalize(HNEring,"L@HNE["+string(i)+"][1]["+ string(j)+","+string(k)+"]",newa); } } setring HNEring; intvec Li2=L@HNE[i][2]; int Li3=L@HNE[i][3]; setring S(ext_0); Laux[1]=Maux; Laux[2]=Li2; Laux[3]=Li3; Laux[4]=paux; BRANCHES=insert(BRANCHES,Laux,N_branches); kill(Laux,Maux,paux,Li2,Li3); PARAMETRIZATIONS[N_branches+1]=param(BRANCHES[N_branches+1],0); N=N+1; intvec iw=ext_0,N_branches+1; Places[N]=iw; if (sing==2) { Conductor[N]=local_conductor(iw[2],S(ext_0)); } kill(iw); PtoPl[i]=N; setring HNEring; update_CURVE[5][ext_0]=list(); update_CURVE[5][ext_0][1]=S(ext_0); update_CURVE[5][ext_0][2]=dgs_points(ext_0); } if (sing==1) { int N_ini=N-n_branches; counter_c=1; for (i=1;i<=n_branches;i=i+1) { dq=delta2[i]; for (j=1;j<=n_geobrs;j=j+1) { dq=dq+I_mult[counter_c,j]; } Conductor[N_ini+i]=dq; counter_c=counter_c+dgs[i]; } } setring base_r; } } update_CURVE[3]=Places; update_CURVE[4]=Conductor; PP[2]=PtoPl; if (Pp[1]==0) { if (Pp[2]==0) { Aff_SPoints[Pp[3]]=PP; } if (Pp[2]==1) { Inf_Points[1][Pp[3]]=PP; } if (Pp[2]==2) { Inf_Points[2][Pp[3]]=PP; } } else { Aff_Points(Pp[2])[Pp[3]]=PP; } update_CURVE[1][1]=base_r; kill(HNEring); kill(local_aux); return(update_CURVE); } /////////////////////////////////////////////////////////////////////////////// static proc local_conductor (int k,SS) { // computes the degree of the local conductor at a place of a plane curve // if the point is non-singular the result will be zero // the computation is carried out with the "Dedekind formula" via // parametrizations int a,b,Cq; def b_ring=basering; setring SS; poly fx=diff(LOC_EQS[k],x); poly fy=diff(LOC_EQS[k],y); int nr=ncols(BRANCHES[k][1]); poly xt=PARAMETRIZATIONS[k][1][1]; poly yt=PARAMETRIZATIONS[k][1][2]; int ordx=PARAMETRIZATIONS[k][2][1]; int ordy=PARAMETRIZATIONS[k][2][2]; map phi_t=basering,xt,yt,1; poly derf; if (fx<>0) { derf=fx; poly tt=diff(yt,t); b=mindeg(tt); if (ordy>-1) { while (b>=ordy) { BRANCHES[k]=extdevelop(BRANCHES[k],2*nr); nr=ncols(BRANCHES[k][1]); PARAMETRIZATIONS[k]=param(BRANCHES[k],0); ordy=PARAMETRIZATIONS[k][2][2]; yt=PARAMETRIZATIONS[k][1][2]; tt=diff(yt,t); b=mindeg(tt); } xt=PARAMETRIZATIONS[k][1][1]; ordx=PARAMETRIZATIONS[k][2][1]; } poly ft=phi_t(derf); } else { derf=fy; poly tt=diff(xt,t); b=mindeg(tt); if (ordx>-1) { while (b>=ordx) { BRANCHES[k]=extdevelop(BRANCHES[k],2*nr); nr=ncols(BRANCHES[k][1]); PARAMETRIZATIONS[k]=param(BRANCHES[k],0); ordx=PARAMETRIZATIONS[k][2][1]; xt=PARAMETRIZATIONS[k][1][1]; tt=diff(xt,t); b=mindeg(tt); } yt=PARAMETRIZATIONS[k][1][2]; ordy=PARAMETRIZATIONS[k][2][2]; } poly ft=phi_t(derf); } a=mindeg(ft); if ( ordx>-1 || ordy>-1 ) { if (ordy==-1) { while (a>ordx) { BRANCHES[k]=extdevelop(BRANCHES[k],2*nr); nr=ncols(BRANCHES[k][1]); PARAMETRIZATIONS[k]=param(BRANCHES[k],0); ordx=PARAMETRIZATIONS[k][2][1]; xt=PARAMETRIZATIONS[k][1][1]; ft=phi_t(derf); a=mindeg(ft); } } else { if (ordx==-1) { while (a>ordy) { BRANCHES[k]=extdevelop(BRANCHES[k],2*nr); nr=ncols(BRANCHES[k][1]); PARAMETRIZATIONS[k]=param(BRANCHES[k],0); ordy=PARAMETRIZATIONS[k][2][2]; yt=PARAMETRIZATIONS[k][1][2]; ft=phi_t(derf); a=mindeg(ft); } } else { int ordf=ordx; if (ordx>ordy) { ordf=ordy; } while (a>ordf) { BRANCHES[k]=extdevelop(BRANCHES[k],2*nr); nr=ncols(BRANCHES[k][1]); PARAMETRIZATIONS[k]=param(BRANCHES[k],0); ordx=PARAMETRIZATIONS[k][2][1]; ordy=PARAMETRIZATIONS[k][2][2]; ordf=ordx; if (ordx>ordy) { ordf=ordy; } xt=PARAMETRIZATIONS[k][1][1]; yt=PARAMETRIZATIONS[k][1][2]; ft=phi_t(derf); a=mindeg(ft); } } } } Cq=a-b; setring b_ring; return(Cq); } /////////////////////////////////////////////////////////////////////////////// static proc max_D (intvec D1,intvec D2) { // computes the maximum of two divisors (intvec) int s1=size(D1); int s2=size(D2); int i; if (s1>s2) { for (i=1;i<=s2;i=i+1) { if (D2[i]>D1[i]) { D1[i]=D2[i]; } } for (i=s2+1;i<=s1;i=i+1) { if (D1[i]<0) { D1[i]=0; } } return(D1); } else { for (i=1;i<=s1;i=i+1) { if (D1[i]>D2[i]) { D2[i]=D1[i]; } } for (i=s1+1;i<=s2;i=i+1) { if (D2[i]<0) { D2[i]=0; } } return(D2); } } /////////////////////////////////////////////////////////////////////////////// static proc deg_D (intvec D,list PP) { // computes the degree of a divisor (intvec or list of integers) int i; int d=0; int s=size(D); for (i=1;i<=s;i=i+1) { d=d+D[i]*PP[i][1]; } return(d); } /////////////////////////////////////////////////////////////////////////////// // ============================================================================ // ******* MAIN PROCEDURES for the "preprocessing" of Brill-Noether ******** // ============================================================================ proc Adj_div (poly f,list #) "USAGE: Adj_div( f [,l] ); f a poly, [l a list] RETURN: list L with the computed data: @format L[1] a list of rings: L[1][1]=aff_r (affine), L[1][2]=Proj_R (projective), L[2] an intvec with 2 entries (degree, genus), L[3] a list of intvec (closed places), L[4] an intvec (conductor), L[5] a list of lists: L[5][d][1] a (local) ring over an extension of degree d, L[5][d][2] an intvec (degrees of base points of places of degree d) @end format NOTE: @code{Adj_div(f);} computes and stores the fundamental data of the plane curve defined by f as needed for AG codes. In the affine ring you can find the following data: @format poly CHI: affine equation of the curve, ideal Aff_SLocus: affine singular locus (std), list Inf_Points: points at infinity Inf_Points[1]: singular points Inf_Points[2]: non-singular points, list Aff_SPoints: affine singular points (if not empty). @end format In the projective ring you can find the projective equation CHI of the curve (poly). In the local rings L[5][d][1] you find: @format list POINTS: base points of the places of degree d, list LOC_EQS: local equations of the curve at the base points, list BRANCHES: Hamburger-Noether developments of the places, list PARAMETRIZATIONS: local parametrizations of the places, @end format Each entry of the list L[3] corresponds to one closed place (i.e., a place and all its conjugates) which is represented by an intvec of size two, the first entry is the degree of the place (in particular, it tells the local ring where to find the data describing one representative of the closed place), and the second one is the position of those data in the lists POINTS, etc., inside this local ring.@* In the intvec L[4] (conductor) the i-th entry corresponds to the i-th entry in the list of places L[3]. With no optional arguments, the conductor is computed by local invariants of the singularities; otherwise it is computed by the Dedekind formula. @* An affine point is represented by a list P where P[1] is std of a prime ideal and P[2] is an intvec containing the position of the places above P in the list of closed places L[3]. @* If the point is at infinity, P[1] is a homogeneous irreducible polynomial in two variables. If @code{printlevel>=0} additional comments are displayed (default: @code{printlevel=0}). KEYWORDS: Hamburger-Noether expansions; adjunction divisor SEE ALSO: closed_points, NSplaces EXAMPLE: example Adj_div; shows an example " { // computes the adjunction divisor and the genus of a (singular) plane curve // as a byproduct, it computes all the singular points with the corresponding // places and the genus of the curve // the adjunction divisor is stored in an intvec // also the non-singular places at infinity are computed // returns a list with all the computed data if (char(basering)==0) { ERROR("Base field not implemented"); } if (npars(basering)>0) { ERROR("Base field not implemented"); } intvec opgt=option(get); option(redSB); def Base_R=basering; list CURVE=curve(f); def aff_r=CURVE[1]; def Proj_R=CURVE[2]; int degX=CURVE[3]; int genusX=(degX-1)*(degX-2); genusX = genusX div 2; intvec iivv=degX,genusX; intvec Conductor; setring aff_r; dbprint(printlevel+1,"Computing affine singular points ... "); list Aff_SPoints=Aff_SL(Aff_SLocus); int s=size(Aff_SPoints); if (s>0) { export(Aff_SPoints); } dbprint(printlevel+1,"Computing all points at infinity ... "); list Inf_Points=inf_P(CHI); export(Inf_Points); list update_CURVE=list(); update_CURVE[1]=list(); update_CURVE[1][1]=aff_r; update_CURVE[1][2]=Proj_R; update_CURVE[2]=iivv; update_CURVE[3]=list(); update_CURVE[4]=Conductor; update_CURVE[5]=list(); int i; intvec pP=0,0,0; if (size(#)==0) { dbprint(printlevel+1,"Computing affine singular places ... "); if (s>0) { for (i=1;i<=s;i=i+1) { pP[3]=i; update_CURVE=place(pP,1,update_CURVE); } } dbprint(printlevel+1,"Computing singular places at infinity ... "); s=size(Inf_Points[1]); if (s>0) { pP[2]=1; for (i=1;i<=s;i=i+1) { pP[3]=i; update_CURVE=place(pP,1,update_CURVE); } } } else { dbprint(printlevel+1,"Computing affine singular places ... "); s=size(Aff_SPoints); if (s>0) { for (i=1;i<=s;i=i+1) { pP[3]=i; update_CURVE=place(pP,2,update_CURVE); } } dbprint(printlevel+1,"Computing singular places at infinity ... "); s=size(Inf_Points[1]); if (s>0) { pP[2]=1; for (i=1;i<=s;i=i+1) { pP[3]=i; update_CURVE=place(pP,2,update_CURVE); } } } dbprint(printlevel+1,"Computing non-singular places at infinity ... "); s=size(Inf_Points[2]); if (s>0) { pP[2]=2; for (i=1;i<=s;i=i+1) { pP[3]=i; update_CURVE=place(pP,0,update_CURVE); } } dbprint(printlevel+1,"Adjunction divisor computed successfully"); list Places=update_CURVE[3]; Conductor=update_CURVE[4]; genusX = genusX - (deg_D(Conductor,Places) div 2); update_CURVE[2][2]=genusX; setring Base_R; dbprint(printlevel+1," "); dbprint(printlevel+2,"The genus of the curve is "+string(genusX)); option(set,opgt); return(update_CURVE); } example { "EXAMPLE:"; echo = 2; int plevel=printlevel; printlevel=-1; ring s=2,(x,y),lp; list C=Adj_div(y9+y8+xy6+x2y3+y2+x3); def aff_R=C[1][1]; // the affine ring setring aff_R; listvar(aff_R); // data in the affine ring CHI; // affine equation of the curve Aff_SLocus; // ideal of the affine singular locus Aff_SPoints[1]; // 1st affine singular point: (1:1:1), no.1 Inf_Points[1]; // singular point(s) at infinity: (1:0:0), no.4 Inf_Points[2]; // list of non-singular points at infinity // pause("press RETURN"); def proj_R=C[1][2]; // the projective ring setring proj_R; CHI; // projective equation of the curve C[2][1]; // degree of the curve C[2][2]; // genus of the curve C[3]; // list of computed places C[4]; // adjunction divisor (all points are singular!) // pause("press RETURN"); // we look at the place(s) of degree 2 by changing to the ring C[5][2][1]; def S(2)=C[5][2][1]; setring S(2); POINTS; // base point(s) of place(s) of degree 2: (1:a:1) LOC_EQS; // local equation(s) PARAMETRIZATIONS; // parametrization(s) and exactness BRANCHES; // Hamburger-Noether development printlevel=plevel; } /////////////////////////////////////////////////////////////////////////////// proc NSplaces (intvec h,list CURVE) "USAGE: NSplaces( h, CURVE ), where h is an intvec and CURVE is a list RETURN: list L with updated data of CURVE after computing all non-singular affine closed places whose degrees are in the intvec h: @* @format in L[1][1]: (affine ring) lists Aff_Points(d) with affine non-singular (closed) points of degree d (if non-empty), in L[3]: the newly computed closed places are added, in L[5]: local rings created/updated to store (repres. of) new places. @end format See @ref{Adj_div} for a description of the entries in L. NOTE: The list_expression should be the output of the procedure Adj_div.@* If @code{printlevel>=0} additional comments are displayed (default: @code{printlevel=0}). SEE ALSO: closed_points, Adj_div EXAMPLE: example NSplaces; shows an example " { // computes all the non-singular affine (closed) places with fixed degrees // creates lists of points and the corresponding places // list CURVE must be the output of the procedure "Adj_div" // notice : if h==0 then nothing is done // warning : non-positive entries are forgotten if (h==0) { return(CURVE); } intvec opgt=option(get); option(redSB); def Base_R=basering; def aff_r=CURVE[1][1]; int M=size(h); list update_CURVE=CURVE; int i,j,l,s; setring aff_r; intvec pP=1,0,0; for (i=1;i<=M;i=i+1) { l=h[i]; if (l>0) { dbprint(printlevel+1,"Computing non-singular affine places of degree " +string(l)+" ... "); if (defined(Aff_Points(l))==0) { list Aff_Points(l)=closed_points_deg(CHI,l,Aff_SLocus); s=size(Aff_Points(l)); if (s>0) { export(Aff_Points(l)); pP[2]=l; for (j=1;j<=s;j=j+1) { pP[3]=j; update_CURVE=place(pP,0,update_CURVE); } } } } } setring Base_R; option(set,opgt); return(update_CURVE); } example { "EXAMPLE:"; echo = 2; int plevel=printlevel; printlevel=-1; ring s=2,(x,y),lp; list C=Adj_div(x3y+y3+x); // The list of computed places: C[3]; // create places up to degree 4 list L=NSplaces(1..4,C); // The list of computed places is now: L[3]; // e.g., affine non-singular points of degree 4 : def aff_r=L[1][1]; setring aff_r; Aff_Points(4); // e.g., base point of the 1st place of degree 4 : def S(4)=L[5][4][1]; setring S(4); POINTS[1]; printlevel=plevel; } /////////////////////////////////////////////////////////////////////////////// // ** SPECIAL PROCEDURES FOR LINEAR ALGEBRA ** static proc Ker (matrix A) { // warning : "lp" ordering is necessary intvec opgt=option(get); option(redSB); matrix M=transpose(syz(A)); option(set,opgt); return(M); } /////////////////////////////////////////////////////////////////////////////// static proc get_NZsol (matrix A) { matrix sol=Ker(A); return(submat(sol,1..1,1..ncols(sol))); } /////////////////////////////////////////////////////////////////////////////// static proc supplement (matrix W,matrix V) "USAGE: supplement(W,V), where W,V are matrices of numbers such that the vector space generated by the rows of W is contained in that generated by the rows of V RETURN: matrix whose rows generate a supplementary vector space of W in V, or a zero row-matrix if == NOTE: W,V must be given with maximal rank " { // W and V represent independent sets of vectors and is assumed to be // contained in // computes matrix S whose rows are l.i. vectors s.t. union is a // basis of // careful : the size of all vectors is assumed to be the same but it is // not checked and neither the linear independence of the vectors is checked // the trivial case W=0 is not covered by this procedure (and thus V<>0) // if = then a zero row-matrix is returned // warning : option(redSB) must be set in advance int n1=nrows(W); int n2=nrows(V); int s=n2-n1; if (s==0) { int n=ncols(W); matrix HH[1][n]; return(HH); } matrix H=transpose(lift(transpose(V),transpose(W))); H=supplem(H); return(H*V); } /////////////////////////////////////////////////////////////////////////////// static proc supplem (matrix M) "USAGE: suplement(M), where M is a matrix of numbers with maximal rank RETURN: matrix whose rows generate a supplementary vector space of in k^n, where k is the base field and n is the number of columns SEE ALSO: supplement NOTE: The rank r is assumed to be 1=0 // exports ideal nFORMS(n) whose generators are ranged with lp order // in Proj_R and returns size(nFORMS(n)) // warning : it is supposed to be called inside Proj_R // if n<=0 then nFORMS(0) is "computed fast" ideal nFORMS(n); int N; if (n>0) { N=(n+1)*(n+2); N=N div 2; N=N+1; int i,j,k; for (i=0;i<=n;i=i+1) { for (j=0;j<=n-i;j=j+1) { k=k+1; nFORMS(n)[N-k]=x^i*y^j*z^(n-i-j); } } export(nFORMS(n)); } else { N=2; nFORMS(0)=1; export(nFORMS(0)); } return(N-1); } /////////////////////////////////////////////////////////////////////////////// static proc nmultiples (int n,int dgX,poly f) { // computes a basis of the space of forms of degree n which are multiple of // CHI // returns a matrix whose rows are the coordinates (related to nFORMS(n)) // of such a basis // warning : it is supposed to be called inside Proj_R // warning : nFORMS(n) is created in the way, together with nFORMS(n-degX) // warning : n must be greater or equal than the degree of the curve if (defined(nFORMS(n))==0) { dbprint(printlevel+1,string(nforms(n))); } int m=n-dgX; if (defined(nFORMS(m))==0) { int k=nforms(m); } else { int k=size(nFORMS(m)); } ideal nmults; int i; for (i=1;i<=k;i=i+1) { nmults[i]=f*nFORMS(m)[i]; } return(transpose(lift(nFORMS(n),nmults))); } /////////////////////////////////////////////////////////////////////////////// static proc interpolating_forms (intvec D,int n,list CURVE) { // computes a vector basis of the space of forms of degree n whose // intersection divisor with the curve is greater or equal than D>=0 // the procedure is supposed to be called inside the ring Proj_R and // assumes that the forms nFORMS(n) are previously computed // the output is a matrix whose rows are the coordinates in nFORMS(n) of // such a basis // remark : the support of D may contain "extra" places def BR=basering; def aff_r=CURVE[1][1]; int N=size(nFORMS(n)); matrix totalM[1][N]; int s=size(D); list Places=CURVE[3]; int NPls=size(Places); int i,j,k,kk,id,ip,RR,ordx,ordy,nr,NR; if (s<=NPls) { for (i=1;i<=s;i=i+1) { if (D[i]>0) { id=Places[i][1]; ip=Places[i][2]; RR=D[i]; def SS=CURVE[5][id][1]; setring SS; poly xt=PARAMETRIZATIONS[ip][1][1]; poly yt=PARAMETRIZATIONS[ip][1][2]; ordx=PARAMETRIZATIONS[ip][2][1]; ordy=PARAMETRIZATIONS[ip][2][2]; nr=ncols(BRANCHES[ip][1]); if ( ordx>-1 || ordy>-1 ) { while ( ( RR>ordx && ordx>-1 ) || ( RR>ordy && ordy>-1 ) ) { BRANCHES[ip]=extdevelop(BRANCHES[ip],2*nr); nr=ncols(BRANCHES[ip][1]); PARAMETRIZATIONS[ip]=param(BRANCHES[ip],0); xt=PARAMETRIZATIONS[ip][1][1]; yt=PARAMETRIZATIONS[ip][1][2]; ordx=PARAMETRIZATIONS[ip][2][1]; ordy=PARAMETRIZATIONS[ip][2][2]; } } if (POINTS[ip][3]==number(1)) { number A=POINTS[ip][1]; number B=POINTS[ip][2]; map Mt=BR,A+xt,B+yt,1; kill(A,B); } else { if (POINTS[ip][2]==number(1)) { number A=POINTS[ip][1]; map Mt=BR,A+xt,1,yt; kill(A); } else { map Mt=BR,1,xt,yt; } } ideal nFORMS(n)=Mt(nFORMS(n)); // rewrite properly the above conditions to obtain the local equations matrix partM[RR][N]; matrix auxMC=coeffs(nFORMS(n),t); NR=nrows(auxMC); if (RR<=NR) { for (j=1;j<=RR;j=j+1) { for (k=1;k<=N;k=k+1) { partM[j,k]=number(auxMC[j,k]); } } } else { for (j=1;j<=NR;j=j+1) { for (k=1;k<=N;k=k+1) { partM[j,k]=number(auxMC[j,k]); } } for (j=NR+1;j<=RR;j=j+1) { for (k=1;k<=N;k=k+1) { partM[j,k]=number(0); } } } matrix localM=partM; matrix conjM=partM; for (j=2;j<=id;j=j+1) { conjM=Frobenius(conjM,1); localM=transpose(concat(transpose(localM),transpose(conjM))); } localM=gauss_row(localM); setring BR; totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM)))); totalM=transpose(compress(transpose(totalM))); setring SS; kill(xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM); setring BR; kill(SS); } } } else { // distinguish between "standard" places and "extra" places for (i=1;i<=NPls;i=i+1) { if (D[i]>0) { id=Places[i][1]; ip=Places[i][2]; RR=D[i]; def SS=CURVE[5][id][1]; setring SS; poly xt=PARAMETRIZATIONS[ip][1][1]; poly yt=PARAMETRIZATIONS[ip][1][2]; ordx=PARAMETRIZATIONS[ip][2][1]; ordy=PARAMETRIZATIONS[ip][2][2]; nr=ncols(BRANCHES[ip][1]); if ( ordx>-1 || ordy>-1 ) { while ( ( RR>ordx && ordx>-1 ) || ( RR>ordy && ordy>-1 ) ) { BRANCHES[ip]=extdevelop(BRANCHES[ip],2*nr); nr=ncols(BRANCHES[ip][1]); PARAMETRIZATIONS[ip]=param(BRANCHES[ip],0); xt=PARAMETRIZATIONS[ip][1][1]; yt=PARAMETRIZATIONS[ip][1][2]; ordx=PARAMETRIZATIONS[ip][2][1]; ordy=PARAMETRIZATIONS[ip][2][2]; } } if (POINTS[ip][3]==number(1)) { number A=POINTS[ip][1]; number B=POINTS[ip][2]; map Mt=BR,A+xt,B+yt,1; kill(A,B); } else { if (POINTS[ip][2]==number(1)) { number A=POINTS[ip][1]; map Mt=BR,A+xt,1,yt; kill(A); } else { map Mt=BR,1,xt,yt; } } ideal nFORMS(n)=Mt(nFORMS(n)); // rewrite properly the above conditions to obtain the local equations matrix partM[RR][N]; matrix auxMC=coeffs(nFORMS(n),t); NR=nrows(auxMC); if (RR<=NR) { for (j=1;j<=RR;j=j+1) { for (k=1;k<=N;k=k+1) { partM[j,k]=number(auxMC[j,k]); } } } else { for (j=1;j<=NR;j=j+1) { for (k=1;k<=N;k=k+1) { partM[j,k]=number(auxMC[j,k]); } } for (j=NR+1;j<=RR;j=j+1) { for (k=1;k<=N;k=k+1) { partM[j,k]=number(0); } } } matrix localM=partM; matrix conjM=partM; for (j=2;j<=id;j=j+1) { conjM=Frobenius(conjM,1); localM=transpose(concat(transpose(localM),transpose(conjM))); } localM=gauss_row(localM); setring BR; totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM)))); totalM=transpose(compress(transpose(totalM))); setring SS; kill(xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM); setring BR; kill(SS); } } k=s-NPls; int l; for (i=1;i<=k;i=i+1) { // in this case D[NPls+i]>0 is assumed to be true during the // Brill-Noether algorithm RR=D[NPls+i]; setring aff_r; def SS=@EXTRA@[2][i]; def extra_dgs=@EXTRA@[3]; setring SS; poly xt=PARAMETRIZATION[1][1]; poly yt=PARAMETRIZATION[1][2]; ordx=PARAMETRIZATION[2][1]; ordy=PARAMETRIZATION[2][2]; nr=ncols(BRANCH[1]); if ( ordx>-1 || ordy>-1 ) { while ( ( RR>ordx && ordx>-1 ) || ( RR>ordy && ordy>-1 ) ) { BRANCH[1]=extdevelop(BRANCH,2*nr); nr=ncols(BRANCH[1]); PARAMETRIZATION=param(BRANCH,0); xt=PARAMETRIZATION[1][1]; yt=PARAMETRIZATION[1][2]; ordx=PARAMETRIZATION[2][1]; ordy=PARAMETRIZATION[2][2]; } } number A=POINT[1]; number B=POINT[2]; map Mt=BR,A+xt,B+yt,1; kill(A,B); ideal nFORMS(n)=Mt(nFORMS(n)); // rewrite properly the above conditions to obtain the local equations matrix partM[RR][N]; matrix auxMC=coeffs(nFORMS(n),t); NR=nrows(auxMC); if (RR<=NR) { for (j=1;j<=RR;j=j+1) { for (kk=1;kk<=N;kk=kk+1) { partM[j,kk]=number(auxMC[j,kk]); } } } else { for (j=1;j<=NR;j=j+1) { for (kk=1;kk<=N;kk=kk+1) { partM[j,kk]=number(auxMC[j,kk]); } } for (j=NR+1;j<=RR;j=j+1) { for (kk=1;kk<=N;kk=kk+1) { partM[j,kk]=number(0); } } } matrix localM=partM; matrix conjM=partM; l=extra_dgs[i]; for (j=2;j<=l;j=j+1) { conjM=Frobenius(conjM,1); localM=transpose(concat(transpose(localM),transpose(conjM))); } localM=gauss_row(localM); setring BR; totalM=transpose(concat(transpose(totalM),transpose(imap(SS,localM)))); totalM=transpose(compress(transpose(totalM))); setring SS; kill(xt,yt,Mt,nFORMS(n),partM,auxMC,conjM,localM); setring BR; kill(SS); } } return(Ker(totalM)); } /////////////////////////////////////////////////////////////////////////////// static proc local_IN (poly h,int m) { // computes the intersection number of h and the curve CHI at a certain place // returns a list with the intersection number and the "leading coefficient" // the procedure must be called inside a local ring, h must be a local // equation with respect to the desired place, and m indicates the // number of place inside that local ring, containing lists // POINT(S)/BRANCH(ES)/PARAMETRIZATION(S) when m=0 an "extra place" is // considered def BR=basering; if (m>0) { int nr=ncols(BRANCHES[m][1]); poly xt=PARAMETRIZATIONS[m][1][1]; poly yt=PARAMETRIZATIONS[m][1][2]; int ordx=PARAMETRIZATIONS[m][2][1]; int ordy=PARAMETRIZATIONS[m][2][2]; map phi=BR,xt,yt,1; poly ht=phi(h); int inum=mindeg(ht); if ( ordx>-1 || ordy>-1 ) { while ( ( inum>ordx && ordx>-1 ) || ( inum>ordy && ordy>-1 ) ) { BRANCHES[m]=extdevelop(BRANCHES[m],2*nr); nr=ncols(BRANCHES[m][1]); PARAMETRIZATIONS[m]=param(BRANCHES[m],0); xt=PARAMETRIZATIONS[m][1][1]; yt=PARAMETRIZATIONS[m][1][2]; ordx=PARAMETRIZATIONS[m][2][1]; ordy=PARAMETRIZATIONS[m][2][2]; ht=phi(h); inum=mindeg(ht); } } } else { int nr=ncols(BRANCH[1]); poly xt=PARAMETRIZATION[1][1]; poly yt=PARAMETRIZATION[1][2]; int ordx=PARAMETRIZATION[2][1]; int ordy=PARAMETRIZATION[2][2]; map phi=basering,xt,yt,1; poly ht=phi(h); int inum=mindeg(ht); if ( ordx>-1 || ordy>-1 ) { while ( ( inum>ordx && ordx>-1 ) || ( inum>ordy && ordy>-1 ) ) { BRANCH=extdevelop(BRANCH,2*nr); nr=ncols(BRANCH[1]); PARAMETRIZATION=param(BRANCH,0); xt=PARAMETRIZATION[1][1]; yt=PARAMETRIZATION[1][2]; ordx=PARAMETRIZATION[2][1]; ordy=PARAMETRIZATION[2][2]; ht=phi(h); inum=mindeg(ht); } } } list answer=list(); answer[1]=inum; number AA=number(coeffs(ht,t)[inum+1,1]); answer[2]=AA; return(answer); } /////////////////////////////////////////////////////////////////////////////// static proc extra_place (ideal P) { // computes the "rational" place which is defined over a (closed) "extra" // point // an "extra" point will be necessarily affine, non-singular and non-rational // creates : a specific local ring to deal with the (unique) place above it // returns : list with the above local ring and the degree of the point/place // warning : the procedure must be called inside the affine ring aff_r def base_r=basering; int ext=deg(P[1]); poly aux=subst(P[2],y,1); ext=ext*deg(aux); // P is assumed to be a std. resp. "(x,y),lp" and thus P[1] depends only // on "y" if (deg(P[1])==1) { // the point is non-rational but the second component needs no field // extension number B=-number(subst(P[1],y,0)); poly aux2=subst(P[2],y,B); // careful : the parameter will be called "a" anyway ring ES=(char(basering),a),(x,y,t),ls; map psi=base_r,a,0; minpoly=number(psi(aux2)); kill(psi); number A=a; number B=imap(base_r,B); } else { if (deg(subst(P[2],y,1))==1) { // the point is non-rational but the needed minpoly is just P[1] // careful : the parameter will be called "a" anyway poly P1=P[1]; poly P2=P[2]; ring ES=(char(basering),a),(x,y,t),ls; map psi=base_r,0,a; minpoly=number(psi(P1)); kill(psi); poly aux1=imap(base_r,P2); poly aux2=subst(aux1,y,a); number A=-number(subst(aux2,x,0)); number B=a; } else { // this is the most complicated case of non-rational point poly P1=P[1]; poly P2=P[2]; int p=char(basering); int Q=p^ext; ring aux_r=(Q,a),(x,y,t),ls; string minpoly_string=string(minpoly); ring ES=(char(basering),a),(x,y,t),ls; execute("minpoly="+minpoly_string+";"); poly P_1=imap(base_r,P1); poly P_2=imap(base_r,P2); ideal factors1=factorize(P_1,1); number B=-number(subst(factors1[1],y,0)); poly P_0=subst(P_2,y,B); ideal factors2=factorize(P_0,1); number A=-number(subst(factors2[1],x,0)); kill(aux_r); } } list POINT=list(); POINT[1]=A; POINT[2]=B; export(POINT); map phi=base_r,x+A,y+B; poly LOC_EQ=phi(CHI); kill(A,B,phi); list L@HNE=essdevelop(LOC_EQ)[1]; export(L@HNE); int m=nrows(L@HNE[1]); int n=ncols(L@HNE[1]); intvec Li2=L@HNE[2]; int Li3=L@HNE[3]; setring ES; string newa=subfield(HNEring); poly paux=importdatum(HNEring,"L@HNE[4]",newa); matrix Maux[m][n]; int i,j; for (i=1;i<=m;i=i+1) { for (j=1;j<=n;j=j+1) { Maux[i,j]=importdatum(HNEring,"L@HNE[1]["+string(i)+","+ string(j)+"]",newa); } } list BRANCH=list(); BRANCH[1]=Maux; BRANCH[2]=Li2; BRANCH[3]=Li3; BRANCH[4]=paux; export(BRANCH); list PARAMETRIZATION=param(BRANCH,0); export(PARAMETRIZATION); kill(HNEring); setring base_r; list answer=list(); answer[1]=ES; answer[2]=ext; kill(ES); return(answer); } /////////////////////////////////////////////////////////////////////////////// static proc intersection_div (poly H,list CURVE) "USAGE: intersection_div(H,CURVE), where H is a homogeneous polynomial in ring Proj_R=p,(x,y,z),lp and CURVE is the list of data for the given curve CREATE: new places which had not been computed in advance if necessary those places are stored each one in a local ring where you find lists POINT,BRANCH,PARAMETRIZATION for the place living in that ring; the degree of the point/place in such a ring is stored in an intvec, and the base points in the remaining list Everything is exported in a list @EXTRA@ inside the ring aff_r=CURVE[1][1] and returned with the updated CURVE RETURN: list with the intersection divisor (intvec) between the underlying curve and H=0, and the list CURVE updated SEE ALSO: Adj_div, NSplaces, closed_points NOTE: The procedure must be called from the ring Proj_R=CURVE[1][2] (projective). If @EXTRA@ already exists, the new places are added to the previous data. " { // computes the intersection divisor of H and the curve CHI // returns a list with (possibly) "extra places" and it must be called // inside Proj_R // in case of extra places, some local rings ES(1) ... ES(m) are created // together with an integer list "extra_dgs" containing the degrees of // those places intvec opgt=option(get); option(redSB); intvec interdiv; def BRing=basering; int Tz1=deg(H); list Places=CURVE[3]; int N=size(Places); def aff_r=CURVE[1][1]; setring aff_r; if (defined(@EXTRA@)==0) { list @EXTRA@=list(); list EP=list(); list ES=list(); list extra_dgs=list(); } else { list EP=@EXTRA@[1]; list ES=@EXTRA@[2]; list extra_dgs=@EXTRA@[3]; } int NN=size(extra_dgs); int counterEPl=0; setring BRing; poly h=subst(H,z,1); int Tz2=deg(h); int Tz3=Tz1-Tz2; int i,j,k,l,m,n,s,np,NP,I_N; if (Tz3==0) { // if this still does not work -> try always with ALL points in // Inf_Points !!!! poly Hinf=subst(H,z,0); setring aff_r; // compute the points at infinity of H and see which of them are in // Inf_Points poly h=imap(BRing,h); poly hinf=imap(BRing,Hinf); ideal pinf=factorize(hinf,1); list TIP=Inf_Points[1]+Inf_Points[2]; s=size(TIP); NP=size(pinf); for (i=1;i<=NP;i=i+1) { for (j=1;j<=s;j=j+1) { if (pinf[i]==TIP[j][1]) { np=size(TIP[j][2]); for (k=1;k<=np;k=k+1) { n=TIP[j][2][k]; l=Places[n][1]; m=Places[n][2]; def SS=CURVE[5][l][1]; setring SS; // local equation h of H if (POINTS[m][2]==number(1)) { number A=POINTS[m][1]; map psi=BRing,x+A,1,y; kill(A); } else { map psi=BRing,1,x,y; } poly h=psi(H); I_N=local_IN(h,m)[1]; interdiv[n]=I_N; kill(h,psi); setring aff_r; kill(SS); } break; } } } kill(hinf,pinf); } else { // H is a multiple of z and hence all the points in Inf_Points intersect // with H setring aff_r; poly h=imap(BRing,h); list TIP=Inf_Points[1]+Inf_Points[2]; s=size(TIP); for (j=1;j<=s;j=j+1) { np=size(TIP[j][2]); for (k=1;k<=np;k=k+1) { n=TIP[j][2][k]; l=Places[n][1]; m=Places[n][2]; def SS=CURVE[5][l][1]; setring SS; // local equation h of H if (POINTS[m][2]==number(1)) { number A=POINTS[m][1]; map psi=BRing,x+A,1,y; kill(A); } else { map psi=BRing,1,x,y; } poly h=psi(H); I_N=local_IN(h,m)[1]; interdiv[n]=I_N; kill(h,psi); setring aff_r; kill(SS); } } } // compute common affine points of H and CHI ideal CAL=h,CHI; CAL=std(CAL); if (CAL<>1) { list TAP=list(); TAP=closed_points(CAL); NP=size(TAP); list auxP=list(); int dP; for (i=1;i<=NP;i=i+1) { if (belongs(TAP[i],Aff_SLocus)==1) { // search the point in the list Aff_SPoints j=isInLP(TAP[i],Aff_SPoints); np=size(Aff_SPoints[j][2]); for (k=1;k<=np;k=k+1) { n=Aff_SPoints[j][2][k]; l=Places[n][1]; m=Places[n][2]; def SS=CURVE[5][l][1]; setring SS; // local equation h of H number A=POINTS[m][1]; number B=POINTS[m][2]; map psi=BRing,x+A,y+B,1; poly h=psi(H); I_N=local_IN(h,m)[1]; interdiv[n]=I_N; kill(A,B,h,psi); setring aff_r; kill(SS); } } else { auxP=list(); auxP[1]=TAP[i]; dP=degree_P(auxP); if (defined(Aff_Points(dP))<>0) { // search the point in the list Aff_Points(dP) j=isInLP(TAP[i],Aff_Points(dP)); n=Aff_Points(dP)[j][2][1]; l=Places[n][1]; m=Places[n][2]; def SS=CURVE[5][l][1]; setring SS; // local equation h of H number A=POINTS[m][1]; number B=POINTS[m][2]; map psi=BRing,x+A,y+B,1; poly h=psi(H); I_N=local_IN(h,m)[1]; interdiv[n]=I_N; kill(A,B,h,psi); setring aff_r; kill(SS); } else { // previously check if it is an existing "extra place" j=isInLP(TAP[i],EP); if (j>0) { def SS=ES[j]; setring SS; // local equation h of H number A=POINT[1]; number B=POINT[2]; map psi=BRing,x+A,y+B,1; poly h=psi(H); I_N=local_IN(h,0)[1]; interdiv[N+j]=I_N; setring aff_r; kill(SS); } else { // then we must create a new "extra place" counterEPl=counterEPl+1; list EXTRA_PLACE=extra_place(TAP[i]); def SS=EXTRA_PLACE[1]; ES[NN+counterEPl]=SS; extra_dgs[NN+counterEPl]=EXTRA_PLACE[2]; EP[NN+counterEPl]=list(); EP[NN+counterEPl][1]=TAP[i]; EP[NN+counterEPl][2]=0; setring SS; // local equation h of H number A=POINT[1]; number B=POINT[2]; map psi=BRing,x+A,y+B,1; poly h=psi(H); I_N=local_IN(h,0)[1]; kill(A,B,h,psi); interdiv[N+NN+counterEPl]=I_N; setring aff_r; kill(SS); } } } } kill(TAP,auxP); } kill(h,CAL,TIP); @EXTRA@[1]=EP; @EXTRA@[2]=ES; @EXTRA@[3]=extra_dgs; kill(EP); list update_CURVE=CURVE; if (size(extra_dgs)>0) { export(@EXTRA@); update_CURVE[1][1]=basering; } else { kill(@EXTRA@); } setring BRing; kill(h); kill(aff_r); list answer=list(); answer[1]=interdiv; answer[2]=update_CURVE; option(set,opgt); return(answer); } /////////////////////////////////////////////////////////////////////////////// static proc local_eq (poly H,SS,int m) { // computes a local equation of poly H in the ring SS related to the place // "m" // list POINT/POINTS is searched depending on wether m=0 or m>0 respectively // warning : the procedure must be called from ring "Proj_R" and returns a // string def BRing=basering; setring SS; if (m>0) { if (POINTS[m][3]==number(1)) { number A=POINTS[m][1]; number B=POINTS[m][2]; map psi=BRing,x+A,y+B,1; } else { if (POINTS[m][2]==number(1)) { number A=POINTS[m][1]; map psi=BRing,x+A,1,y; } else { map psi=BRing,1,x,y; } } } else { number A=POINT[1]; number B=POINT[2]; map psi=BRing,x+A,y+B,1; } poly h=psi(H); string str_h=string(h); setring BRing; return(str_h); } /////////////////////////////////////////////////////////////////////////////// static proc min_wt_rmat (matrix M) { // finds the row of M with minimum non-zero entries, i.e. minimum // "Hamming-weight" int m=nrows(M); int n=ncols(M); int i,j; int Hwt=0; for (j=1;j<=n;j=j+1) { if (M[1,j]<>0) { Hwt=Hwt+1; } } int minHwt=Hwt; int k=1; for (i=2;i<=m;i=i+1) { Hwt=0; for (j=1;j<=n;j=j+1) { if (M[i,j]<>0) { Hwt=Hwt+1; } } if (Hwt0) { kill(@EXTRA@); } setring BRing; update_CURVE[1][1]=aux_RING; kill(aux_RING); matrix B0=supplement(W,V2); if (Hamming_wt(B0)==0) { return(list()); } int ld=nrows(B0); list Bres=list(); ideal HH; poly H; for (j=1;j<=ld;j=j+1) { H=0; for (i=1;i<=N;i=i+1) { H=H+(B0[j,i]*nFORMS(n)[i]); } HH=H,Ho; Bres[j]=simplifyRF(HH); } kill(nFORMS(n)); dbprint(printlevel+1," "); dbprint(printlevel+2,"Vector basis successfully computed "); dbprint(printlevel+1," "); return(Bres); } example { "EXAMPLE:"; echo = 2; int plevel=printlevel; printlevel=-1; ring s=2,(x,y),lp; list C=Adj_div(x3y+y3+x); C=NSplaces(1..4,C); // the first 3 Places in C[3] are of degree 1. // we define the rational divisor G = 4*C[3][1]+4*C[3][3] (of degree 8): intvec G=4,0,4; def R=C[1][2]; setring R; list LG=BrillNoether(G,C); // here is the vector basis of L(G): LG; printlevel=plevel; } /////////////////////////////////////////////////////////////////////////////// // *** procedures for dealing with "RATIONAL FUNCTIONS" over a plane curve *** // a rational function F may be given by (homogeneous) ideal or (affine) poly // (or number) static proc polytoRF (F) { // converts a poly (or number) into a "rational function" of type "ideal" // warning : it must be called inside "R" and poly should be affine ideal RF; RF[1]=homog(F,z); RF[2]=z^(deg(F)); return(RF); } /////////////////////////////////////////////////////////////////////////////// static proc simplifyRF (ideal F) { // simplifies a rational function f/g extracting the gcd(f,g) // maybe add a "restriction" to the curve "CHI" but it is not easy to // programm poly auxp=gcd(F[1],F[2]); return(ideal(division(F,auxp)[1])); } /////////////////////////////////////////////////////////////////////////////// static proc sumRF (F,G) { // sum of two "rational functions" F,G given by either a pair // numerator/denominator or a poly if ( typeof(F)=="ideal" && typeof(G)=="ideal" ) { ideal FG; FG[1]=F[1]*G[2]+F[2]*G[1]; FG[2]=F[2]*G[2]; return(simplifyRF(FG)); } else { if (typeof(F)=="ideal") { ideal GG=polytoRF(G); ideal FG; FG[1]=F[1]*GG[2]+F[2]*GG[1]; FG[2]=F[2]*GG[2]; return(simplifyRF(FG)); } else { if (typeof(G)=="ideal") { ideal FF=polytoRF(F); ideal FG; FG[1]=FF[1]*G[2]+FF[2]*G[1]; FG[2]=FF[2]*G[2]; return(simplifyRF(FG)); } else { return(F+G); } } } } /////////////////////////////////////////////////////////////////////////////// static proc negRF (F) { // returns -F as rational function if (typeof(F)=="ideal") { ideal FF=F; FF[1]=-F[1]; return(FF); } else { return(-F); } } /////////////////////////////////////////////////////////////////////////////// static proc escprodRF (l,F) { // computes l*F as rational function // l should be either a number or a poly of degree zero if (typeof(F)=="ideal") { ideal lF=F; lF[1]=l*F[1]; return(lF); } else { return(l*F); } } /////////////////////////////////////////////////////////////////////////////// // ******** procedures to compute Weierstrass semigroups ******** static proc orderRF (ideal F,SS,int m) "USAGE: orderRF(F,SS,m); F an ideal, SS a ring and m an integer RETURN: list with the order (int) and the leading coefficient (number) NOTE: F represents a rational function, thus the procedure must be called from global ring R or R(d). SS contains the name of a local ring where rational places are stored, and then we take that which is in position m in the corresponding lists of data. The order of F at the place given by SS,m is returned together with the coefficient of minimum degree in the corresponding power series. " { // computes the order of a rational function F at a RATIONAL place given by // a local ring SS and a position "m" inside SS // warning : the procedure must be called from global projective ring "R" or // "R(i)" // returns a list with the order (int) and the "leading coefficient" (number) def BR=basering; poly f=F[1]; string sf=local_eq(f,SS,m); poly g=F[2]; string sg=local_eq(g,SS,m); setring SS; execute("poly ff="+sf+";"); execute("poly gg="+sg+";"); list o1=local_IN(ff,m); list o2=local_IN(gg,m); int oo=o1[1]-o2[1]; number lc=o1[2]/o2[2]; setring BR; number LC=number(imap(SS,lc)); return(list(oo,LC)); } /////////////////////////////////////////////////////////////////////////////// proc Weierstrass (int P,int m,list CURVE) "USAGE: Weierstrass( i, m, CURVE ); i,m integers and CURVE a list RETURN: list WS of two lists: @format WS[1] list of integers (Weierstr. semigroup of the curve at place i up to m) WS[2] list of ideals (the associated rational functions) @end format NOTE: The procedure must be called from the ring CURVE[1][2], where CURVE is the output of the procedure @code{NSplaces}. @* i represents the place CURVE[3][i]. @* Rational functions are represented by numerator/denominator in form of ideals with two homogeneous generators. WARNING: The place must be rational, i.e., necessarily CURVE[3][i][1]=1. @* SEE ALSO: Adj_div, NSplaces, BrillNoether EXAMPLE: example Weierstrass; shows an example " { // computes the Weierstrass semigroup at a RATIONAL place P up to a bound "m" // together with the functions achieving each value up to m, via // Brill-Noether // returns 2 lists : the first consists of all the poles up to m in // increasing order and the second consists of the corresponging rational // functions list Places=CURVE[3]; intvec pl=Places[P]; int dP=pl[1]; int nP=pl[2]; if (dP<>1) { ERROR("the given place is not defined over the prime field"); } if (m<=0) { if (m==0) { list semig=list(); int auxint=0; semig[1]=auxint; list funcs=list(); ideal auxF; auxF[1]=1; auxF[2]=1; funcs[1]=auxF; return(list(semig,funcs)); } else { ERROR("second argument must be non-negative"); } } int auxint=0; ideal auxF; auxF[1]=1; auxF[2]=1; // Brill-Noether algorithm intvec mP; mP[P]=m; list LmP=BrillNoether(mP,CURVE); int lmP=size(LmP); if (lmP==1) { list semig=list(); semig[1]=auxint; list funcs=list(); funcs[1]=auxF; return(list(semig,funcs)); } def SS=CURVE[5][dP][1]; list ordLmP=list(); list polLmP=list(); list sortpol=list(); int maxpol; int i,j,k; for (i=1;i<=lmP-1;i=i+1) { for (j=1;j<=lmP-i+1;j=j+1) { ordLmP[j]=orderRF(LmP[j],SS,nP); polLmP[j]=-ordLmP[j][1]; } sortpol=sort(polLmP); polLmP=sortpol[1]; maxpol=polLmP[lmP-i+1]; LmP=permute_L(LmP,sortpol[2]); ordLmP=permute_L(ordLmP,sortpol[2]); // triangulate LmP for (k=1;k<=lmP-i;k=k+1) { if (polLmP[lmP-i+1-k]==maxpol) { LmP[lmP-i+1-k]=sumRF(LmP[lmP-i+1-k],negRF(escprodRF( ordLmP[lmP-i+1-k][2]/ordLmP[lmP-i+1][2],LmP[lmP-i+1]))); } else { break; } } } polLmP[1]=auxint; LmP[1]=auxF; return(list(polLmP,LmP)); } example { "EXAMPLE:"; echo = 2; int plevel=printlevel; printlevel=-1; ring s=2,(x,y),lp; list C=Adj_div(x3y+y3+x); C=NSplaces(1..4,C); def R=C[1][2]; setring R; // Place C[3][1] has degree 1 (i.e it is rational); list WS=Weierstrass(1,7,C); // the first part of the list is the Weierstrass semigroup up to 7 : WS[1]; // and the second part are the corresponding functions : WS[2]; printlevel=plevel; } /////////////////////////////////////////////////////////////////////////////// // axiliary procedure for permuting a list or intvec proc permute_L (L,P) "USAGE: permute_L( L, P ); L,P either intvecs or lists RETURN: list obtained from L by applying the permutation given by P. NOTE: If P is a list, all entries must be integers. SEE ALSO: sys_code, AGcode_Omega, prepSV EXAMPLE: example permute_L; shows an example " { // permutes the list L according to the permutation P (both intvecs or // lists of integers) int s=size(L); int n=size(P); int i; if (s0) { return(number(0)); } else { ERROR("the function is not well-defined at the given place"); } } } /////////////////////////////////////////////////////////////////////////////// // // ******** procedures for constructing AG codes ******** // /////////////////////////////////////////////////////////////////////////////// static proc gen_mat (list LF,intvec LP,RP) "USAGE: gen_mat(LF,LP,RP); LF list of rational functions, LP intvec of rational places and RP a local ring RETURN: a generator matrix of the evaluation code given by LF and LP SEE ALSO: extcurve KEYWORDS: evaluation codes NOTE: Rational places are searched in local ring RP The procedure must be called from R or R(d) fromlist CURVE after having executed extcurve(d,CURVE) " { // computes a generator matrix (with numbers) of the evaluation code given // by a list of rational functions LF and a list of RATIONAL places LP int m=size(LF); int n=size(LP); matrix GM[m][n]; int i,j; for (i=1;i<=m;i=i+1) { for (j=1;j<=n;j=j+1) { GM[i,j]=evalRF(LF[i],RP,LP[j]); } } return(GM); } /////////////////////////////////////////////////////////////////////////////// proc dual_code (matrix G) "USAGE: dual_code(G); G a matrix of numbers RETURN: a generator matrix of the dual code generated by G NOTE: The input should be a matrix G of numbers. @* The output is also a parity check matrix for the code defined by G KEYWORDS: linear code, dual EXAMPLE: example dual_code; shows an example " { // computes the dual code of C given by a generator matrix G // i.e. computes a parity-check matrix H of C // conversely : computes also G if H is given return(Ker(G)); } example { "EXAMPLE:"; echo = 2; ring s=2,T,lp; // here is the Hamming code of length 7 and dimension 3 matrix G[3][7]=1,0,1,0,1,0,1,0,1,1,0,0,1,1,0,0,0,1,1,1,1; print(G); matrix H=dual_code(G); print(H); } /////////////////////////////////////////////////////////////////////////////// // ====================================================================== // *********** initial test for disjointness *************** // ====================================================================== static proc disj_divs (intvec H,intvec P,list EC) { int s1=size(H); int s2=size(P); list PLACES=EC[3]; def BRing=basering; def auxR=EC[1][5]; setring auxR; int s=res_deg(); setring BRing; kill(auxR); int i,j,k,d,l,N,M; intvec auxIV,iw; for (i=1;i<=s;i=i+1) { if ( (s mod i) == 0 ) { if (typeof(EC[5][i])=="list") { def auxR=EC[5][i][1]; setring auxR; auxIV[i]=size(POINTS); setring BRing; kill(auxR); } else { auxIV[i]=0; } } else { auxIV[i]=0; } } for (i=1;i<=s1;i=i+1) { if (H[i]<>0) { iw=PLACES[i]; d=iw[1]; if ( (s mod d) == 0 ) { l=iw[2]; // check that this place(s) is/are not in sup(D) if (d==1) { for (j=1;j<=s2;j=j+1) { if (l==P[j]) { return(0); } } } else { N=0; for (j=1;j1) { def R=EC[1][2]; setring R; list LG=BrillNoether(G,EC); setring BR; list LG=imap(R,LG); setring R; kill(LG); setring BR; kill(R); } else { list LG=BrillNoether(G,EC); } def RP=EC[1][5]; matrix M=gen_mat(LG,D,RP); kill(LG); return(M); } example { "EXAMPLE:"; echo = 2; int plevel=printlevel; printlevel=-1; ring s=2,(x,y),lp; list HC=Adj_div(x3+y2+y); HC=NSplaces(1..2,HC); HC=extcurve(2,HC); def ER=HC[1][4]; setring ER; intvec G=5; // the rational divisor G = 5*HC[3][1] intvec D=2..9; // D = sum of the rational places no. 2..9 over F_4 // let us construct the corresponding evaluation AG code : matrix C=AGcode_L(G,D,HC); // here is a linear code of type [8,5,>=3] over F_4 print(C); printlevel=plevel; } /////////////////////////////////////////////////////////////////////////////// proc AGcode_Omega (intvec G,intvec D,list EC) "USAGE: AGcode_Omega( G, D, EC ); G,D intvec, EC a list RETURN: a generator matrix for the residual AG code defined by the divisors G and D. NOTE: The procedure must be called within the ring EC[1][4], where EC is the output of @code{extcurve(d)} (or within the ring EC[1][2] if d=1). @* The entry i in the intvec D refers to the i-th rational place in EC[1][5] (i.e., to POINTS[i], etc., see @ref{extcurve}).@* The intvec G represents a rational divisor (see @ref{BrillNoether} for more details).@* The code computes the residues of a vector space basis of @math{\Omega(G-D)} at the rational places given by D. WARNINGS: G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is not checked by the algorithm. G and D should have disjoint supports (checked by the algorithm). SEE ALSO: Adj_div, BrillNoether, extcurve, AGcode_L EXAMPLE: example AGcode_Omega; shows an example " { // returns a generator matrix for the residual AG code given by G and D // G must be a divisor defined over the prime field and D an intvec or // "rational places" // it must be called inside R or R(d) and requires previously "extcurve(d)" return(dual_code(AGcode_L(G,D,EC))); } example { "EXAMPLE:"; echo = 2; int plevel=printlevel; printlevel=-1; ring s=2,(x,y),lp; list HC=Adj_div(x3+y2+y); HC=NSplaces(1..2,HC); HC=extcurve(2,HC); def ER=HC[1][4]; setring ER; intvec G=5; // the rational divisor G = 5*HC[3][1] intvec D=2..9; // D = sum of the rational places no. 2..9 over F_4 // let us construct the corresponding residual AG code : matrix C=AGcode_Omega(G,D,HC); // here is a linear code of type [8,3,>=5] over F_4 print(C); printlevel=plevel; } /////////////////////////////////////////////////////////////////////////////// // ============================================================================ // ******* auxiliary procedure to define AG codes over extensions ******** // ============================================================================ proc extcurve (int d,list CURVE) "USAGE: extcurve( d, CURVE ); d an integer, CURVE a list RETURN: list L which is the update of the list CURVE with additional entries @format L[1][3]: ring (p,a),(x,y),lp (affine), L[1][4]: ring (p,a),(x,y,z),lp (projective), L[1][5]: ring (p,a),(x,y,t),ls (local), L[2][3]: int (the number of rational places), @end format the rings being defined over a field extension of degree d. @* If d<2 then @code{extcurve(d,CURVE);} creates a list L which is the update of the list CURVE with additional entries @format L[1][5]: ring p,(x,y,t),ls, L[2][3]: int (the number of computed places over the base field). @end format In both cases, in the ring L[1][5] lists with the data for all the computed rational places (after a field extension of degree d) are created (see @ref{Adj_div}): @format lists POINTS, LOC_EQS, BRANCHES, PARAMETRIZATIONS. @end format NOTE: The list CURVE should be the output of @code{NSplaces}, and must contain (at least) one place of degree d. @* You actually need all the places with degree dividing d. Otherwise, not all the places are computed, but only part of them. @* This procedure must be executed before constructing AG codes, even if no extension is needed. The ring L[1][4] must be active when constructing codes over the field extension.@* SEE ALSO: closed_points, Adj_div, NSplaces, AGcode_L, AGcode_Omega EXAMPLE: example extcurve; shows an example " { // extends the underlying curve and all its associated objects to a larger // base field in order to evaluate points over such a extension // if d<2 then the only change is that a local ring "RatPl" (which is a // copy of "S(1)") is created in order to store the rational places // where we can do evaluations // otherwise, such a ring contains all places which are rational over the // extension // warning : list Places does not change so that only divisors which are // "rational over the prime field" are allowed; this probably will // change in the future // warning : the places in RatPl are ranged in increasing degree, respecting // the order from list Places and placing the conjugate branches all // together def BR=basering; list ext_CURVE=CURVE; if (d<2) { def SS=CURVE[5][1][1]; ring RatPl=char(basering),(x,y,t),ls; list POINTS=imap(SS,POINTS); list LOC_EQS=imap(SS,LOC_EQS); list BRANCHES=imap(SS,BRANCHES); list PARAMETRIZATIONS=imap(SS,PARAMETRIZATIONS); export(POINTS); export(LOC_EQS); export(BRANCHES); export(PARAMETRIZATIONS); int NrRatPl=size(POINTS); ext_CURVE[2][3]=NrRatPl; setring BR; ext_CURVE[1][5]=RatPl; dbprint(printlevel+1,""); dbprint(printlevel+2,"Total number of rational places : "+string(NrRatPl)); dbprint(printlevel+1,""); kill(RatPl); return(ext_CURVE); } else { // exclude the case when no place of degree d was previously computed/found int dd=size(CURVE[5]); if (dd"list") { ERROR("you did not compute/find any place of degree "+string(d)); } // Define the ring "RatPl" : def S(d)=CURVE[5][d][1]; setring S(d); string smin=string(minpoly); setring BR; ring RatPl=(char(basering),a),(x,y,t),ls; execute("minpoly="+smin+";"); // import data from local ring S(d) list POINTS=imap(S(d),POINTS); list LOC_EQS=imap(S(d),LOC_EQS); list BRANCHES=imap(S(d),BRANCHES); list PARAMETRIZATIONS=imap(S(d),PARAMETRIZATIONS); kill(S(d)); // conjugate data from S(d) int s=size(POINTS); int counter=0; int piv=0; int i,j,k; for (j=1;j<=s;j=j+1) { counter=counter+1; piv=counter; for (k=1;k=2;i=i-1) { if ( (d mod i)==0 ) { if (typeof(CURVE[5][i])=="list") { def S(i)=CURVE[5][i][1]; // embedd S(i) inside basering == RatPl " ==S(d) " olda=subfield(S(i)); setring S(i); s=size(POINTS); setring RatPl; // import data from S(i) for (j=s;j>=1;j=j-1) { counter=0; POINTS=insert(POINTS,list(),0); POINTS[1][1]=number(importdatum(S(i),"POINTS["+string(j) +"][1]",olda)); POINTS[1][2]=number(importdatum(S(i),"POINTS["+string(j) +"][2]",olda)); POINTS[1][3]=number(importdatum(S(i),"POINTS["+string(j) +"][3]",olda)); LOC_EQS=insert(LOC_EQS,importdatum(S(i),"LOC_EQS["+string(j) +"]",olda),0); BRANCHES=insert(BRANCHES,list(),0); setring S(i); m=nrows(BRANCHES[j][1]); n=ncols(BRANCHES[j][1]); iv=BRANCHES[j][2]; kk=BRANCHES[j][3]; poly par@1=subst(PARAMETRIZATIONS[j][1][1],t,x); poly par@2=subst(PARAMETRIZATIONS[j][1][2],t,x); export(par@1); export(par@2); iw=PARAMETRIZATIONS[j][2]; setring RatPl; paux=importdatum(S(i),"BRANCHES["+string(j)+"][4]",olda); matrix Maux[m][n]; for (ii=1;ii<=m;ii=ii+1) { for (jj=1;jj<=n;jj=jj+1) { Maux[ii,jj]=importdatum(S(i),"BRANCHES["+string(j) +"][1]["+string(ii)+","+string(jj)+"]",olda); } } BRANCHES[1][1]=Maux; BRANCHES[1][2]=iv; BRANCHES[1][3]=kk; BRANCHES[1][4]=paux; kill(Maux); PARAMETRIZATIONS=insert(PARAMETRIZATIONS,list(),0); PARAMETRIZATIONS[1][1]=ideal(0); PARAMETRIZATIONS[1][1][1]=importdatum(S(i),"par@1",olda); PARAMETRIZATIONS[1][1][2]=importdatum(S(i),"par@2",olda); PARAMETRIZATIONS[1][1][1]=subst(PARAMETRIZATIONS[1][1][1],x,t); PARAMETRIZATIONS[1][1][2]=subst(PARAMETRIZATIONS[1][1][2],x,t); PARAMETRIZATIONS[1][2]=iw; for (k=1;k0) { w=w+1; } } } return(w); } /////////////////////////////////////////////////////////////////////////////// // Basic Algorithm of Skorobogatov and Vladut for decoding AG codes // warning : the user must choose carefully the parameters for the code and // the decoding since they will never be checked by the procedures proc prepSV (intvec G,intvec D,intvec F,list EC) "USAGE: prepSV( G, D, F, EC ); G,D,F intvecs and EC a list RETURN: list E of size n+3, where n=size(D). All its entries but E[n+3] are matrices: @format E[1]: parity check matrix for the current AG code E[2] ... E[n+2]: matrices used in the procedure decodeSV E[n+3]: intvec with E[n+3][1]: correction capacity @math{epsilon} of the algorithm E[n+3][2]: designed Goppa distance @math{delta} of the current AG code @end format NOTE: Computes the preprocessing for the basic (Skorobogatov-Vladut) decoding algorithm.@* The procedure must be called within the ring EC[1][4], where EC is the output of @code{extcurve(d)} (or in the ring EC[1][2] if d=1) @* The intvec G and F represent rational divisors (see @ref{BrillNoether} for more details).@* The intvec D refers to rational places (see @ref{AGcode_Omega} for more details.). The current AG code is @code{AGcode_Omega(G,D,EC)}.@* If you know the exact minimum distance d and you want to use it in @code{decodeSV} instead of @math{delta}, you can change the value of E[n+3][2] to d before applying decodeSV. If you have a systematic encoding for the current code and want to keep it during the decoding, you must previously permute D (using @code{permute_L(D,P);}), e.g., according to the permutation P=L[3], L being the output of @code{sys_code}. WARNINGS: F must be a divisor with support disjoint from the support of D and with degree @math{epsilon + genus}, where @math{epsilon:=[(deg(G)-3*genus+1)/2]}.@* G should satisfy @math{ 2*genus-2 < deg(G) < size(D) }, which is not checked by the algorithm. G and D should also have disjoint supports (checked by the algorithm). KEYWORDS: SV-decoding algorithm, preprocessing SEE ALSO: extcurve, AGcode_Omega, decodeSV, sys_code, permute_L EXAMPLE: example prepSV; shows an example " { if (disj_divs(F,D,EC)==0) { dbprint(printlevel+3,"? the divisors do not have disjoint supports, empty list returned ?"); return(list()); } list E=list(); list Places=EC[3]; int m=deg_D(G,Places); int genusX=EC[2][2]; int e=(m+1-3*genusX)/2; if (e<1) { dbprint(printlevel+3,"? the correction capacity of the basic algorithm is zero, empty list returned ?"); return(list()); } // deg(F)==e+genusX should be satisfied, and sup(D),sup(F) should be // disjoint !!!! int n=size(D); // 2*genusX-2=5] over F_4: matrix C=AGcode_Omega(G,D,HC); // we can correct 1 error and the genus is 1, thus F must have degree 2 // and support disjoint from that of D; intvec F=2; list SV=prepSV(G,D,F,HC); // now everything is prepared to decode with the basic algorithm; // for example, here is a parity check matrix to compute the syndrome : print(SV[1]); // and here you have the correction capacity of the algorithm : int epsilon=SV[size(D)+3][1]; epsilon; printlevel=plevel; } /////////////////////////////////////////////////////////////////////////////// proc decodeSV (matrix y,list K) "USAGE: decodeSV( y, K ); y a row-matrix and K a list RETURN: a codeword (row-matrix) if possible, resp. the 0-matrix (of size 1) if decoding is impossible. For decoding the basic (Skorobogatov-Vladut) decoding algorithm is applied. NOTE: The list_expression should be the output K of the procedure @code{prepSV}.@* The matrix_expression should be a (1 x n)-matrix, where n = ncols(K[1]).@* The decoding may fail if the number of errors is greater than the correction capacity of the algorithm. KEYWORDS: SV-decoding algorithm SEE ALSO: extcurve, AGcode_Omega, prepSV EXAMPLE: example decodeSV; shows an example " { // decodes y with the "basic decoding algorithm", if possible // requires the preprocessing done by the procedure "prepSV" // the procedure must be called from ring R or R(d) // returns either a codeword (matrix) of none (in case of too many errors) matrix syndr=K[1]*transpose(y); if (Hamming_wt(syndr)==0) { return(y); } matrix Ey=y[1,1]*K[2]; int n=ncols(y); int i; for (i=2;i<=n;i=i+1) { Ey=Ey+y[1,i]*K[i+1]; } matrix Ky=get_NZsol(Ey); if (Hamming_wt(Ky)==0) { dbprint(printlevel+3,"? no error-locator found ?"); dbprint(printlevel+3,"? too many errors occur, 0-matrix returned ?"); matrix answer; return(answer); } int r=nrows(K[n+2]); matrix ErrLoc[1][n]; list Z=list(); list NZ=list(); int j; for (j=1;j<=n;j=j+1) { for (i=1;i<=r;i=i+1) { ErrLoc[1,j]=ErrLoc[1,j]+K[n+2][i,j]*Ky[1,i]; } if (ErrLoc[1,j]==0) { Z=insert(Z,j,size(Z)); } else { NZ=insert(NZ,j,size(NZ)); } } int k=size(NZ); int l=nrows(K[1]); int s=l+k; matrix A[s][n]; matrix b[s][1]; for (i=1;i<=l;i=i+1) { for (j=1;j<=n;j=j+1) { A[i,j]=K[1][i,j]; } b[i,1]=syndr[i,1]; } for (i=1;i<=k;i=i+1) { A[l+i,NZ[i]]=number(1); } intvec opgt=option(get); option(redSB); matrix L=transpose(syz(concat(A,-b))); if (nrows(L)==1) { if (L[1,n+1]<>0) { poly pivote=L[1,n+1]; matrix sol=submat(L,1..1,1..n); if (pivote<>1) { sol=(number(1)/number(pivote))*sol; } // check at least that the number of comitted errors is less than half // the Goppa distance // imposing Hamming_wt(sol)<=K[n+3][1] would be more correct, but maybe // is too strong // on the other hand, if Hamming_wt(sol) is too large the decoding may // not be acceptable if ( Hamming_wt(sol) <= (K[n+3][2]-1)/2 ) { option(set,opgt); return(y-sol); } else { dbprint(printlevel+3,"? non-acceptable decoding ?"); } } else { dbprint(printlevel+3,"? no solution found ?"); } } else { dbprint(printlevel+3,"? non-unique solution ?"); } option(set,opgt); dbprint(printlevel+3,"? too many errors occur, 0-matrix returned ?"); matrix answer; return(answer); } example { "EXAMPLE:"; echo = 2; int plevel=printlevel; printlevel=-1; ring s=2,(x,y),lp; list HC=Adj_div(x3+y2+y); HC=NSplaces(1..2,HC); HC=extcurve(2,HC); def ER=HC[1][4]; setring ER; intvec G=5; // the rational divisor G = 5*HC[3][1] intvec D=2..9; // D = sum of the rational places no. 2..9 over F_4 // construct the corresp. residual AG code of type [8,3,>=5] over F_4: matrix C=AGcode_Omega(G,D,HC); // we can correct 1 error and the genus is 1, thus F must have degree 2 // and support disjoint from that of D intvec F=2; list SV=prepSV(G,D,F,HC); // now we produce 1 error on the zero-codeword : matrix y[1][8]; y[1,3]=a; // and then we decode : print(decodeSV(y,SV)); printlevel=plevel; } /////////////////////////////////////////////////////////////////////////////// proc sys_code (matrix C) "USAGE: sys_code(C); C is a matrix of constants RETURN: list L with: @format L[1] is the generator matrix in standard form of an equivalent code, L[2] is the parity check matrix in standard form of such code, L[3] is an intvec which represents the needed permutation. @end format NOTE: Computes a systematic code which is equivalent to the given one.@* The input should be a matrix of numbers.@* The output has to be interpreted as follows: if the input was the generator matrix of an AG code then one should apply the permutation L[3] to the divisor D of rational points by means of @code{permute_L(D,L[3]);} before continuing to work with the code (for instance, if you want to use the systematic encoding together with a decoding algorithm). KEYWORDS: linear code, systematic SEE ALSO: permute_L, AGcode_Omega, prepSV EXAMPLE: example sys_code; shows an example " { // computes a systematic code which is equivalent to that given by C int i,j,k,l,h,r; int m=nrows(C); int n=ncols(C); int mr=m; matrix A[m][n]=C; poly c,p; list corners=list(); if(m>n) { mr=n; } // first the matrix A will be reduced with elementary operations by rows for(i=1;i<=mr;i=i+1) { if((i+l)>n) { // the process is finished break; } // look for a non-zero element if(A[i,i+l]==0) { h=i; p=0; for (j=i+1;j<=m;j=j+1) { c=A[j,i+l]; if (c!=0) { p=c; h=j; break; } } if (h!=i) { // permutation of rows i and h for (j=1;j<=n;j=j+1) { c=A[i,j]; A[i,j]=A[h,j]; A[h,j]=c; } } if(p==0) { // non-zero element not found in the current column l=l+1; continue; } } // non-zero element was found in "strategic" position corners[i]=i+l; // make zeros below that position for(j=i+1;j<=m;j=j+1) { c=A[j,i+l]/A[i,i+l]; for(k=i+l+1;k<=n;k=k+1) { A[j,k]=A[j,k]-A[i,k]*c; } A[j,i+l]=0; A[j,i]=0; } // the rank is at least r // when the process stops the last r is actually the true rank of A=a r=i; } if (rj) { // interchange columns j and corners[j] for (i=1;i<=m;i=i+1) { c=A[i,j]; A[i,j]=A[i,corners[j]]; A[i,corners[j]]=c; } k=PCols[j]; PCols[j]=PCols[corners[j]]; PCols[corners[j]]=k; } } // convert the diagonal into ones for (i=1;i<=m;i=i+1) { for (j=i+1;j<=n;j=j+1) { A[i,j]=A[i,j]/A[i,i]; } A[i,i]=1; } // complete a block with the identity matrix for (k=1;k1; else use R instead setring ER; Now you have to decide the divisor G and the sequence of rational points D to use for constructing the codes; first notice that the syntax for G and D is different: (a) G is a divisor supported on the closed places of CURVE[3], and you must say just the coefficient of each such a point; for example, G=2,0,-1 means 2 times the place 1 minus 1 times the place 3. (b) D is a sequence of rational points (all different and counted 1 time each one), whose data are read from the lists inside CURVE[1][5] and now you must say just the order how you range the chosen point; for example, D=2,4,1 means that you choose the rational places 1,2,4 and you range them as 2,4,1. (3.1) matrix C=AGcode_L(divisor,places,CURVE); (3.2) AGcode_Omega(divisor,places,CURVE); In the same way as for defining the codes, the auxiliary divisor F must have disjoint support to that of D, and its degree has to be given by a formula (see help prepSV). (3.3) list SV=prepSV(divisor,places,auxiliary_divisor,CURVE); (3.4) decodeSV(word,SV); Special Issues : (I) AG codes with systematic encoding : matrix C=AGcode_Omega(G,D,CURVE); list CODE=sys_code(G); C=CODE[1]; // generator matrix in standard form D=permute_L(D,CODE[3]); // suitable permutation of coordinates list SV=prepSV(G,D,F,CURVE); SV[1]=CODE[2]; // parity-check matrix in standard form (II) Using the true minimum distance d for decoding : matrix C=AGcode_Omega(G,D,CURVE); int n=size(D); list SV=prepSV(G,D,F,CURVE); SV[n+3][2]=d; // then use decodeSV(y,SV); // ============================================================================ // *** Some "macros" with typical examples of curves in Coding Theory **** // ============================================================================ proc Klein () { list KLEIN=Adj_div(x3y+y3+x); KLEIN=NSplaces(1..3,KLEIN); KLEIN=extcurve(3,KLEIN); dbprint(printlevel+1,"Klein quartic over F_8 successfully constructed"); return(KLEIN); } proc Hermite (int m) { int p=char(basering); int r=p^m; list HERMITE=Adj_div(y^r+y-x^(r+1)); HERMITE=NSplaces(1..2*m,HERMITE); HERMITE=extcurve(2*m,HERMITE); dbprint(printlevel+1,"Hermitian curve over F_("+string(r)+"^2) successfully constructed"); return(HERMITE); } proc Suzuki () { list SUZUKI=Adj_div(x10+x3+y8+y); SUZUKI=NSplaces(1..3,SUZUKI); SUZUKI=extcurve(3,SUZUKI); dbprint(printlevel+1,"Suzuki curve over F_8 successfully constructed"); return(SUZUKI); } // **** Other interesting examples : // A hyperelliptic curve with 33 rational points over F_16 list CURVE=Adj_div(x5+y2+y); CURVE=NSplaces(1..4,CURVE); CURVE=extcurve(4,CURVE); // A nice curve with 113 rational points over F_64 list CURVE=Adj_div(y9+y8+xy6+x2y3+y2+x3); intvec ww=1,2,3,6; CURVE=NSplaces(ww,CURVE); CURVE=extcurve(6,CURVE); */