source: git/Singular/LIB/paraplanecurves.lib @ e8eba7

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