[1a96a4] | 1 | /////////////////////////////////////////////////////////////////////////////// |
---|
[e8eba7] | 2 | version="$Id: paraplanecurves.lib 14461 2011-12-14 11:00:17Z motsak $"; |
---|
[d2ea299] | 3 | category="Algebraic Geometry"; |
---|
[1a96a4] | 4 | info=" |
---|
[d2ea299] | 5 | LIBRARY: paraplanecurves.lib Rational parametrization of rational plane curves |
---|
[1a96a4] | 6 | |
---|
[ad71e36] | 7 | AUTHORS: J. Boehm, j.boehm at mx.uni-saarland.de @* |
---|
| 8 | W. Decker, decker at mathematik.uni-kl.de> @* |
---|
| 9 | S. Laplagne, slaplagn at dm.uba.ar @* |
---|
| 10 | F. Seelisch, seelisch at mathematik.uni-kl.de |
---|
[1a96a4] | 11 | |
---|
[6c265b] | 12 | OVERVIEW: |
---|
[1a96a4] | 13 | |
---|
[5e8ee4c] | 14 | Suppose C = {f(x,y,z)=0} is a rational plane curve, where f is homogeneous |
---|
| 15 | of degree n with coefficients in Q and absolutely irreducible (these |
---|
[a212fe] | 16 | conditions are checked automatically.) @* |
---|
[1a96a4] | 17 | After a first step, realized by a projective automorphism in the procedure |
---|
| 18 | adjointIdeal, C satisfies: @* |
---|
| 19 | - C does not have singularities at infinity z=0. @* |
---|
| 20 | - C does not contain the point (0:1:0) (that is, the dehomogenization of f |
---|
[a212fe] | 21 | with respect to z is monic as a polynomial in y). @* |
---|
[1a96a4] | 22 | Considering C in the chart z<>0, the algorithm regards x as transcendental |
---|
| 23 | and y as algebraic and computes an integral basis in C(x)[y] of the integral |
---|
| 24 | closure of C[x] in C(x,y) using the normalization algorithm from |
---|
[5e8ee4c] | 25 | @ref{normal_lib}: see @ref{integralbasis_lib}. In a future edition of the |
---|
| 26 | library, also van Hoeij's algorithm for computing the integral basis will |
---|
[a212fe] | 27 | be available. @* |
---|
[e8eba7] | 28 | >From the integral basis, the adjoint ideal is obtained by linear algebra. |
---|
[1a96a4] | 29 | Alternatively, the algorithm starts with a local analysis of the singular |
---|
| 30 | locus of C. Then, for each primary component of the singular locus which |
---|
[6c265b] | 31 | does not correspond to ordinary multiple points or cusps, the integral |
---|
[1a96a4] | 32 | basis algorithm is applied separately. The ordinary multiple points and |
---|
[a212fe] | 33 | cusps, in turn, are addressed by a straightforward direct algorithm. The |
---|
[6c265b] | 34 | adjoint ideal is obtained by intersecting all ideals obtained locally. |
---|
[1a96a4] | 35 | The local variant of the algorithm is used by default. @* |
---|
| 36 | The linear system corresponding to the adjoint ideal maps the curve |
---|
| 37 | birationally to a rational normal curve in P^(n-2). @* |
---|
| 38 | Iterating the anticanonical map, the algorithm projects the rational normal |
---|
| 39 | curve to PP1 for n odd resp. to a conic C2 in PP2 for n even. @* |
---|
| 40 | In case n is even, the algorithm tests whether there is a rational point on |
---|
| 41 | C2 and if so gives a parametrization of C2 which is defined over Q. Otherwise, |
---|
| 42 | the parametrization given is defined over a quadratic field extension of Q. @* |
---|
| 43 | By inverting the birational map of C to PP1 resp. to C2, a parametrization |
---|
| 44 | of C is obtained (defined over Q or the quadratic field extension). |
---|
| 45 | |
---|
[2797244] | 46 | REFERENCES: |
---|
[1a96a4] | 47 | |
---|
[a212fe] | 48 | Janko Boehm: Parametrisierung rationaler Kurven, Diploma Thesis, |
---|
| 49 | http://www.math.uni-sb.de/ag/schreyer/jb/diplom%20janko%20boehm.pdf |
---|
| 50 | |
---|
[1a96a4] | 51 | Theo de Jong: An algorithm for computing the integral closure, |
---|
| 52 | Journal of Symbolic Computation 26 (3) (1998), p. 273-277 |
---|
| 53 | |
---|
| 54 | Gert-Martin Greuel, Santiago Laplagne, Frank Seelisch: Normalization of Rings, |
---|
| 55 | Journal of Symbolic Computation 9 (2010), p. 887-901 |
---|
| 56 | |
---|
| 57 | Mark van Hoeij: An Algorithm for Computing an Integral Basis in an Algebraic |
---|
| 58 | Function Field, Journal of Symbolic Computation 18 (1994), p. 353-363, |
---|
| 59 | http://www.math.fsu.edu/~hoeij/papers/comments/jsc1994.html |
---|
| 60 | |
---|
[6c265b] | 61 | KEYWORDS: |
---|
[1a96a4] | 62 | Curves; Parametrization; Rational curves; Adjoint ideal; Geometric genus |
---|
| 63 | |
---|
| 64 | |
---|
| 65 | PROCEDURES: |
---|
| 66 | |
---|
| 67 | adjointIdeal(poly, [...]); Adjoint ideal of a plane curve |
---|
| 68 | invertBirMap(ideal,ideal); Invert a birational map of algebraic |
---|
| 69 | varieties |
---|
| 70 | paraPlaneCurve(poly, [...]); Compute a rational parametrization of a |
---|
| 71 | rational plane curve |
---|
[d2ea299] | 72 | rncAntiCanonicalMap(ideal); Anticanonical map of a rational normal curve |
---|
[1a96a4] | 73 | rationalPointConic(poly); Finds a point on the conic. This point has |
---|
| 74 | either coefficients in Q or in a quadratic |
---|
| 75 | extension field of Q |
---|
| 76 | mapToRatNormCurve(poly,ideal); Map a plane rational curve to a rational |
---|
| 77 | normal curve (RNC) |
---|
| 78 | rncItProjOdd(ideal); Map a RNC via successive anticanonical maps |
---|
| 79 | to PP1 |
---|
| 80 | rncItProjEven(ideal); Map a RNC via successive anticanonical maps |
---|
| 81 | to a conic in PP2 |
---|
| 82 | paraConic(poly); Compute a rational parametrization of a conic |
---|
| 83 | testParametrization(poly,ring); Checks whether a given curve is parametrized |
---|
[6c265b] | 84 | by a given rational map (defined in the |
---|
[1a96a4] | 85 | given ring) |
---|
[6c265b] | 86 | testPointConic(poly,ring); Checks whether a given point (defined in the |
---|
[1a96a4] | 87 | given ring) lies on the given conic. |
---|
| 88 | "; |
---|
| 89 | |
---|
| 90 | LIB "elim.lib"; |
---|
| 91 | LIB "general.lib"; |
---|
| 92 | LIB "primdec.lib"; |
---|
[d2ea299] | 93 | LIB "absfact.lib"; |
---|
[1a96a4] | 94 | LIB "matrix.lib"; |
---|
| 95 | LIB "random.lib"; |
---|
[d2ea299] | 96 | LIB "homolog.lib"; |
---|
[6c265b] | 97 | LIB "integralbasis.lib"; |
---|
[d2ea299] | 98 | LIB "normal.lib"; |
---|
[1a96a4] | 99 | |
---|
| 100 | |
---|
| 101 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 102 | proc invertBirMap(ideal phi, ideal I) |
---|
| 103 | "USAGE: invertBirMap(phi, I); phi ideal, I ideal |
---|
| 104 | ASSUME: The ideal phi in the basering R represents a birational map of the |
---|
| 105 | variety given by the ideal I in R to its image in projective space |
---|
| 106 | P = PP^(size(phi)-1). |
---|
| 107 | NOTE: The procedure might fail or give a wrong output if phi does |
---|
| 108 | not define a birational map. |
---|
[5e8ee4c] | 109 | RETURN: ring, the coordinate ring of P, with an ideal named J and an ideal |
---|
[a212fe] | 110 | named psi.@* |
---|
| 111 | The ideal J defines the image of phi.@* |
---|
[1a96a4] | 112 | The ideal psi gives the inverse of phi.@* |
---|
| 113 | Note that the entries of psi should be considered as representatives |
---|
| 114 | of classes in the quotient ring R/J.@* |
---|
| 115 | THEORY: We compute the ideal I(G) in R**S of the graph G of phi.@* |
---|
| 116 | The ideal J is given by the intersection of I(G) with S.@* |
---|
| 117 | The map psi is given by a relation mod J of those relations |
---|
| 118 | in I(G) which are linear in the variables of R.@* |
---|
| 119 | KEYWORDS: birational map, image, inverse. |
---|
| 120 | EXAMPLE: example invertBirMap; shows an example |
---|
| 121 | " |
---|
| 122 | { |
---|
| 123 | def Roriginal = basering; |
---|
| 124 | int n = nvars(Roriginal); |
---|
[6c265b] | 125 | int m = size(phi); |
---|
[1a96a4] | 126 | /*phi: P^(n-1) --> P^(m-1)*/ |
---|
| 127 | list rl = ringlist(Roriginal); |
---|
| 128 | int k; |
---|
| 129 | for(k = 1; k <= n; k++) |
---|
| 130 | { |
---|
[6c265b] | 131 | rl[2][k] = "x("+string(k)+")"; |
---|
| 132 | } |
---|
[1a96a4] | 133 | for(k = 1; k <= m; k++) |
---|
| 134 | { |
---|
[6c265b] | 135 | rl[2][k+n] = "y("+string(k)+")"; |
---|
| 136 | } |
---|
[1a96a4] | 137 | rl[3]= list(list("dp",1:(n+m)),list("C",0)); |
---|
| 138 | /*Use Hilbert driven Buchberger*/ |
---|
| 139 | def Rbig0 = ring(rl); |
---|
| 140 | setring Rbig0; |
---|
| 141 | ideal I = fetch(Roriginal,I); |
---|
| 142 | ideal phi = fetch(Roriginal,phi); |
---|
| 143 | ideal mi = maxideal(1); |
---|
| 144 | ideal xv = mi[1..n]; |
---|
| 145 | ideal yv = mi[n+1..n+m]; |
---|
| 146 | matrix HM[2][m] = concat(transpose(yv),transpose(phi)); |
---|
| 147 | ideal graph = sat(I+minor(HM,2),phi)[1]; |
---|
| 148 | graph = sat(graph,xv)[1]; |
---|
| 149 | intvec Hgraph = hilb(graph,1); |
---|
| 150 | setring Roriginal; |
---|
[6c265b] | 151 | rl[3]= list(list("dp",1:n),list("dp",1:m),list("C",0)); |
---|
[1a96a4] | 152 | def Rbig = ring(rl); |
---|
| 153 | setring Rbig; |
---|
| 154 | ideal graph = imap(Rbig0,graph); |
---|
[6c265b] | 155 | graph = std(graph,Hgraph); |
---|
[1a96a4] | 156 | ideal xv = imap(Rbig0,xv); |
---|
| 157 | /*The ideal J defines the image of phi*/ |
---|
| 158 | ideal J = graph; |
---|
| 159 | for(k = 1; k <= n; k++) |
---|
| 160 | { |
---|
[6c265b] | 161 | J = subst(J,xv[k],0); |
---|
| 162 | } |
---|
[1a96a4] | 163 | J = compress(J); |
---|
| 164 | /*now we start inverting phi to psi*/ |
---|
| 165 | matrix relpsi = diff(xv,graph); |
---|
| 166 | for(k = 1; k <= n; k++) |
---|
| 167 | { |
---|
[6c265b] | 168 | relpsi = subst(relpsi,xv[k],0); |
---|
| 169 | } |
---|
[1a96a4] | 170 | relpsi = compress(relpsi); |
---|
| 171 | list rl = ringlist(Rbig); |
---|
| 172 | list rl2 = rl[2]; |
---|
| 173 | rl[2] = list(rl2[n+1..n+m]); |
---|
| 174 | rl[3]= list(list("dp",1:m),list("C",0)); |
---|
| 175 | def Rtarget = ring(rl); |
---|
| 176 | setring Rtarget; |
---|
| 177 | ideal J = imap(Rbig,J); |
---|
| 178 | qring RtargetmodJ = std(J); |
---|
| 179 | matrix relpsi = imap(Rbig,relpsi); |
---|
| 180 | relpsi = syz(transpose(relpsi)); |
---|
| 181 | ideal psi = submat(relpsi,1..nrows(relpsi),1); |
---|
| 182 | setring Rtarget; |
---|
| 183 | ideal psi = imap(RtargetmodJ,psi); |
---|
| 184 | export(J,psi); |
---|
| 185 | int p = printlevel - voice + 3; |
---|
| 186 | dbprint(p,"// 'invertBirMap' created a ring together with two ideals J and psi."); |
---|
| 187 | dbprint(p,"// Supposing you typed, say, def RPn = invertBirMap(phi,I);"); |
---|
| 188 | dbprint(p,"// you may access the ideals by typing"); |
---|
| 189 | dbprint(p,"// setring RPn; J; psi;"); |
---|
| 190 | return(Rtarget); |
---|
| 191 | } |
---|
| 192 | |
---|
| 193 | example |
---|
[75e82e] | 194 | { "EXAMPLE:"; echo=2; |
---|
[1a96a4] | 195 | ring R = 0,(x,y,z),dp; |
---|
| 196 | poly f = y^8-x^3*(z+x)^5; |
---|
| 197 | ideal adj = adjointIdeal(f); |
---|
| 198 | def Rn = invertBirMap(adj,ideal(f)); |
---|
| 199 | setring(Rn); |
---|
| 200 | J; |
---|
| 201 | psi; |
---|
| 202 | } |
---|
| 203 | |
---|
| 204 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 205 | static proc checkAssumptions(poly f) |
---|
| 206 | "USAGE: checkAssumptions(f); f poly |
---|
| 207 | RETURN: 1 if assumptions are satisfied, 0 otherwise.@* |
---|
[6c265b] | 208 | Assumptions checked are: basering is polynomial ring in 3 variables |
---|
[1a96a4] | 209 | with coefficients in Q, f is homogeneous and absolutely irreducible |
---|
| 210 | " |
---|
| 211 | { |
---|
| 212 | def Roriginal = basering; |
---|
| 213 | list rl = ringlist(Roriginal); |
---|
[6c265b] | 214 | rl[3] = list(list("dp",1:3),list("C",0)); |
---|
[1a96a4] | 215 | if(size(rl[1])>1){ERROR("ground field is not Q");} |
---|
| 216 | if(rl[1]!=0){ERROR("ground field is not Q");} |
---|
| 217 | if(nvars(Roriginal)!=3) |
---|
| 218 | { ERROR("not a projective plane curve: wrong number of variables"); } |
---|
| 219 | if(homog(f)==0) |
---|
| 220 | { ERROR("not a projective plane curve: polynomial is not homogeneous"); } |
---|
| 221 | def RP2 = ring(rl); |
---|
[6c265b] | 222 | setring RP2; |
---|
[1a96a4] | 223 | poly f = fetch(Roriginal,f); |
---|
[6c265b] | 224 | if(isIrreducible(f)==0){ERROR("curve is not absolutely irreducible");} |
---|
[1a96a4] | 225 | setring Roriginal; |
---|
| 226 | } |
---|
| 227 | |
---|
| 228 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 229 | proc paraPlaneCurve(poly f, list #) |
---|
[2cdfc8] | 230 | "USAGE: paraPlaneCurve(f [, s]); f poly , s optional string@* |
---|
| 231 | optional string s can be: @* |
---|
| 232 | 'normal': compute integral basis via normalization. @* |
---|
| 233 | 'local': make local analysis of singularities first and apply |
---|
| 234 | normalization separately. @* |
---|
[1a96a4] | 235 | The default is 2. @* |
---|
| 236 | ASSUME: The basering must be a polynomial ring in three variables, say x,y,z, |
---|
| 237 | with coefficients in Q. @* |
---|
| 238 | The polynomial f must be homogeneous and absolutely irreducible. @* |
---|
| 239 | The curve C = {f = 0} must be rational, i.e., have geometric genus 0 |
---|
[a212fe] | 240 | (see @ref{genus}). @* |
---|
| 241 | These conditions will be checked automatically. |
---|
[2cdfc8] | 242 | RETURN: ring with an ideal PARA which contains a rational parametrization of |
---|
| 243 | the rational plane curve given by f; the ground field of the returned |
---|
| 244 | polynomial ring is either Q or some algebraic extension Q(a); PARA |
---|
| 245 | consists of three generators that parametrize the three coordinates |
---|
| 246 | of the rational curve |
---|
[1a96a4] | 247 | THEORY: After a first step, realized by a projective automorphism in the |
---|
| 248 | procedure adjointIdeal, C satisfies: @* |
---|
| 249 | - C does not have singularities at infinity z=0. @* |
---|
| 250 | - C does not contain the point (0:1:0) (that is, the dehomogenization of f |
---|
[a212fe] | 251 | with respect to z is monic as a polynomial in y). @* |
---|
[1a96a4] | 252 | Considering C in the chart z<>0, the algorithm regards x as transcendental |
---|
| 253 | and y as algebraic and computes an integral basis in C(x)[y] of the integral |
---|
[7f92483] | 254 | closure of C[x] in C(x,y) using the normalization algorithm from @ref{normal_lib}: |
---|
| 255 | see @ref{integralbasis_lib}. In a future edition of the library, also van Hoeij's |
---|
[1a96a4] | 256 | algorithm for computing the integral basis will be available. @* |
---|
[e8eba7] | 257 | >From the integral basis, the adjoint ideal is obtained by linear algebra. |
---|
[1a96a4] | 258 | Alternatively, the algorithm starts with a local analysis of the singular |
---|
| 259 | locus of C. Then, for each primary component of the singular locus which |
---|
[6c265b] | 260 | does not correspond to ordinary multiple points or cusps, the integral |
---|
[1a96a4] | 261 | basis algorithm is applied separately. The ordinary multiple points and |
---|
[a212fe] | 262 | cusps, in turn, are addressed by a straightforward direct algorithm. The |
---|
[6c265b] | 263 | adjoint ideal is obtained by intersecting all ideals obtained locally. |
---|
[1a96a4] | 264 | The local variant of the algorithm is used by default. @* |
---|
| 265 | The linear system corresponding to the adjoint ideal maps the curve |
---|
| 266 | birationally to a rational normal curve in P^(n-2). @* |
---|
| 267 | Iterating the anticanonical map, the algorithm projects the rational normal |
---|
| 268 | curve to PP1 for n odd resp. to a conic C2 in PP2 for n even. @* |
---|
| 269 | In case n is even, the algorithm tests whether there is a rational point on C2 |
---|
| 270 | and if so gives a parametrization of C2 which is defined over Q. Otherwise the |
---|
| 271 | parametrization is defined over a quadratic field extension of Q. @* |
---|
| 272 | By inverting the birational map of C to PP1 resp. to C2, a parametrization of |
---|
| 273 | C is obtained (defined over Q or the quadratic field extension). |
---|
| 274 | KEYWORDS: rational curves, rational parametrization of rational curves. |
---|
| 275 | EXAMPLE: example paraPlaneCurve; shows an example |
---|
| 276 | " |
---|
[6c265b] | 277 | { |
---|
[2cdfc8] | 278 | int choice = 2; |
---|
| 279 | if (size(#) != 0) |
---|
| 280 | { |
---|
| 281 | if (typeof(#[1]) == "string") |
---|
| 282 | { |
---|
| 283 | string s = string(#[1]); |
---|
| 284 | if (s == "normal") { choice = 1; } |
---|
| 285 | else |
---|
| 286 | { |
---|
| 287 | if (s == "local") { choice = 2; } |
---|
| 288 | else { ERROR("expected optional argument to be either" |
---|
| 289 | + " 'local' or 'normal'"); } |
---|
| 290 | } |
---|
| 291 | } |
---|
| 292 | else { ERROR("expected optional argument to be a string"); } |
---|
| 293 | } |
---|
[6c265b] | 294 | def Roriginal = basering; |
---|
[1a96a4] | 295 | /*checking assumptions and handling the conic case*/ |
---|
| 296 | checkAssumptions(f); |
---|
| 297 | list rl = ringlist(Roriginal); |
---|
| 298 | rl[2] = list("x","y","z"); |
---|
| 299 | rl[3] = list(list("dp",1:3),list("C",0)); |
---|
| 300 | def RP2 = ring(rl); |
---|
[6c265b] | 301 | setring RP2; |
---|
[1a96a4] | 302 | poly f = fetch(Roriginal,f); |
---|
| 303 | int d = deg(f); |
---|
| 304 | if(d==2) |
---|
| 305 | { |
---|
| 306 | def RP1 = paraConic(f); // (ring, PARACONIC) |
---|
| 307 | setring RP1; |
---|
| 308 | ideal PARA = PARACONIC; |
---|
| 309 | export(PARA); |
---|
| 310 | "// 'paraPlaneCurve' created a ring together with an ideal PARA."; |
---|
| 311 | "// Supposing you typed, say, def RP1 = paraPlaneCurve(f);"; |
---|
| 312 | "// you may access the ideal by typing"; |
---|
| 313 | "// setring RP1; PARA;"; |
---|
| 314 | return(RP1); |
---|
[6c265b] | 315 | } |
---|
[1a96a4] | 316 | int k; |
---|
[6c265b] | 317 | /*the adjoint ideal*/ |
---|
| 318 | ideal AI = adjointIdeal(f,list(choice,"rattestyes/firstchecksdone")); |
---|
[1a96a4] | 319 | /*rattestyes -> causes error message if curve is not rational*/ |
---|
| 320 | /*firstchecksdone -> prevents that check of assumptions will be done again*/ |
---|
| 321 | /*mapping the curve to a rational normal curve V(RNC) in Proj(Rrnc) via AI*/ |
---|
| 322 | def Rrnc = mapToRatNormCurve(f, AI); // ring containing ideal RNC |
---|
| 323 | setring Rrnc; |
---|
| 324 | int m = d-1; // the size of AI |
---|
| 325 | /*the odd dimensional case*/ |
---|
| 326 | if((d mod 2) == 1) |
---|
| 327 | { |
---|
| 328 | /*mapping the rational normal curve to P^1 creating PHI*/ |
---|
| 329 | ideal PHI = rncItProjOdd(RNC); |
---|
| 330 | /*composing the maps AI and PHI*/ |
---|
| 331 | def Rbig = Rrnc + RP2; |
---|
| 332 | setring Rbig; |
---|
| 333 | ideal AI = imap(RP2,AI); |
---|
| 334 | ideal PROJ = imap(Rrnc,PHI); |
---|
| 335 | for(k = 1; k <= m; k++) |
---|
| 336 | { |
---|
[6c265b] | 337 | PROJ = subst(PROJ,var(k),AI[k]); |
---|
[1a96a4] | 338 | } |
---|
| 339 | setring RP2; |
---|
| 340 | ideal If = ideal(fetch(Roriginal,f)); |
---|
| 341 | ideal PROJ = imap(Rbig,PROJ); |
---|
| 342 | /*inverting the composed map to psi*/ |
---|
| 343 | def rp1 = invertBirMap(PROJ,If); // ring containing ideal psi |
---|
| 344 | setring rp1; |
---|
| 345 | list rl1 = ringlist(rp1); |
---|
| 346 | rl1[2] = list("s","t"); |
---|
| 347 | def RP1 = ring(rl1); |
---|
| 348 | setring RP1; |
---|
| 349 | ideal PARA = fetch(rp1,psi); |
---|
| 350 | export(PARA); |
---|
| 351 | "// 'paraPlaneCurve' created a ring together with an ideal PARA."; |
---|
| 352 | "// Supposing you typed, say, def RP1 = paraPlaneCurve(f);"; |
---|
| 353 | "// you may access the ideal by typing"; |
---|
| 354 | "// setring RP1; PARA;"; |
---|
| 355 | return(RP1); |
---|
[6c265b] | 356 | } |
---|
[1a96a4] | 357 | /*the even dimensional case*/ |
---|
| 358 | /*mapping the rational normal curve to a CONIC in P^2* creating PHI*/ |
---|
| 359 | def RP2conic = rncItProjEven(RNC); // exports PHI, returns ring |
---|
| 360 | // containing CONIC |
---|
[6c265b] | 361 | setring RP2conic; |
---|
[1a96a4] | 362 | /*mapping the conic to P^1 via pencil defined by Ipoint*/ |
---|
| 363 | def RP2conicNew = projConic(CONIC); // ring containing ideal Ipoint |
---|
| 364 | // possibly defined over algebraic number field |
---|
| 365 | // variables u,v,w |
---|
| 366 | /*composing the maps AI and PHI and Ipoint in two steps*/ |
---|
| 367 | def Rbig = RP2conicNew + Rrnc; |
---|
| 368 | setring Rbig; |
---|
| 369 | ideal PHI = imap(Rrnc,PHI); |
---|
| 370 | ideal PROJ = imap(RP2conicNew,Ipoint); |
---|
| 371 | for(k = 1; k <= 3; k++) |
---|
| 372 | { |
---|
[6c265b] | 373 | PROJ = subst(PROJ,var(k),PHI[k]); |
---|
| 374 | } |
---|
[1a96a4] | 375 | ideal AI = fetch(RP2,AI); |
---|
| 376 | for(k = 1; k <= m; k++) |
---|
| 377 | { |
---|
[6c265b] | 378 | PROJ = subst(PROJ,var(k+3),AI[k]); |
---|
[1a96a4] | 379 | } |
---|
| 380 | setring RP2conicNew; |
---|
| 381 | ideal If = ideal(fetch(Roriginal,f)); |
---|
| 382 | ideal PROJ = imap(Rbig,PROJ); |
---|
| 383 | /*inverting the composed map to psi*/ |
---|
[6c265b] | 384 | def rp1 = invertBirMap(PROJ,If); // (ring, (J,psi)) |
---|
[1a96a4] | 385 | setring rp1; |
---|
| 386 | list rl1 = ringlist(rp1); |
---|
| 387 | rl1[2] = list("s","t"); |
---|
| 388 | def RP1 = ring(rl1); |
---|
| 389 | setring RP1; |
---|
| 390 | ideal PARA = fetch(rp1,psi); |
---|
| 391 | export(PARA); |
---|
| 392 | "// 'paraPlaneCurve' created a ring together with an ideal PARA."; |
---|
| 393 | "// Supposing you typed, say, def RP1 = paraPlaneCurve(f);"; |
---|
| 394 | "// you may access the ideal by typing"; |
---|
| 395 | "// setring RP1; PARA;"; |
---|
| 396 | return(RP1); |
---|
| 397 | } |
---|
| 398 | |
---|
| 399 | example |
---|
[75e82e] | 400 | { "EXAMPLE:"; echo=2; |
---|
[1a96a4] | 401 | ring R = 0,(x,y,z),dp; |
---|
| 402 | poly f1 = 1/2*x^5+x^2*y*z^2+x^3*y*z+1/2*x*y^2*z^2-2*x*y^3*z+y^5; |
---|
| 403 | def Rp1 = paraPlaneCurve(f1); |
---|
| 404 | setring Rp1; |
---|
| 405 | PARA; |
---|
| 406 | setring R; |
---|
| 407 | poly f2 = x6+3x4y2+3x2y4+y6-4x4z2-34x3yz2-7x2y2z2+12xy3z2+6y4z2; |
---|
| 408 | f2 = f2+36x2z4+36xyz4+9y2z4; |
---|
| 409 | def Rp2 = paraPlaneCurve(f2); |
---|
| 410 | setring Rp2; |
---|
[6c265b] | 411 | PARA; |
---|
[1a96a4] | 412 | } |
---|
| 413 | |
---|
[2ab830] | 414 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 415 | // compute the weighted degree of p; |
---|
| 416 | // this code is an exact copy of the proc in paraplanecurves.lib |
---|
| 417 | // (since we do not want to make it non-static) |
---|
| 418 | static proc w_deg(poly p, intvec v) |
---|
| 419 | { |
---|
| 420 | if(p==0){return(-1);} |
---|
| 421 | int d=0; |
---|
| 422 | while(jet(p,d,v)==0){d++;} |
---|
| 423 | d=(transpose(leadexp(jet(p,d,v)))*v)[1]; |
---|
| 424 | return(d); |
---|
| 425 | } |
---|
| 426 | |
---|
[1a96a4] | 427 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 428 | static proc findCoordChange(poly f, ideal JAC) |
---|
| 429 | "USAGE: findCoordChange(f, JAC); f poly, JAC ideal. |
---|
| 430 | ASSUME: The polynomial f is homogeneous in three variables, JAC is |
---|
| 431 | the Jacobi ideal of f. |
---|
[6c265b] | 432 | RETURN: intvec, say a,b,c. After the coordinate change |
---|
[1a96a4] | 433 | var(3) --> a*var(1)+b*var(2)+c*var(3), the curve {f=0} |
---|
| 434 | has no singularities at infinity {var(3)=0}. |
---|
| 435 | " |
---|
| 436 | { |
---|
| 437 | int h = 2; |
---|
| 438 | int a,b,c; |
---|
| 439 | ideal Jfinfty; |
---|
| 440 | while(h) |
---|
| 441 | { |
---|
| 442 | c = 1; |
---|
| 443 | while(c<=h) |
---|
| 444 | { |
---|
| 445 | b = 0; |
---|
| 446 | while(b<=(h-c)) |
---|
| 447 | { |
---|
| 448 | a = h-b-c; |
---|
| 449 | Jfinfty = JAC,a*var(1)+b*var(2)+c*var(3); |
---|
| 450 | if(dim(std(Jfinfty)) == 0) |
---|
[6c265b] | 451 | { |
---|
[1a96a4] | 452 | return(a,b,c); |
---|
| 453 | } |
---|
| 454 | b = b+1; |
---|
| 455 | } |
---|
| 456 | c = c+1; |
---|
| 457 | } |
---|
| 458 | h = h+1; |
---|
[6c265b] | 459 | } |
---|
[1a96a4] | 460 | } |
---|
| 461 | |
---|
| 462 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 463 | proc adjointIdeal(poly f, list #) |
---|
[a212fe] | 464 | "USAGE: adjointIdeal(f [, choices]); f polynomial in three variables, choices |
---|
[1a96a4] | 465 | optional list consisting of one integer or of one string or of one |
---|
| 466 | integer followed by one string. @* |
---|
[a212fe] | 467 | Optional integer can be: @* |
---|
[1a96a4] | 468 | 1: compute integral basis via normalization. @* |
---|
| 469 | 2: make local analysis of singularities first and apply normalization |
---|
| 470 | separately. @* |
---|
[76bbd3] | 471 | 3: normalization via ideal quotient. @* |
---|
| 472 | 4: normalization via local ideal quotient. @* |
---|
[1a96a4] | 473 | The default is 2. @* |
---|
[a212fe] | 474 | Optional string may contain substrings: @* |
---|
[1a96a4] | 475 | - rattestyes -> causes error message if curve is not rational. @* |
---|
| 476 | - firstchecksdone -> prevents that check of assumptions will be done |
---|
| 477 | more than once. |
---|
| 478 | ASSUME: The basering must be a polynomial ring in three variables, say x,y,z, |
---|
| 479 | with coefficients in Q. @* |
---|
| 480 | The polynomial f must be homogeneous and absolutely irreducible.@* |
---|
| 481 | All these conditions will be checked automatically.@* |
---|
| 482 | RETURN: ideal, the adjoint ideal of the curve defined by f. |
---|
| 483 | THEORY: Considering C in the chart z<>0, the algorithm regards x as |
---|
| 484 | transcendental and y as algebraic and computes an integral basis in C(x)[y] of |
---|
| 485 | the integral closure of C[x] in C(x,y) using the normalization algorithm |
---|
[7f92483] | 486 | from @ref{normal_lib}: see @ref{integralbasis_lib}. In a future edition of the library, |
---|
[3576f6] | 487 | also van Hoeij's algorithm for computing the integral basis will be available.@* |
---|
[e8eba7] | 488 | >From the integral basis, the adjoint ideal is obtained by linear algebra. |
---|
[1a96a4] | 489 | Alternatively, the algorithm starts with a local analysis of the singular |
---|
| 490 | locus of C. Then, for each primary component of the singular locus which |
---|
[6c265b] | 491 | does not correspond to ordinary multiple points or cusps, the integral |
---|
[1a96a4] | 492 | basis algorithm is applied separately. The ordinary multiple points and |
---|
[a212fe] | 493 | cusps, in turn, are addressed by a straightforward direct algorithm. The |
---|
[6c265b] | 494 | adjoint ideal is obtained by intersecting all ideals obtained locally. |
---|
[1a96a4] | 495 | The local variant of the algorithm is used by default. @* |
---|
| 496 | KEYWORDS: integral basis; normalization. |
---|
| 497 | EXAMPLE: example adjointIdeal; shows an example |
---|
| 498 | " |
---|
[6c265b] | 499 | { |
---|
[1a96a4] | 500 | list choices = #; |
---|
| 501 | if(size(#)==0) |
---|
| 502 | { |
---|
| 503 | checkAssumptions(f); |
---|
| 504 | choices = list(2, "rattestno"); |
---|
| 505 | } |
---|
| 506 | if(size(#)==1) |
---|
| 507 | { |
---|
| 508 | if(typeof(choices[1])=="int") |
---|
| 509 | { |
---|
| 510 | checkAssumptions(f); |
---|
| 511 | choices = list(choices[1], "rattestno"); |
---|
| 512 | } |
---|
| 513 | else |
---|
| 514 | { |
---|
| 515 | if(not(find(choices[1], "firstchecksdone"))) |
---|
| 516 | { |
---|
| 517 | checkAssumptions(f); |
---|
| 518 | } |
---|
| 519 | else |
---|
| 520 | { |
---|
| 521 | choices = list(2, choices[1]); |
---|
| 522 | } |
---|
[6c265b] | 523 | } |
---|
| 524 | } |
---|
[1a96a4] | 525 | if(size(#) == 2) |
---|
[6c265b] | 526 | { |
---|
[1a96a4] | 527 | if(not(find(choices[2],"firstchecksdone"))) |
---|
| 528 | { |
---|
| 529 | checkAssumptions(f); |
---|
[6c265b] | 530 | } |
---|
[1a96a4] | 531 | } |
---|
| 532 | ideal JAC = diff(maxideal(1),ideal(f)); |
---|
| 533 | ideal Jfinfty = JAC,var(3); |
---|
| 534 | /*applying a change of coordinates if (f=0) has singularities at infinity*/ |
---|
| 535 | int bb1; |
---|
| 536 | if(dim(std(Jfinfty)) >= 1) |
---|
| 537 | { |
---|
| 538 | bb1 = 1; |
---|
[6c265b] | 539 | int a,b,c = findCoordChange(f,JAC); |
---|
[62c2b0] | 540 | f = subst(f,var(3),var(3)-number(a)/c*var(1)-number(b)/c*var(2)); |
---|
[1a96a4] | 541 | } |
---|
| 542 | /*applying a change of coordinates if the point (0:1:0) lies on the curve*/ |
---|
| 543 | matrix co = coeffs(f,var(2)); |
---|
| 544 | int bb2 = ((size(co)-1) != deg(f)); |
---|
| 545 | if(bb2) |
---|
[6c265b] | 546 | { |
---|
[1a96a4] | 547 | co = coeffs(f,var(1)); |
---|
| 548 | int bb2x = ((size(co)-1) == deg(f)); |
---|
[6c265b] | 549 | if(bb2x) |
---|
[1a96a4] | 550 | { |
---|
| 551 | map perm = basering, var(2), var(1), var(3); |
---|
| 552 | f = perm(f); |
---|
| 553 | } |
---|
| 554 | else |
---|
| 555 | { |
---|
[6c265b] | 556 | f = subst(f,var(1),var(1)+var(2)); |
---|
| 557 | } |
---|
[1a96a4] | 558 | } |
---|
| 559 | co = coeffs(f,var(2)); |
---|
| 560 | f = f/co[size(co),1]; |
---|
| 561 | /*the actual computation*/ |
---|
| 562 | ideal AI = adjointIdealAtWork(f,choices); |
---|
| 563 | /*reversing the changes of coordinates if needed*/ |
---|
| 564 | if(bb2) |
---|
| 565 | { |
---|
[6c265b] | 566 | if(bb2x) |
---|
[1a96a4] | 567 | { |
---|
| 568 | map perm = basering, var(2), var(1), var(3); |
---|
| 569 | AI = mstd(perm(AI))[2]; |
---|
| 570 | } |
---|
| 571 | else |
---|
| 572 | { |
---|
[6c265b] | 573 | AI = mstd(substitute(AI,var(1),var(1)-var(2)))[2]; |
---|
| 574 | } |
---|
| 575 | |
---|
[1a96a4] | 576 | } |
---|
| 577 | if(bb1==1) |
---|
| 578 | { |
---|
[62c2b0] | 579 | AI = mstd(substitute(AI,var(3),var(3)+number(a)/c*var(1)+number(b)/c*var(2)))[2]; |
---|
[1a96a4] | 580 | } |
---|
| 581 | return(AI); |
---|
| 582 | } |
---|
| 583 | |
---|
| 584 | example |
---|
[75e82e] | 585 | { "EXAMPLE:"; echo=2; |
---|
[1a96a4] | 586 | ring R = 0,(x,y,z),dp; |
---|
| 587 | poly f = y^8-x^3*(z+x)^5; |
---|
| 588 | adjointIdeal(f); |
---|
| 589 | } |
---|
| 590 | |
---|
| 591 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 592 | static proc adjointIdealAtWork(poly f, list choices) |
---|
| 593 | "USAGE: adjointIdealAtWork(f, choices); f polynomial in three variables, |
---|
| 594 | choices list consisting of one integer followed by one string. @* |
---|
| 595 | integer can be: @* |
---|
| 596 | 1: compute integral basis via normalization. @* |
---|
| 597 | 2: make local analysis of singularities first and apply normalization |
---|
| 598 | separately. @* |
---|
[76bbd3] | 599 | 3: normalization via ideal quotient. @* |
---|
| 600 | 4: normalization via local ideal quotient. @* |
---|
[6c265b] | 601 | The default is 2. @* |
---|
[1a96a4] | 602 | string may contain substring: @* |
---|
| 603 | - rattestyes -> causes error message if curve is not rational. @* |
---|
| 604 | ASSUME: The basering must be a polynomial ring in three variables, say x,y,z, |
---|
| 605 | with coefficients in Q. @* |
---|
| 606 | The polynomial f must be homogeneous and absolutely irreducible. @* |
---|
| 607 | Its dehomogenization with respect to the third variable must be monic |
---|
| 608 | as a polynomial in the second variable (that is, C does not contain |
---|
| 609 | the point (0:1:0)).@* |
---|
[6c265b] | 610 | The curve C is not allowed to have singularities |
---|
[1a96a4] | 611 | at infinity (z = 0). @* |
---|
| 612 | RETURN: ideal, the adjoint ideal of the curve defined by f. |
---|
| 613 | " |
---|
[6c265b] | 614 | { |
---|
| 615 | def Roriginal = basering; |
---|
[1a96a4] | 616 | list rl = ringlist(Roriginal); |
---|
| 617 | rl[3] = list(list("dp",1:nvars(Roriginal)),list("C",0)); |
---|
[6c265b] | 618 | def RP2 = ring(rl); |
---|
[1a96a4] | 619 | setring RP2; |
---|
| 620 | poly f = imap(Roriginal,f); |
---|
| 621 | poly dhf = subst(f,var(3),1); |
---|
| 622 | int n = deg(f); |
---|
| 623 | if((choices[1]==1) || (choices[1]==3)) // no local analysis of singularities |
---|
| 624 | { |
---|
[76bbd3] | 625 | ideal AI; |
---|
| 626 | if (choices[1]==1) |
---|
| 627 | { AI = adjointIdealIB(f,insert(choices,ideal(0),size(choices))); } |
---|
| 628 | else |
---|
| 629 | { AI = adjointIdealIQ(f,insert(choices,ideal(0),size(choices))); } |
---|
[1a96a4] | 630 | AI = homog(std(AI),var(3)); |
---|
| 631 | AI = sat(AI, maxideal(1))[1]; |
---|
| 632 | AI = minbase(AI); |
---|
| 633 | setring Roriginal; |
---|
| 634 | return(imap(RP2,AI)); |
---|
| 635 | } |
---|
| 636 | list LL = geomGenusLA(f); // local analysis of singularities |
---|
| 637 | int sizeLL2 = size(LL[2]); |
---|
| 638 | int sizeLL3 = size(LL[3]); |
---|
| 639 | list LL3 = LL[3]; |
---|
| 640 | ideal LL4 = LL[4]; |
---|
| 641 | if((LL[1]!=0) && (find(choices[2],"rattestyes"))) |
---|
| 642 | { ERROR("not a rational curve"); } |
---|
| 643 | if((LL[2]==0) && (sizeLL3==0) && (LL4==1)) // smooth case |
---|
| 644 | { |
---|
| 645 | setring Roriginal; |
---|
| 646 | return(ideal(1)); |
---|
| 647 | } |
---|
[ce136a] | 648 | int j,k; |
---|
[1a96a4] | 649 | list rl = ringlist(RP2); |
---|
| 650 | rl[2] = list(var(1), var(2)); |
---|
| 651 | rl[3] = list(list("dp",1:2),list("C",0)); |
---|
| 652 | def Rdummy = ring(rl); |
---|
| 653 | ideal B; |
---|
| 654 | if(sizeLL3==0){B = 1;} // no ordinary multiple points |
---|
| 655 | // (except possibly nodes) |
---|
| 656 | else // there are ordinary multiple points |
---|
| 657 | // (other than nodes) |
---|
| 658 | { |
---|
[ce136a] | 659 | setring Rdummy; |
---|
[1a96a4] | 660 | list OMP = imap(RP2,LL3); |
---|
| 661 | int ub; |
---|
| 662 | for(k=1;k<=size(OMP);k++) |
---|
| 663 | { |
---|
| 664 | if(OMP[k][1]>ub) |
---|
| 665 | { |
---|
| 666 | ub = OMP[k][1]; |
---|
| 667 | } |
---|
| 668 | } |
---|
| 669 | int lb = ub; |
---|
| 670 | for(k=1;k<=size(OMP);k++) |
---|
| 671 | { |
---|
| 672 | if(OMP[k][1]<lb) |
---|
| 673 | { |
---|
| 674 | lb = OMP[k][1]; |
---|
| 675 | } |
---|
| 676 | } |
---|
| 677 | for(k=lb;k<=ub;k++) |
---|
| 678 | { |
---|
| 679 | ideal A(k) = 1; |
---|
| 680 | } |
---|
| 681 | for(k=1;k<=size(OMP);k++) |
---|
| 682 | { |
---|
| 683 | A(OMP[k][1]) = intersect(A(OMP[k][1]), OMP[k][2]); |
---|
| 684 | } |
---|
[6c265b] | 685 | int i = ub; |
---|
[1a96a4] | 686 | setring RP2; |
---|
[ce136a] | 687 | for(k=lb;k<=ub;k++) |
---|
| 688 | { |
---|
| 689 | ideal A(k) = homog(std(fetch(Rdummy,A(k))),var(3)); |
---|
| 690 | } |
---|
[1a96a4] | 691 | B = maxideal(n-i); |
---|
[6c265b] | 692 | ideal A; |
---|
[1a96a4] | 693 | while(i>=lb) |
---|
[6c265b] | 694 | { |
---|
[1a96a4] | 695 | A = A(i)**(i-1); |
---|
[ce136a] | 696 | j=1; |
---|
| 697 | while(j<=ncols(A)) |
---|
| 698 | { |
---|
| 699 | if(deg(A[j]>(n-2))) |
---|
| 700 | { |
---|
| 701 | A = sat(A, maxideal(1))[1]; |
---|
| 702 | break; |
---|
| 703 | } |
---|
| 704 | j = j+1; |
---|
[5e8ee4c] | 705 | } |
---|
[1a96a4] | 706 | B = intersect(B,A); |
---|
| 707 | i = i-1; |
---|
| 708 | } |
---|
| 709 | } //end else |
---|
[5e8ee4c] | 710 | B = intersect(B,homog(std(LL4),var(3))); // add nodes and cusps |
---|
[1a96a4] | 711 | if(sizeLL2==0) // ordinary multiple points plus cusps only |
---|
| 712 | { |
---|
[ce136a] | 713 | ideal AI = sat(B, maxideal(1))[1]; |
---|
[1a96a4] | 714 | AI = minbase(AI); |
---|
| 715 | setring Roriginal; |
---|
| 716 | return(imap(RP2,AI)); |
---|
| 717 | } |
---|
[ce136a] | 718 | setring Rdummy; |
---|
[1a96a4] | 719 | poly f = imap(RP2,dhf); |
---|
| 720 | ideal SL = jacob(f),f; |
---|
[ce136a] | 721 | SL = sat(SL, fetch(RP2,LL4))[1]; |
---|
[1a96a4] | 722 | if(sizeLL3!=0) |
---|
| 723 | { |
---|
| 724 | for(k=lb;k<=ub;k++) |
---|
| 725 | { |
---|
| 726 | SL = sat(SL, A(k))[1]; |
---|
| 727 | } |
---|
| 728 | } |
---|
[6c265b] | 729 | list PD = primdecGTZ(SL); // could be made faster -- see minAssGTZ |
---|
[1a96a4] | 730 | // in deltaLocMod -- only one PD needed |
---|
| 731 | int pd = size(PD); |
---|
| 732 | setring RP2; |
---|
[6c265b] | 733 | list PD = imap(Rdummy,PD); |
---|
[1a96a4] | 734 | ideal AI = 1; |
---|
| 735 | for(k=1;k<=pd;k++) |
---|
[6c265b] | 736 | { |
---|
[76bbd3] | 737 | if (choices[1]==2) |
---|
| 738 | { AI = intersect(AI,adjointIdealIB(f,insert(choices,PD[k][1], |
---|
| 739 | size(choices)))); } |
---|
| 740 | else |
---|
| 741 | { AI = intersect(AI,adjointIdealIQ(f,insert(choices,PD[k][1], |
---|
| 742 | size(choices)))); } |
---|
[1a96a4] | 743 | } |
---|
| 744 | AI = homog(std(AI),var(3)); |
---|
[ce136a] | 745 | AI = intersect(AI,B); |
---|
[1a96a4] | 746 | AI = sat(AI, maxideal(1))[1]; |
---|
| 747 | AI = minbase(AI); |
---|
| 748 | setring Roriginal; |
---|
| 749 | return(imap(RP2,AI)); |
---|
| 750 | } |
---|
| 751 | |
---|
| 752 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 753 | static proc adjointIdealIB(poly f, list choices) |
---|
[6c265b] | 754 | "USAGE: adjointIdealIB(f, choices); f polynomial in three variables, choices |
---|
[1a96a4] | 755 | list consisting of one integer followed by one string followed by one |
---|
| 756 | ideal. @* |
---|
| 757 | integer can be: @* |
---|
| 758 | 1, 2 : compute integral basis via normalization @* |
---|
| 759 | The default is 2. @* |
---|
| 760 | string may contain substring: @* |
---|
| 761 | - rattestyes -> causes error message if curve is not rational. @* |
---|
| 762 | ideal serves as input for @ref integralBasis. |
---|
| 763 | ASSUME: The basering must be a polynomial ring in three variables, say x,y,z, |
---|
| 764 | with coefficients in Q. @* |
---|
| 765 | The polynomial f must be homogeneous and absolutely irreducible.@* |
---|
[6c265b] | 766 | Its dehomogenization with respect to the third variable must be monic |
---|
| 767 | as a polynomial in the second variable (that is, C does not contain |
---|
[1a96a4] | 768 | the point (0:1:0)).@* |
---|
[6c265b] | 769 | The curve C is not allowed to have singularities |
---|
[1a96a4] | 770 | at infinity (z = 0). @* |
---|
| 771 | RETURN: ideal containing the adjoint ideal of the curve defined by f. @* |
---|
| 772 | " |
---|
[6c265b] | 773 | { |
---|
[1a96a4] | 774 | poly dhf = subst(f,var(3),1); |
---|
| 775 | def Roriginal = basering; |
---|
| 776 | list rl = ringlist(Roriginal); |
---|
| 777 | rl[2] = list(var(1), var(2)); |
---|
| 778 | rl[3] = list(list("dp",1:2),list("C",0)); |
---|
| 779 | def Rdummy = ring(rl); |
---|
| 780 | setring Rdummy; |
---|
| 781 | poly f = imap(Roriginal,dhf); |
---|
| 782 | poly d2f = diff(f,var(2)); |
---|
| 783 | list DATA = imap(Roriginal,choices); |
---|
| 784 | /* Creating rings for later use */ |
---|
| 785 | list rl = ringlist(Rdummy); |
---|
| 786 | rl[2] = list(var(2), var(1)); |
---|
| 787 | rl[3] = list(list("lp",1:2),list("C",0)); |
---|
| 788 | def Rred = ring(rl); // make var(2) > var(1) |
---|
| 789 | rl = ringlist(Rdummy); |
---|
[6c265b] | 790 | rl[1] = list(0,list(var(1)),list(list("dp",1)),ideal(0)); |
---|
[1a96a4] | 791 | rl[2] = list(var(2)); |
---|
| 792 | rl[3] = list(list("dp",1),list("C",0)); |
---|
| 793 | def QF = ring(rl); // make var(1) transcendental |
---|
| 794 | list LIntB; |
---|
[76bbd3] | 795 | if(DATA[1] <= 4) // use normalization algorithm |
---|
[1a96a4] | 796 | { |
---|
| 797 | LIntB = integralBasis(f, 2, list(list("inputC", DATA[3]),"isIrred")); |
---|
[6c265b] | 798 | } |
---|
[1a96a4] | 799 | else // use van Hoeij's algorithm |
---|
[6c265b] | 800 | { |
---|
[1a96a4] | 801 | LIntB = integralBasisVH(f,DATA[3],2); // van Hoeij in future version |
---|
[76bbd3] | 802 | // used when DATA[1] = 5 |
---|
[1a96a4] | 803 | } |
---|
| 804 | if(find(DATA[2],"rattestyes") && (DATA[3]==0)) |
---|
| 805 | { |
---|
| 806 | setring Roriginal; |
---|
| 807 | int gg = geomGenusIB(f,imap(Rdummy, LIntB)); |
---|
[6c265b] | 808 | if(gg!=0){ERROR("not a rational curve");} |
---|
[1a96a4] | 809 | setring Rdummy; |
---|
| 810 | } |
---|
| 811 | int i,j,k,l; |
---|
| 812 | ideal IB = LIntB[1]; |
---|
| 813 | poly d = LIntB[2]; |
---|
| 814 | int sL=size(IB); |
---|
| 815 | setring Rred; |
---|
| 816 | ideal IB = imap(Rdummy,IB); |
---|
| 817 | ideal fred = std(imap(Rdummy,f)); |
---|
| 818 | IB = reduce(IB,fred); |
---|
| 819 | matrix M = coeffs(IB,var(1)); |
---|
| 820 | setring QF; |
---|
| 821 | matrix M = imap(Rred,M); |
---|
| 822 | poly d = imap(Rdummy,d); |
---|
| 823 | M=1/d*M; |
---|
[6c265b] | 824 | list LUM = ludecomp(M); |
---|
| 825 | list LS; |
---|
[1a96a4] | 826 | matrix dummyvector[sL][1]; |
---|
| 827 | matrix Gij[sL][sL]; |
---|
| 828 | matrix Tr[sL][sL]; |
---|
| 829 | setring Rred; |
---|
| 830 | poly eiej; |
---|
| 831 | list Iij, empty; |
---|
| 832 | matrix Gij[sL][sL]; |
---|
| 833 | for(i = 1; i <= sL; i++) |
---|
| 834 | { |
---|
| 835 | for(j = i; j <= sL; j++) |
---|
[6c265b] | 836 | { |
---|
[1a96a4] | 837 | setring Rred; |
---|
[6c265b] | 838 | Gij = 0; |
---|
| 839 | eiej = IB[i]*IB[j]; |
---|
[1a96a4] | 840 | Iij=empty; |
---|
| 841 | for(k = 1; k <= sL; k++) |
---|
| 842 | { |
---|
| 843 | Iij[k] = reduce(eiej*IB[k],fred); |
---|
[6c265b] | 844 | } |
---|
| 845 | Gij = coeffs(ideal(Iij[1..sL]),var(1)); |
---|
[1a96a4] | 846 | setring QF; |
---|
[6c265b] | 847 | Gij = imap (Rred, Gij); |
---|
| 848 | for(k = 1; k <= sL; k++) |
---|
[1a96a4] | 849 | { |
---|
| 850 | dummyvector = Gij[1..sL,k]; |
---|
[6c265b] | 851 | LS = lusolve(LUM[1], LUM[2], LUM[3], dummyvector); |
---|
[1a96a4] | 852 | Tr[i,j] = Tr[i,j] + 1/d^3*LS[2][k,1]; |
---|
[6c265b] | 853 | } |
---|
[1a96a4] | 854 | } |
---|
| 855 | } |
---|
| 856 | for(i = 1; i <= sL; i++) |
---|
| 857 | { |
---|
| 858 | for(j = 1; j < i; j++) |
---|
[6c265b] | 859 | { |
---|
[1a96a4] | 860 | Tr[i,j] = Tr[j,i]; |
---|
[6c265b] | 861 | } |
---|
| 862 | } |
---|
| 863 | LUM = ludecomp(Tr); |
---|
[1a96a4] | 864 | setring Rred; |
---|
| 865 | poly d2f = imap(Rdummy,d2f); |
---|
| 866 | IB = d2f*IB; |
---|
| 867 | IB = reduce(IB,fred); |
---|
| 868 | setring QF; |
---|
| 869 | matrix IB = transpose(matrix(imap(Rred,IB))); |
---|
| 870 | IB = 1/d*IB; |
---|
| 871 | LS = lusolve(LUM[1], LUM[2], LUM[3], IB); |
---|
| 872 | ideal LL = ideal(LS[2]); |
---|
| 873 | setring Roriginal; |
---|
| 874 | ideal AI = imap(QF,LL); |
---|
| 875 | return(AI); |
---|
| 876 | } |
---|
| 877 | |
---|
[76bbd3] | 878 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 879 | static proc adjointIdealIQ(poly f, list choices) |
---|
| 880 | { |
---|
| 881 | poly dhf = subst(f,var(3),1); |
---|
| 882 | def Roriginal = basering; |
---|
| 883 | list rl = ringlist(Roriginal); |
---|
| 884 | rl[2] = list(var(1), var(2)); |
---|
| 885 | rl[3] = list(list("dp",1:2),list("C",0)); |
---|
| 886 | def Rdummy = ring(rl); |
---|
| 887 | setring Rdummy; |
---|
| 888 | list DATA = imap(Roriginal,choices); |
---|
| 889 | poly f = imap(Roriginal,dhf); |
---|
| 890 | list LIntB; |
---|
| 891 | if(DATA[1] <= 4) // use normalization algorithm |
---|
| 892 | { |
---|
| 893 | LIntB = integralBasis(f, 2, list(list("inputC", DATA[3]),"isIrred")); |
---|
| 894 | } |
---|
| 895 | else // use van Hoeij's algorithm |
---|
| 896 | { |
---|
| 897 | LIntB = integralBasisVH(f,DATA[3],2); // van Hoeij in future version |
---|
| 898 | // used when DATA[1] = 5 |
---|
| 899 | } |
---|
| 900 | if(find(DATA[2],"rattestyes") && (DATA[3]==0)) |
---|
| 901 | { |
---|
| 902 | setring Roriginal; |
---|
| 903 | int gg = geomGenusIB(f,imap(Rdummy, LIntB)); |
---|
| 904 | if(gg!=0){ERROR("not a rational curve");} |
---|
| 905 | setring Rdummy; |
---|
| 906 | } |
---|
| 907 | int i,j,k,l; |
---|
| 908 | ideal IB = LIntB[1]; |
---|
| 909 | poly d = LIntB[2]; |
---|
| 910 | ideal fd = f, d; |
---|
| 911 | ideal IBf = IB, f; |
---|
| 912 | ideal AI = quotient(fd, IBf); |
---|
| 913 | //"#### IB:"; IB; |
---|
| 914 | //"#### d:", d; |
---|
| 915 | setring Roriginal; |
---|
| 916 | ideal AI = imap(Rdummy,AI); |
---|
| 917 | //"#### AI:"; AI; |
---|
[64cf44] | 918 | return(AI); |
---|
[76bbd3] | 919 | } |
---|
| 920 | |
---|
[1a96a4] | 921 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 922 | proc mapToRatNormCurve(poly f, ideal AI) |
---|
| 923 | "USAGE: mapToRatNormCurve(f, AI); f polynomial, AI ideal |
---|
| 924 | ASSUME: The polynomial f is homogeneous in three variables and absolutely |
---|
| 925 | irreducible. |
---|
| 926 | The plane curve C defined by f is rational. |
---|
[a212fe] | 927 | The ideal AI is the adjoint ideal of C. |
---|
[1a96a4] | 928 | RETURN: ring with an ideal RNC. |
---|
| 929 | EXAMPLE: example mapToRatNormCurve; shows an example |
---|
| 930 | " |
---|
[6c265b] | 931 | { |
---|
[1a96a4] | 932 | int n = size(AI); |
---|
| 933 | int k; |
---|
| 934 | //if(n!=deg(f)-1){ERROR("not a rational curve");} |
---|
| 935 | def Roriginal = basering; |
---|
| 936 | ideal IC = f; |
---|
[6c265b] | 937 | list rl = ringlist(Roriginal); |
---|
[1a96a4] | 938 | /* begin workaround elimination*/ |
---|
| 939 | for(k = 1; k <= 3; k++) |
---|
| 940 | { |
---|
[6c265b] | 941 | rl[2][k] = "x("+string(k)+")"; |
---|
| 942 | } |
---|
[1a96a4] | 943 | for(k = 1; k <= n; k++) |
---|
| 944 | { |
---|
[6c265b] | 945 | rl[2][k+3] = "y("+string(k)+")"; |
---|
| 946 | } |
---|
[1a96a4] | 947 | rl[3]= list(list("dp",1:(3+n)),list("C",0)); |
---|
| 948 | def Relim = ring(rl); |
---|
| 949 | setring Relim; |
---|
| 950 | ideal IC = fetch(Roriginal,IC); |
---|
| 951 | ideal AI = fetch(Roriginal,AI); |
---|
[6c265b] | 952 | ideal J; |
---|
[1a96a4] | 953 | J = IC; |
---|
| 954 | for(k=1;k<=n;k++) |
---|
| 955 | { |
---|
| 956 | J=J,var(k+3)-AI[k]; |
---|
| 957 | } |
---|
| 958 | ideal SJ = std(J); |
---|
| 959 | intvec HJ = hilb(SJ,1); |
---|
| 960 | ideal RNC = eliminate(J,x(1)*x(2)*x(3),HJ); |
---|
| 961 | list rl = ringlist(Relim); |
---|
| 962 | list rl2 = rl[2]; |
---|
| 963 | rl[2] = list(rl2[4..n+3]); |
---|
| 964 | rl[3]= list(list("dp",1:n),list("C",0)); |
---|
| 965 | def Rtarget = ring(rl); |
---|
| 966 | setring Rtarget; |
---|
| 967 | ideal RNC = imap(Relim,RNC); |
---|
| 968 | /* end workaround elimination*/ |
---|
| 969 | export(RNC); |
---|
| 970 | int p = printlevel - voice + 3; |
---|
| 971 | dbprint(p,"//'mapToRatNorm' created a ring together with an ideal RNC."); |
---|
| 972 | dbprint(p,"// Supposing you typed, say, def RPn = mapToRatNorm(f,AI);"); |
---|
| 973 | dbprint(p,"// you may access the ideal by typing"); |
---|
| 974 | dbprint(p,"// setring RPn; RNC;"); |
---|
| 975 | return(Rtarget); |
---|
| 976 | } |
---|
| 977 | |
---|
| 978 | example |
---|
[75e82e] | 979 | { "EXAMPLE:"; echo=2; |
---|
[1a96a4] | 980 | ring R = 0,(x,y,z),dp; |
---|
| 981 | poly f = y^8-x^3*(z+x)^5; |
---|
| 982 | ideal adj = adjointIdeal(f); |
---|
| 983 | def Rn = mapToRatNormCurve(f,adj); |
---|
| 984 | setring(Rn); |
---|
| 985 | RNC; |
---|
| 986 | } |
---|
| 987 | |
---|
| 988 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 989 | proc rncAntiCanonicalMap(ideal I) |
---|
[6c265b] | 990 | "USAGE: rncAntiCanonicalMap(I); I ideal |
---|
[5e8ee4c] | 991 | ASSUME: I is a homogeneous ideal in the basering |
---|
[1a96a4] | 992 | defining a rational normal curve C in PP^n. |
---|
[6c265b] | 993 | NOTE: The procedure will fail or give a wrong output if I is not the |
---|
[1a96a4] | 994 | ideal of a rational normal curve. |
---|
| 995 | RETURN: ideal defining the anticanonical map C --> PP^(n-2). @* |
---|
[6c265b] | 996 | Note that the entries of the ideal should be considered as |
---|
[1a96a4] | 997 | representatives of elements in R/I, where R is the basering. |
---|
[6c265b] | 998 | THEORY: The anti-canonical map of a rational normal curve |
---|
[1a96a4] | 999 | maps C isomorpically to a rational normal curve in PP^(n-2). |
---|
| 1000 | KEYWORDS: rational normal curve, projection. |
---|
| 1001 | EXAMPLE: example rncAntiCanonicalMap; shows an example |
---|
| 1002 | " |
---|
[6c265b] | 1003 | { |
---|
[1a96a4] | 1004 | def Roriginal = basering; |
---|
| 1005 | list rl = ringlist(Roriginal); |
---|
| 1006 | rl[3] = list(list("dp",1:nvars(Roriginal)),list("C",0)); |
---|
[6c265b] | 1007 | def RoriginalDP = ring(rl); |
---|
[1a96a4] | 1008 | setring RoriginalDP; |
---|
| 1009 | ideal I = imap(Roriginal,I); |
---|
| 1010 | int cc = nvars(RoriginalDP)-2; |
---|
| 1011 | module AKD = Ext_R(cc,I); |
---|
| 1012 | qring qI = std(I); |
---|
| 1013 | matrix AKD = imap(RoriginalDP,AKD); |
---|
| 1014 | AKD = syz(transpose(AKD)); |
---|
| 1015 | ideal PR = submat(AKD,1..nrows(AKD),1); |
---|
| 1016 | setring Roriginal; |
---|
| 1017 | return(imap(qI,PR)); |
---|
| 1018 | } |
---|
| 1019 | |
---|
| 1020 | example |
---|
[75e82e] | 1021 | { "EXAMPLE:"; echo=2; |
---|
[1a96a4] | 1022 | ring R = 0,(x,y,z),dp; |
---|
| 1023 | poly f = y^8-x^3*(z+x)^5; |
---|
| 1024 | ideal adj = adjointIdeal(f); |
---|
| 1025 | def Rn = mapToRatNormCurve(f,adj); |
---|
| 1026 | setring(Rn); |
---|
| 1027 | RNC; |
---|
| 1028 | rncAntiCanonicalMap(RNC); |
---|
| 1029 | } |
---|
| 1030 | |
---|
| 1031 | |
---|
| 1032 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1033 | proc rncItProjOdd(ideal I) |
---|
[6c265b] | 1034 | "USAGE: rncItProjOdd(I); I ideal |
---|
[5e8ee4c] | 1035 | ASSUME: I is a homogeneous ideal in the basering with n+1 variables |
---|
[a212fe] | 1036 | defining a rational normal curve C in PP^n with n odd. |
---|
[1a96a4] | 1037 | NOTE: The procedure will fail or give a wrong output if I is not the |
---|
| 1038 | ideal of a rational normal curve. It will test whether n is odd. |
---|
| 1039 | RETURN: ideal PHI defining an isomorphic projection of C to PP^1.@* |
---|
[6c265b] | 1040 | Note that the entries of PHI should be considered as |
---|
[1a96a4] | 1041 | representatives of elements in R/I, where R is the basering. |
---|
[a212fe] | 1042 | THEORY: We iterate the procedure @ref{rncAntiCanonicalMap} to obtain PHI. |
---|
[1a96a4] | 1043 | KEYWORDS: rational normal curve, projection. |
---|
| 1044 | SEE ALSO: rncItProjEven. |
---|
| 1045 | EXAMPLE: example rncItProjOdd; shows an example |
---|
| 1046 | " |
---|
[6c265b] | 1047 | { |
---|
[1a96a4] | 1048 | int n = nvars(basering); |
---|
| 1049 | if((n mod 2) == 1){ERROR("Pn has even dimension");} |
---|
[6c265b] | 1050 | def Roriginal = basering; |
---|
[1a96a4] | 1051 | list rlo = ringlist(Roriginal); |
---|
| 1052 | rlo[3]= list(list("dp",1:n),list("C",0)); |
---|
| 1053 | int k; |
---|
| 1054 | for(k = 1; k <= n; k++) |
---|
| 1055 | { |
---|
[6c265b] | 1056 | rlo[2][k] = "z("+string(k)+")"; |
---|
| 1057 | } |
---|
[1a96a4] | 1058 | def RoriginalCopy = ring(rlo); |
---|
| 1059 | for(k = 1; k <= n; k++) |
---|
| 1060 | { |
---|
[6c265b] | 1061 | rlo[2][k] = "y("+string(k)+")"; |
---|
| 1062 | } |
---|
[1a96a4] | 1063 | def Rold = ring(rlo); |
---|
| 1064 | setring RoriginalCopy; |
---|
| 1065 | ideal PHI = maxideal(1); |
---|
| 1066 | setring Rold; |
---|
| 1067 | ideal J = fetch(Roriginal,I); |
---|
| 1068 | list rl2; |
---|
| 1069 | def Rnew; |
---|
| 1070 | def Rbig; |
---|
| 1071 | def Relim; |
---|
| 1072 | intvec HJJ; |
---|
| 1073 | while(n>2) |
---|
[6c265b] | 1074 | { |
---|
[1a96a4] | 1075 | ideal PR = rncAntiCanonicalMap(J); |
---|
| 1076 | list rl = ringlist(Rold); |
---|
| 1077 | Rbig = Rold + RoriginalCopy; |
---|
| 1078 | setring Rbig; |
---|
| 1079 | ideal PHI = imap(RoriginalCopy,PHI); |
---|
| 1080 | ideal dummy = imap(Rold,PR); |
---|
| 1081 | for(k = 1; k <= n; k++) |
---|
| 1082 | { |
---|
[6c265b] | 1083 | dummy = subst(dummy,var(k),PHI[k]); |
---|
| 1084 | } |
---|
[1a96a4] | 1085 | setring RoriginalCopy; |
---|
| 1086 | PHI = imap(Rbig,dummy); |
---|
| 1087 | /* begin workaround elimination*/ |
---|
| 1088 | setring Rold; |
---|
| 1089 | for(k = 1; k <= n; k++) |
---|
| 1090 | { |
---|
[6c265b] | 1091 | rl[2][k] = "x("+string(k)+")"; |
---|
| 1092 | } |
---|
[1a96a4] | 1093 | for(k = 1; k <= n-2; k++) |
---|
| 1094 | { |
---|
[6c265b] | 1095 | rl[2][k+n] = "y("+string(k)+")"; |
---|
| 1096 | } |
---|
[1a96a4] | 1097 | rl[3]= list(list("dp",1:(2*n-2)),list("C",0)); |
---|
| 1098 | Relim = ring(rl); |
---|
| 1099 | setring Relim; |
---|
| 1100 | ideal J = fetch(Rold,J); |
---|
| 1101 | ideal PR = fetch(Rold,PR); |
---|
| 1102 | ideal JJ = J; |
---|
| 1103 | poly pvar=1; |
---|
| 1104 | for(k = 1; k <= n; k++) |
---|
| 1105 | { |
---|
[6c265b] | 1106 | pvar = pvar*var(k); |
---|
| 1107 | } |
---|
[1a96a4] | 1108 | for(k=1;k<=n-2;k++) |
---|
| 1109 | { |
---|
| 1110 | JJ=JJ,var(k+n)-PR[k]; |
---|
| 1111 | } |
---|
| 1112 | ideal SJJ = std(JJ); |
---|
| 1113 | HJJ = hilb(SJJ,1); |
---|
| 1114 | J = eliminate(JJ,pvar,HJJ); |
---|
| 1115 | list rl = ringlist(Relim); |
---|
| 1116 | rl2 = rl[2]; |
---|
| 1117 | rl[2] = list(rl2[n+1..2*n-2]); |
---|
| 1118 | rl[3]= list(list("dp",1:(n-2)),list("C",0)); |
---|
| 1119 | Rnew = ring(rl); |
---|
| 1120 | setring Rnew; |
---|
| 1121 | ideal J = imap(Relim,J); |
---|
| 1122 | /* end workaround elimination*/ |
---|
| 1123 | Rold = Rnew; |
---|
| 1124 | setring Rold; |
---|
| 1125 | n = n-2; |
---|
| 1126 | } |
---|
| 1127 | setring Roriginal; |
---|
| 1128 | return(fetch(RoriginalCopy,PHI)); |
---|
| 1129 | } |
---|
| 1130 | |
---|
| 1131 | example |
---|
[75e82e] | 1132 | { "EXAMPLE:"; echo=2; |
---|
[1a96a4] | 1133 | ring R = 0,(x,y,z),dp; |
---|
| 1134 | poly f = -x7-10x5y2-10x4y3-3x3y4+8x2y5+7xy6+11y7+3x6+10x5y +30x4y2 |
---|
| 1135 | +26x3y3-13x2y4-29xy5-33y6-3x5-20x4y-33x3y2-8x2y3+37xy4+33y5 |
---|
| 1136 | +x4+10x3y+13x2y2-15xy3-11y4; |
---|
| 1137 | f = homog(f,z); |
---|
| 1138 | ideal adj = adjointIdeal(f); |
---|
| 1139 | def Rn = mapToRatNormCurve(f,adj); |
---|
| 1140 | setring(Rn); |
---|
| 1141 | RNC; |
---|
| 1142 | rncItProjOdd(RNC); |
---|
| 1143 | } |
---|
| 1144 | |
---|
| 1145 | |
---|
| 1146 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1147 | proc rncItProjEven(ideal I) |
---|
[6c265b] | 1148 | "USAGE: rncItProjEven(I); I ideal |
---|
[5e8ee4c] | 1149 | ASSUME: I is a homogeneous ideal in the basering with n+1 variables |
---|
[a212fe] | 1150 | defining a rational normal curve C in PP^n with n even. |
---|
[6c265b] | 1151 | NOTE: The procedure will fail or give a wrong output if I is not the |
---|
[1a96a4] | 1152 | ideal of a rational normal curve. It will test whether n is odd. |
---|
[a212fe] | 1153 | RETURN: ring with an ideal CONIC defining a conic C2 in PP^2.@* |
---|
[5e8ee4c] | 1154 | In addition, an ideal PHI in the basering defining an isomorphic |
---|
[a212fe] | 1155 | projection of C to C2 will be exported.@* |
---|
[6c265b] | 1156 | Note that the entries of PHI should be considered as |
---|
[1a96a4] | 1157 | representatives of elements in R/I, where R is the basering. |
---|
[a212fe] | 1158 | THEORY: We iterate the procedure @ref{rncAntiCanonicalMap} to obtain PHI. |
---|
[1a96a4] | 1159 | KEYWORDS: rational normal curve, projection. |
---|
| 1160 | SEE ALSO: rncItProjOdd. |
---|
| 1161 | EXAMPLE: example rncItProjEven; shows an example |
---|
| 1162 | " |
---|
[6c265b] | 1163 | { |
---|
[1a96a4] | 1164 | int n = nvars(basering); |
---|
| 1165 | if((n mod 2) == 0){ERROR("Pn has odd dimension");} |
---|
[6c265b] | 1166 | def Roriginal = basering; |
---|
[1a96a4] | 1167 | list rlo = ringlist(Roriginal); |
---|
| 1168 | rlo[3]= list(list("dp",1:n),list("C",0)); |
---|
| 1169 | int k; |
---|
| 1170 | for(k = 1; k <= n; k++) |
---|
| 1171 | { |
---|
[6c265b] | 1172 | rlo[2][k] = "z("+string(k)+")"; |
---|
| 1173 | } |
---|
[1a96a4] | 1174 | def RoriginalCopy = ring(rlo); |
---|
| 1175 | for(k = 1; k <= n; k++) |
---|
| 1176 | { |
---|
[6c265b] | 1177 | rlo[2][k] = "y("+string(k)+")"; |
---|
| 1178 | } |
---|
[1a96a4] | 1179 | def Rold = ring(rlo); |
---|
| 1180 | setring RoriginalCopy; |
---|
| 1181 | ideal PHI = maxideal(1); |
---|
| 1182 | setring Rold; |
---|
| 1183 | ideal J = fetch(Roriginal,I); |
---|
| 1184 | list rl2; |
---|
| 1185 | def Rnew; |
---|
| 1186 | def Rbig; |
---|
| 1187 | def Relim; |
---|
| 1188 | intvec HJJ; |
---|
| 1189 | while(n>3) |
---|
[6c265b] | 1190 | { |
---|
[1a96a4] | 1191 | ideal PR = rncAntiCanonicalMap(J); |
---|
| 1192 | list rl = ringlist(Rold); |
---|
| 1193 | Rbig = Rold + RoriginalCopy; |
---|
| 1194 | setring Rbig; |
---|
| 1195 | ideal PHI = imap(RoriginalCopy,PHI); |
---|
| 1196 | ideal dummy = imap(Rold,PR); |
---|
| 1197 | for(k = 1; k <= n; k++) |
---|
| 1198 | { |
---|
[6c265b] | 1199 | dummy = subst(dummy,var(k),PHI[k]); |
---|
| 1200 | } |
---|
[1a96a4] | 1201 | setring RoriginalCopy; |
---|
| 1202 | PHI = imap(Rbig,dummy); |
---|
| 1203 | /* begin workaround elimination*/ |
---|
| 1204 | setring Rold; |
---|
| 1205 | for(k = 1; k <= n; k++) |
---|
| 1206 | { |
---|
[6c265b] | 1207 | rl[2][k] = "x("+string(k)+")"; |
---|
| 1208 | } |
---|
[1a96a4] | 1209 | for(k = 1; k <= n-2; k++) |
---|
| 1210 | { |
---|
[6c265b] | 1211 | rl[2][k+n] = "y("+string(k)+")"; |
---|
| 1212 | } |
---|
[1a96a4] | 1213 | rl[3]= list(list("dp",1:(2*n-2)),list("C",0)); |
---|
| 1214 | Relim = ring(rl); |
---|
| 1215 | setring Relim; |
---|
| 1216 | ideal J = fetch(Rold,J); |
---|
| 1217 | ideal PR = fetch(Rold,PR); |
---|
| 1218 | ideal JJ = J; |
---|
| 1219 | poly pvar=1; |
---|
| 1220 | for(k = 1; k <= n; k++) |
---|
| 1221 | { |
---|
[6c265b] | 1222 | pvar = pvar*var(k); |
---|
| 1223 | } |
---|
[1a96a4] | 1224 | for(k=1;k<=n-2;k++) |
---|
| 1225 | { |
---|
| 1226 | JJ=JJ,var(k+n)-PR[k]; |
---|
| 1227 | } |
---|
| 1228 | ideal SJJ = std(JJ); |
---|
| 1229 | HJJ = hilb(SJJ,1); |
---|
| 1230 | J = eliminate(JJ,pvar,HJJ); |
---|
| 1231 | list rl = ringlist(Relim); |
---|
| 1232 | rl2 = rl[2]; |
---|
| 1233 | rl[2] = list(rl2[n+1..2*n-2]); |
---|
| 1234 | rl[3]= list(list("dp",1:(n-2)),list("C",0)); |
---|
| 1235 | Rnew = ring(rl); |
---|
| 1236 | setring Rnew; |
---|
| 1237 | ideal J = imap(Relim,J); |
---|
| 1238 | /* end workaround elimination*/ |
---|
| 1239 | Rold = Rnew; |
---|
| 1240 | setring Rold; |
---|
[6c265b] | 1241 | n = n-2; |
---|
[1a96a4] | 1242 | } |
---|
| 1243 | poly CONIC = J[1]; |
---|
| 1244 | export(CONIC); |
---|
| 1245 | setring Roriginal; |
---|
| 1246 | ideal PHI = fetch(RoriginalCopy,PHI); |
---|
| 1247 | export(PHI); |
---|
| 1248 | return(Rold); |
---|
| 1249 | } |
---|
| 1250 | |
---|
| 1251 | example |
---|
[75e82e] | 1252 | { "EXAMPLE:"; echo=2; |
---|
[1a96a4] | 1253 | ring R = 0,(x,y,z),dp; |
---|
| 1254 | poly f = y^8-x^3*(z+x)^5; |
---|
| 1255 | ideal adj = adjointIdeal(f); |
---|
| 1256 | def Rn = mapToRatNormCurve(f,adj); |
---|
| 1257 | setring(Rn); |
---|
| 1258 | RNC; |
---|
| 1259 | def Rc = rncItProjEven(RNC); |
---|
| 1260 | PHI; |
---|
| 1261 | setring Rc; |
---|
| 1262 | CONIC; |
---|
| 1263 | } |
---|
| 1264 | |
---|
| 1265 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1266 | static proc geomGenusIB(poly f, list #) |
---|
| 1267 | "USAGE: geomGenusIB(f [, L]); f poly, L optional list representing the |
---|
| 1268 | integral basis as returned by @ref integralBasisJ. |
---|
| 1269 | ASSUME: The basering must be a polynomial ring in three variables, say x,y,z, |
---|
| 1270 | with coefficients in Q. @* |
---|
| 1271 | The polynomial f must be homogeneous and absolutely irreducible.@* |
---|
[6c265b] | 1272 | Its dehomogenization with respect to the third variable must be monic |
---|
[1a96a4] | 1273 | as a polynomial in the second variable (that is, the curve C = {f = 0} |
---|
| 1274 | does not contain the point (0:1:0)).@* |
---|
[6c265b] | 1275 | The curve C is not allowed to have singularities |
---|
[1a96a4] | 1276 | at infinity (z = 0). @* |
---|
| 1277 | NOTE: The last two conditions can be met by a suitable change of coordinates in PGL(3) |
---|
| 1278 | as applied in the procedure @ref adjointIdeal. The other conditions |
---|
| 1279 | can be tested using @ref checkAssumptions.@* |
---|
| 1280 | RETURN: int, the geometric genus of C. |
---|
| 1281 | THEORY: We compute an integral basis of the integral closure of the coordinate |
---|
| 1282 | ring of C and from that the geometric genus.@* |
---|
| 1283 | KEYWORDS: geometric genus, plane curves. |
---|
| 1284 | SEE ALSO: genus. |
---|
| 1285 | " |
---|
[6c265b] | 1286 | { |
---|
[1a96a4] | 1287 | int bb = size(#); |
---|
| 1288 | poly dhf = subst(f,var(3),1); |
---|
| 1289 | def Roriginal = basering; |
---|
| 1290 | list rl = ringlist(Roriginal); |
---|
| 1291 | rl[2] = list(var(1), var(2)); |
---|
| 1292 | rl[3] = list(list("dp",1:2),list("C",0)); |
---|
| 1293 | def Rdummy = ring(rl); |
---|
| 1294 | setring Rdummy; |
---|
| 1295 | poly f = imap(Roriginal,dhf); |
---|
| 1296 | list LIntB; |
---|
| 1297 | if(bb == 0) |
---|
| 1298 | { |
---|
[6c265b] | 1299 | LIntB = integralBasis(f,2,"isIrred"); |
---|
[1a96a4] | 1300 | } |
---|
| 1301 | else |
---|
| 1302 | { |
---|
| 1303 | LIntB = imap(Roriginal,#); |
---|
| 1304 | } |
---|
| 1305 | ideal IB = LIntB[1]; |
---|
| 1306 | poly d = LIntB[2]; |
---|
| 1307 | int ud = deg(d); |
---|
| 1308 | int sL = size(IB); |
---|
| 1309 | int k; |
---|
[62c2b0] | 1310 | int gg = (sL-1)*(sL-2) div 2-sL*ud; |
---|
[1a96a4] | 1311 | for(k = 1; k <= sL; k++) |
---|
| 1312 | { |
---|
[6c265b] | 1313 | gg = gg + deg(gcd(d,IB[k])); |
---|
[1a96a4] | 1314 | } |
---|
| 1315 | setring Roriginal; |
---|
| 1316 | return(gg); |
---|
| 1317 | } |
---|
| 1318 | |
---|
| 1319 | |
---|
| 1320 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1321 | static proc geomGenusLA(poly F) |
---|
| 1322 | "USAGE: geomGenusLA(F); F polynomial |
---|
| 1323 | ASSUME: The basering must be a polynomial ring in three variables. @* |
---|
| 1324 | The polynomial F must be homogeneous.@* |
---|
| 1325 | RETURN: list L: |
---|
| 1326 | @texinfo |
---|
| 1327 | @table @asis |
---|
| 1328 | @item @code{L[1]}; int: |
---|
| 1329 | the geometric genus p_g = p_a - delta of the projective |
---|
| 1330 | curve C defined by F, where p_a is the arithmetic genus. |
---|
| 1331 | @item @code{L[2]}; int: |
---|
| 1332 | is positive if C has singularities other |
---|
| 1333 | than ordinary multiple points.@* |
---|
| 1334 | @item @code{L[3]}; list: |
---|
| 1335 | consists of one list for each primary component |
---|
| 1336 | of the singular locus of C which correponds to a set of conjugated |
---|
| 1337 | ordinary multiple points. Each list consists of an int, the |
---|
| 1338 | multiplicity of the points, and an ideal, the primary component. |
---|
| 1339 | @end table |
---|
| 1340 | @end texinfo |
---|
| 1341 | NOTE: delta is the sum of all local delta-invariants of the singularities, |
---|
| 1342 | i.e. dim(R'/R), R' the normalization of the local ring R of the |
---|
| 1343 | singularity. @* |
---|
| 1344 | SEE ALSO: genus |
---|
| 1345 | " |
---|
[6c265b] | 1346 | { |
---|
[1a96a4] | 1347 | int w = printlevel-voice+2; // w=printlevel (default: w=0) |
---|
| 1348 | int d = deg(F); |
---|
| 1349 | def R = basering; |
---|
| 1350 | execute("ring S=("+charstr(R)+"),(x,y,t),dp;"); |
---|
| 1351 | execute("ring C=("+charstr(R)+"),(x,y),ds;"); |
---|
[62c2b0] | 1352 | int genus=(d-1)*(d-2) div 2; |
---|
[1a96a4] | 1353 | if(w>=1){"the arithmetic genus of the plane curve:";genus;pause();} |
---|
| 1354 | |
---|
| 1355 | int delt,deltaloc,deltainf,tau,tauinf,cusps,iloc,iglob,l,nsing, |
---|
| 1356 | tauloc,tausing,k,rat,nbranchinf,nbranch,nodes,cuspsinf,nodesinf; |
---|
[6c265b] | 1357 | list inv; |
---|
[1a96a4] | 1358 | execute("ring newR=("+charstr(R)+"),(x,y),dp;"); |
---|
| 1359 | //the singularities at the affine part |
---|
| 1360 | map sigma=R,var(1),var(2),1; |
---|
| 1361 | ideal I=sigma(F); |
---|
| 1362 | |
---|
[6c265b] | 1363 | list OMPButNodes; |
---|
| 1364 | int sizeOMPButNodes; |
---|
| 1365 | int NotOnlyOMPPlusCusps; |
---|
[1a96a4] | 1366 | |
---|
| 1367 | ideal I1=jacob(I); |
---|
| 1368 | matrix Hess[2][2]=jacob(I1); |
---|
| 1369 | ideal ID=I+I1+ideal(det(Hess));//singular locus of I+I1 |
---|
| 1370 | ideal radID=std(radical(ID));//the non-nodal locus |
---|
| 1371 | if(w>=1){"the non-nodal locus:";"";radID;pause();"";} |
---|
| 1372 | if(deg(radID[1])==0) |
---|
| 1373 | { |
---|
| 1374 | ideal IDsing=1; |
---|
| 1375 | } |
---|
| 1376 | else |
---|
| 1377 | { |
---|
| 1378 | ideal IDsing=minor(jacob(ID),2)+radID;//singular locus of ID |
---|
| 1379 | } |
---|
| 1380 | |
---|
| 1381 | iglob=vdim(std(IDsing)); |
---|
| 1382 | |
---|
| 1383 | ideal radIDsing = 1; |
---|
| 1384 | |
---|
| 1385 | if(iglob!=0)//computation of the radical of IDsing |
---|
| 1386 | { |
---|
| 1387 | radIDsing=reduce(IDsing,radID); |
---|
| 1388 | if(size(radIDsing)==0) |
---|
| 1389 | { |
---|
| 1390 | radIDsing=radID; |
---|
| 1391 | attrib(radIDsing,"isSB",1); |
---|
| 1392 | } |
---|
| 1393 | else |
---|
| 1394 | { |
---|
| 1395 | radIDsing=std(radical(IDsing)); |
---|
| 1396 | } |
---|
| 1397 | iglob=vdim(radIDsing); |
---|
| 1398 | if((w>=1)&&(iglob)) |
---|
| 1399 | {"the non-nodal-cuspidal locus:";radIDsing;pause();"";} |
---|
| 1400 | } |
---|
| 1401 | cusps=vdim(radID)-iglob; |
---|
| 1402 | |
---|
| 1403 | ideal NodesPlusCusps = radical(sat(I+I1, radIDsing)[1]); |
---|
[6c265b] | 1404 | |
---|
[1a96a4] | 1405 | nsing=nsing+cusps; |
---|
| 1406 | |
---|
| 1407 | if(iglob==0) |
---|
| 1408 | { |
---|
| 1409 | if(w>=1){" there are only cusps and nodes";"";} |
---|
| 1410 | tau=vdim(std(I+jacob(I))); |
---|
| 1411 | tauinf=tauinf+tau; |
---|
| 1412 | nodes=tau-2*cusps; |
---|
| 1413 | delt=nodes+cusps; |
---|
| 1414 | nbranch=2*tau-3*cusps; |
---|
| 1415 | nsing=nsing+nodes; |
---|
| 1416 | } |
---|
| 1417 | else |
---|
| 1418 | { |
---|
| 1419 | if(w>=1){"the non-nodal-cuspidal singularities";"";} |
---|
| 1420 | setring C; |
---|
[fa932ac] | 1421 | ideal I1=imap(newR,radIDsing); |
---|
[1a96a4] | 1422 | iloc=vdim(std(I1)); |
---|
| 1423 | if(iglob==iloc) |
---|
| 1424 | { |
---|
| 1425 | if(w>=1){"only cusps and nodes outside (0,0,1)";} |
---|
| 1426 | setring newR; |
---|
| 1427 | tau=vdim(std(I+jacob(I))); |
---|
| 1428 | tauinf=tauinf+tau; |
---|
| 1429 | inv=deltaLocMod(I[1],maxideal(1)); |
---|
| 1430 | delt=inv[1]; |
---|
| 1431 | tauloc=inv[2]; |
---|
| 1432 | nodes=tau-tauloc-2*cusps; |
---|
| 1433 | nsing=nsing+nodes; |
---|
[fa932ac] | 1434 | if(inv[2]!=0) |
---|
| 1435 | { |
---|
| 1436 | nsing=nsing+1; |
---|
| 1437 | } |
---|
[1a96a4] | 1438 | nbranch=inv[3]+ 2*nodes+cusps; |
---|
| 1439 | delt=delt+nodes+cusps; |
---|
| 1440 | if((w>=1)&&(inv[2]==0)){"smooth at (0,0,1)";} |
---|
[6c265b] | 1441 | if(inv[4]!=0) |
---|
[1a96a4] | 1442 | { |
---|
| 1443 | OMPButNodes = insert(OMPButNodes,list(inv[4],maxideal(1)), |
---|
| 1444 | sizeOMPButNodes); |
---|
| 1445 | sizeOMPButNodes = size(OMPButNodes); // new |
---|
| 1446 | } |
---|
| 1447 | else |
---|
| 1448 | { |
---|
[6c265b] | 1449 | NotOnlyOMPPlusCusps = NotOnlyOMPPlusCusps + 1; |
---|
[1a96a4] | 1450 | } |
---|
| 1451 | } |
---|
| 1452 | else |
---|
| 1453 | { |
---|
| 1454 | setring newR; |
---|
[fa932ac] | 1455 | list pr=minAssGTZ(radIDsing); |
---|
[1a96a4] | 1456 | if(w>=1){pr;} |
---|
| 1457 | |
---|
| 1458 | for(k=1;k<=size(pr);k++) |
---|
| 1459 | { |
---|
| 1460 | if(w>=1){nsing=nsing+vdim(std(pr[k]));} |
---|
| 1461 | inv=deltaLocMod(I[1],pr[k]); |
---|
| 1462 | delt=delt+inv[1]; |
---|
| 1463 | tausing=tausing+inv[2]; |
---|
| 1464 | nbranch=nbranch+inv[3]; |
---|
[6c265b] | 1465 | if(inv[4]!=0) |
---|
[1a96a4] | 1466 | { |
---|
| 1467 | OMPButNodes = insert(OMPButNodes,list(inv[4],pr[k]), |
---|
| 1468 | sizeOMPButNodes); |
---|
[6c265b] | 1469 | sizeOMPButNodes = size(OMPButNodes); |
---|
[1a96a4] | 1470 | } |
---|
| 1471 | else |
---|
| 1472 | { |
---|
[6c265b] | 1473 | NotOnlyOMPPlusCusps = NotOnlyOMPPlusCusps + 1; |
---|
[1a96a4] | 1474 | } |
---|
| 1475 | } |
---|
| 1476 | tau=vdim(std(I+jacob(I))); |
---|
| 1477 | tauinf=tauinf+tau; |
---|
| 1478 | nodes=tau-tausing-2*cusps; |
---|
| 1479 | nsing=nsing+nodes; |
---|
| 1480 | delt=delt+nodes+cusps; |
---|
| 1481 | nbranch=nbranch+2*nodes+cusps; |
---|
| 1482 | } |
---|
| 1483 | } |
---|
| 1484 | genus=genus-delt-deltainf; |
---|
| 1485 | if(w>=1) |
---|
| 1486 | { |
---|
| 1487 | "The projected plane curve has locally:";""; |
---|
| 1488 | "singularities:";nsing; |
---|
| 1489 | "branches:";nbranch+nbranchinf; |
---|
| 1490 | "nodes:"; nodes+nodesinf; |
---|
| 1491 | "cusps:";cusps+cuspsinf; |
---|
| 1492 | "Tjurina number:";tauinf; |
---|
| 1493 | "Milnor number:";2*(delt+deltainf)-nbranch-nbranchinf+nsing; |
---|
| 1494 | "delta of the projected curve:";delt+deltainf; |
---|
| 1495 | //"delta of the curve:";p_a-genus; |
---|
| 1496 | "genus:";genus; |
---|
| 1497 | "===================================================="; |
---|
| 1498 | ""; |
---|
| 1499 | } |
---|
| 1500 | setring R; |
---|
[6c265b] | 1501 | if(sizeOMPButNodes>0) |
---|
[1a96a4] | 1502 | { |
---|
[6c265b] | 1503 | list OMPButNodes = fetch(newR,OMPButNodes); |
---|
[1a96a4] | 1504 | } |
---|
| 1505 | return(list(genus,NotOnlyOMPPlusCusps,OMPButNodes, |
---|
| 1506 | fetch(newR,NodesPlusCusps))); |
---|
| 1507 | } |
---|
| 1508 | |
---|
| 1509 | |
---|
| 1510 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1511 | static proc deltaLocMod(poly f,ideal singL) |
---|
| 1512 | "USAGE: deltaLoc(f,J); f poly, J ideal |
---|
| 1513 | ASSUME: f is reduced bivariate polynomial; basering has exactly two variables; |
---|
| 1514 | J is irreducible prime component of the singular locus of f (e.g., one |
---|
| 1515 | entry of the output of @code{minAssGTZ(I);}, I = <f,jacob(f)>). |
---|
| 1516 | RETURN: list L: |
---|
| 1517 | @texinfo |
---|
| 1518 | @table @asis |
---|
| 1519 | @item @code{L[1]}; int: |
---|
| 1520 | the sum of (local) delta invariants of f at the (conjugated) singular |
---|
| 1521 | points given by J. |
---|
| 1522 | @item @code{L[2]}; int: |
---|
| 1523 | the sum of (local) Tjurina numbers of f at the (conjugated) singular |
---|
| 1524 | points given by J. |
---|
| 1525 | @item @code{L[3]}; int: |
---|
| 1526 | the sum of (local) number of branches of f at the (conjugated) |
---|
| 1527 | singular points given by J. |
---|
| 1528 | @item @code{L[3]}; int: |
---|
| 1529 | the multiplicity of f at the (conjugated) singular points given by J, |
---|
| 1530 | if these are ordinary multiple points, and 0 otherwise. |
---|
| 1531 | @end table |
---|
| 1532 | @end texinfo |
---|
| 1533 | NOTE: procedure makes use of @code{execute}; increasing printlevel displays |
---|
| 1534 | more comments (default: printlevel=0). |
---|
| 1535 | SEE ALSO: deltaLoc, delta, tjurina |
---|
| 1536 | KEYWORDS: delta invariant; Tjurina number |
---|
| 1537 | " |
---|
| 1538 | { |
---|
| 1539 | option(redSB); |
---|
| 1540 | def R=basering; |
---|
| 1541 | execute("ring S=("+charstr(R)+"),(x,y),lp;"); |
---|
| 1542 | map phi=R,x,y; |
---|
| 1543 | ideal singL=phi(singL); |
---|
| 1544 | singL=simplify(std(singL),1); |
---|
| 1545 | attrib(singL,"isSB",1); |
---|
| 1546 | int d=vdim(singL); |
---|
| 1547 | poly f=phi(f); |
---|
| 1548 | int i; |
---|
| 1549 | int w = printlevel-voice+2; // w=printlevel (default: w=0) |
---|
| 1550 | if(d==1) |
---|
| 1551 | { |
---|
| 1552 | map alpha=S,var(1)-singL[2][2],var(2)-singL[1][2]; |
---|
| 1553 | f=alpha(f); |
---|
| 1554 | execute("ring C=("+charstr(S)+"),("+varstr(S)+"),ds;"); |
---|
| 1555 | poly f=imap(S,f); |
---|
| 1556 | ideal singL=imap(S,singL); |
---|
| 1557 | if((w>=1)&&(ord(f)>=2)) |
---|
| 1558 | { |
---|
| 1559 | "local analysis of the singularities";""; |
---|
| 1560 | basering; |
---|
| 1561 | singL; |
---|
| 1562 | f; |
---|
| 1563 | pause(); |
---|
| 1564 | } |
---|
| 1565 | } |
---|
| 1566 | else |
---|
| 1567 | { |
---|
| 1568 | poly p; |
---|
| 1569 | poly c; |
---|
| 1570 | map psi; |
---|
| 1571 | number co; |
---|
| 1572 | |
---|
| 1573 | while((deg(lead(singL[1]))>1)&&(deg(lead(singL[2]))>1)) |
---|
| 1574 | { |
---|
| 1575 | psi=S,x,y+random(-100,100)*x; |
---|
| 1576 | singL=psi(singL); |
---|
| 1577 | singL=std(singL); |
---|
| 1578 | f=psi(f); |
---|
| 1579 | } |
---|
| 1580 | |
---|
| 1581 | if(deg(lead(singL[2]))==1) |
---|
| 1582 | { |
---|
| 1583 | p=singL[1]; |
---|
| 1584 | c=singL[2]-lead(singL[2]); |
---|
| 1585 | co=leadcoef(singL[2]); |
---|
| 1586 | } |
---|
| 1587 | if(deg(lead(singL[1]))==1) |
---|
| 1588 | { |
---|
| 1589 | psi=S,y,x; |
---|
| 1590 | f=psi(f); |
---|
| 1591 | singL=psi(singL); |
---|
| 1592 | p=singL[2]; |
---|
| 1593 | c=singL[1]-lead(singL[1]); |
---|
| 1594 | co=leadcoef(singL[1]); |
---|
| 1595 | } |
---|
| 1596 | |
---|
| 1597 | execute("ring B=("+charstr(S)+"),a,dp;"); |
---|
| 1598 | map beta=S,a,a; |
---|
| 1599 | poly p=beta(p); |
---|
| 1600 | |
---|
| 1601 | execute("ring C=("+charstr(S)+",a),("+varstr(S)+"),ds;"); |
---|
| 1602 | number p=number(imap(B,p)); |
---|
| 1603 | minpoly=p; |
---|
| 1604 | |
---|
| 1605 | map iota=S,a,a; |
---|
| 1606 | number c=number(iota(c)); |
---|
| 1607 | number co=iota(co); |
---|
| 1608 | |
---|
| 1609 | map alpha=S,x-c/co,y+a; |
---|
| 1610 | poly f=alpha(f); |
---|
| 1611 | f=cleardenom(f); |
---|
| 1612 | if((w>=1)&&(ord(f)>=2)) |
---|
| 1613 | { |
---|
| 1614 | "local analysis of the singularities";""; |
---|
| 1615 | basering; |
---|
| 1616 | alpha; |
---|
| 1617 | f; |
---|
| 1618 | pause(); |
---|
| 1619 | ""; |
---|
| 1620 | } |
---|
| 1621 | } |
---|
[6c265b] | 1622 | int intMult = deg(lead(f)); |
---|
| 1623 | poly fdummy = f; |
---|
| 1624 | poly gdummy = lead(f); |
---|
| 1625 | int ivr = 1; |
---|
| 1626 | while(ivr) |
---|
[1a96a4] | 1627 | { |
---|
| 1628 | fdummy = fdummy - lead(fdummy); |
---|
| 1629 | if((fdummy ==0) || (deg(lead(fdummy))>intMult)){break;} |
---|
[6c265b] | 1630 | gdummy = gdummy + lead(fdummy); |
---|
[1a96a4] | 1631 | } |
---|
[fee24e] | 1632 | poly SQRpart = sqrfree(gdummy, 3); |
---|
[d17cfd] | 1633 | int IntForRet; |
---|
[fee24e] | 1634 | if(deg(SQRpart)==intMult) |
---|
[1a96a4] | 1635 | { |
---|
[6c265b] | 1636 | IntForRet = intMult; |
---|
[1a96a4] | 1637 | } |
---|
| 1638 | option(noredSB); |
---|
| 1639 | ideal fstd=std(ideal(f)+jacob(f)); |
---|
| 1640 | poly hc=highcorner(fstd); |
---|
| 1641 | int tau=vdim(fstd); |
---|
| 1642 | int o=ord(f); |
---|
| 1643 | int delt,nb; |
---|
| 1644 | |
---|
| 1645 | if(tau==0) //smooth case |
---|
| 1646 | { |
---|
| 1647 | setring R; |
---|
| 1648 | return(list(0,0,1,0)); |
---|
| 1649 | } |
---|
| 1650 | if((char(basering)>=181)||(char(basering)==0)) |
---|
| 1651 | { |
---|
| 1652 | if(o==2) //A_k-singularity |
---|
| 1653 | { |
---|
| 1654 | if(w>=1){"A_k-singularity";"";} |
---|
| 1655 | setring R; |
---|
[62c2b0] | 1656 | delt=(tau+1) div 2; |
---|
[1a96a4] | 1657 | return(list(d*delt,d*tau,d*(2*delt-tau+1),IntForRet)); |
---|
| 1658 | } |
---|
| 1659 | if((lead(f)==var(1)*var(2)^2)||(lead(f)==var(1)^2*var(2))) |
---|
| 1660 | { |
---|
| 1661 | if(w>=1){"D_k- singularity";"";} |
---|
| 1662 | |
---|
| 1663 | setring R; |
---|
[62c2b0] | 1664 | delt=(tau+2) div 2; |
---|
[1a96a4] | 1665 | return(list(d*delt,d*tau,d*(2*delt-tau+1),IntForRet)); |
---|
| 1666 | } |
---|
| 1667 | |
---|
| 1668 | int mu=vdim(std(jacob(f))); |
---|
| 1669 | poly g=f+var(1)^mu+var(2)^mu; //to obtain a convenient Newton-polygon |
---|
| 1670 | |
---|
| 1671 | list NP=newtonpoly(g); |
---|
| 1672 | if(w>=1){"Newton-Polygon:";NP;"";} |
---|
| 1673 | int s=size(NP); |
---|
| 1674 | |
---|
| 1675 | if(is_NND(f,mu,NP)) |
---|
| 1676 | { // the Newton-polygon is non-degenerate |
---|
| 1677 | // compute nb, the number of branches |
---|
| 1678 | for(i=1;i<=s-1;i++) |
---|
| 1679 | { |
---|
| 1680 | nb=nb+gcd(NP[i][2]-NP[i+1][2],NP[i][1]-NP[i+1][1]); |
---|
| 1681 | } |
---|
| 1682 | if(w>=1){"Newton-Polygon is non-degenerated";"";} |
---|
[62c2b0] | 1683 | return(list(d*(mu+nb-1) div 2,d*tau,d*nb,IntForRet)); |
---|
[1a96a4] | 1684 | } |
---|
| 1685 | |
---|
| 1686 | if(w>=1){"Newton-Polygon is degenerated";"";} |
---|
| 1687 | |
---|
| 1688 | // the following can certainly be made more efficient when replacing |
---|
| 1689 | // 'hnexpansion' (used only for computing number of branches) by |
---|
| 1690 | // successive blowing-up + test if Newton polygon degenerate: |
---|
| 1691 | if(s>2) // splitting of f |
---|
| 1692 | { |
---|
| 1693 | if(w>=1){"Newton polygon can be used for splitting";"";} |
---|
| 1694 | intvec v=NP[1][2]-NP[2][2],NP[2][1]; |
---|
| 1695 | int de=w_deg(g,v); |
---|
| 1696 | int st=w_deg(hc,v)+v[1]+v[2]; |
---|
| 1697 | poly f1=var(2)^NP[2][2]; |
---|
| 1698 | poly f2=jet(g,de,v)/var(2)^NP[2][2]; |
---|
| 1699 | poly h=g-f1*f2; |
---|
| 1700 | de=w_deg(h,v); |
---|
| 1701 | poly k; |
---|
| 1702 | ideal wi=var(2)^NP[2][2],f2; |
---|
| 1703 | matrix li; |
---|
| 1704 | while(de<st) |
---|
| 1705 | { |
---|
| 1706 | k=jet(h,de,v); |
---|
| 1707 | li=lift(wi,k); |
---|
| 1708 | f1=f1+li[2,1]; |
---|
| 1709 | f2=f2+li[1,1]; |
---|
| 1710 | h=g-f1*f2; |
---|
| 1711 | de=w_deg(h,v); |
---|
| 1712 | } |
---|
| 1713 | nb=deltaLocMod(f1,maxideal(1))[3]+deltaLocMod(f2,maxideal(1))[3]; |
---|
| 1714 | setring R; |
---|
[62c2b0] | 1715 | return(list(d*(mu+nb-1) div 2,d*tau,d*nb,IntForRet)); |
---|
[1a96a4] | 1716 | } |
---|
| 1717 | |
---|
| 1718 | f=jet(f,deg(hc)+2); |
---|
| 1719 | if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";} |
---|
| 1720 | ideal fac=factorize(f,1); |
---|
| 1721 | if(size(fac)>1) |
---|
| 1722 | { |
---|
| 1723 | nb=0; |
---|
| 1724 | for(i=1;i<=size(fac);i++) |
---|
| 1725 | { |
---|
| 1726 | nb=nb+deltaLocMod(fac[i],maxideal(1))[3]; |
---|
| 1727 | } |
---|
| 1728 | setring R; |
---|
[62c2b0] | 1729 | return(list(d*(mu+nb-1) div 2,d*tau,d*nb,IntForRet)); |
---|
[1a96a4] | 1730 | } |
---|
| 1731 | list HNEXP=hnexpansion(f); |
---|
| 1732 | if (typeof(HNEXP[1])=="ring") { |
---|
| 1733 | def altring = basering; |
---|
| 1734 | def HNEring = HNEXP[1]; setring HNEring; |
---|
| 1735 | nb=size(hne); |
---|
| 1736 | setring R; |
---|
| 1737 | kill HNEring; |
---|
| 1738 | } |
---|
| 1739 | else |
---|
| 1740 | { |
---|
| 1741 | nb=size(HNEXP); |
---|
| 1742 | } |
---|
[62c2b0] | 1743 | return(list(d*(mu+nb-1) div 2,d*tau,d*nb,IntForRet)); |
---|
[1a96a4] | 1744 | } |
---|
| 1745 | else //the case of small characteristic |
---|
| 1746 | { |
---|
| 1747 | f=jet(f,deg(hc)+2); |
---|
| 1748 | if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";} |
---|
| 1749 | delt=delta(f); |
---|
| 1750 | return(list(d*delt,d*tau,d,IntForRet)); |
---|
| 1751 | } |
---|
| 1752 | } |
---|
| 1753 | |
---|
| 1754 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1755 | proc paraConic(poly q) |
---|
| 1756 | "USAGE: paraConic(q); q poly |
---|
[6c265b] | 1757 | ASSUME: The basering must be a polynomial ring in three variables with |
---|
[1a96a4] | 1758 | coefficients in Q. @* |
---|
| 1759 | The polynomial q must be homogeneous of degree 2 and absolutely |
---|
| 1760 | irreducible. @* |
---|
| 1761 | NOTE: The procedure might fail or give a wrong output if the assumptions |
---|
| 1762 | do not hold. |
---|
| 1763 | |
---|
| 1764 | RETURN: ring with an ideal PARACONIC. The ring should be considered as the |
---|
| 1765 | homogeneous coordinate ring of PP^1, the ideal defines a rational |
---|
| 1766 | parametrization PP^1 --> C2 = {q=0}. |
---|
| 1767 | |
---|
[a212fe] | 1768 | THEORY: We compute a point on C2 via @ref{rationalPointConic}. The pencil of |
---|
[1a96a4] | 1769 | lines through this point projects C2 birationally to PP^1. Inverting |
---|
| 1770 | the projection gives the result. |
---|
| 1771 | KEYWORDS: conic, parametrization, rational point. |
---|
| 1772 | SEE ALSO: rationalPointConic. |
---|
| 1773 | EXAMPLE: example paraConic; shows an example |
---|
| 1774 | " |
---|
[6c265b] | 1775 | { |
---|
[1a96a4] | 1776 | def Roriginal = basering; |
---|
| 1777 | def RP2 = projConic(q); // ring with ideal Ipoint |
---|
| 1778 | // possibly defined over algebraic number field |
---|
| 1779 | setring RP2; |
---|
| 1780 | def rp1 = invertBirMap(Ipoint, ideal(fetch(Roriginal,q))); |
---|
| 1781 | setring rp1; |
---|
| 1782 | list rl = ringlist(rp1); |
---|
| 1783 | rl[2] = list("s","t"); |
---|
| 1784 | def RP1 = ring(rl); |
---|
| 1785 | setring RP1; |
---|
[6c265b] | 1786 | ideal PARACONIC = fetch(rp1,psi); |
---|
[1a96a4] | 1787 | export(PARACONIC); |
---|
| 1788 | "// 'paraConic' created a ring together with an ideal RNC."; |
---|
| 1789 | "// Supposing you typed, say, def RP1 = paraConic(q);"; |
---|
| 1790 | "// you may access the ideal by typing"; |
---|
| 1791 | "// setring RP1; PARACONIC;"; |
---|
| 1792 | return(RP1); |
---|
| 1793 | } |
---|
| 1794 | example |
---|
[75e82e] | 1795 | { "EXAMPLE:"; echo=2; |
---|
[1a96a4] | 1796 | ring R = 0,(x,y,z),dp; |
---|
| 1797 | poly f = y^8-x^3*(z+x)^5; |
---|
| 1798 | ideal adj = adjointIdeal(f); |
---|
| 1799 | def Rn = invertBirMap(adj,ideal(f)); |
---|
| 1800 | setring(Rn); |
---|
| 1801 | J; |
---|
| 1802 | def Rc = rncItProjEven(J); |
---|
| 1803 | PHI; |
---|
| 1804 | setring Rc; |
---|
| 1805 | CONIC; |
---|
| 1806 | def RPc = paraConic(CONIC); |
---|
| 1807 | setring RPc; |
---|
| 1808 | PARACONIC; |
---|
| 1809 | } |
---|
| 1810 | |
---|
| 1811 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1812 | static proc projConic(poly q) |
---|
| 1813 | "USAGE: projConic(q); q poly |
---|
[6c265b] | 1814 | ASSUME: The basering must be a polynomial ring in three variables with |
---|
[1a96a4] | 1815 | coefficients in Q. @* |
---|
| 1816 | The polynomial q must be homogeneous of degree 2 and absolutely |
---|
| 1817 | irreducible. @* |
---|
| 1818 | NOTE: The procedure might fail or give a wrong output if the assumptions |
---|
| 1819 | do not hold. |
---|
| 1820 | RETURN: ring with an ideal Ipoint defining a pencil of lines through a point |
---|
| 1821 | on the conic C2 = {q=0}. This point has either coefficients in Q or |
---|
| 1822 | in a quadratic extension field of Q. |
---|
| 1823 | THEORY: We compute the point on C2 via @ref rationalPointConic. |
---|
| 1824 | KEYWORDS: conic, parametrization, rational point. |
---|
| 1825 | SEE ALSO: rationalPointConic. |
---|
| 1826 | " |
---|
[6c265b] | 1827 | { |
---|
[1a96a4] | 1828 | def Roriginal = basering; |
---|
| 1829 | list rl = ringlist(Roriginal); |
---|
| 1830 | rl[3] = list(list("dp",1:3),list("C",0)); |
---|
| 1831 | def RP20 = ring(rl); |
---|
| 1832 | setring RP20; |
---|
| 1833 | poly q = imap(Roriginal,q); |
---|
| 1834 | def RP21 = rationalPointConic(q); // ring with ideal point representing |
---|
| 1835 | // point on conic |
---|
| 1836 | // possibly defined over algebraic number |
---|
| 1837 | // field |
---|
[6c265b] | 1838 | setring RP21; |
---|
[1a96a4] | 1839 | list rl1 = ringlist(RP21); |
---|
[6c265b] | 1840 | rl1[2] = list("u","v","w"); |
---|
[1a96a4] | 1841 | rl1[3] = list(list("dp",1:3),list("C",0)); |
---|
[6c265b] | 1842 | def RP2 = ring(rl1); |
---|
[1a96a4] | 1843 | setring RP2; |
---|
| 1844 | ideal point = fetch(RP21,point); |
---|
| 1845 | matrix bP = syz(point); |
---|
| 1846 | ideal Ipoint = matrix(maxideal(1))*bP; // defines pencil of lines through |
---|
| 1847 | // point |
---|
| 1848 | export(Ipoint); |
---|
| 1849 | return(RP2); |
---|
| 1850 | } |
---|
| 1851 | |
---|
| 1852 | |
---|
| 1853 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1854 | static proc isIrreducible(poly f) |
---|
| 1855 | "USAGE: isIrreducible(f); f poly |
---|
| 1856 | RETURN: 1 iff the given polynomial f is irreducible; 0 otherwise. |
---|
| 1857 | THEORY: This test is performed by computing the absolute factorization of f. |
---|
| 1858 | KEYWORDS: irreducible. |
---|
| 1859 | " |
---|
| 1860 | { |
---|
| 1861 | def r = basering; |
---|
| 1862 | def s = absFactorize(f); |
---|
| 1863 | setring s; |
---|
| 1864 | list L = absolute_factors; |
---|
| 1865 | int result = 0; |
---|
| 1866 | if (L[4] == 1){result = 1;} |
---|
| 1867 | setring r; |
---|
| 1868 | kill s; |
---|
| 1869 | return (result); |
---|
| 1870 | } |
---|
| 1871 | |
---|
| 1872 | |
---|
| 1873 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1874 | static proc isQuadratic(poly p) |
---|
| 1875 | "USAGE: isQuadratic(p); p poly |
---|
| 1876 | RETURN: checks whether p is a homogeneous, quadratic polynomial in the |
---|
| 1877 | first three ring variables, x, y, z say; |
---|
| 1878 | If so, the method extracs the coefficients a, b, ..., f such that |
---|
| 1879 | p = a*x2 + b*xy + c * y2 + d * xz + e * yz + f * z2 |
---|
| 1880 | and returns them as a list of seven entries, [1, a, b, c, d, e, f]; |
---|
| 1881 | otherwise, a list with the single entry [0] is returned |
---|
| 1882 | " |
---|
| 1883 | { |
---|
| 1884 | bigint a = bigint(leadcoef(subst(p, var(2), 0, var(3), 0))); |
---|
| 1885 | bigint c = bigint(leadcoef(subst(p, var(1), 0, var(3), 0))); |
---|
| 1886 | bigint f = bigint(leadcoef(subst(p, var(1), 0, var(2), 0))); |
---|
| 1887 | poly h = p - a * var(1)^2 - c * var(2)^2 - f * var(3)^2; |
---|
| 1888 | bigint b = bigint(leadcoef(subst(h, var(3), 0))); |
---|
| 1889 | bigint d = bigint(leadcoef(subst(h, var(2), 0))); |
---|
| 1890 | bigint e = bigint(leadcoef(subst(h, var(1), 0))); |
---|
| 1891 | list L = 0; |
---|
| 1892 | if (h - b * var(1) * var(2) - d * var(1) * var(3) |
---|
| 1893 | - e * var(2) * var(3) != 0) { return (L); } |
---|
| 1894 | L = 1, a, b, c, d, e, f; |
---|
| 1895 | return (L); |
---|
| 1896 | } |
---|
| 1897 | |
---|
| 1898 | |
---|
| 1899 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1900 | static proc largestSquare(bigint n) |
---|
| 1901 | "USAGE: largestSquare(n); n bigint |
---|
| 1902 | ASSUME: n <> 0 |
---|
| 1903 | RETURN: returns the largest positive number m (as bigint) such that m^2 |
---|
| 1904 | divides n. |
---|
| 1905 | THEORY: This computation is done by prime factorization of n. |
---|
| 1906 | KEYWORDS: prime factorization. |
---|
| 1907 | " |
---|
| 1908 | { |
---|
| 1909 | if (n == 0) { "ERROR: largestSquare(0) had been invoked"; } |
---|
| 1910 | |
---|
| 1911 | bigint nn = n; if (nn < 0) { nn = -n; } |
---|
| 1912 | list L = primefactors(nn); |
---|
[797a9ff] | 1913 | if (L[3] != 1) |
---|
[1a96a4] | 1914 | { "WARNING: command 'primefactors(.)' did not find all prime factors"; } |
---|
| 1915 | int i; bigint m = bigint(1); int e; int j; |
---|
[797a9ff] | 1916 | for (i = 1; i <= size(L[1]); i++) |
---|
[1a96a4] | 1917 | { |
---|
[62c2b0] | 1918 | e = L[2][i] div 2; |
---|
[797a9ff] | 1919 | for (j = 1; j <= e; j++) { m = m * bigint(L[1][i]); } |
---|
[1a96a4] | 1920 | } |
---|
| 1921 | return (m); |
---|
| 1922 | } |
---|
| 1923 | |
---|
| 1924 | |
---|
| 1925 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1926 | static proc jIndex(bigint a, bigint b, bigint c) |
---|
| 1927 | "USAGE: jIndex(a, b, c); a, b, c bigint's |
---|
| 1928 | RETURN: returns the middle of the three numbers |ab|, |bc|, and |ca|. |
---|
| 1929 | " |
---|
| 1930 | { |
---|
| 1931 | bigint n1 = a*b; if (n1 < 0) { n1 = -n1; } |
---|
| 1932 | bigint n2 = b*c; if (n2 < 0) { n2 = -n2; } |
---|
| 1933 | bigint n3 = c*a; if (n3 < 0) { n3 = -n3; } |
---|
| 1934 | if ((n1 <= n2) && (n2 <= n3)) { return (n2); } |
---|
| 1935 | if ((n1 >= n2) && (n2 >= n3)) { return (n2); } |
---|
| 1936 | if ((n2 <= n1) && (n1 <= n3)) { return (n1); } |
---|
| 1937 | if ((n2 >= n1) && (n1 >= n3)) { return (n1); } |
---|
| 1938 | return (n3); |
---|
| 1939 | } |
---|
| 1940 | |
---|
| 1941 | |
---|
[e8eba7] | 1942 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1943 | static proc aMod(bigint a, bigint b) |
---|
| 1944 | "USAGE: aMod(a,b); a, b bigint |
---|
| 1945 | RETURN: r bigint |
---|
| 1946 | THEORY: The asymmetric residue r of the division with remainder a mod b. |
---|
| 1947 | " |
---|
| 1948 | { |
---|
| 1949 | bigint c = a mod b; |
---|
| 1950 | if (c<0) |
---|
| 1951 | { |
---|
| 1952 | return(c+b); |
---|
| 1953 | } |
---|
| 1954 | return(c); |
---|
| 1955 | } |
---|
| 1956 | |
---|
| 1957 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1958 | static proc aDiv(bigint a, bigint b) |
---|
| 1959 | "USAGE: aDiv(a,b); a, b bigint |
---|
| 1960 | RETURN: q bigint |
---|
| 1961 | THEORY: Quotient with remainder q = a div b with asymmetric residue. |
---|
| 1962 | " |
---|
| 1963 | { |
---|
| 1964 | bigint q = a div b; |
---|
| 1965 | if ((a mod b)<0) |
---|
| 1966 | { |
---|
| 1967 | return(q-1); |
---|
| 1968 | } |
---|
| 1969 | return(q); |
---|
| 1970 | } |
---|
| 1971 | |
---|
| 1972 | |
---|
| 1973 | |
---|
[1a96a4] | 1974 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1975 | static proc polyModP(poly q, bigint p) |
---|
| 1976 | "USAGE: polyModP(q, p); q poly, p bigint |
---|
| 1977 | RETURN: takes each coefficient of q modulo p and returns the resulting poly |
---|
| 1978 | " |
---|
| 1979 | { |
---|
| 1980 | poly qq = q; poly res = 0; |
---|
| 1981 | bigint c; |
---|
| 1982 | while (qq != 0) |
---|
| 1983 | { |
---|
| 1984 | c = bigint(leadcoef(qq) mod p); |
---|
| 1985 | res = res + c * leadmonom(qq); |
---|
| 1986 | qq = qq - lead(qq); |
---|
| 1987 | } |
---|
| 1988 | return (res); |
---|
| 1989 | } |
---|
| 1990 | |
---|
| 1991 | |
---|
| 1992 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1993 | static proc rootModP(bigint r, bigint p) |
---|
| 1994 | "USAGE: rootModP(r, p); r, p bigint's |
---|
| 1995 | ASSUME: 0 <= r < p, and p prime; |
---|
| 1996 | Furthermore it is assumes that there is some x in {0, 1, ..., p-1} |
---|
| 1997 | such that x^2 = r mod p; |
---|
| 1998 | RETURN: an x in {0, 1, ..., p-1} such that x^2 = r mod p; |
---|
| 1999 | THEORY: For p larger than 32003, this computation is done using Cantor- |
---|
| 2000 | Zassenhaus' algorithm. Otherwise a brute force approach is used. |
---|
| 2001 | KEYWORDS: Cantor-Zassenhaus algorithm. |
---|
| 2002 | " |
---|
| 2003 | { |
---|
| 2004 | if (r == 0) { return (0); } |
---|
| 2005 | if (r == 1) { return (1); } |
---|
| 2006 | if (p <= 32003) |
---|
| 2007 | { |
---|
| 2008 | /* For small p, we use a brute force approach: */ |
---|
| 2009 | int i; |
---|
[797a9ff] | 2010 | for (i = 2; i < p; i++) |
---|
[1a96a4] | 2011 | { |
---|
| 2012 | if (((i*i) mod p) == r) { return (i); } |
---|
| 2013 | } |
---|
| 2014 | /* should never be reached: */ |
---|
| 2015 | return (-1); |
---|
| 2016 | } |
---|
| 2017 | |
---|
| 2018 | /* For p > 32003, we use Cantor-Zassenhaus' algorithm: */ |
---|
| 2019 | def br = basering; |
---|
| 2020 | ring rTemp = 0, x, dp; |
---|
| 2021 | bigint b; bigint exponent; poly factor; |
---|
| 2022 | poly h = x^2 - r; |
---|
| 2023 | ideal redI = h; redI = std(redI); |
---|
| 2024 | poly q = x^2; bigint root = 0; |
---|
| 2025 | while (root == 0) |
---|
| 2026 | { |
---|
| 2027 | b = bigint(random(1, 2^30)); |
---|
[62c2b0] | 2028 | exponent = bigint((p - 1) div 2); |
---|
[1a96a4] | 2029 | /* We need to compute q^exponent mod (x^2 - a) and mod p: */ |
---|
| 2030 | factor = x + b; q = 1; |
---|
| 2031 | while (exponent > 0) |
---|
| 2032 | { |
---|
| 2033 | if ((exponent mod 2) == 1) |
---|
| 2034 | { |
---|
| 2035 | q = q * factor; |
---|
| 2036 | q = reduce(q, redI); |
---|
| 2037 | q = polyModP(q, p); |
---|
| 2038 | } |
---|
[62c2b0] | 2039 | exponent = bigint(exponent div 2); |
---|
[1a96a4] | 2040 | factor = factor * factor; |
---|
| 2041 | factor = reduce(factor, redI); |
---|
| 2042 | factor = polyModP(factor, p); |
---|
| 2043 | } |
---|
| 2044 | if (deg(q) == 1) |
---|
| 2045 | { |
---|
| 2046 | q = q - 1; |
---|
| 2047 | b = inverseModP(bigint(leadcoef(q)), p); |
---|
| 2048 | q = q - lead(q); |
---|
[e8eba7] | 2049 | root = aMod((bigint(q) * b),p); |
---|
[1a96a4] | 2050 | if (((root * root - r) mod p) != 0) { root = 0; } |
---|
| 2051 | } |
---|
| 2052 | } |
---|
| 2053 | setring br; kill rTemp; |
---|
| 2054 | return (root); |
---|
| 2055 | } |
---|
| 2056 | |
---|
| 2057 | |
---|
| 2058 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2059 | static proc inverseModP(bigint r, bigint p) |
---|
| 2060 | "USAGE: inverseModP(r, p); r, p bigint's |
---|
| 2061 | ASSUME: 0 <= r < p, and r and p coprime; |
---|
| 2062 | RETURN: returns the inverse of r in Z/p represented by an element in |
---|
| 2063 | {1, 2, ..., p-1} |
---|
| 2064 | THEORY: This uses Euclid's extended gcd algorithm. |
---|
| 2065 | " |
---|
| 2066 | { |
---|
| 2067 | list L = extendedEuclid(r, p); |
---|
| 2068 | if (L[1] != 1) { "ERROR: GCD of", r, "and", p, "should be 1."; } |
---|
[e8eba7] | 2069 | L[2] = aMod(L[2],p); |
---|
[1a96a4] | 2070 | return (L[2]); |
---|
| 2071 | } |
---|
| 2072 | |
---|
| 2073 | |
---|
| 2074 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2075 | static proc squareRoot(bigint r, bigint m, int justCheck) |
---|
| 2076 | "USAGE: squareRoot(r, m, j); r, m bigint's, j int |
---|
| 2077 | RETURN: checks whether r is a square modulo m, i.e., checks whether there is |
---|
| 2078 | some x such that x^2 = r mod m; |
---|
| 2079 | If justCheck is 1, then the method will terminate after the check |
---|
| 2080 | and return 1 if r is a square and -1 otherwise. |
---|
| 2081 | If justCheck is 0 and r is a square, then the method continues and |
---|
| 2082 | computes a solution x in {0, 1, m-1} with x^2 = r mod m, which will |
---|
| 2083 | then be returned |
---|
| 2084 | THEORY: This algorithm checks solvability by computing the Legendre symbols |
---|
| 2085 | modulo all primes in m. Then, individual roots will be computed and |
---|
| 2086 | lifted to the desired square root modulo m using Chinese |
---|
| 2087 | remaindering. |
---|
| 2088 | " |
---|
| 2089 | { |
---|
| 2090 | if (m == 0) { "ERROR: squareRoot had been invoked with m = 0"; } |
---|
| 2091 | |
---|
| 2092 | list L = primefactors(m); |
---|
[2a33e6] | 2093 | if ((L[3] != 1) && (L[3] != -1)) |
---|
| 2094 | { "WARNING: command 'primefactors(.)' did not find all prime factors"; } |
---|
[1a96a4] | 2095 | int i; |
---|
[797a9ff] | 2096 | for (i = 1; i <= size(L[2]); i++) |
---|
[1a96a4] | 2097 | { |
---|
[797a9ff] | 2098 | if (legendreSymbol(r, L[1][i]) == -1) { return (-1); } |
---|
[1a96a4] | 2099 | } |
---|
| 2100 | /* now we know that there is some x in {0, 1, m-1} with |
---|
| 2101 | x^2 = r mod m */ |
---|
| 2102 | if (justCheck == 1) { return (1); } |
---|
| 2103 | else |
---|
| 2104 | { |
---|
| 2105 | // now we need to compute x; this works in two stages: |
---|
| 2106 | // 1) write m = p1^e1 * ... * pk^ek (we already have that), |
---|
| 2107 | // 2) for each i in {1, 2, ..., k} |
---|
| 2108 | // 2.1) compute a yi such that yi^2 = r mod pi, |
---|
| 2109 | // 2.2) lift yi to an xi such that xi^2 = r mod (pi^ei), |
---|
| 2110 | // 3) lift (x1, x2, ..., xk) in Z/p1^e1 * ... * Z/pk^ek |
---|
| 2111 | // to x in Z/m via Chinese remainder theorem |
---|
| 2112 | |
---|
| 2113 | list roots; |
---|
| 2114 | // 2.1): |
---|
[797a9ff] | 2115 | for (i = 1; i <= size(L[1]); i++) |
---|
[1a96a4] | 2116 | { |
---|
[e8eba7] | 2117 | roots = insert(roots, rootModP(aMod(r,L[1][i]), L[1][i]), size(roots)); |
---|
[1a96a4] | 2118 | } |
---|
| 2119 | |
---|
| 2120 | // 2.2): |
---|
| 2121 | bigint c; bigint l; bigint temp; bigint pPower; int e; |
---|
[797a9ff] | 2122 | for (i = 1; i <= size(roots); i++) |
---|
[1a96a4] | 2123 | { |
---|
[797a9ff] | 2124 | pPower = bigint(L[1][i]); |
---|
| 2125 | for (e = 2; e <= L[2][i]; e++) |
---|
[1a96a4] | 2126 | { |
---|
| 2127 | c = bigint(roots[i]); l = pPower; |
---|
| 2128 | temp = r - c * c; l = bigint(2) * c * l; c = temp; |
---|
[e8eba7] | 2129 | c = aDiv(c,pPower); l = aDiv(l,pPower); |
---|
| 2130 | c = aMod(c,L[1][i]); l = aMod(l,L[1][i]); |
---|
| 2131 | c = aMod((c * bigint(inverseModP(l, L[1][i]))), L[1][i]); |
---|
[1a96a4] | 2132 | c = bigint(roots[i]) + c * pPower; |
---|
[797a9ff] | 2133 | pPower = pPower * L[1][i]; roots[i] = c; |
---|
[1a96a4] | 2134 | } |
---|
| 2135 | } |
---|
| 2136 | |
---|
| 2137 | // 2.3): |
---|
| 2138 | list mm; bigint z; int j; |
---|
[797a9ff] | 2139 | for (i = 1; i <= size(L[1]); i++) |
---|
[1a96a4] | 2140 | { |
---|
[797a9ff] | 2141 | z = bigint(L[1][i]); |
---|
| 2142 | for (j = 2; j <= L[2][i]; j++) |
---|
[1a96a4] | 2143 | { |
---|
[797a9ff] | 2144 | z = z * bigint(L[1][i]); |
---|
[1a96a4] | 2145 | } |
---|
| 2146 | mm = insert(mm, z, size(mm)); |
---|
| 2147 | } |
---|
[e8eba7] | 2148 | return (aMod(chineseRemainder(roots, mm) , m)); |
---|
[1a96a4] | 2149 | } |
---|
| 2150 | } |
---|
| 2151 | |
---|
| 2152 | |
---|
| 2153 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2154 | static proc chineseRemainder(list rr, list mm) |
---|
| 2155 | "USAGE: chineseRemainder(rr, mm); rr, mm lists of bigint's |
---|
| 2156 | ASSUME: lists rr and mm must have same sizes; |
---|
| 2157 | Furthermore the entries of mm must be mutually coprime. |
---|
| 2158 | RETURN: an x which fulfills the simultaneous remainder conditions |
---|
| 2159 | x = rr[i] mod mm[i], 1 <= i <= size(rr) |
---|
| 2160 | KEYWORDS: Chinese remainder. |
---|
| 2161 | " |
---|
| 2162 | { |
---|
| 2163 | bigint x = bigint(0); int i; bigint N; list l; |
---|
| 2164 | bigint M = bigint(mm[1]); |
---|
[797a9ff] | 2165 | for (i = 2; i <= size(mm); i++) { M = M * bigint(mm[i]); } |
---|
| 2166 | for (i = 1; i <= size(mm); i++) |
---|
[1a96a4] | 2167 | { |
---|
[e8eba7] | 2168 | N = aDiv(M,mm[i]); |
---|
[1a96a4] | 2169 | l = extendedEuclid(mm[i], N); |
---|
| 2170 | x = x + rr[i]*l[3]*N; |
---|
| 2171 | } |
---|
| 2172 | return (x); |
---|
| 2173 | } |
---|
| 2174 | |
---|
| 2175 | |
---|
| 2176 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2177 | static proc rationalPointSpecial(bigint b1, bigint c1) |
---|
| 2178 | "USAGE: rationalPointSpecial(b1, c1); b1, c1 bigint's |
---|
| 2179 | ASSUME: b1 <> 0 and c1 <> 0; |
---|
| 2180 | RETURN: with poly p = var(1)^2 + b1 * var(2)^2 + c1 * var(3)^2, the method |
---|
| 2181 | returns a list L with either one entry or four entries: |
---|
| 2182 | case 'three entries': |
---|
| 2183 | L[1] = 0 signaling that there is no rational point on V(p), |
---|
| 2184 | L[2] the largest number b such that b^2 divides b1 |
---|
| 2185 | (for subsequent use by the caller of this method), |
---|
| 2186 | L[3] the largest number c such that c^2 divides c1 |
---|
| 2187 | (for subsequent use by the caller of this method); |
---|
| 2188 | case 'four entries': |
---|
| 2189 | L[1] = 1 signaling that there is a rational point on V(p), |
---|
| 2190 | L[2], L[3], L[4] rational numbers such that the tuple |
---|
| 2191 | (L[2], L[3], L[4]) is on V(p) |
---|
| 2192 | " |
---|
| 2193 | { |
---|
| 2194 | if (b1 == 0) { "ERROR: rationalPointSpecial(0, c1) had been invoked"; } |
---|
| 2195 | if (c1 == 0) { "ERROR: rationalPointSpecial(b1, 0) had been invoked"; } |
---|
| 2196 | |
---|
| 2197 | bigint b_s = largestSquare(b1); bigint b_r = b1/b_s/b_s; |
---|
| 2198 | bigint c_s = largestSquare(c1); bigint c_r = c1/c_s/c_s; |
---|
| 2199 | bigint g = gcd(b_r, c_r); |
---|
[74270f] | 2200 | def S=basering; |
---|
[1a96a4] | 2201 | ideal mi = maxideal(1); |
---|
| 2202 | map mm = basering, mi; map mTemp; |
---|
| 2203 | mm[1] = var(1); mm[2] = var(2)/b_s/g; mm[3] = var(3)/c_s/g; |
---|
| 2204 | bigint a = g; bigint aa = a; if (aa <= 0) { aa = -aa; } |
---|
| 2205 | bigint b = b_r/g; bigint bb = b; if (bb <= 0) { bb = -bb; } |
---|
| 2206 | bigint c = c_r/g; bigint cc = c; if (cc <= 0) { cc = -cc; } |
---|
| 2207 | bigint R1 = squareRoot(-a*b, cc, 1); |
---|
| 2208 | if (R1 == -1) { list L = 0, b_s, c_s; return (L); } |
---|
| 2209 | bigint R2 = squareRoot(-a*c, bb, 1); |
---|
| 2210 | if (R2 == -1) { list L = 0, b_s, c_s; return (L); } |
---|
| 2211 | bigint R3 = squareRoot(-b*c, aa, 1); |
---|
| 2212 | if (R3 == -1) { list L = 0, b_s, c_s; return (L); } |
---|
| 2213 | bigint t; bigint r1; bigint Q; bigint A; bigint B; bigint C; |
---|
| 2214 | bigint alpha; bigint beta; bigint gamma; |
---|
| 2215 | while (jIndex(a, b, c) > 1) |
---|
| 2216 | { |
---|
| 2217 | mTemp = basering, mi; |
---|
| 2218 | if (aa > cc) |
---|
| 2219 | { |
---|
| 2220 | t = a; a = c; c = t; |
---|
| 2221 | t = aa; aa = cc; cc = t; |
---|
| 2222 | mTemp = basering, mi; |
---|
| 2223 | mTemp[1] = var(3); mTemp[3] = var(1); mm = mTemp(mm); |
---|
| 2224 | } |
---|
| 2225 | if (bb > cc) |
---|
| 2226 | { |
---|
| 2227 | t = b; b = c; c = t; |
---|
| 2228 | t = bb; bb = cc; cc = t; |
---|
| 2229 | mTemp = basering, mi; |
---|
| 2230 | mTemp[2] = var(3); mTemp[3] = var(2); mm = mTemp(mm); |
---|
| 2231 | } |
---|
| 2232 | if (bb < aa) |
---|
| 2233 | { |
---|
| 2234 | t = b; b = a; a = t; |
---|
| 2235 | t = bb; bb = aa; aa = t; |
---|
| 2236 | mTemp = basering, mi; |
---|
| 2237 | mTemp[1] = var(2); mTemp[2] = var(1); mm = mTemp(mm); |
---|
| 2238 | } |
---|
| 2239 | /* now, we have established |a| <= |b| <= |c|; and permuted |
---|
| 2240 | the map mm, accordingly */ |
---|
| 2241 | cc = c; if (cc <= 0) { cc = -cc; } |
---|
| 2242 | R1 = squareRoot(-a*b, cc, 0); |
---|
[e8eba7] | 2243 | r1 = aMod((R1 * inverseModP(a, cc)), cc); |
---|
[1a96a4] | 2244 | if (r1*bigint(2) > cc) { r1 = r1 - cc; } |
---|
| 2245 | Q = (a*r1*r1 + b)/c; |
---|
| 2246 | if (Q == 0) |
---|
| 2247 | { |
---|
| 2248 | list L = 1, subst(mm[1], var(1), 1, var(2), 1, var(3), 0), |
---|
| 2249 | subst(mm[2], var(1), 1, var(2), 1, var(3), 0), |
---|
| 2250 | subst(mm[3], var(1), 1, var(2), 1, var(3), 0); |
---|
| 2251 | return (L); |
---|
| 2252 | } |
---|
| 2253 | A = gcd(gcd(a*r1*r1, b), c*Q); |
---|
| 2254 | alpha = r1/A; beta = b/A; |
---|
| 2255 | B = a*beta; |
---|
| 2256 | gamma = largestSquare(Q/A); |
---|
| 2257 | C = Q/A/gamma/gamma; |
---|
| 2258 | mTemp = basering, mi; |
---|
| 2259 | mTemp[1] = A*alpha*var(1) - beta*var(2); |
---|
| 2260 | mTemp[2] = var(1) + a*alpha*var(2); |
---|
| 2261 | mTemp[3] = C*gamma*var(3); |
---|
| 2262 | mm = mTemp(mm); |
---|
| 2263 | a = A; b = B; c = C; |
---|
| 2264 | aa = a; if (aa <= 0) { aa = -aa; } |
---|
| 2265 | bb = b; if (bb <= 0) { bb = -bb; } |
---|
| 2266 | cc = c; if (cc <= 0) { cc = -cc; } |
---|
| 2267 | } |
---|
| 2268 | if (a*b < 0) |
---|
| 2269 | { |
---|
| 2270 | list L = 1, subst(mm[1], var(1), 1, var(2), 1, var(3), 0), |
---|
| 2271 | subst(mm[2], var(1), 1, var(2), 1, var(3), 0), |
---|
| 2272 | subst(mm[3], var(1), 1, var(2), 1, var(3), 0); |
---|
| 2273 | return (L); |
---|
| 2274 | } |
---|
| 2275 | if (a*c < 0) |
---|
| 2276 | { |
---|
| 2277 | list L = 1, subst(mm[1], var(1), 1, var(2), 0, var(3), 1), |
---|
| 2278 | subst(mm[2], var(1), 1, var(2), 0, var(3), 1), |
---|
| 2279 | subst(mm[3], var(1), 1, var(2), 0, var(3), 1); |
---|
| 2280 | return (L); |
---|
| 2281 | } |
---|
| 2282 | if (b*c < 0) |
---|
| 2283 | { |
---|
| 2284 | list L = 1, subst(mm[1], var(1), 0, var(2), 1, var(3), 1), |
---|
| 2285 | subst(mm[2], var(1), 0, var(2), 1, var(3), 1), |
---|
| 2286 | subst(mm[3], var(1), 0, var(2), 1, var(3), 1); |
---|
| 2287 | return (L); |
---|
| 2288 | } |
---|
| 2289 | list L = 0, b_s, c_s; return (L); |
---|
| 2290 | } |
---|
| 2291 | |
---|
| 2292 | |
---|
| 2293 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2294 | static proc extendedEuclid(bigint a, bigint b) |
---|
| 2295 | "USAGE: extendedEuclid(a, b); a, b bigint's |
---|
| 2296 | ASSUME: a <> 0 or b <> 0; |
---|
| 2297 | RETURN: returns a list with three entries: |
---|
| 2298 | _[1]: gcd(a,b) > 0, |
---|
| 2299 | _[2], _[3]: s, t, such that s*a + t*b = gcd(a,b) |
---|
| 2300 | KEYWORDS: extended Euclidean algorithm. |
---|
| 2301 | " |
---|
| 2302 | { |
---|
| 2303 | list l = 0; bigint temp; |
---|
[62dc18e] | 2304 | if (a == 0) { l = b, 0, 1; if (b < 0) { l = -b, 0, -1; } } |
---|
| 2305 | if (b == 0) { l = a, 1, 0; if (a < 0) { l = -a, -1, 0; } } |
---|
[e8eba7] | 2306 | if (aMod(a , b) == 0) { l = b, 0, 1; if (b < 0) { l = -b, 0, -1; } } |
---|
| 2307 | if (aMod(b , a) == 0) { l = a, 1, 0; if (a < 0) { l = -a, -1, 0; } } |
---|
[1a96a4] | 2308 | if (size(l) > 1) { return (l); } |
---|
[6c265b] | 2309 | |
---|
[e8eba7] | 2310 | temp = aMod(a , b); |
---|
[1a96a4] | 2311 | l = extendedEuclid(b, temp); |
---|
| 2312 | temp = (a - temp) / b; |
---|
| 2313 | temp = bigint(l[2]) - temp * bigint(l[3]); |
---|
| 2314 | l = l[1], l[3], temp; |
---|
| 2315 | return (l); |
---|
| 2316 | } |
---|
| 2317 | |
---|
| 2318 | static proc legendreSymbol(bigint r, bigint p) |
---|
| 2319 | "assumes p prime; |
---|
| 2320 | returns the Legendre symbol (r/p), that is |
---|
| 2321 | 1 if r appears as residue modulo p of a square, |
---|
| 2322 | -1 if not, |
---|
| 2323 | 0 if r is a multiple of p |
---|
| 2324 | " |
---|
| 2325 | { |
---|
[e8eba7] | 2326 | bigint rr = aMod(r , p); |
---|
[62dc18e] | 2327 | if (rr == 0) { return (0) } |
---|
| 2328 | if (rr == 1) { return (1) } |
---|
[1a96a4] | 2329 | /* now, p must be at least 3 */ |
---|
| 2330 | bigint e = (p - 1) / bigint(2); |
---|
| 2331 | bigint result = 1; |
---|
| 2332 | bigint power = rr; |
---|
| 2333 | while (e > 0) |
---|
| 2334 | { |
---|
[e8eba7] | 2335 | if ((e mod 2) == 1) { result = aMod((result * power), p); } |
---|
[1a96a4] | 2336 | e = e / bigint(2); |
---|
[e8eba7] | 2337 | power = aMod((power * power), p); |
---|
[1a96a4] | 2338 | } |
---|
| 2339 | if (result > 1) { result = result - p; /* should be -1 */ } |
---|
| 2340 | return (result); |
---|
| 2341 | } |
---|
| 2342 | |
---|
| 2343 | |
---|
| 2344 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2345 | static proc buildExtension(bigint b, bigint c, bigint bs, bigint cs) |
---|
| 2346 | "USAGE: buildExtension(b, c, bs, cs); b, c, bs, cs bigint's |
---|
| 2347 | ASSUME: assumes that bs is the largest positive number such that bs^2 |
---|
| 2348 | divides b; analogously for cs regarding c; |
---|
| 2349 | Assumes furthermore that there is no rational point on the conic |
---|
| 2350 | X^2 + b*Y^2 + c*Z^2 = 0. |
---|
| 2351 | Assumes that the ground field of the basering is Q. |
---|
| 2352 | RETURN: builds an appropriate quadratic field extension Q(a) in which a |
---|
| 2353 | point exists that lies on the given conic. This point is stored in |
---|
| 2354 | a newly defined and exported (1x3) matrix named 'point'. |
---|
| 2355 | The method returns the resulting polynomial ring over Q(a). |
---|
| 2356 | " |
---|
| 2357 | { |
---|
| 2358 | bigint br = b/bs/bs; |
---|
| 2359 | bigint cr = c/cs/cs; |
---|
| 2360 | /* X^2 + br*bs^2*Y^2 + cr*cs^2*Z^2 = 0 */ |
---|
| 2361 | def bRing = basering; |
---|
| 2362 | list L = ringlist(bRing); |
---|
| 2363 | |
---|
| 2364 | if (b != 0) |
---|
| 2365 | { |
---|
| 2366 | L[1] = list(0, list("a"), list(list("lp", 1)), ideal(0)); |
---|
| 2367 | def RTemp = ring(L); |
---|
| 2368 | setring RTemp; list L = ringlist(RTemp); |
---|
| 2369 | L[1][4] = ideal(a^2 + br); |
---|
| 2370 | def R = ring(L); |
---|
| 2371 | setring R; kill RTemp; |
---|
| 2372 | matrix point[1][3]; |
---|
| 2373 | point[1, 1] = a * bs; point[1, 2] = 1; point[1, 3] = 0; |
---|
| 2374 | export point; |
---|
| 2375 | setring bRing; |
---|
| 2376 | return (R); |
---|
| 2377 | } |
---|
| 2378 | if (c != 0) |
---|
| 2379 | { |
---|
| 2380 | L[1] = list(0, list("a"), list(list("lp", 1)), ideal(0)); |
---|
| 2381 | def RTemp = ring(L); |
---|
| 2382 | setring RTemp; list L = ringlist(RTemp); |
---|
| 2383 | L[1][4] = ideal(a^2 + cr); |
---|
| 2384 | def R = ring(L); |
---|
| 2385 | setring R; kill RTemp; |
---|
| 2386 | matrix point[1][3]; |
---|
| 2387 | point[1, 1] = a * cs; point[1, 2] = 0; point[1, 3] = 1; |
---|
| 2388 | export point; |
---|
| 2389 | setring bRing; |
---|
| 2390 | return (R); |
---|
| 2391 | } |
---|
| 2392 | |
---|
| 2393 | "ERROR: unexpectedly encountered conic X^2 + 0*Y^2 + 0*Z^2 = 0"; |
---|
| 2394 | return (bRing); |
---|
| 2395 | } |
---|
| 2396 | |
---|
| 2397 | |
---|
| 2398 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2399 | static proc testRationalPointConic(poly pp) |
---|
| 2400 | "USAGE: testRationalPointConic(pp); pp poly |
---|
| 2401 | RETURN: returns 0 in case of unexpected input (e.g. non-quadratic, |
---|
| 2402 | reducible); 1 otherwise |
---|
| 2403 | NOTE: This method calles rationalPointConic, measures time consumption |
---|
| 2404 | and checks whether the computed point lies indeed on the conic pp. |
---|
| 2405 | The results will be printed to standard output. |
---|
| 2406 | " |
---|
| 2407 | { |
---|
| 2408 | "testing rationalPointConic(poly) for poly p:"; |
---|
| 2409 | "p =", pp; |
---|
| 2410 | if (isQuadratic(pp)[1] == 1) { "p is quadratic."; } |
---|
| 2411 | else { "p is not quadratic."; return (0); } |
---|
| 2412 | if (isIrreducible(pp) == 1) { "p is irreducible."; } |
---|
| 2413 | else { "p is not irreducible."; return (0); } |
---|
| 2414 | def rOrig = basering; |
---|
| 2415 | int t = rtimer; |
---|
| 2416 | def rNew = rationalPointConic(pp); |
---|
| 2417 | t = rtimer - t; |
---|
[22d1b1b] | 2418 | "time for finding a point on the conic [sec] =", t; |
---|
[1a96a4] | 2419 | setring rNew; |
---|
| 2420 | poly ff = fetch(rOrig, pp); |
---|
| 2421 | if (minpoly == 0) |
---|
| 2422 | { "there is a rational point on the conic p"; |
---|
| 2423 | "x =", point[1,1], " y =", point[1,2], " z =", point[1,3]; |
---|
| 2424 | "check (should be zero):", subst(ff, var(1), point[1,1], |
---|
| 2425 | var(2), point[1,2], |
---|
| 2426 | var(3), point[1,3]); |
---|
| 2427 | } |
---|
| 2428 | else |
---|
| 2429 | { |
---|
| 2430 | "there is no rational point on the conic p"; |
---|
| 2431 | "but there is a point on the conic in the field extension Q(a),"; |
---|
| 2432 | "with minpoly =", minpoly; |
---|
| 2433 | "x =", point[1,1], " y =", point[1,2], " z =", point[1,3]; |
---|
| 2434 | "check (should be zero):", subst(ff, var(1), point[1,1], |
---|
| 2435 | var(2), point[1,2], |
---|
| 2436 | var(3), point[1,3]); |
---|
| 2437 | } |
---|
| 2438 | setring rOrig; |
---|
| 2439 | } |
---|
| 2440 | |
---|
| 2441 | example |
---|
[75e82e] | 2442 | { "EXAMPLE:"; echo=2; |
---|
[1a96a4] | 2443 | ring r = 0, (x,y,z, u, v, w), dp; |
---|
| 2444 | poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z; |
---|
| 2445 | testRationalPointConic(p); |
---|
| 2446 | } |
---|
| 2447 | |
---|
| 2448 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2449 | proc rationalPointConic(poly p) |
---|
| 2450 | "USAGE: rationalPointConic(p); p poly |
---|
| 2451 | ASSUME: assumes that p is an irreducible quadratic polynomial in the first |
---|
| 2452 | three ring variables; |
---|
[a212fe] | 2453 | ground field is expected to be Q. |
---|
[1a96a4] | 2454 | RETURN: The method finds a point on the given conic. There are two |
---|
| 2455 | possibilities: |
---|
| 2456 | 1) There is a rational point on the curve. |
---|
| 2457 | 2) There is no rational point on the curve. |
---|
| 2458 | In the second case, the method creates a modification of the current |
---|
| 2459 | basering which is a polynomial ring over some quadratic field |
---|
| 2460 | extension Q(a) of Q. Apart from the replacement of Q by Q(a), the |
---|
| 2461 | new polynomial ring, R say, is the same as the original basering. |
---|
| 2462 | (In the first case, R is identical with the basering.) |
---|
| 2463 | In both cases, the method will then define a (1x3) matrix named |
---|
| 2464 | 'point' which lives in R and which contains the coordinates of the |
---|
| 2465 | desired point on q. |
---|
| 2466 | Finally, the method returns the ring R (which will in the 1st case |
---|
| 2467 | be the original base ring). |
---|
| 2468 | EXAMPLE: example rationalPointConic; shows an example |
---|
| 2469 | " |
---|
| 2470 | { |
---|
| 2471 | list L = isQuadratic(p); |
---|
| 2472 | bigint a = bigint(L[2]); bigint b = bigint(L[3]); bigint c = bigint(L[4]); |
---|
| 2473 | bigint d = bigint(L[5]); bigint e = bigint(L[6]); bigint f = bigint(L[7]); |
---|
| 2474 | bigint x; bigint y; bigint z; bigint nn; |
---|
| 2475 | def R = basering; |
---|
| 2476 | |
---|
| 2477 | if (b^2 == 4*a*c) |
---|
| 2478 | { |
---|
| 2479 | if (c == 0) |
---|
| 2480 | { |
---|
| 2481 | x = -2*d*e; y = d^2-4*a*f; z = e*4*a; |
---|
| 2482 | nn = gcd(gcd(absValue(x), absValue(y)), absValue(z)); |
---|
| 2483 | matrix point[1][3]; |
---|
| 2484 | point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn; |
---|
| 2485 | export point; return (R); |
---|
| 2486 | } |
---|
| 2487 | else |
---|
| 2488 | { |
---|
| 2489 | bigint fs = 4*c*f - e^2; |
---|
| 2490 | bigint ds = 4*c*d - 2*b*e; |
---|
| 2491 | x = -fs*2*c; y = b*fs-e*ds; z = ds*2*c; |
---|
| 2492 | nn = gcd(gcd(absValue(x), absValue(y)), absValue(z)); |
---|
| 2493 | matrix point[1][3]; |
---|
| 2494 | point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn; |
---|
| 2495 | export point; return (R); |
---|
| 2496 | } |
---|
| 2497 | } |
---|
| 2498 | |
---|
| 2499 | if (d^2 == 4*a*f) |
---|
| 2500 | { |
---|
| 2501 | if (f == 0) |
---|
| 2502 | { |
---|
| 2503 | x = -b*e*2; y = e*4*a; z = b^2-4*a*c; |
---|
| 2504 | nn = gcd(gcd(absValue(x), absValue(y)), absValue(z)); |
---|
| 2505 | matrix point[1][3]; |
---|
| 2506 | point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn; |
---|
| 2507 | export point; return (R); |
---|
| 2508 | } |
---|
| 2509 | else |
---|
| 2510 | { |
---|
| 2511 | bigint c_s = 4*c*f - e^2; |
---|
| 2512 | bigint b_s = 4*f*b - 2*d*e; |
---|
| 2513 | x = -c_s*2*f; y = b_s*2*f; z = d*c_s-e*b_s; |
---|
| 2514 | nn = gcd(gcd(absValue(x), absValue(y)), absValue(z)); |
---|
| 2515 | matrix point[1][3]; |
---|
| 2516 | point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn; |
---|
| 2517 | export point; return (R); |
---|
| 2518 | } |
---|
| 2519 | } |
---|
| 2520 | |
---|
| 2521 | if (e^2 == 4*c*f) |
---|
| 2522 | { |
---|
| 2523 | if (c == 0) |
---|
| 2524 | { |
---|
| 2525 | x = b*4*f; y = d^2-4*a*f; z = -b*d*2; |
---|
| 2526 | nn = gcd(gcd(absValue(x), absValue(y)), absValue(z)); |
---|
| 2527 | matrix point[1][3]; |
---|
| 2528 | point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn; |
---|
| 2529 | export point; return (R); |
---|
| 2530 | } |
---|
| 2531 | else |
---|
| 2532 | { |
---|
| 2533 | bigint as = 4*c*a - b^2; |
---|
| 2534 | bigint ds = 4*c*d - 2*b*e; |
---|
| 2535 | x = ds*2*c; y = e*as-b*ds; z = -as*2*c; |
---|
| 2536 | nn = gcd(gcd(absValue(x), absValue(y)), absValue(z)); |
---|
| 2537 | matrix point[1][3]; |
---|
| 2538 | point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn; |
---|
| 2539 | export point; return (R); |
---|
| 2540 | } |
---|
| 2541 | } |
---|
| 2542 | |
---|
| 2543 | ideal mi = maxideal(1); |
---|
| 2544 | map mm = R, mi; |
---|
| 2545 | bigint B; bigint C; bigint D; |
---|
| 2546 | |
---|
| 2547 | if ((a == 0) && (c == 0)) |
---|
| 2548 | { |
---|
| 2549 | B = -1; C = 4*b*f - 4*d*e; |
---|
| 2550 | /* now, b <> 0 since otherwise p would have the factor z, |
---|
| 2551 | and hence not be irreducible */ |
---|
| 2552 | mm[1] = (var(1)+var(2)-2*e*var(3))/(2*b); |
---|
| 2553 | mm[2] = (var(1)-var(2)-2*d*var(3))/(2*b); |
---|
| 2554 | } |
---|
| 2555 | if ((a != 0) && (c == 0)) |
---|
| 2556 | { |
---|
| 2557 | mm[1] = var(2); |
---|
| 2558 | mm[2] = var(1); |
---|
| 2559 | bigint t = a; a = c; c = t; |
---|
| 2560 | t = e; e = d; d = t; |
---|
| 2561 | } |
---|
| 2562 | if (c != 0) |
---|
| 2563 | { |
---|
| 2564 | D = 4*a*c-b^2; |
---|
| 2565 | mm[2] = (var(2)-e*var(3)-b*var(1))/(2*c); |
---|
| 2566 | map mTemp = basering, mi; |
---|
| 2567 | mTemp[1] = (var(1)-2*d*c*var(3)+b*e*var(3))/D; |
---|
| 2568 | mm = mTemp(mm); |
---|
| 2569 | B = D; |
---|
| 2570 | C = 16*a*c^2*f-4*a*c*e^2-4*b^2*c*f+4*b*c*d*e-4*c^2*d^2; |
---|
| 2571 | } |
---|
| 2572 | list K; |
---|
| 2573 | if ((B > 0) && (C >= 0)) { K = 0; } |
---|
| 2574 | if ((B >= 0) && (C > 0)) { K = 0; } |
---|
| 2575 | if (B == 0) |
---|
| 2576 | { |
---|
| 2577 | /* looking for a point on X^2 = |C| * Z^2 */ |
---|
| 2578 | bigint root = largestSquare(absValue(C)); |
---|
| 2579 | if (absValue(C)/root/root == 1) { K = 1, root, 0, 1; } |
---|
| 2580 | else { K = 0; } |
---|
| 2581 | } |
---|
| 2582 | if (C == 0) |
---|
| 2583 | { |
---|
| 2584 | /* looking for a point on X^2 = |B| * Y^2 */ |
---|
| 2585 | bigint root = largestSquare(absValue(B)); |
---|
| 2586 | if (absValue(B)/root/root == 1) { K = 1, root, 1, 0; } |
---|
| 2587 | else { K = 0; } |
---|
| 2588 | } |
---|
| 2589 | else { K = rationalPointSpecial(B, C); } |
---|
| 2590 | if (K[1] == 0) |
---|
| 2591 | { |
---|
| 2592 | /* no rational point on conic; |
---|
| 2593 | we need to move to an appropriate field extension Q(a) */ |
---|
| 2594 | poly h1 = mm[1]; poly h2 = mm[2]; poly h3 = mm[3]; |
---|
| 2595 | def extendedR = buildExtension(B, C, K[2], K[3]); |
---|
| 2596 | setring extendedR; |
---|
| 2597 | poly g1 = fetch(R, h1); |
---|
| 2598 | poly g2 = fetch(R, h2); |
---|
| 2599 | poly g3 = fetch(R, h3); |
---|
| 2600 | matrix temp[1][3]; |
---|
| 2601 | temp[1, 1] = subst(g1, var(1), point[1, 1], var(2), point[1, 2], |
---|
| 2602 | var(3), point[1, 3]); |
---|
| 2603 | temp[1, 2] = subst(g2, var(1), point[1, 1], var(2), point[1, 2], |
---|
| 2604 | var(3), point[1, 3]); |
---|
| 2605 | temp[1, 3] = subst(g3, var(1), point[1, 1], var(2), point[1, 2], |
---|
| 2606 | var(3), point[1, 3]); |
---|
| 2607 | point[1, 1] = temp[1, 1]; point[1, 2] = temp[1, 2]; |
---|
| 2608 | point[1, 3] = temp[1, 3]; |
---|
| 2609 | setring R; |
---|
| 2610 | return (extendedR); |
---|
| 2611 | } |
---|
| 2612 | else |
---|
| 2613 | { |
---|
[e9259f] | 2614 | string dummyString = string(K); // without this useless line, we |
---|
| 2615 | // sometimes get a seg fault because |
---|
| 2616 | // mm is corrupted; strange!?!?!?!? |
---|
[1a96a4] | 2617 | number nx = number(subst(mm[1], var(1), K[2], var(2), K[3], var(3), K[4])); |
---|
| 2618 | number ny = number(subst(mm[2], var(1), K[2], var(2), K[3], var(3), K[4])); |
---|
| 2619 | number nz = number(subst(mm[3], var(1), K[2], var(2), K[3], var(3), K[4])); |
---|
| 2620 | /* the point (nx, ny, nz) is already a solution; |
---|
| 2621 | the following lines will just remove denominators and reduce |
---|
| 2622 | numerators in order to return a nice tuple from Z^3 */ |
---|
| 2623 | bigint nxd = bigint(denominator(absValue(nx))); |
---|
| 2624 | bigint nyd = bigint(denominator(absValue(ny))); |
---|
| 2625 | bigint nzd = bigint(denominator(absValue(nz))); |
---|
| 2626 | nn = nxd * nyd / gcd(nxd, nyd); |
---|
| 2627 | nn = nn * nzd / gcd(nn, nzd); |
---|
| 2628 | x = bigint(nx*nn); y = bigint(ny*nn); z = bigint(nz*nn); |
---|
| 2629 | nn = gcd(gcd(absValue(x), absValue(y)), absValue(z)); |
---|
| 2630 | matrix point[1][3]; |
---|
| 2631 | point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn; |
---|
| 2632 | export point; |
---|
| 2633 | return (R); |
---|
| 2634 | } |
---|
| 2635 | } |
---|
| 2636 | |
---|
[d2ea299] | 2637 | example |
---|
[75e82e] | 2638 | { "EXAMPLE:"; echo=2; |
---|
[d2ea299] | 2639 | ring R = 0, (x,y,z), dp; |
---|
| 2640 | system("random", 4711); |
---|
| 2641 | poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z; |
---|
| 2642 | def S = rationalPointConic(p); // quadratic field extension, |
---|
| 2643 | // minpoly = a^2 - 2 |
---|
| 2644 | testPointConic(p, S); |
---|
| 2645 | setring R; |
---|
| 2646 | p = x^2 - 1857669520 * y^2 + 86709575222179747132487270400 * z^2; |
---|
| 2647 | S = rationalPointConic(p); // same as current basering, |
---|
| 2648 | // no extension needed |
---|
| 2649 | testPointConic(p, S); |
---|
| 2650 | } |
---|
[1a96a4] | 2651 | /////////////////////////////////////////////////////////////////////////////// |
---|
[6c265b] | 2652 | proc testParametrization(poly f, def rTT) |
---|
[a212fe] | 2653 | "USAGE: testParametrization(f, rTT); f poly, rTT ring |
---|
[6c265b] | 2654 | ASSUME: The assumptions on the basering and the polynomial f are as required |
---|
[a212fe] | 2655 | by @ref{paraPlaneCurve}. The ring rTT has two variables and contains |
---|
[1a96a4] | 2656 | an ideal PARA (such as the ring obtained by applying |
---|
[a212fe] | 2657 | @ref{paraPlaneCurve} to f). |
---|
[6c265b] | 2658 | RETURN: int which is 1 if PARA defines a parametrization of the curve |
---|
[1a96a4] | 2659 | {f=0} and 0, otherwise. |
---|
[5e8ee4c] | 2660 | THEORY: We compute the polynomial defining the image of PARA |
---|
[a212fe] | 2661 | and compare it with f. |
---|
[1a96a4] | 2662 | KEYWORDS: Parametrization, image. |
---|
| 2663 | EXAMPLE: example testParametrization; shows an example |
---|
| 2664 | " |
---|
| 2665 | { |
---|
| 2666 | def Roriginal = basering; |
---|
| 2667 | setring rTT; |
---|
| 2668 | /* begin workaround elimination*/ |
---|
| 2669 | int k; |
---|
| 2670 | list rl = ringlist(rTT); |
---|
| 2671 | rl[2] = list("s","t","x","y","z"); |
---|
| 2672 | rl[3]= list(list("dp",1:5),list("C",0)); |
---|
| 2673 | def Relim = ring(rl); |
---|
| 2674 | setring Relim; |
---|
| 2675 | ideal PARA = fetch(rTT,PARA); |
---|
[6c265b] | 2676 | ideal JJ; |
---|
[1a96a4] | 2677 | for(k=1;k<=3;k++) |
---|
| 2678 | { |
---|
| 2679 | JJ=JJ,var(k+2)-PARA[k]; |
---|
| 2680 | } |
---|
| 2681 | ideal SJJ = std(JJ); |
---|
| 2682 | intvec HJJ = hilb(SJJ,1); |
---|
| 2683 | ideal J = eliminate(JJ,var(1)*var(2),HJJ); |
---|
| 2684 | setring rTT; |
---|
| 2685 | /*end workaround elimination*/ |
---|
| 2686 | rl[2] = list("x","y","z"); |
---|
| 2687 | rl[3] = list(list("dp",1:3),list("C",0)); |
---|
| 2688 | def RP2 = ring(rl); |
---|
| 2689 | setring RP2; |
---|
| 2690 | ideal f = fetch(Roriginal,f); |
---|
| 2691 | ideal ftest = imap(Relim,J); |
---|
| 2692 | poly g = reduce(f[1],std(ftest)); |
---|
[62dc18e] | 2693 | if(g!=0){return(0)} |
---|
[1a96a4] | 2694 | g = reduce(ftest[1],std(ideal(f))); |
---|
[62dc18e] | 2695 | if(g!=0){return(0)} |
---|
[6c265b] | 2696 | return (1); |
---|
[1a96a4] | 2697 | } |
---|
| 2698 | |
---|
| 2699 | example |
---|
[75e82e] | 2700 | { "EXAMPLE:"; echo=2; |
---|
[1a96a4] | 2701 | ring R = 0,(x,y,z),dp; |
---|
| 2702 | poly f = y^8-x^3*(z+x)^5; |
---|
| 2703 | def RP1 = paraPlaneCurve(f); |
---|
| 2704 | testParametrization(f, RP1); |
---|
| 2705 | } |
---|
| 2706 | |
---|
| 2707 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 2708 | proc testPointConic(poly p, def r) |
---|
| 2709 | "USAGE: testPointConic(p, r); p poly, r ring |
---|
| 2710 | ASSUME: assumes that p is a homogeneous quadratic polynomial in the |
---|
| 2711 | first three ring variables of the current basering; |
---|
| 2712 | Assumes that there is a (1x3) matrix named 'point' in r with |
---|
| 2713 | entries from the ground field of r. |
---|
| 2714 | RETURN: returns 1 iff the point named 'point', residing in r, lies on |
---|
| 2715 | the conic given by p; 0 otherwise |
---|
| 2716 | NOTE: This method temporarily changes the basering to r. Afterwards, |
---|
| 2717 | the basering will be the same as before. |
---|
| 2718 | EXAMPLE: example testPointConic; shows an example |
---|
| 2719 | " |
---|
| 2720 | { |
---|
| 2721 | def rOrig = basering; |
---|
| 2722 | "conic:", p; |
---|
| 2723 | setring r; |
---|
| 2724 | string s = "point: " + string(point[1,1]) + ", " + string(point[1,2]); |
---|
| 2725 | s = s + ", " + string(point[1,3]); |
---|
| 2726 | s; |
---|
| 2727 | if (minpoly != 0) { "minpoly:", minpoly; } |
---|
| 2728 | poly f = fetch(rOrig, p); |
---|
| 2729 | poly g = subst(f, var(1), point[1,1], |
---|
| 2730 | var(2), point[1,2], |
---|
| 2731 | var(3), point[1,3]); |
---|
| 2732 | int result = 0; if (g == 0) { result = 1; } |
---|
| 2733 | setring rOrig; |
---|
| 2734 | return (result); |
---|
| 2735 | } |
---|
| 2736 | |
---|
| 2737 | example |
---|
[75e82e] | 2738 | { "EXAMPLE:"; echo=2; |
---|
[1a96a4] | 2739 | ring R = 0, (x,y,z), dp; |
---|
| 2740 | system("random", 4711); |
---|
| 2741 | poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z; |
---|
| 2742 | def S = rationalPointConic(p); |
---|
| 2743 | if (testPointConic(p, S) == 1) |
---|
| 2744 | { "point lies on conic"; } |
---|
| 2745 | else |
---|
| 2746 | { "point does not lie on conic"; } |
---|
| 2747 | } |
---|
| 2748 | |
---|
| 2749 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 2750 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 2751 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 2752 | /* |
---|
| 2753 | ///////////////////////////////////////////////////////////////////////////// |
---|
| 2754 | /// Further examples for testing the main procedures |
---|
| 2755 | /// Timings on wawa Sept 29 |
---|
| 2756 | ///////////////////////////////////////////////////////////////////////////// |
---|
[d2ea299] | 2757 | LIB"paraplanecurves.lib"; |
---|
[1a96a4] | 2758 | // ------------------------------------------------------- |
---|
| 2759 | // Example 1 |
---|
| 2760 | // ------------------------------------------------------- |
---|
| 2761 | ring RR = 0, (x,y,z), dp; |
---|
| 2762 | poly f = 7*z2+11*y2+13*z*y+17*x2+19*x*y; // conic |
---|
| 2763 | def RP1 = paraConic(f); |
---|
| 2764 | setring RP1; PARACONIC; |
---|
| 2765 | setring RR; |
---|
| 2766 | RP1 = paraPlaneCurve(f); |
---|
| 2767 | testParametrization(f,RP1); |
---|
| 2768 | setring RP1; PARA; |
---|
| 2769 | kill RR;kill RP1; |
---|
| 2770 | // ------------------------------------------------------- |
---|
| 2771 | // Example 2 |
---|
| 2772 | // ------------------------------------------------------- |
---|
| 2773 | ring RR = 0, (x,y,z), dp; |
---|
| 2774 | poly f = y3-x2z; // cusp at origin |
---|
| 2775 | adjointIdeal(f,1); |
---|
| 2776 | adjointIdeal(f,2); |
---|
| 2777 | def RP1 = paraPlaneCurve(f); // time 0 |
---|
| 2778 | testParametrization(f,RP1); |
---|
| 2779 | setring RP1; PARA; |
---|
| 2780 | kill RR;kill RP1; |
---|
| 2781 | // ------------------------------------------------------- |
---|
| 2782 | // Example 3 |
---|
| 2783 | // ------------------------------------------------------- |
---|
| 2784 | ring RR = 0, (x,y,z), dp; |
---|
| 2785 | poly f=(xz-y^2)^2-x*y^3; // 1 sing at origin, 1 cusp, no OMPs |
---|
| 2786 | adjointIdeal(f,1); |
---|
| 2787 | adjointIdeal(f,2); |
---|
| 2788 | def RP1 = paraPlaneCurve(f); // time 0 |
---|
| 2789 | testParametrization(f,RP1); |
---|
| 2790 | setring RP1; PARA; |
---|
| 2791 | kill RR;kill RP1; |
---|
| 2792 | // ------------------------------------------------------- |
---|
| 2793 | // Example 4 |
---|
| 2794 | // ------------------------------------------------------- |
---|
| 2795 | ring RR = 0, (x,y,z), dp; |
---|
| 2796 | poly f = y5-y4x+4y2x2z-x4z; // 1 sing at origin, no OMPs, no cusps |
---|
| 2797 | adjointIdeal(f,1); |
---|
| 2798 | adjointIdeal(f,2); |
---|
| 2799 | def RP1 = paraPlaneCurve(f); // time 0 |
---|
| 2800 | testParametrization(f,RP1); |
---|
| 2801 | setring RP1; PARA; |
---|
| 2802 | kill RR;kill RP1; |
---|
| 2803 | // ------------------------------------------------------- |
---|
| 2804 | // Example 5 |
---|
| 2805 | // ------------------------------------------------------- |
---|
| 2806 | ring RR = 0, (x,y,z), dp; |
---|
| 2807 | poly f = 259x5-31913x4y+939551x3y2+2871542x2y3+2845801xy4; |
---|
| 2808 | f = f+914489y5+32068x4z-1884547x3yz-8472623x2y2z-11118524xy3z; |
---|
| 2809 | f = f-4589347y4z+944585x3z2+8563304x2yz2+16549772xy2z2+9033035y3z2; |
---|
| 2810 | f = f-2962425x2z3-11214315xyz3-8951744y2z3+2937420xz4+4547571yz4-953955z5; |
---|
| 2811 | // 6 nodes |
---|
| 2812 | adjointIdeal(f,1); |
---|
| 2813 | adjointIdeal(f,2); |
---|
| 2814 | def RP1 = paraPlaneCurve(f); // time 0 |
---|
| 2815 | testParametrization(f,RP1); |
---|
| 2816 | setring RP1; PARA; |
---|
| 2817 | kill RR;kill RP1; |
---|
| 2818 | // ------------------------------------------------------- |
---|
| 2819 | // Example 7 |
---|
| 2820 | // ------------------------------------------------------- |
---|
| 2821 | ring RR = 0, (x,y,z), dp; |
---|
| 2822 | poly f = y^8-x^3*(z+x)^5; // 1 sing at origin, 1 further sing, no OMPs, |
---|
| 2823 | // no cusps |
---|
| 2824 | adjointIdeal(f,1); |
---|
| 2825 | adjointIdeal(f,2); |
---|
| 2826 | def RP1 = paraPlaneCurve(f); // time 0 |
---|
| 2827 | testParametrization(f,RP1); |
---|
| 2828 | setring RP1; PARA; |
---|
| 2829 | kill RR;kill RP1; |
---|
| 2830 | // ------------------------------------------------------- |
---|
| 2831 | // Example 8 |
---|
| 2832 | // ------------------------------------------------------- |
---|
| 2833 | ring RR = 0, (x,y,z), dp; |
---|
| 2834 | poly f = 11y7+7y6x+8y5x2-3y4x3-10y3x4-10y2x5-x7-33y6-29y5x-13y4x2+26y3x3; |
---|
| 2835 | f = f+30y2x4+10yx5+3x6+33y5+37y4x-8y3x2-33y2x3-20yx4-3x5-11y4-15y3x; |
---|
| 2836 | f = f+13y2x2+10yx3+x4; // 3 OMPs of mult 3, 1 OMP of mult 4 |
---|
| 2837 | f = homog(f,z); |
---|
| 2838 | adjointIdeal(f,1); |
---|
| 2839 | adjointIdeal(f,2); |
---|
| 2840 | def RP1 = paraPlaneCurve(f); // time 0 |
---|
| 2841 | testParametrization(f,RP1); |
---|
| 2842 | setring RP1; PARA; |
---|
| 2843 | kill RR;kill RP1; |
---|
| 2844 | // ------------------------------------------------------- |
---|
| 2845 | // Example 9 |
---|
| 2846 | // ------------------------------------------------------- |
---|
| 2847 | ring RR = 0, (x,y,z), dp; |
---|
| 2848 | poly f = y^8-x^3*(z+x)^5; // 1 sing at origin, 1 further sing, no OMPs, |
---|
| 2849 | // no cusps |
---|
| 2850 | adjointIdeal(f,1); |
---|
| 2851 | adjointIdeal(f,2); |
---|
| 2852 | def RP1 = paraPlaneCurve(f); // time 0 |
---|
[6c265b] | 2853 | testParametrization(f,RP1); |
---|
[1a96a4] | 2854 | setring RP1; PARA; |
---|
| 2855 | kill RR;kill RP1; |
---|
| 2856 | // ------------------------------------------------------- |
---|
| 2857 | // Example 10 |
---|
| 2858 | // ------------------------------------------------------- |
---|
| 2859 | ring SS = 0, (u,v,z), dp; |
---|
| 2860 | poly f = u^4-14*u^2*v^2+v^4+8*u^2*v*z+8*v^3*z; // 1 OMP of mult 3 at orgin |
---|
| 2861 | adjointIdeal(f,1); |
---|
| 2862 | adjointIdeal(f,2); |
---|
| 2863 | def RP1 = paraPlaneCurve(f); // time 0 |
---|
| 2864 | testParametrization(f,RP1); |
---|
| 2865 | setring RP1; PARA; |
---|
| 2866 | kill SS;kill RP1; |
---|
| 2867 | // ------------------------------------------------------- |
---|
| 2868 | // Example 11 |
---|
| 2869 | // ------------------------------------------------------- |
---|
| 2870 | ring SS = 0, (u,v,z), dp; |
---|
| 2871 | poly f = 14440*u^5-16227*u^4*v+10812*u^3*v^2-13533*u^2*v^3+3610*u*v^4; |
---|
| 2872 | f = f+1805*v^5+14440*u^4*z-18032*u^3*v*z+16218*u^2*v^2*z-12626*u*v^3*z; |
---|
| 2873 | f = f+3610*v^4*z+3610*u^3*z^2-4508*u^2*v*z^2+5406*u*v^2*z^2-2703*v^3*z^2; |
---|
| 2874 | // 1 OMP of mult 3 at origin, 2 nodes |
---|
| 2875 | adjointIdeal(f,1); |
---|
| 2876 | adjointIdeal(f,2); |
---|
| 2877 | def RP1 = paraPlaneCurve(f); // time 0 |
---|
| 2878 | testParametrization(f,RP1); |
---|
| 2879 | setring RP1; PARA; |
---|
| 2880 | kill SS;kill RP1; |
---|
| 2881 | // ------------------------------------------------------- |
---|
| 2882 | // Example 12 |
---|
| 2883 | // ------------------------------------------------------- |
---|
| 2884 | ring SS = 0, (u,v,z), dp; |
---|
| 2885 | poly f = u^6+3*u^4*v^2+3*u^2*v^4+v^6-4*u^4*z^2-34*u^3*v*z^2-7*u^2*v^2*z^2; |
---|
| 2886 | f = f+12*u*v^3*z^2+6*v^4*z^2+36*u^2*z^4+36*u*v*z^4+9*v^2*z^4; |
---|
| 2887 | // needs field extension *** 6 nodes, 2 cusps, 1 sing at 0 |
---|
| 2888 | adjointIdeal(f,1); |
---|
| 2889 | adjointIdeal(f,2); |
---|
| 2890 | def RP1 = paraPlaneCurve(f); // time 0 |
---|
| 2891 | testParametrization(f,RP1); |
---|
| 2892 | setring RP1; PARA; |
---|
| 2893 | kill SS;kill RP1; |
---|
| 2894 | // ------------------------------------------------------- |
---|
| 2895 | // Example 13 |
---|
| 2896 | // ------------------------------------------------------- |
---|
| 2897 | ring SS = 0, (u,v,z), dp; |
---|
| 2898 | poly f = -24135/322*u^6-532037/6440*u^5*v+139459/560*u^4*v^2; |
---|
| 2899 | f = f-1464887/12880*u^3*v^3+72187/25760*u^2*v^4+9/8*u*v^5+1/8*v^6; |
---|
| 2900 | f = f-403511/3220*u^5*z-40817/920*u^4*v*z+10059/80*u^3*v^2*z; |
---|
| 2901 | f = f-35445/1288*u^2*v^3*z+19/4*u*v^4*z+3/4*v^5*z-20743/805*u^4*z^2; |
---|
| 2902 | f = f+126379/3220*u^3*v*z^2-423417/6440*u^2*v^2*z^2+11/2*u*v^3*z^2; |
---|
| 2903 | f = f+3/2*v^4*z^2+3443/140*u^3*z^3+u^2*v*z^3+u*v^2*z^3+v^3*z^3; |
---|
| 2904 | // 2 OMPs of mult 3 (1 at origin), 4 nodes |
---|
| 2905 | adjointIdeal(f,1); |
---|
| 2906 | adjointIdeal(f,2); |
---|
[ce136a] | 2907 | def RP1 = paraPlaneCurve(f); // time 14 |
---|
[1a96a4] | 2908 | testParametrization(f,RP1); |
---|
| 2909 | setring RP1; PARA; |
---|
| 2910 | kill SS;kill RP1; |
---|
| 2911 | // ------------------------------------------------------- |
---|
| 2912 | // Example 14 |
---|
| 2913 | // ------------------------------------------------------- |
---|
| 2914 | ring SS = 0, (u,v,z), dp; |
---|
[6c265b] | 2915 | poly f = |
---|
[1a96a4] | 2916 | 2*u^7+u^6*v+3*u^5*v^2+u^4*v^3+2*u^3*v^4+u^2*v^5+2*u*v^6+v^7 |
---|
| 2917 | -7780247/995328*u^6*z-78641/9216*u^5*v*z-10892131/995328*u^4*v^2*z |
---|
| 2918 | -329821/31104*u^3*v^3*z-953807/331776*u^2*v^4*z-712429/248832*u*v^5*z |
---|
| 2919 | +1537741/331776*v^6*z+2340431/248832*u^5*z^2+5154337/248832*u^4*v*z^2 |
---|
| 2920 | +658981/41472*u^3*v^2*z^2+1737757/124416*u^2*v^3*z^2 |
---|
| 2921 | -1234733/248832*u*v^4*z^2-1328329/82944*v^5*z^2-818747/248832*u^4*z^3 |
---|
| 2922 | -1822879/124416*u^3*v*z^3-415337/31104*u^2*v^2*z^3 |
---|
| 2923 | +1002655/124416*u*v^3*z^3+849025/82944*v^4*z^3; |
---|
| 2924 | // 3 OMPs of mult 3, 1 OMP of mult 4 at origin |
---|
| 2925 | adjointIdeal(f,2); |
---|
| 2926 | def RP1 = paraPlaneCurve(f); // time 1 |
---|
| 2927 | testParametrization(f,RP1); |
---|
| 2928 | setring RP1; PARA; |
---|
| 2929 | kill SS;kill RP1; |
---|
| 2930 | // ------------------------------------------------------- |
---|
| 2931 | // Example 15 |
---|
| 2932 | // ------------------------------------------------------- |
---|
| 2933 | ring SS = 0, (u,v,z), dp; |
---|
| 2934 | poly f = 590819418867856650536224u7-147693905508217596067968u6v; |
---|
| 2935 | f = f+229117518934972047619978u5v2-174050799674982973889542u4v3; |
---|
| 2936 | f = f-92645796479789150855110u3v4-65477418713685583062704u2v5; |
---|
| 2937 | f = f+4529961835917468460168uv6+7715404057796585983136v7; |
---|
| 2938 | f = f-413640780091141905428104u6z+571836835577486968144618u5vz; |
---|
| 2939 | f = f-551807810327826605739444u4v2z-488556410340789283359926u3v3z; |
---|
| 2940 | f = f-473466023008413178155962u2v4z+48556741573432247323608uv5z; |
---|
| 2941 | f = f+77647371229172269259528v6z+340450118906560552282893u5z2; |
---|
| 2942 | f = f-433598825064368371610344u4vz2-937281070591684636591672u3v2z2; |
---|
| 2943 | f = f-1388949843915129934647751u2v3z2+204081793110898617103998uv4z2; |
---|
| 2944 | f = f+335789953068251652554308v5z2+6485661002496681852577u4z3; |
---|
| 2945 | f = f-772700266516318390630202u3vz3-2068348417248100329533330u2v2z3; |
---|
| 2946 | f = f+440320154612359641806108uv3z3+808932515589210854581618v4z3; |
---|
| 2947 | f = f-229384307132237615286548u3z4-1564303565658228216055227u2vz4; |
---|
| 2948 | f = f+520778334468674798322974uv2z4+1172483905704993294097655v3z4; |
---|
| 2949 | f = f-480789741398016816562100u2z5+322662751598958620410786uvz5; |
---|
| 2950 | f = f+1022525576391791616258310v2z5+82293493608853837667471uz6; |
---|
| 2951 | f = f+496839109904761426785889vz6+103766136235628614937587z7; // 15 nodes |
---|
| 2952 | adjointIdeal(f,2); |
---|
| 2953 | def RP1 = paraPlaneCurve(f); // time 72 |
---|
[6c265b] | 2954 | testParametrization(f,RP1); |
---|
[1a96a4] | 2955 | setring RP1; PARA; |
---|
| 2956 | kill SS;kill RP1; |
---|
| 2957 | |
---|
| 2958 | // ------------------------------------------------------- |
---|
| 2959 | // Example 16 |
---|
| 2960 | // ------------------------------------------------------- |
---|
| 2961 | ring SS = 0, (u,v,z), dp; |
---|
| 2962 | poly f = 25*u^8+184*u^7*v+518*u^6*v^2+720*u^5*v^3+576*u^4*v^4+282*u^3*v^5; |
---|
| 2963 | f = f+84*u^2*v^6+14*u*v^7+v^8+244*u^7*z+1326*u^6*v*z+2646*u^5*v^2*z; |
---|
| 2964 | f = f+2706*u^4*v^3*z+1590*u^3*v^4*z+546*u^2*v^5*z+102*u*v^6*z+8*v^7*z; |
---|
| 2965 | f = f+854*u^6*z^2+3252*u^5*v*z^2+4770*u^4*v^2*z^2+3582*u^3*v^3*z^2; |
---|
| 2966 | f = f+1476*u^2*v^4*z^2+318*u*v^5*z^2+28*v^6*z^2+1338*u^5*z^3+3740*u^4*v*z^3; |
---|
| 2967 | f = f+4030*u^3*v^2*z^3+2124*u^2*v^3*z^3+550*u*v^4*z^3+56*v^5*z^3+1101*u^4*z^4; |
---|
| 2968 | f = f+2264*u^3*v*z^4+1716*u^2*v^2*z^4+570*u*v^3*z^4+70*v^4*z^4+508*u^3*z^5; |
---|
| 2969 | f = f+738*u^2*v*z^5+354*u*v^2*z^5+56*v^3*z^5+132*u^2*z^6+122*u*v*z^6; |
---|
| 2970 | f = f+28*v^2*z^6+18*u*z^7+8*v*z^7+z^8; // 3 nodes, 1 sing |
---|
| 2971 | adjointIdeal(f,1); |
---|
| 2972 | adjointIdeal(f,2); |
---|
| 2973 | def RP1 = paraPlaneCurve(f); // time 20 |
---|
| 2974 | testParametrization(f,RP1); |
---|
| 2975 | setring RP1; PARA; |
---|
| 2976 | kill SS;kill RP1; |
---|
| 2977 | // ------------------------------------------------------- |
---|
| 2978 | // Example 17 |
---|
| 2979 | // ------------------------------------------------------- |
---|
| 2980 | ring SS = 0, (u,v,z), dp; |
---|
| 2981 | poly f = -2*u*v^4*z^4+u^4*v^5+12*u^4*v^3*z^2+12*u^2*v^4*z^3-u^3*v*z^5; |
---|
| 2982 | f = f+11*u^3*v^2*z^4-21*u^3*v^3*z^3-4*u^4*v*z^4+2*u^4*v^2*z^3-6*u^4*v^4*z; |
---|
| 2983 | f = f+u^5*z^4-3*u^5*v^2*z^2+u^5*v^3*z-3*u*v^5*z^3-2*u^2*v^3*z^4+u^3*v^4*z^2; |
---|
[fa932ac] | 2984 | f = f+v^5*z^4; // 2 OMPs of mult 4, 1 OMP of mult 5, 1 sing at origin |
---|
[1a96a4] | 2985 | f = subst(f,z,u+v+z); |
---|
| 2986 | adjointIdeal(f,2); |
---|
| 2987 | def RP1 = paraPlaneCurve(f); // time 5 |
---|
| 2988 | testParametrization(f,RP1); |
---|
| 2989 | setring RP1; PARA; |
---|
| 2990 | kill SS;kill RP1; |
---|
| 2991 | // ------------------------------------------------------- |
---|
| 2992 | // Example 18 |
---|
| 2993 | // ------------------------------------------------------- |
---|
| 2994 | ring SS = 0, (u,v,z), dp; |
---|
| 2995 | poly f = u^5*v^5+21*u^5*v^4*z-36*u^4*v^5*z-19*u^5*v^3*z^2+12*u^4*v^4*z^2; |
---|
| 2996 | f = f+57*u^3*v^5*z^2+u^5*v^2*z^3+u^4*v^3*z^3-53*u^3*v^4*z^3-19*u^2*v^5*z^3; |
---|
| 2997 | f = f+u^5*v*z^4+43*u^3*v^3*z^4+u*v^5*z^4+u^5*z^5-15*u^3*v^2*z^5+u^2*v^3*z^5; |
---|
| 2998 | f = f+u*v^4*z^5+v^5*z^5; // 1 OMP of mult 4, 3 OMPs of mult 5 (1 at origin) |
---|
| 2999 | adjointIdeal(f,2); |
---|
| 3000 | def RP1 = paraPlaneCurve(f); // time 8 |
---|
| 3001 | testParametrization(f,RP1); |
---|
| 3002 | setring RP1; PARA; |
---|
| 3003 | kill SS;kill RP1; |
---|
| 3004 | // ------------------------------------------------------- |
---|
| 3005 | // Example 19 |
---|
| 3006 | // ------------------------------------------------------- |
---|
| 3007 | ring SS = 0, (u,v,z), dp; |
---|
| 3008 | poly f = u^10+6*u^9*v-30*u^7*v^3-15*u^6*v^4+u^5*v^5+u^4*v^6+6*u^3*v^7; |
---|
| 3009 | f = f+u^2*v^8+7*u*v^9+v^10+5*u^9*z+24*u^8*v*z-30*u^7*v^2*z-120*u^6*v^3*z; |
---|
| 3010 | f = f-43*u^5*v^4*z+5*u^4*v^5*z+20*u^3*v^6*z+10*u^2*v^7*z+29*u*v^8*z+5*v^9*z; |
---|
| 3011 | f = f+10*u^8*z^2+36*u^7*v*z^2-105*u^6*v^2*z^2-179*u^5*v^3*z^2-38*u^4*v^4*z^2; |
---|
| 3012 | f = f+25*u^3*v^5*z^2+25*u^2*v^6*z^2+46*u*v^7*z^2+10*v^8*z^2+10*u^7*z^3; |
---|
| 3013 | f = f+24*u^6*v*z^3-135*u^5*v^2*z^3-117*u^4*v^3*z^3-u^3*v^4*z^3+25*u^2*v^5*z^3; |
---|
| 3014 | f = f+34*u*v^6*z^3+10*v^7*z^3+5*u^6*z^4+6*u^5*v*z^4-75*u^4*v^2*z^4; |
---|
| 3015 | f = f-27*u^3*v^3*z^4+10*u^2*v^4*z^4+11*u*v^5*z^4+5*v^6*z^4+u^5*z^5; |
---|
| 3016 | f = f-15*u^3*v^2*z^5+u^2*v^3*z^5+u*v^4*z^5+v^5*z^5; |
---|
| 3017 | // 1 OMP of mult 4, 3 OMPs of mult 5 (1 at origin) |
---|
| 3018 | adjointIdeal(f,2); |
---|
[d2ea299] | 3019 | def RP1 = paraPlaneCurve(f); // time 2 // see Ex. 18 |
---|
[1a96a4] | 3020 | testParametrization(f,RP1); |
---|
| 3021 | setring RP1; PARA; |
---|
| 3022 | kill SS;kill RP1; |
---|
| 3023 | // ------------------------------------------------------- |
---|
| 3024 | // Example 20 |
---|
| 3025 | // ------------------------------------------------------- |
---|
| 3026 | ring R = 0, (x,y,z), dp; |
---|
| 3027 | system("random", 4711); |
---|
| 3028 | poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z; |
---|
| 3029 | def S = rationalPointConic(p); // quadratic field extension, |
---|
| 3030 | // minpoly = a^2 - 2 |
---|
| 3031 | if (testPointConic(p, S) == 1) |
---|
| 3032 | { "point lies on conic"; } |
---|
| 3033 | else |
---|
| 3034 | { "point does not lie on conic"; } |
---|
[ce136a] | 3035 | kill R;kill S; |
---|
[1a96a4] | 3036 | // ------------------------------------------------------- |
---|
| 3037 | // Example 21 |
---|
| 3038 | // ------------------------------------------------------- |
---|
| 3039 | ring R = 0, (x,y,z), dp; |
---|
| 3040 | system("random", 4711); |
---|
| 3041 | poly p = x^2 - 1857669520 * y^2 + 86709575222179747132487270400 * z^2; |
---|
| 3042 | def S = rationalPointConic(p); // same as current basering, |
---|
| 3043 | // no extension needed |
---|
| 3044 | if (testPointConic(p, S) == 1) |
---|
| 3045 | { "point lies on conic"; } |
---|
| 3046 | else |
---|
| 3047 | { "point does not lie on conic"; } |
---|
[ce136a] | 3048 | kill R;kill S; |
---|
| 3049 | // ------------------------------------------------------- |
---|
| 3050 | // Example 21 |
---|
| 3051 | // ------------------------------------------------------- |
---|
| 3052 | ring RR = 0, (x,y,z), dp; |
---|
| 3053 | poly f = -1965466244509920x5y+34871245546721380061760x4y2; |
---|
| 3054 | f = f+104613747941595046117320x3y3+113331564241941002407560x2y4; |
---|
| 3055 | f = f+52306876673313609259800xy5+8717812860780028397880y6; |
---|
| 3056 | f = f+1040297748510024x5z+4468147845634872x4yz; |
---|
| 3057 | f = f-22398508728211453743258x3y2z-33223996581074443306854x2y3z; |
---|
| 3058 | f = f-10638598235041298082366xy4z+186886189971594356382y5z; |
---|
| 3059 | f = f-1385078844909312x4z2-34893092731637052532683x3yz2; |
---|
| 3060 | f = f-98591463214095439056609x2y2z2-92339459334829609336485xy3z2; |
---|
| 3061 | f = f-24923289542522905755711y4z2+472440640471377x3z3; |
---|
| 3062 | f = f+33821511925664516716011x2yz3+49745237303968344397437xy2z3; |
---|
| 3063 | f = f+11040465960074786720475y3z3+8728735735878837099404x2z4; |
---|
| 3064 | f = f+17676785754519678518537xyz4+17935885079051421934609y2z4; |
---|
| 3065 | f = f-11314701999743172607075xz5-16164284825803158969425yz5; |
---|
| 3066 | f = f+3666695988537425618750z6; |
---|
| 3067 | // 4 nodes, 1 OMP of mult 4 |
---|
| 3068 | adjointIdeal(f,2); |
---|
| 3069 | kill RR; |
---|
[1a96a4] | 3070 | */ |
---|