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