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

jengelh-datetimespielwiese
Last change on this file since fee24e was fee24e, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: sqrfree command now returns also multiplicities by default
  • Property mode set to 100644
File size: 100.1 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id$";
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. @*
28From 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. @*
257From 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.@*
488From 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 polyModP(poly q, bigint p)
1944"USAGE:  polyModP(q, p); q poly, p bigint
1945RETURN:  takes each coefficient of q modulo p and returns the resulting poly
1946"
1947{
1948  poly qq = q; poly res = 0;
1949  bigint c;
1950  while (qq != 0)
1951  {
1952    c = bigint(leadcoef(qq) mod p);
1953    res = res + c * leadmonom(qq);
1954    qq = qq - lead(qq);
1955  }
1956  return (res);
1957}
1958
1959
1960///////////////////////////////////////////////////////////////////////////////
1961static proc rootModP(bigint r, bigint p)
1962"USAGE:  rootModP(r, p); r, p bigint's
1963ASSUME:  0 <= r < p, and p prime;
1964         Furthermore it is assumes that there is some x in {0, 1, ..., p-1}
1965         such that x^2 = r mod p;
1966RETURN:  an x in {0, 1, ..., p-1} such that x^2 = r mod p;
1967THEORY:  For p larger than 32003, this computation is done using Cantor-
1968         Zassenhaus' algorithm. Otherwise a brute force approach is used.
1969KEYWORDS: Cantor-Zassenhaus algorithm.
1970"
1971{
1972  if (r == 0) { return (0); }
1973  if (r == 1) { return (1); }
1974  if (p <= 32003)
1975  {
1976    /* For small p, we use a brute force approach: */
1977    int i;
1978    for (i = 2; i < p; i++)
1979    {
1980      if (((i*i) mod p) == r) { return (i); }
1981    }
1982    /* should never be reached: */
1983    return (-1);
1984  }
1985
1986  /* For p > 32003, we use Cantor-Zassenhaus' algorithm: */
1987  def br = basering;
1988  ring rTemp = 0, x, dp;
1989  bigint b; bigint exponent; poly factor;
1990  poly h = x^2 - r;
1991  ideal redI = h; redI = std(redI);
1992  poly q = x^2; bigint root = 0;
1993  while (root == 0)
1994  {
1995    b = bigint(random(1, 2^30));
1996    exponent = bigint((p - 1) div 2);
1997    /* We need to compute q^exponent mod (x^2 - a) and mod p: */
1998    factor = x + b; q = 1;
1999    while (exponent > 0)
2000    {
2001      if ((exponent mod 2) == 1)
2002      {
2003        q = q * factor;
2004        q = reduce(q, redI);
2005        q = polyModP(q, p);
2006      }
2007      exponent = bigint(exponent div 2);
2008      factor = factor * factor;
2009      factor = reduce(factor, redI);
2010      factor = polyModP(factor, p);
2011    }
2012    if (deg(q) == 1)
2013    {
2014      q = q - 1;
2015      b = inverseModP(bigint(leadcoef(q)), p);
2016      q = q - lead(q);
2017      root = (bigint(q) * b) mod p;
2018      if (((root * root - r) mod p) != 0) { root = 0; }
2019    }
2020  }
2021  setring br; kill rTemp;
2022  return (root);
2023}
2024
2025
2026///////////////////////////////////////////////////////////////////////////////
2027static proc inverseModP(bigint r, bigint p)
2028"USAGE:  inverseModP(r, p); r, p bigint's
2029ASSUME:  0 <= r < p, and r and p coprime;
2030RETURN:  returns the inverse of r in Z/p represented by an element in
2031         {1, 2, ..., p-1}
2032THEORY:  This uses Euclid's extended gcd algorithm.
2033"
2034{
2035  list L = extendedEuclid(r, p);
2036  if (L[1] != 1) { "ERROR: GCD of", r, "and", p, "should be 1."; }
2037  L[2] = L[2] mod p;
2038  return (L[2]);
2039}
2040
2041
2042///////////////////////////////////////////////////////////////////////////////
2043static proc squareRoot(bigint r, bigint m, int justCheck)
2044"USAGE:  squareRoot(r, m, j); r, m bigint's, j int
2045RETURN:  checks whether r is a square modulo m, i.e., checks whether there is
2046         some x such that x^2 = r mod m;
2047         If justCheck is 1, then the method will terminate after the check
2048         and return 1 if r is a square and -1 otherwise.
2049         If justCheck is 0 and r is a square, then the method continues and
2050         computes a solution x in {0, 1, m-1} with x^2 = r mod m, which will
2051         then be returned
2052THEORY:  This algorithm checks solvability by computing the Legendre symbols
2053         modulo all primes in m. Then, individual roots will be computed and
2054         lifted to the desired square root modulo m using Chinese
2055         remaindering.
2056"
2057{
2058  if (m == 0) { "ERROR: squareRoot had been invoked with m = 0"; }
2059
2060  list L = primefactors(m);
2061  if ((L[3] != 1) && (L[3] != -1))
2062  { "WARNING: command 'primefactors(.)' did not find all prime factors"; }
2063  int i;
2064  for (i = 1; i <= size(L[2]); i++)
2065  {
2066    if (legendreSymbol(r, L[1][i]) == -1) { return (-1); }
2067  }
2068  /* now we know that there is some x in {0, 1, m-1} with
2069     x^2 = r mod m */
2070  if (justCheck == 1) { return (1); }
2071  else
2072  {
2073    // now we need to compute x; this works in two stages:
2074    // 1) write m = p1^e1 * ... * pk^ek (we already have that),
2075    // 2) for each i in {1, 2, ..., k}
2076    //    2.1) compute a yi such that yi^2 = r mod pi,
2077    //    2.2) lift yi to an xi such that xi^2 = r mod (pi^ei),
2078    // 3) lift (x1, x2, ..., xk) in Z/p1^e1 * ... * Z/pk^ek
2079    //    to x in Z/m via Chinese remainder theorem
2080
2081    list roots;
2082    // 2.1):
2083    for (i = 1; i <= size(L[1]); i++)
2084    {
2085      roots = insert(roots, rootModP(r mod L[1][i], L[1][i]), size(roots));
2086    }
2087
2088    // 2.2):
2089    bigint c; bigint l; bigint temp; bigint pPower; int e;
2090    for (i = 1; i <= size(roots); i++)
2091    {
2092      pPower = bigint(L[1][i]);
2093      for (e = 2; e <= L[2][i]; e++)
2094      {
2095        c = bigint(roots[i]); l = pPower;
2096        temp = r - c * c; l = bigint(2) * c * l; c = temp;
2097        c = c div pPower; l = l div pPower;
2098        c = c mod L[1][i]; l = l mod L[1][i];
2099        c = (c * bigint(inverseModP(l, L[1][i]))) mod L[1][i];
2100        c = bigint(roots[i]) + c * pPower;
2101        pPower = pPower * L[1][i]; roots[i] = c;
2102      }
2103    }
2104
2105    // 2.3):
2106    list mm; bigint z; int j;
2107    for (i = 1; i <= size(L[1]); i++)
2108    {
2109      z = bigint(L[1][i]);
2110      for (j = 2; j <= L[2][i]; j++)
2111      {
2112        z = z * bigint(L[1][i]);
2113      }
2114      mm = insert(mm, z, size(mm));
2115    }
2116    return (chineseRemainder(roots, mm) mod m);
2117  }
2118}
2119
2120
2121///////////////////////////////////////////////////////////////////////////////
2122static proc chineseRemainder(list rr, list mm)
2123"USAGE:  chineseRemainder(rr, mm); rr, mm lists of bigint's
2124ASSUME:  lists rr and mm must have same sizes;
2125         Furthermore the entries of mm must be mutually coprime.
2126RETURN:  an x which fulfills the simultaneous remainder conditions
2127         x = rr[i] mod mm[i], 1 <= i <= size(rr)
2128KEYWORDS: Chinese remainder.
2129"
2130{
2131  bigint x = bigint(0); int i; bigint N; list l;
2132  bigint M = bigint(mm[1]);
2133  for (i = 2; i <= size(mm); i++) { M = M * bigint(mm[i]); }
2134  for (i = 1; i <= size(mm); i++)
2135  {
2136    N = M div mm[i];
2137    l = extendedEuclid(mm[i], N);
2138    x = x + rr[i]*l[3]*N;
2139  }
2140  return (x);
2141}
2142
2143
2144///////////////////////////////////////////////////////////////////////////////
2145static proc rationalPointSpecial(bigint b1, bigint c1)
2146"USAGE:  rationalPointSpecial(b1, c1); b1, c1 bigint's
2147ASSUME:  b1 <> 0 and c1 <> 0;
2148RETURN:  with poly p = var(1)^2 + b1 * var(2)^2 + c1 * var(3)^2, the method
2149         returns a list L with either one entry or four entries:
2150         case 'three entries':
2151           L[1] = 0 signaling that there is no rational point on V(p),
2152           L[2] the largest number b such that b^2 divides b1
2153                (for subsequent use by the caller of this method),
2154           L[3] the largest number c such that c^2 divides c1
2155                (for subsequent use by the caller of this method);
2156         case 'four entries':
2157           L[1] = 1 signaling that there is a rational point on V(p),
2158           L[2], L[3], L[4] rational numbers such that the tuple
2159                (L[2], L[3], L[4]) is on V(p)
2160"
2161{
2162  if (b1 == 0) { "ERROR: rationalPointSpecial(0, c1) had been invoked"; }
2163  if (c1 == 0) { "ERROR: rationalPointSpecial(b1, 0) had been invoked"; }
2164
2165  bigint b_s = largestSquare(b1); bigint b_r = b1/b_s/b_s;
2166  bigint c_s = largestSquare(c1); bigint c_r = c1/c_s/c_s;
2167  bigint g = gcd(b_r, c_r);
2168  def S=basering;
2169  ideal mi = maxideal(1);
2170  map mm = basering, mi; map mTemp;
2171  mm[1] = var(1); mm[2] = var(2)/b_s/g; mm[3] = var(3)/c_s/g;
2172  bigint a = g;     bigint aa = a; if (aa <= 0) { aa = -aa; }
2173  bigint b = b_r/g; bigint bb = b; if (bb <= 0) { bb = -bb; }
2174  bigint c = c_r/g; bigint cc = c; if (cc <= 0) { cc = -cc; }
2175  bigint R1 = squareRoot(-a*b, cc, 1);
2176  if (R1 == -1) { list L = 0, b_s, c_s; return (L); }
2177  bigint R2 = squareRoot(-a*c, bb, 1);
2178  if (R2 == -1) { list L = 0, b_s, c_s; return (L); }
2179  bigint R3 = squareRoot(-b*c, aa, 1);
2180  if (R3 == -1) { list L = 0, b_s, c_s; return (L); }
2181  bigint t; bigint r1; bigint Q; bigint A; bigint B; bigint C;
2182  bigint alpha; bigint beta; bigint gamma;
2183  while (jIndex(a, b, c) > 1)
2184  {
2185    mTemp = basering, mi;
2186    if (aa > cc)
2187    {
2188      t = a; a = c; c = t;
2189      t = aa; aa = cc; cc = t;
2190      mTemp = basering, mi;
2191      mTemp[1] = var(3); mTemp[3] = var(1); mm = mTemp(mm);
2192    }
2193    if (bb > cc)
2194    {
2195      t = b; b = c; c = t;
2196      t = bb; bb = cc; cc = t;
2197      mTemp = basering, mi;
2198      mTemp[2] = var(3); mTemp[3] = var(2); mm = mTemp(mm);
2199    }
2200    if (bb < aa)
2201    {
2202      t = b; b = a; a = t;
2203      t = bb; bb = aa; aa = t;
2204      mTemp = basering, mi;
2205      mTemp[1] = var(2); mTemp[2] = var(1); mm = mTemp(mm);
2206    }
2207    /* now, we have established |a| <= |b| <= |c|; and permuted
2208       the map mm, accordingly */
2209    cc = c; if (cc <= 0) { cc = -cc; }
2210    R1 = squareRoot(-a*b, cc, 0);
2211    r1 = (R1 * inverseModP(a, cc)) mod cc;
2212    if (r1*bigint(2) > cc) { r1 = r1 - cc; }
2213    Q = (a*r1*r1 + b)/c;
2214    if (Q == 0)
2215    {
2216      list L = 1, subst(mm[1], var(1), 1, var(2), 1, var(3), 0),
2217                  subst(mm[2], var(1), 1, var(2), 1, var(3), 0),
2218                  subst(mm[3], var(1), 1, var(2), 1, var(3), 0);
2219      return (L);
2220    }
2221    A = gcd(gcd(a*r1*r1, b), c*Q);
2222    alpha = r1/A; beta = b/A;
2223    B = a*beta;
2224    gamma = largestSquare(Q/A);
2225    C = Q/A/gamma/gamma;
2226    mTemp = basering, mi;
2227    mTemp[1] = A*alpha*var(1) - beta*var(2);
2228    mTemp[2] = var(1) + a*alpha*var(2);
2229    mTemp[3] = C*gamma*var(3);
2230    mm = mTemp(mm);
2231    a = A; b = B; c = C;
2232    aa = a; if (aa <= 0) { aa = -aa; }
2233    bb = b; if (bb <= 0) { bb = -bb; }
2234    cc = c; if (cc <= 0) { cc = -cc; }
2235  }
2236  if (a*b < 0)
2237  {
2238    list L = 1, subst(mm[1], var(1), 1, var(2), 1, var(3), 0),
2239                subst(mm[2], var(1), 1, var(2), 1, var(3), 0),
2240                subst(mm[3], var(1), 1, var(2), 1, var(3), 0);
2241    return (L);
2242  }
2243  if (a*c < 0)
2244  {
2245    list L = 1, subst(mm[1], var(1), 1, var(2), 0, var(3), 1),
2246                subst(mm[2], var(1), 1, var(2), 0, var(3), 1),
2247                subst(mm[3], var(1), 1, var(2), 0, var(3), 1);
2248    return (L);
2249  }
2250  if (b*c < 0)
2251  {
2252    list L = 1, subst(mm[1], var(1), 0, var(2), 1, var(3), 1),
2253                subst(mm[2], var(1), 0, var(2), 1, var(3), 1),
2254                subst(mm[3], var(1), 0, var(2), 1, var(3), 1);
2255    return (L);
2256  }
2257  list L = 0, b_s, c_s; return (L);
2258}
2259
2260
2261///////////////////////////////////////////////////////////////////////////////
2262static proc extendedEuclid(bigint a, bigint b)
2263"USAGE:  extendedEuclid(a, b); a, b bigint's
2264ASSUME:  a <> 0 or b <> 0;
2265RETURN:  returns a list with three entries:
2266         _[1]: gcd(a,b) > 0,
2267         _[2], _[3]: s, t, such that s*a + t*b = gcd(a,b)
2268KEYWORDS: extended Euclidean algorithm.
2269"
2270{
2271  list l = 0; bigint temp;
2272  if (a == 0)         { l = b, 0, 1; if (b < 0) { l = -b, 0, -1; } }
2273  if (b == 0)         { l = a, 1, 0; if (a < 0) { l = -a, -1, 0; } }
2274  if ((a mod b) == 0) { l = b, 0, 1; if (b < 0) { l = -b, 0, -1; } }
2275  if ((b mod a) == 0) { l = a, 1, 0; if (a < 0) { l = -a, -1, 0; } }
2276  if (size(l) > 1) { return (l); }
2277
2278  temp = a mod b;
2279  l = extendedEuclid(b, temp);
2280  temp = (a - temp) / b;
2281  temp = bigint(l[2]) - temp * bigint(l[3]);
2282  l = l[1], l[3], temp;
2283  return (l);
2284}
2285
2286static proc legendreSymbol(bigint r, bigint p)
2287"assumes p prime;
2288returns the Legendre symbol (r/p), that is
2289 1 if r appears as residue modulo p of a square,
2290-1 if not,
2291 0 if r is a multiple of p
2292"
2293{
2294  bigint rr = r mod p;
2295  if (rr == 0) { return (0) }
2296  if (rr == 1) { return (1) }
2297  /* now, p must be at least 3 */
2298  bigint e = (p - 1) / bigint(2);
2299  bigint result = 1;
2300  bigint power = rr;
2301  while (e > 0)
2302  {
2303    if ((e mod 2) == 1) { result = (result * power) mod p; }
2304    e = e / bigint(2);
2305    power = (power * power) mod p;
2306  }
2307  if (result > 1) { result = result - p; /* should be -1 */ }
2308  return (result);
2309}
2310
2311
2312///////////////////////////////////////////////////////////////////////////////
2313static proc buildExtension(bigint b, bigint c, bigint bs, bigint cs)
2314"USAGE:  buildExtension(b, c, bs, cs); b, c, bs, cs bigint's
2315ASSUME:  assumes that bs is the largest positive number such that bs^2
2316         divides b; analogously for cs regarding c;
2317         Assumes furthermore that there is no rational point on the conic
2318         X^2 + b*Y^2 + c*Z^2 = 0.
2319         Assumes that the ground field of the basering is Q.
2320RETURN:  builds an appropriate quadratic field extension Q(a) in which a
2321         point exists that lies on the given conic. This point is stored in
2322         a newly defined and exported (1x3) matrix named 'point'.
2323         The method returns the resulting polynomial ring over Q(a).
2324"
2325{
2326  bigint br = b/bs/bs;
2327  bigint cr = c/cs/cs;
2328  /* X^2 + br*bs^2*Y^2 + cr*cs^2*Z^2 = 0 */
2329  def bRing = basering;
2330  list L = ringlist(bRing);
2331
2332  if (b != 0)
2333  {
2334    L[1] = list(0, list("a"), list(list("lp", 1)), ideal(0));
2335    def RTemp = ring(L);
2336    setring RTemp; list L = ringlist(RTemp);
2337    L[1][4] = ideal(a^2 + br);
2338    def R = ring(L);
2339    setring R; kill RTemp;
2340    matrix point[1][3];
2341    point[1, 1] = a * bs; point[1, 2] = 1; point[1, 3] = 0;
2342    export point;
2343    setring bRing;
2344    return (R);
2345  }
2346  if (c != 0)
2347  {
2348    L[1] = list(0, list("a"), list(list("lp", 1)), ideal(0));
2349    def RTemp = ring(L);
2350    setring RTemp; list L = ringlist(RTemp);
2351    L[1][4] = ideal(a^2 + cr);
2352    def R = ring(L);
2353    setring R; kill RTemp;
2354    matrix point[1][3];
2355    point[1, 1] = a * cs; point[1, 2] = 0; point[1, 3] = 1;
2356    export point;
2357    setring bRing;
2358    return (R);
2359  }
2360
2361  "ERROR: unexpectedly encountered conic X^2 + 0*Y^2 + 0*Z^2 = 0";
2362  return (bRing);
2363}
2364
2365
2366///////////////////////////////////////////////////////////////////////////////
2367static proc testRationalPointConic(poly pp)
2368"USAGE:  testRationalPointConic(pp); pp poly
2369RETURN:  returns 0 in case of unexpected input (e.g. non-quadratic,
2370         reducible); 1 otherwise
2371NOTE:    This method calles rationalPointConic, measures time consumption
2372         and checks whether the computed point lies indeed on the conic pp.
2373         The results will be printed to standard output.
2374"
2375{
2376  "testing rationalPointConic(poly) for poly p:";
2377  "p =", pp;
2378  if (isQuadratic(pp)[1] == 1) { "p is quadratic."; }
2379  else                         { "p is not quadratic.";   return (0); }
2380  if (isIrreducible(pp) == 1)  { "p is irreducible."; }
2381  else                         { "p is not irreducible."; return (0); }
2382  def rOrig = basering;
2383  int t = rtimer;
2384  def rNew = rationalPointConic(pp);
2385  t = rtimer - t;
2386  "time for finding a point on the conic [sec] =", t;
2387  setring rNew;
2388  poly ff = fetch(rOrig, pp);
2389  if (minpoly == 0)
2390  { "there is a rational point on the conic p";
2391    "x =", point[1,1], "  y =", point[1,2], "  z =", point[1,3];
2392    "check (should be zero):", subst(ff, var(1), point[1,1],
2393                                         var(2), point[1,2],
2394                                         var(3), point[1,3]);
2395  }
2396  else
2397  {
2398    "there is no rational point on the conic p";
2399    "but there is a point on the conic in the field extension Q(a),";
2400    "with minpoly =", minpoly;
2401    "x =", point[1,1], "  y =", point[1,2], "  z =", point[1,3];
2402    "check (should be zero):", subst(ff, var(1), point[1,1],
2403                                         var(2), point[1,2],
2404                                         var(3), point[1,3]);
2405  }
2406  setring rOrig;
2407}
2408
2409example
2410{ "EXAMPLE:"; echo=2;
2411  ring r = 0, (x,y,z, u, v, w), dp;
2412  poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z;
2413  testRationalPointConic(p);
2414}
2415
2416///////////////////////////////////////////////////////////////////////////////
2417proc rationalPointConic(poly p)
2418"USAGE:  rationalPointConic(p); p poly
2419ASSUME:  assumes that p is an irreducible quadratic polynomial in the first
2420         three ring variables;
2421         ground field is expected to be Q.
2422RETURN:  The method finds a point on the given conic. There are two
2423         possibilities:
2424         1) There is a rational point on the curve.
2425         2) There is no rational point on the curve.
2426         In the second case, the method creates a modification of the current
2427         basering which is a polynomial ring over some quadratic field
2428         extension Q(a) of Q. Apart from the replacement of Q by Q(a), the
2429         new polynomial ring, R say, is the same as the original basering.
2430         (In the first case, R is identical with the basering.)
2431         In both cases, the method will then define a (1x3) matrix named
2432         'point' which lives in R and which contains the coordinates of the
2433         desired point on q.
2434         Finally, the method returns the ring R (which will in the 1st case
2435         be the original base ring).
2436EXAMPLE: example rationalPointConic; shows an example
2437"
2438{
2439  list L = isQuadratic(p);
2440  bigint a = bigint(L[2]); bigint b = bigint(L[3]); bigint c = bigint(L[4]);
2441  bigint d = bigint(L[5]); bigint e = bigint(L[6]); bigint f = bigint(L[7]);
2442  bigint x; bigint y; bigint z; bigint nn;
2443  def R = basering;
2444
2445  if (b^2 == 4*a*c)
2446  {
2447    if (c == 0)
2448    {
2449      x = -2*d*e; y = d^2-4*a*f; z = e*4*a;
2450      nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
2451      matrix point[1][3];
2452      point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
2453      export point; return (R);
2454    }
2455    else
2456    {
2457      bigint fs = 4*c*f - e^2;
2458      bigint ds = 4*c*d - 2*b*e;
2459      x = -fs*2*c; y = b*fs-e*ds; z = ds*2*c;
2460      nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
2461      matrix point[1][3];
2462      point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
2463      export point; return (R);
2464    }
2465  }
2466
2467  if (d^2 == 4*a*f)
2468  {
2469    if (f == 0)
2470    {
2471      x = -b*e*2; y = e*4*a; z = b^2-4*a*c;
2472      nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
2473      matrix point[1][3];
2474      point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
2475      export point; return (R);
2476    }
2477    else
2478    {
2479      bigint c_s = 4*c*f - e^2;
2480      bigint b_s = 4*f*b - 2*d*e;
2481      x = -c_s*2*f; y = b_s*2*f; z = d*c_s-e*b_s;
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  }
2488
2489  if (e^2 == 4*c*f)
2490  {
2491    if (c == 0)
2492    {
2493      x = b*4*f; y = d^2-4*a*f; z = -b*d*2;
2494      nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
2495      matrix point[1][3];
2496      point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
2497      export point; return (R);
2498    }
2499    else
2500    {
2501      bigint as = 4*c*a - b^2;
2502      bigint ds = 4*c*d - 2*b*e;
2503      x = ds*2*c; y = e*as-b*ds; z = -as*2*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  }
2510
2511  ideal mi = maxideal(1);
2512  map mm = R, mi;
2513  bigint B; bigint C; bigint D;
2514
2515  if ((a == 0) && (c == 0))
2516  {
2517    B = -1; C = 4*b*f - 4*d*e;
2518    /* now, b <> 0 since otherwise p would have the factor z,
2519       and hence not be irreducible */
2520    mm[1] = (var(1)+var(2)-2*e*var(3))/(2*b);
2521    mm[2] = (var(1)-var(2)-2*d*var(3))/(2*b);
2522  }
2523  if ((a != 0) && (c == 0))
2524  {
2525    mm[1] = var(2);
2526    mm[2] = var(1);
2527    bigint t = a; a = c; c = t;
2528    t = e; e = d; d = t;
2529  }
2530  if (c != 0)
2531  {
2532    D = 4*a*c-b^2;
2533    mm[2] = (var(2)-e*var(3)-b*var(1))/(2*c);
2534    map mTemp = basering, mi;
2535    mTemp[1] = (var(1)-2*d*c*var(3)+b*e*var(3))/D;
2536    mm = mTemp(mm);
2537    B = D;
2538    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;
2539  }
2540  list K;
2541  if ((B > 0) && (C >= 0)) { K = 0; }
2542  if ((B >= 0) && (C > 0)) { K = 0; }
2543  if (B == 0)
2544  {
2545    /* looking for a point on X^2 = |C| * Z^2 */
2546    bigint root = largestSquare(absValue(C));
2547    if (absValue(C)/root/root == 1) { K = 1, root, 0, 1; }
2548    else                            { K = 0; }
2549  }
2550  if (C == 0)
2551  {
2552    /* looking for a point on X^2 = |B| * Y^2 */
2553    bigint root = largestSquare(absValue(B));
2554    if (absValue(B)/root/root == 1) { K = 1, root, 1, 0; }
2555    else                            { K = 0; }
2556  }
2557  else { K = rationalPointSpecial(B, C); }
2558  if (K[1] == 0)
2559  {
2560    /* no rational point on conic;
2561       we need to move to an appropriate field extension Q(a) */
2562    poly h1 = mm[1]; poly h2 = mm[2]; poly h3 = mm[3];
2563    def extendedR = buildExtension(B, C, K[2], K[3]);
2564    setring extendedR;
2565    poly g1 = fetch(R, h1);
2566    poly g2 = fetch(R, h2);
2567    poly g3 = fetch(R, h3);
2568    matrix temp[1][3];
2569    temp[1, 1] = subst(g1, var(1), point[1, 1], var(2), point[1, 2],
2570                           var(3), point[1, 3]);
2571    temp[1, 2] = subst(g2, var(1), point[1, 1], var(2), point[1, 2],
2572                           var(3), point[1, 3]);
2573    temp[1, 3] = subst(g3, var(1), point[1, 1], var(2), point[1, 2],
2574                           var(3), point[1, 3]);
2575    point[1, 1] = temp[1, 1]; point[1, 2] = temp[1, 2];
2576    point[1, 3] = temp[1, 3];
2577    setring R;
2578    return (extendedR);
2579  }
2580  else
2581  {
2582    string dummyString = string(K); // without this useless line, we
2583                                    // sometimes get a seg fault because
2584                                    // mm is corrupted; strange!?!?!?!?
2585    number nx = number(subst(mm[1], var(1), K[2], var(2), K[3], var(3), K[4]));
2586    number ny = number(subst(mm[2], var(1), K[2], var(2), K[3], var(3), K[4]));
2587    number nz = number(subst(mm[3], var(1), K[2], var(2), K[3], var(3), K[4]));
2588    /* the point (nx, ny, nz) is already a solution;
2589       the following lines will just remove denominators and reduce
2590       numerators in order to return a nice tuple from Z^3 */
2591    bigint nxd = bigint(denominator(absValue(nx)));
2592    bigint nyd = bigint(denominator(absValue(ny)));
2593    bigint nzd = bigint(denominator(absValue(nz)));
2594    nn = nxd * nyd / gcd(nxd, nyd);
2595    nn =  nn * nzd / gcd(nn, nzd);
2596    x = bigint(nx*nn); y = bigint(ny*nn); z = bigint(nz*nn);
2597    nn = gcd(gcd(absValue(x), absValue(y)), absValue(z));
2598    matrix point[1][3];
2599    point[1, 1] = x/nn; point[1, 2] = y/nn; point[1, 3] = z/nn;
2600    export point;
2601    return (R);
2602  }
2603}
2604
2605example
2606{ "EXAMPLE:"; echo=2;
2607ring R = 0, (x,y,z), dp;
2608system("random", 4711);
2609poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z;
2610def S = rationalPointConic(p); // quadratic field extension,
2611                               // minpoly = a^2 - 2
2612testPointConic(p, S);
2613setring R;
2614p = x^2 - 1857669520 * y^2 + 86709575222179747132487270400 * z^2;
2615S = rationalPointConic(p); // same as current basering,
2616                           // no extension needed
2617testPointConic(p, S);
2618}
2619///////////////////////////////////////////////////////////////////////////////
2620proc testParametrization(poly f, def rTT)
2621"USAGE:  testParametrization(f, rTT); f poly, rTT ring
2622ASSUME:  The assumptions on the basering and the polynomial f are as required
2623         by @ref{paraPlaneCurve}. The ring rTT has two variables and contains
2624         an ideal PARA (such as the ring obtained by applying
2625         @ref{paraPlaneCurve} to f).
2626RETURN: int which is 1 if PARA defines a parametrization of the curve
2627        {f=0} and 0, otherwise.
2628THEORY: We compute the polynomial defining the image of PARA
2629        and compare it with f.
2630KEYWORDS: Parametrization, image.
2631EXAMPLE: example testParametrization; shows an example
2632"
2633{
2634  def Roriginal = basering;
2635  setring rTT;
2636  /* begin workaround elimination*/
2637  int k;
2638  list rl = ringlist(rTT);
2639  rl[2] = list("s","t","x","y","z");
2640  rl[3]= list(list("dp",1:5),list("C",0));
2641  def Relim = ring(rl);
2642  setring Relim;
2643  ideal PARA = fetch(rTT,PARA);
2644  ideal JJ;
2645  for(k=1;k<=3;k++)
2646     {
2647       JJ=JJ,var(k+2)-PARA[k];
2648     }
2649  ideal SJJ = std(JJ);
2650  intvec HJJ = hilb(SJJ,1);
2651  ideal J = eliminate(JJ,var(1)*var(2),HJJ);
2652  setring rTT;
2653  /*end workaround elimination*/
2654  rl[2] = list("x","y","z");
2655  rl[3] = list(list("dp",1:3),list("C",0));
2656  def RP2 = ring(rl);
2657  setring RP2;
2658  ideal f = fetch(Roriginal,f);
2659  ideal ftest = imap(Relim,J);
2660  poly g = reduce(f[1],std(ftest));
2661  if(g!=0){return(0)}
2662  g = reduce(ftest[1],std(ideal(f)));
2663  if(g!=0){return(0)}
2664  return (1);
2665}
2666
2667example
2668{ "EXAMPLE:"; echo=2;
2669  ring R = 0,(x,y,z),dp;
2670  poly f = y^8-x^3*(z+x)^5;
2671  def RP1 = paraPlaneCurve(f);
2672  testParametrization(f, RP1);
2673}
2674
2675///////////////////////////////////////////////////////////////////////////////
2676proc testPointConic(poly p, def r)
2677"USAGE:  testPointConic(p, r); p poly, r ring
2678ASSUME:  assumes that p is a homogeneous quadratic polynomial in the
2679        first three ring variables of the current basering;
2680        Assumes that there is a (1x3) matrix named 'point' in r with
2681        entries from the ground field of r.
2682RETURN:  returns 1 iff the point named 'point', residing in r, lies on
2683        the conic given by p; 0 otherwise
2684NOTE:    This method temporarily changes the basering to r. Afterwards,
2685        the basering will be the same as before.
2686EXAMPLE: example testPointConic; shows an example
2687"
2688{
2689 def rOrig = basering;
2690 "conic:", p;
2691 setring r;
2692 string s = "point: " + string(point[1,1]) + ", " + string(point[1,2]);
2693 s = s + ", " + string(point[1,3]);
2694 s;
2695 if (minpoly != 0) { "minpoly:", minpoly; }
2696 poly f = fetch(rOrig, p);
2697 poly g = subst(f, var(1), point[1,1],
2698                   var(2), point[1,2],
2699                   var(3), point[1,3]);
2700 int result = 0; if (g == 0) { result = 1; }
2701 setring rOrig;
2702 return (result);
2703}
2704
2705example
2706{ "EXAMPLE:"; echo=2;
2707 ring R = 0, (x,y,z), dp;
2708 system("random", 4711);
2709 poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z;
2710 def S = rationalPointConic(p);
2711 if (testPointConic(p, S) == 1)
2712 { "point lies on conic"; }
2713 else
2714 { "point does not lie on conic"; }
2715}
2716
2717/////////////////////////////////////////////////////////////////////////////
2718/////////////////////////////////////////////////////////////////////////////
2719/////////////////////////////////////////////////////////////////////////////
2720/*
2721/////////////////////////////////////////////////////////////////////////////
2722/// Further examples for testing the main procedures
2723/// Timings on wawa Sept 29
2724/////////////////////////////////////////////////////////////////////////////
2725LIB"paraplanecurves.lib";
2726// -------------------------------------------------------
2727// Example 1
2728// -------------------------------------------------------
2729ring RR = 0, (x,y,z), dp;
2730poly f = 7*z2+11*y2+13*z*y+17*x2+19*x*y; // conic
2731def RP1 = paraConic(f);
2732setring RP1; PARACONIC;
2733setring RR;
2734RP1 = paraPlaneCurve(f);
2735testParametrization(f,RP1);
2736setring RP1; PARA;
2737kill RR;kill RP1;
2738// -------------------------------------------------------
2739// Example 2
2740// -------------------------------------------------------
2741ring RR = 0, (x,y,z), dp;
2742poly f = y3-x2z;  // cusp at origin
2743adjointIdeal(f,1);
2744adjointIdeal(f,2);
2745def RP1 = paraPlaneCurve(f);  // time 0
2746testParametrization(f,RP1);
2747setring RP1; PARA;
2748kill RR;kill RP1;
2749// -------------------------------------------------------
2750// Example 3
2751// -------------------------------------------------------
2752ring RR = 0, (x,y,z), dp;
2753poly f=(xz-y^2)^2-x*y^3; // 1 sing at origin, 1 cusp, no OMPs
2754adjointIdeal(f,1);
2755adjointIdeal(f,2);
2756def RP1 = paraPlaneCurve(f); // time 0
2757testParametrization(f,RP1);
2758setring RP1; PARA;
2759kill RR;kill RP1;
2760// -------------------------------------------------------
2761// Example 4
2762// -------------------------------------------------------
2763ring RR = 0, (x,y,z), dp;
2764poly f = y5-y4x+4y2x2z-x4z;  // 1 sing at origin, no OMPs, no cusps
2765adjointIdeal(f,1);
2766adjointIdeal(f,2);
2767def RP1 = paraPlaneCurve(f);  // time 0
2768testParametrization(f,RP1);
2769setring RP1; PARA;
2770kill RR;kill RP1;
2771// -------------------------------------------------------
2772// Example 5
2773// -------------------------------------------------------
2774ring RR = 0, (x,y,z), dp;
2775poly f = 259x5-31913x4y+939551x3y2+2871542x2y3+2845801xy4;
2776f = f+914489y5+32068x4z-1884547x3yz-8472623x2y2z-11118524xy3z;
2777f = f-4589347y4z+944585x3z2+8563304x2yz2+16549772xy2z2+9033035y3z2;
2778f = f-2962425x2z3-11214315xyz3-8951744y2z3+2937420xz4+4547571yz4-953955z5;
2779// 6 nodes
2780adjointIdeal(f,1);
2781adjointIdeal(f,2);
2782def RP1 = paraPlaneCurve(f);  // time 0
2783testParametrization(f,RP1);
2784setring RP1; PARA;
2785kill RR;kill RP1;
2786// -------------------------------------------------------
2787// Example 7
2788// -------------------------------------------------------
2789ring RR = 0, (x,y,z), dp;
2790poly f = y^8-x^3*(z+x)^5;  // 1 sing at origin, 1 further sing, no OMPs,
2791                           // no cusps
2792adjointIdeal(f,1);
2793adjointIdeal(f,2);
2794def RP1 = paraPlaneCurve(f);  // time 0
2795testParametrization(f,RP1);
2796setring RP1; PARA;
2797kill RR;kill RP1;
2798// -------------------------------------------------------
2799// Example 8
2800// -------------------------------------------------------
2801ring RR = 0, (x,y,z), dp;
2802poly f = 11y7+7y6x+8y5x2-3y4x3-10y3x4-10y2x5-x7-33y6-29y5x-13y4x2+26y3x3;
2803f = f+30y2x4+10yx5+3x6+33y5+37y4x-8y3x2-33y2x3-20yx4-3x5-11y4-15y3x;
2804f = f+13y2x2+10yx3+x4; // 3 OMPs of mult 3, 1 OMP of mult 4
2805f = homog(f,z);
2806adjointIdeal(f,1);
2807adjointIdeal(f,2);
2808def RP1 = paraPlaneCurve(f);  // time 0
2809testParametrization(f,RP1);
2810setring RP1; PARA;
2811kill RR;kill RP1;
2812// -------------------------------------------------------
2813// Example 9
2814// -------------------------------------------------------
2815ring RR = 0, (x,y,z), dp;
2816poly f = y^8-x^3*(z+x)^5;  // 1 sing at origin, 1 further sing, no OMPs,
2817                           // no cusps
2818adjointIdeal(f,1);
2819adjointIdeal(f,2);
2820def RP1 = paraPlaneCurve(f);  // time 0
2821testParametrization(f,RP1);
2822setring RP1; PARA;
2823kill RR;kill RP1;
2824// -------------------------------------------------------
2825// Example 10
2826// -------------------------------------------------------
2827ring SS = 0, (u,v,z), dp;
2828poly 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
2829adjointIdeal(f,1);
2830adjointIdeal(f,2);
2831def RP1 = paraPlaneCurve(f);  // time 0
2832testParametrization(f,RP1);
2833setring RP1; PARA;
2834kill SS;kill RP1;
2835// -------------------------------------------------------
2836// Example 11
2837// -------------------------------------------------------
2838ring SS = 0, (u,v,z), dp;
2839poly f = 14440*u^5-16227*u^4*v+10812*u^3*v^2-13533*u^2*v^3+3610*u*v^4;
2840f = f+1805*v^5+14440*u^4*z-18032*u^3*v*z+16218*u^2*v^2*z-12626*u*v^3*z;
2841f = 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;
2842// 1 OMP of mult 3 at origin, 2 nodes
2843adjointIdeal(f,1);
2844adjointIdeal(f,2);
2845def RP1 = paraPlaneCurve(f);  // time 0
2846testParametrization(f,RP1);
2847setring RP1; PARA;
2848kill SS;kill RP1;
2849// -------------------------------------------------------
2850// Example 12
2851// -------------------------------------------------------
2852ring SS = 0, (u,v,z), dp;
2853poly 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;
2854f = 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;
2855// needs field extension *** 6 nodes, 2 cusps, 1 sing at 0
2856adjointIdeal(f,1);
2857adjointIdeal(f,2);
2858def RP1 = paraPlaneCurve(f);  // time 0
2859testParametrization(f,RP1);
2860setring RP1; PARA;
2861kill SS;kill RP1;
2862// -------------------------------------------------------
2863// Example 13
2864// -------------------------------------------------------
2865ring SS = 0, (u,v,z), dp;
2866poly f = -24135/322*u^6-532037/6440*u^5*v+139459/560*u^4*v^2;
2867f = f-1464887/12880*u^3*v^3+72187/25760*u^2*v^4+9/8*u*v^5+1/8*v^6;
2868f = f-403511/3220*u^5*z-40817/920*u^4*v*z+10059/80*u^3*v^2*z;
2869f = 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;
2870f = f+126379/3220*u^3*v*z^2-423417/6440*u^2*v^2*z^2+11/2*u*v^3*z^2;
2871f = 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;
2872// 2 OMPs of mult 3 (1 at origin), 4 nodes
2873adjointIdeal(f,1);
2874adjointIdeal(f,2);
2875def RP1 = paraPlaneCurve(f);  // time 14
2876testParametrization(f,RP1);
2877setring RP1; PARA;
2878kill SS;kill RP1;
2879// -------------------------------------------------------
2880// Example 14
2881// -------------------------------------------------------
2882ring SS = 0, (u,v,z), dp;
2883poly f =
28842*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
2885-7780247/995328*u^6*z-78641/9216*u^5*v*z-10892131/995328*u^4*v^2*z
2886-329821/31104*u^3*v^3*z-953807/331776*u^2*v^4*z-712429/248832*u*v^5*z
2887+1537741/331776*v^6*z+2340431/248832*u^5*z^2+5154337/248832*u^4*v*z^2
2888+658981/41472*u^3*v^2*z^2+1737757/124416*u^2*v^3*z^2
2889-1234733/248832*u*v^4*z^2-1328329/82944*v^5*z^2-818747/248832*u^4*z^3
2890-1822879/124416*u^3*v*z^3-415337/31104*u^2*v^2*z^3
2891+1002655/124416*u*v^3*z^3+849025/82944*v^4*z^3;
2892// 3 OMPs of mult 3, 1 OMP of mult 4 at origin
2893adjointIdeal(f,2);
2894def RP1 = paraPlaneCurve(f);  // time 1
2895testParametrization(f,RP1);
2896setring RP1; PARA;
2897kill SS;kill RP1;
2898// -------------------------------------------------------
2899// Example 15
2900// -------------------------------------------------------
2901ring SS = 0, (u,v,z), dp;
2902poly f = 590819418867856650536224u7-147693905508217596067968u6v;
2903f = f+229117518934972047619978u5v2-174050799674982973889542u4v3;
2904f = f-92645796479789150855110u3v4-65477418713685583062704u2v5;
2905f = f+4529961835917468460168uv6+7715404057796585983136v7;
2906f = f-413640780091141905428104u6z+571836835577486968144618u5vz;
2907f = f-551807810327826605739444u4v2z-488556410340789283359926u3v3z;
2908f = f-473466023008413178155962u2v4z+48556741573432247323608uv5z;
2909f = f+77647371229172269259528v6z+340450118906560552282893u5z2;
2910f = f-433598825064368371610344u4vz2-937281070591684636591672u3v2z2;
2911f = f-1388949843915129934647751u2v3z2+204081793110898617103998uv4z2;
2912f = f+335789953068251652554308v5z2+6485661002496681852577u4z3;
2913f = f-772700266516318390630202u3vz3-2068348417248100329533330u2v2z3;
2914f = f+440320154612359641806108uv3z3+808932515589210854581618v4z3;
2915f = f-229384307132237615286548u3z4-1564303565658228216055227u2vz4;
2916f = f+520778334468674798322974uv2z4+1172483905704993294097655v3z4;
2917f = f-480789741398016816562100u2z5+322662751598958620410786uvz5;
2918f = f+1022525576391791616258310v2z5+82293493608853837667471uz6;
2919f = f+496839109904761426785889vz6+103766136235628614937587z7; // 15 nodes
2920adjointIdeal(f,2);
2921def RP1 = paraPlaneCurve(f);  // time 72
2922testParametrization(f,RP1);
2923setring RP1; PARA;
2924kill SS;kill RP1;
2925
2926// -------------------------------------------------------
2927// Example 16
2928// -------------------------------------------------------
2929ring SS = 0, (u,v,z), dp;
2930poly 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;
2931f = 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;
2932f = 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;
2933f = 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;
2934f = 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;
2935f = 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;
2936f = 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;
2937f = 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;
2938f = f+28*v^2*z^6+18*u*z^7+8*v*z^7+z^8; // 3 nodes, 1 sing
2939adjointIdeal(f,1);
2940adjointIdeal(f,2);
2941def RP1 = paraPlaneCurve(f);  // time 20
2942testParametrization(f,RP1);
2943setring RP1; PARA;
2944kill SS;kill RP1;
2945// -------------------------------------------------------
2946// Example 17
2947// -------------------------------------------------------
2948ring SS = 0, (u,v,z), dp;
2949poly 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;
2950f = 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;
2951f = 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;
2952f = f+v^5*z^4; // 2 OMPs of mult 4, 1 OMP of mult 5, 1 sing at origin
2953f = subst(f,z,u+v+z);
2954adjointIdeal(f,2);
2955def RP1 = paraPlaneCurve(f);  // time 5
2956testParametrization(f,RP1);
2957setring RP1; PARA;
2958kill SS;kill RP1;
2959// -------------------------------------------------------
2960// Example 18
2961// -------------------------------------------------------
2962ring SS = 0, (u,v,z), dp;
2963poly 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;
2964f = 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;
2965f = 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;
2966f = f+u*v^4*z^5+v^5*z^5; // 1 OMP of mult 4, 3 OMPs of mult 5 (1 at origin)
2967adjointIdeal(f,2);
2968def RP1 = paraPlaneCurve(f);  // time 8
2969testParametrization(f,RP1);
2970setring RP1; PARA;
2971kill SS;kill RP1;
2972// -------------------------------------------------------
2973// Example 19
2974// -------------------------------------------------------
2975ring SS = 0, (u,v,z), dp;
2976poly 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;
2977f = 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;
2978f = 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;
2979f = 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;
2980f = 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;
2981f = 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;
2982f = 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;
2983f = 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;
2984f = f-15*u^3*v^2*z^5+u^2*v^3*z^5+u*v^4*z^5+v^5*z^5;
2985// 1 OMP of mult 4, 3 OMPs of mult 5 (1 at origin)
2986adjointIdeal(f,2);
2987def RP1 = paraPlaneCurve(f);  // time 2 // see Ex. 18
2988testParametrization(f,RP1);
2989setring RP1; PARA;
2990kill SS;kill RP1;
2991// -------------------------------------------------------
2992// Example 20
2993// -------------------------------------------------------
2994ring R = 0, (x,y,z), dp;
2995system("random", 4711);
2996poly p = x^2 + 2*y^2 + 5*z^2 - 4*x*y + 3*x*z + 17*y*z;
2997def S = rationalPointConic(p); // quadratic field extension,
2998                              // minpoly = a^2 - 2
2999if (testPointConic(p, S) == 1)
3000{ "point lies on conic"; }
3001else
3002{ "point does not lie on conic"; }
3003kill R;kill S;
3004// -------------------------------------------------------
3005// Example 21
3006// -------------------------------------------------------
3007ring R = 0, (x,y,z), dp;
3008system("random", 4711);
3009poly p = x^2 - 1857669520 * y^2 + 86709575222179747132487270400 * z^2;
3010def S = rationalPointConic(p); // same as current basering,
3011                              // no extension needed
3012if (testPointConic(p, S) == 1)
3013{ "point lies on conic"; }
3014else
3015{ "point does not lie on conic"; }
3016kill R;kill S;
3017// -------------------------------------------------------
3018// Example 21
3019// -------------------------------------------------------
3020ring RR = 0, (x,y,z), dp;
3021poly f = -1965466244509920x5y+34871245546721380061760x4y2;
3022f = f+104613747941595046117320x3y3+113331564241941002407560x2y4;
3023f = f+52306876673313609259800xy5+8717812860780028397880y6;
3024f = f+1040297748510024x5z+4468147845634872x4yz;
3025f = f-22398508728211453743258x3y2z-33223996581074443306854x2y3z;
3026f = f-10638598235041298082366xy4z+186886189971594356382y5z;
3027f = f-1385078844909312x4z2-34893092731637052532683x3yz2;
3028f = f-98591463214095439056609x2y2z2-92339459334829609336485xy3z2;
3029f = f-24923289542522905755711y4z2+472440640471377x3z3;
3030f = f+33821511925664516716011x2yz3+49745237303968344397437xy2z3;
3031f = f+11040465960074786720475y3z3+8728735735878837099404x2z4;
3032f = f+17676785754519678518537xyz4+17935885079051421934609y2z4;
3033f = f-11314701999743172607075xz5-16164284825803158969425yz5;
3034f = f+3666695988537425618750z6;
3035// 4 nodes, 1 OMP of mult 4
3036adjointIdeal(f,2);
3037kill RR;
3038*/
Note: See TracBrowser for help on using the repository browser.