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