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