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