source: git/Singular/LIB/hnoether.lib @ bf42f0

spielwiese
Last change on this file since bf42f0 was bf42f0, checked in by Hans Schoenemann <hannes@…>, 13 years ago
format git-svn-id: file:///usr/local/Singular/svn/trunk@13483 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 153.1 KB
Line 
1version="$Id$";
2// (GMG, last modified 08.06.08: comments to Newton non-degenerate added)
3///////////////////////////////////////////////////////////////////////////////
4
5category="Singularities";
6info="
7LIBRARY:  hnoether.lib   Hamburger-Noether (Puiseux) Expansion
8AUTHORS:   Martin Lamm,      lamm@mathematik.uni-kl.de
9           Christoph Lossen, lossen@mathematik.uni-kl.de
10
11OVERVIEW:
12 A library for computing the Hamburger-Noether expansion (analogue of
13 Puiseux expansion over fields of arbitrary characteristic) of a reduced
14 plane curve singularity following [Campillo, A.: Algebroid curves in
15 positive characteristic, Springer LNM 813 (1980)]. @*
16 The library contains also procedures for computing the (topological)
17 numerical invariants of plane curve singularities.
18
19PROCEDURES:
20 hnexpansion(f [,\"ess\"]); Hamburger-Noether (HN) expansion of f
21 develop(f [,n]);           HN expansion of irreducible plane curve germs
22 extdevelop(hne,n);         extension of the H-N expansion hne of f
23 param(hne [,s]);           parametrization of branches described by HN data
24 displayHNE(hne);           display HN expansion as an ideal
25 invariants(hne);           invariants of f, e.g. the characteristic exponents
26 displayInvariants(hne);    display invariants of f
27 multsequence(hne);         sequence of multiplicities
28 displayMultsequence(hne);  display sequence of multiplicities
29 intersection(hne1,hne2);   intersection multiplicity of two local branches
30 is_irred(f);               test whether f is irreducible as power series
31 delta(f);                  delta invariant of f
32 newtonpoly(f);             (local) Newton polygon of f
33 is_NND(f);                 test whether f is Newton non-degenerate
34
35
36 stripHNE(hne);             reduce amount of memory consumed by hne
37 puiseux2generators(m,n);   convert Puiseux pairs to generators of semigroup
38 separateHNE(hne1,hne2);    number of quadratic transf. needed for separation
39 squarefree(f);             a squarefree divisor of the polynomial f
40 allsquarefree(f,l);        the maximal squarefree divisor of the polynomial f
41 further_hn_proc();         show further procedures useful for interactive use
42
43KEYWORDS: Hamburger-Noether expansion; Puiseux expansion; curve singularities
44";
45
46// essdevelop(f);             HN expansion of essential branches
47// multiplicities(hne);       multiplicities of blowed up curves
48
49///////////////////////////////////////////////////////////////////////////////
50LIB "primitiv.lib";
51LIB "inout.lib";
52LIB "sing.lib";
53
54///////////////////////////////////////////////////////////////////////////////
55
56proc further_hn_proc()
57"USAGE: further_hn_proc();
58NOTE:  The library @code{hnoether.lib} contains some more procedures which
59       are not shown when typing @code{help hnoether.lib;}. They may be useful
60       for interactive use (e.g. if you want to do the calculation of an HN
61       development \"by hand\" to see the intermediate results), and they
62       can be enumerated by calling @code{further_hn_proc()}. @*
63       Use @code{help <procedure>;} for detailed information about each of
64       them.
65"
66{
67 "
68 The following procedures are also part of `hnoether.lib':
69
70 getnm(f);           intersection pts. of Newton polygon with axes
71 T_Transform(f,Q,N); returns f(y,xy^Q)/y^NQ (f: poly, Q,N: int)
72 T1_Transform(f,d,M); returns f(x,y+d*x^M)  (f: poly,d:number,M:int)
73 T2_Transform(f,d,M,N,ref);   a composition of T1 & T
74 koeff(f,I,J);       gets coefficient of indicated monomial of polynomial f
75 redleit(f,S,E);     restriction of monomials of f to line (S-E)
76 leit(f,n,m);        special case of redleit (for irred. polynomials)
77 testreducible(f,n,m); tests whether f is reducible
78 charPoly(f,M,N);    characteristic polynomial of f
79 find_in_list(L,p);  find int p in list L
80 get_last_divisor(M,N); last divisor in Euclid's algorithm
81 factorfirst(f,M,N); try to factor f without `factorize'
82 factorlist(L);      factorize a list L of polynomials
83 referencepoly(D);   a polynomial f s.t. D is the Newton diagram of f";
84
85//       static procedures not useful for interactive use:
86// polytest(f);        tests coefficients and exponents of polynomial f
87// extractHNEs(H,t);   extracts output H of HN to output of hnexpansion
88// HN(f,grenze);       recursive subroutine for hnexpansion
89// constructHNEs(...); subroutine for HN
90}
91example
92{ echo=2;
93  further_hn_proc();
94}
95///////////////////////////////////////////////////////////////////////////////
96
97proc getnm (poly f)
98"USAGE:   getnm(f); f bivariate polynomial
99RETURN:  intvec(n,m) : (0,n) is the intersection point of the Newton
100         polygon of f with the y-axis, n=-1 if it doesn't exist
101         (m,0) is its intersection point with the x-axis,
102         m=-1 if this point doesn't exist
103ASSUME:  ring has ordering `ls' or `ds'
104EXAMPLE: example getnm; shows an example
105"
106{
107 // assume being called by develop ==> ring ordering is ls (ds would also work)
108 return(ord(subst(f,var(1),0)),ord(subst(f,var(2),0)));
109}
110example
111{ "EXAMPLE:"; echo = 2;
112   ring r = 0,(x,y),ds;
113   poly f = x5+x4y3-y2+y4;
114   getnm(f);
115}
116///////////////////////////////////////////////////////////////////////////////
117
118proc leit (poly f, int n, int m)
119"USAGE:   leit(f,n,m);  poly f, int n,m
120RETURN:  all monomials on the line from (0,n) to (m,0) in the Newton diagram
121EXAMPLE: example leit;  shows an example
122"
123{
124 return(jet(f,m*n,intvec(n,m))-jet(f,m*n-1,intvec(n,m)))
125}
126example
127{ "EXAMPLE:"; echo = 2;
128   ring r = 0,(x,y),ds;
129   poly f = x5+x4y3-y2+y4;
130   leit(f,2,5);
131}
132///////////////////////////////////////////////////////////////////////////////
133proc testreducible (poly f, int n, int m)
134"USAGE:   testreducible(f,n,m);  f poly, n,m int
135RETURN:  1 if there are points in the Newton diagram below the line (0,n)-(m,0)
136         0 else
137EXAMPLE: example testreducible;  shows an example
138"
139{
140 return(size(jet(f,m*n-1,intvec(n,m))) != 0)
141}
142example
143{ "EXAMPLE:"; echo = 2;
144  ring rg=0,(x,y),ls;
145  testreducible(x2+y3-xy4,3,2);
146}
147///////////////////////////////////////////////////////////////////////////////
148proc T_Transform (poly f, int Q, int N)
149"USAGE:   T_Transform(f,Q,N);  f poly, Q,N int
150RETURN:  f(y,xy^Q)/y^NQ   if x,y are the ring variables
151NOTE:    this is intended for irreducible power series f
152EXAMPLE: example T_Transform;  shows an example
153"
154{
155 map T = basering,var(2),var(1)*var(2)^Q;
156 return(T(f)/var(2)^(N*Q));
157}
158example
159{ "EXAMPLE:"; echo = 2;
160  ring exrg=0,(x,y),ls;
161  export exrg;
162  T_Transform(x3+y2-xy3,1,2);
163  kill exrg;
164}
165///////////////////////////////////////////////////////////////////////////////
166proc T1_Transform (poly f, number d, int Q)
167"USAGE:   T1_Transform(f,d,Q);  f poly, d number, Q int
168RETURN:  f(x,y+d*x^Q)   if x,y are the ring variables
169EXAMPLE: example T1_Transform;  shows an example
170"
171{
172 map T1 = basering,var(1),var(2)+d*var(1)^Q;
173 return(T1(f));
174}
175example
176{ "EXAMPLE:"; echo = 2;
177  ring exrg=0,(x,y),ls;
178  export exrg;
179  T1_Transform(y2-2xy+x2+x2y,1,1);
180  kill exrg;
181}
182///////////////////////////////////////////////////////////////////////////////
183
184proc T2_Transform (poly f_neu, number d, int M, int N, poly refpoly)
185"USAGE:   T2_Transform(f,d,M,N,ref); f poly, d number; M,N int; ref poly
186RETURN:  list: poly T2(f,d',M,N), number d' in \{ d, 1/d \}
187ASSUME:  ref has the same Newton polygon as f (but can be simpler)
188         for this you can e.g. use the proc `referencepoly' or simply f again
189COMMENT: T2 is a composition of T_Transform and T1_Transform; the exact
190         definition can be found in  Rybowicz: `Sur le calcul des places ...'
191         or in  Lamm: `Hamburger-Noether-Entwicklung von Kurvensingularitaeten'
192SEE ALSO: T_Transform, T1_Transform, referencepoly
193EXAMPLE: example T2_Transform;  shows an example
194"
195{
196 //---------------------- compute gcd and extgcd of N,M -----------------------
197  int ggt=gcd(M,N);
198  M=M/ggt; N=N/ggt;
199  list ts=extgcd(M,N);
200  int tau,sigma=ts[2],-ts[3];
201  int s,t;
202  poly xp=var(1);
203  poly yp=var(2);
204  poly hilf;
205  if (sigma<0) { tau=-tau; sigma=-sigma;}
206 // es gilt: 0<=tau<=N, 0<=sigma<=M, |N*sigma-M*tau| = 1 = ggT(M,N)
207  if (N*sigma < M*tau) { d = 1/d; }
208 //--------------------------- euklid. Algorithmus ----------------------------
209  int R;
210  int M1,N1=M,N;
211  for ( R=M1%N1; R!=0; ) { M1=N1; N1=R; R=M1%N1;}
212  int Q=M1 / N1;
213  map T1 = basering,xp,yp+d*xp^Q;
214  map Tstar=basering,xp^(N-Q*tau)*yp^tau,xp^(M-sigma*Q)*yp^sigma;
215  if (defined(HNDebugOn)) {
216   "Trafo. T2: x->x^"+string(N-Q*tau)+"*y^"+string(tau)+", y->x^"
217    +string(M-sigma*Q)+"*y^"+string(sigma);
218   "delt =",d,"Q =",Q,"tau,sigma =",tau,sigma;
219  }
220 //------------------- Durchfuehrung der Transformation T2 --------------------
221  f_neu=Tstar(f_neu);
222  refpoly=Tstar(refpoly);
223  //--- dividiere f_neu so lange durch x & y, wie die Division aufgeht,
224  //    benutze ein Referenzpolynom mit gleichem Newtonpolynom wie f_neu zur
225  //    Beschleunigung: ---
226  for (hilf=refpoly/xp; hilf*xp==refpoly; hilf=refpoly/xp) {refpoly=hilf; s++;}
227  for (hilf=refpoly/yp; hilf*yp==refpoly; hilf=refpoly/yp) {refpoly=hilf; t++;}
228  f_neu=f_neu/(xp^s*yp^t);
229  return(list(T1(f_neu),d));
230}
231example
232{ "EXAMPLE:"; echo = 2;
233  ring exrg=0,(x,y),ds;
234  export exrg;
235  poly f=y2-2x2y+x6-x5y+x4y2;
236  T2_Transform(f,1/2,4,1,f);
237  T2_Transform(f,1/2,4,1,referencepoly(newtonpoly(f,1)));
238  // if  size(referencepoly) << size(f)  the 2nd example would be faster
239  referencepoly(newtonpoly(f,1));
240  kill exrg;
241}
242///////////////////////////////////////////////////////////////////////////////
243
244proc koeff (poly f, int I, int J)
245"USAGE:   koeff(f,I,J); f bivariate polynomial, I,J integers
246RETURN:  if f = sum(a(i,j)*x^i*y^j), then koeff(f,I,J)= a(I,J) (of type number)
247NOTE:    J must be in the range of the exponents of the 2nd ring variable
248EXAMPLE: example koeff;  shows an example
249"
250{
251  matrix mat = coeffs(coeffs(f,var(2))[J+1,1],var(1));
252  if (size(mat) <= I) { return(0);}
253  else { return(leadcoef(mat[I+1,1]));}
254}
255example
256{ "EXAMPLE:"; echo = 2;
257  ring r=0,(x,y),dp;
258  koeff(x2+2xy+3xy2-x2y-2y3,1,2);
259}
260///////////////////////////////////////////////////////////////////////////////
261
262proc squarefree (poly f)
263"USAGE:  squarefree(f);  f poly
264ASSUME:  f is a bivariate polynomial (in the first 2 ring variables).
265RETURN:  poly, a squarefree divisor of f.
266NOTE:    Usually, the return value is the greatest squarefree divisor, but
267         there is one exception: factors with a p-th root, p the
268         characteristic of the basering, are lost.
269SEE ALSO: allsquarefree
270EXAMPLE: example squarefree; shows some examples.
271"
272{
273 //----------------- Wechsel in geeigneten Ring & Variablendefinition ---------
274  if (nvars(basering)!=2)
275  { ERROR("basering must have exactly 2 variables for Hnoether::squarefree"); }
276  def altring = basering;
277  int e;
278  int gcd_ok=1;
279  string mipl="0";
280  if (size(parstr(altring))==1) { mipl=string(minpoly); }
281 //---- test: char = (p^k,a) (-> gcd not implemented) or (p,a) (gcd works) ----
282  if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) {
283    string tststr=charstr(basering);
284    tststr=tststr[1..find(tststr,",")-1];           //-> "p^k" bzw. "p"
285    gcd_ok=(tststr==string(char(basering)));
286  }
287  execute("ring rsqrf = ("+charstr(altring)+"),(x,y),dp;");
288  if ((gcd_ok!=0) && (mipl!="0")) { execute("minpoly="+mipl+";"); }
289  poly f=fetch(altring,f);
290  poly dif,g,l;
291  if ((char(basering)==0) and (charstr(basering)!=string(char(basering)))
292      and (mipl!="0")) {
293    gcd_ok=0;                   // since Singular 1.2 gcd no longer implemented
294  }
295  if (gcd_ok!=0) {
296 //--------------------- Berechne f/ggT(f,df/dx,df/dy) ------------------------
297    dif=diff(f,x);
298    if (dif==0) { g=f; }        // zur Beschleunigung
299    else { g=gcd(f,dif); }
300    if (g!=1) {                 // sonst schon sicher, dass f quadratfrei
301     dif=diff(f,y);
302     if (dif!=0) { g=gcd(g,dif); }
303    }
304    if (g!=1) {
305     e=0;
306     if (g==f) { l=1; }         // zur Beschleunigung
307     else {
308       module m=syz(ideal(g,f));
309       if (deg(m[2,1])>0) {
310         "!! The Singular command 'syz' has returned a wrong result !!";
311         l=1;                   // Division f/g muss aufgehen
312       }
313       else { l=m[1,1]; }
314     }
315    }
316    else { e=1; }
317  }
318  else {
319 //------------------- Berechne syz(f,df/dx) oder syz(f,df/dy) ----------------
320 //-- Achtung: Ist f reduzibel, koennen Faktoren mit Ableitung Null verloren --
321 //-- gehen! Ist aber nicht weiter schlimm, weil char (p^k,a) nur im irred.  --
322 //-- Fall vorkommen kann. Wenn f nicht g^p ist, wird auf jeden Fall         --
323 //------------------------ ein Faktor gefunden. ------------------------------
324    dif=diff(f,x);
325    if (dif == 0) {
326     dif=diff(f,y);
327     if (dif==0) { e=2; l=1; } // f is of power divisible by char of basefield
328     else { l=syz(ideal(dif,f))[1,1];  // x^p+y^(p-1) abgedeckt
329            if (subst(f,x,0)==0) { l=l*x; }
330            if (deg(l)==deg(f))  { e=1;}
331            else {e=0;}
332     }
333    }
334    else { l=syz(ideal(dif,f))[1,1];
335           if (subst(f,y,0)==0) { l=l*y; }
336           if (deg(l)==deg(f))  { e=1;}
337           else {e=0;}
338    }
339  }
340 //--------------- Wechsel in alten Ring und Rueckgabe des Ergebnisses --------
341  setring altring;
342  if (e==1) { return(f); }    // zur Beschleunigung
343  else {
344   poly l=fetch(rsqrf,l);
345   return(l);
346  }
347}
348example
349{ "EXAMPLE:"; echo = 2;
350 ring exring=3,(x,y),dp;
351 squarefree((x3+y)^2);
352 squarefree((x+y)^3*(x-y)^2); // Warning: (x+y)^3 is lost
353 squarefree((x+y)^4*(x-y)^2); // result is (x+y)*(x-y)
354}
355///////////////////////////////////////////////////////////////////////////////
356
357proc allsquarefree (poly f, poly l)
358"USAGE : allsquarefree(f,g);  f,g poly
359ASSUME: g is the output of @code{squarefree(f)}.
360RETURN: the greatest squarefree divisor of f.
361NOTE  : This proc uses factorize to get the missing factors of f not in g and,
362        therefore, may be slow.
363SEE ALSO: squarefree
364EXAMPLE: example allsquarefree;  shows an example
365"
366{
367 //------------------------ Wechsel in geeigneten Ring ------------------------
368 def altring = basering;
369 string mipl="0";
370 if (size(parstr(altring))==1) { mipl=string(minpoly); }
371 if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) {
372   string tststr=charstr(basering);
373   tststr=tststr[1..find(tststr,",")-1];           //-> "p^k" bzw. "p"
374   if (tststr!=string(char(basering))) {
375     " Sorry -- not implemented for this ring (gcd doesn't work)";
376     return(l);
377   }
378 }
379 execute("ring rsqrf = ("+charstr(altring)+"),(x,y),dp;");
380 if (mipl!="0") { execute("minpoly="+mipl+";"); }
381 poly f=fetch(altring,f);
382 poly l=fetch(altring,l);
383 //---------- eliminiere bereits mit squarefree gefundene Faktoren ------------
384 poly g=l;
385 while (deg(g)!=0) {
386   f=syz(ideal(g,f))[1,1];                         // f=f/g;
387   g=gcd(f,l);
388 }                                                 // jetzt f=h^p
389 //--------------- Berechne uebrige Faktoren mit factorize --------------------
390 if (deg(f)>0) {
391  g=1;
392//*CL old:  ideal factf=factorize(f,1);
393//*         for (int i=1; i<=size(factf); i++) { g=g*factf[i]; }
394  ideal factf=factorize(f)[1];
395  for (int i=2; i<=size(factf); i++) { g=g*factf[i]; }
396  poly testp=squarefree(g);
397  if (deg(testp)<deg(g)) {
398    "!! factorize has not worked correctly !!";
399    if (testp==1) {" We cannot proceed ..."; g=1;}
400    else {" But we could recover some factors..."; g=testp;}
401  }
402  l=l*g;
403 }
404 //--------------- Wechsel in alten Ring und Rueckgabe des Ergebnisses --------
405 setring altring;
406 l=fetch(rsqrf,l);
407 return(l);
408}
409example
410{ "EXAMPLE:"; echo = 2;
411  ring exring=7,(x,y),dp;
412  poly f=(x+y)^7*(x-y)^8;
413  poly g=squarefree(f);
414  g;                      // factor x+y lost, since characteristic=7
415  allsquarefree(f,g);     // all factors (x+y)*(x-y) found
416}
417///////////////////////////////////////////////////////////////////////////////
418
419proc is_irred (poly f)
420"USAGE:   is_irred(f); f poly
421ASSUME:  f is a squarefree bivariate polynomial (in the first 2 ring
422         variables).
423RETURN:  int (0 or 1): @*
424         - @code{is_irred(f)=1} if f is irreducible as a formal power
425         series over the algebraic closure of its coefficient field (f
426         defines an analytically irreducible curve at zero), @*
427         - @code{is_irred(f)=0} otherwise.
428NOTE:    0 and units in the ring of formal power series are considered to be
429         not irreducible.
430KEYWORDS: irreducible power series
431EXAMPLE: example is_irred;  shows an example
432"
433{
434  int pl=printlevel;
435  printlevel=-1;
436  list hnl=develop(f,-1);
437  printlevel=pl;
438  return(hnl[5]);
439}
440example
441{ "EXAMPLE:"; echo = 2;
442  ring exring=0,(x,y),ls;
443  is_irred(x2+y3);
444  is_irred(x2+y2);
445  is_irred(x2+y3+1);
446}
447///////////////////////////////////////////////////////////////////////////////
448
449static proc polytest(poly f)
450"USAGE : polytest(f); f poly in x and y
451RETURN: a monomial of f with |coefficient| > 16001
452          or exponent divisible by 32003, if there is one
453        0 else (in this case computing a squarefree divisor
454                in characteristic 32003 could make sense)
455NOTE:   this procedure is only useful in characteristic zero, because otherwise
456        there is no appropriate ordering of the leading coefficients
457"
458{
459 poly verbrecher=0;
460 intvec leitexp;
461 for (; (f<>0) and (verbrecher==0); f=f-lead(f)) {
462  if ((leadcoef(f)<-16001) or (leadcoef(f)>16001)) {verbrecher=lead(f);}
463  leitexp=leadexp(f);
464  if (( ((leitexp[1] % 32003) == 0)   and (leitexp[1]<>0))
465     or ( ((leitexp[2] % 32003) == 0) and (leitexp[2]<>0)) )
466       {verbrecher=lead(f);}
467 }
468 return(verbrecher);
469}
470
471//////////////////////////////////////////////////////////////////////////////
472
473
474proc develop
475"USAGE:   develop(f [,n]); f poly, n int
476ASSUME:  f is a bivariate polynomial (in the first 2 ring variables) and
477         irreducible as power series (for reducible f use @code{hnexpansion}).
478RETURN:  list @code{L} with:
479@texinfo
480@table @asis
481@item @code{L[1]}; matrix:
482         Each row contains the coefficients of the corresponding line of the
483         Hamburger-Noether expansion (HNE). The end of the line is marked in
484         the matrix by the first ring variable (usually x).
485@item @code{L[2]}; intvec:
486         indicating the length of lines of the HNE
487@item @code{L[3]}; int:
488         0  if the 1st ring variable was transversal (with respect to f), @*
489         1  if the variables were changed at the beginning of the
490            computation, @*
491        -1  if an error has occurred.
492@item @code{L[4]}; poly:
493         the transformed polynomial of f to make it possible to extend the
494         Hamburger-Noether development a posteriori without having to do
495         all the previous calculation once again (0 if not needed)
496@item @code{L[5]}; int:
497         1  if the curve has exactly one branch (i.e., is irreducible), @*
498         0  else (i.e., the curve has more than one HNE, or f is not valid).
499@end table
500@end texinfo
501DISPLAY: The (non zero) elements of the HNE (if not called by another proc).
502NOTE:    The optional parameter @code{n} affects only the computation of
503         the LAST line of the HNE. If it is given, the HN-matrix @code{L[1]}
504         will have at least @code{n} columns. @*
505         Otherwise, the number of columns will be chosen minimal such that the
506         matrix contains all necessary information (i.e., all lines of the HNE
507         but the last (which is in general infinite) have place). @*
508         If @code{n} is negative, the algorithm is stopped as soon as the
509         computed information is sufficient for @code{invariants(L)}, but the
510         HN-matrix @code{L[1]} may still contain undetermined elements, which
511         are marked with the 2nd variable (of the basering). @*
512         For time critical computations it is recommended to use
513         @code{ring ...,(x,y),ls} as basering - it increases the algorithm's
514         speed. @*
515         If @code{printlevel>=0} comments are displayed (default is
516         @code{printlevel=0}).
517SEE ALSO: hnexpansion, extdevelop, displayHNE
518EXAMPLES: example develop;         shows an example
519          example parametrize;     shows an example for using the 2nd parameter
520"
521{
522 //--------- Abfangen unzulaessiger Ringe: 1) nur eine Unbestimmte ------------
523 poly f=#[1];
524 if (size(#) > 1) {int maxspalte=#[2];}
525 else             {int maxspalte= 1 ; }
526 if (nvars(basering) < 2) {
527   " Sorry. I need two variables in the ring.";
528   return(list(matrix(maxideal(1)[1]),intvec(0),-1,poly(0),0));}
529 if (nvars(basering) > 2) {
530   dbprint(printlevel-voice+2,
531   " Warning! You have defined too many variables!
532 All variables except the first two will be ignored!"
533           );
534 }
535
536 string namex=varstr(1); string namey=varstr(2);
537 list return_error=matrix(maxideal(1)[2]),intvec(0),int(-1),poly(0),int(0);
538
539 //------------- 2) mehrere Unbestimmte, weitere unzulaessige Ringe -----------
540 // Wir koennen einheitlichen Rueckgabewert benutzen, aus dem ersichtlich ist,
541 // dass ein Fehler aufgetreten ist: return_error.
542 //----------------------------------------------------------------------------
543
544 if (charstr(basering)=="real") {
545  " The algorithm doesn't work with 'real' as coefficient field.";
546                     // denn : map from characteristic -1 to -1 not implemented
547  return(return_error);
548 }
549 if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) {
550 //-- teste, ob char = (p^k,a) (-> a primitiv; erlaubt) oder (p,a[,b,...]) ----
551    string tststr=charstr(basering);
552    tststr=tststr[1..find(tststr,",")-1];           //-> "p^k" bzw. "p"
553    int primit=(tststr==string(char(basering)));
554    if (primit!=0) {
555      " Such extensions of Z/p are not implemented.";
556      " Please try (p^k,a) as ground field or use `hnexpansion'.";
557      return(return_error);
558    }
559 }
560 //---- Ende der unzulaessigen Ringe; Ringwechsel in einen guenstigen Ring: ---
561
562 int ringwechsel=(varstr(basering)!="x,y") or (ordstr(basering)!="ls(2),C");
563
564 def altring = basering;
565 if (ringwechsel) {
566   string mipl=string(minpoly);
567   execute("ring guenstig = ("+charstr(altring)+"),(x,y),ls;");
568   if ((char(basering)==0) && (mipl!="0")) {
569     execute("minpoly="+mipl+";");
570   }}
571 else { def guenstig=basering; }
572 export guenstig;
573
574 //-------------------------- Initialisierungen -------------------------------
575 map m=altring,x,y;
576 if (ringwechsel) { poly f=m(f); }
577 if (defined(HNDebugOn))
578 {"received polynomial: ",f,", where x =",namex,", y =",namey;}
579 kill m;
580 int M,N,Q,R,l,e,hilf,eps,getauscht,Abbruch,zeile,exponent,Ausgabe;
581
582 // Werte von Ausgabe: 0 : normale HNE-Matrix,
583 // 1 : Fehler aufgetreten - Matrix (namey) zurueck
584 // 2 : Die HNE ist eine Nullzeile - Matrix (0) zurueck
585 // int maxspalte=1; geaendert: wird jetzt am Anfang gesetzt
586
587 int minimalHNE=0;          // Flag fuer minimale HNE-Berechnung
588 int einzweig=1;            // Flag fuer Irreduzibilit"at
589 intvec hqs;                // erhaelt die Werte von h(zeile)=Q;
590
591 if (maxspalte<0) {
592   minimalHNE=1;
593   maxspalte=1;
594 }
595
596 number c,delt;
597 int p = char(basering);
598 string ringchar=charstr(basering);
599 map xytausch = basering,y,x;
600 if ((p!=0) and (ringchar != string(p))) {
601                            // coefficient field is extension of Z/pZ
602   execute("int n_elements="+
603           ringchar[1,size(ringchar)-size(parstr(basering))-1]+";");
604                            // number of elements of actual ring
605   number generat=par(1);   // generator of the coefficient field of the ring
606 }
607
608
609 //========= Abfangen von unzulaessigen oder trivialen Eingaben ===============
610 //------------ Nullpolynom oder Einheit im Potenzreihenring: -----------------
611 if (f == 0) {
612   dbprint(printlevel+1,"The given polynomial is the zero-polynomial !");
613   Abbruch=1; Ausgabe=1;
614 }
615 else {
616   intvec nm = getnm(f);
617   N = nm[1]; M = nm[2]; // Berechne Schnittpunkte Newtonpolygon mit Achsen
618   if (N == 0) {
619     dbprint(printlevel+1,"The given polynomial is a unit as power series !");
620     Abbruch=1; Ausgabe=1;
621   }
622   else {
623    if (N == -1) {
624      if ((voice==2) && (printlevel > -1)) { "The HNE is x = 0"; }
625      Abbruch=1; Ausgabe=2; getauscht=1;
626      if (M <> 1) { einzweig=0; }
627    }
628    else {
629     if (M == -1) {
630       if ((voice==2) && (printlevel > -1)) { "The HNE is y = 0"; }
631       Abbruch=1; Ausgabe=2;
632       if (N <> 1) { einzweig=0; }
633   }}}
634 }
635 //--------------------- Test auf Quadratfreiheit -----------------------------
636 if (Abbruch==0) {
637
638 //-------- Fall basering==0,... : Wechsel in Ring mit char >0 ----------------
639 // weil squarefree eine Standardbasis berechnen muss (verwendet Syzygien)
640 // -- wenn f in diesem Ring quadratfrei ist, dann erst recht im Ring guenstig
641 //----------------------------------------------------------------------------
642
643  if ((p==0) and (size(charstr(basering))==1)) {
644   int testerg=(polytest(f)==0);
645   ring zweitring = 32003,(x,y),dp;
646   map polyhinueber=guenstig,x,y;     // fetch geht nicht
647   poly f=polyhinueber(f);
648   poly test_sqr=squarefree(f);
649   if (test_sqr != f) {
650    if (printlevel>0) {
651      "Most probably the given polynomial is not squarefree. But the test was";
652      "made in characteristic 32003 and not 0 to improve speed. You can";
653      "(r) redo the test in char 0 (but this may take some time)";
654      "(c) continue the development, if you're sure that the polynomial",
655      "IS squarefree";
656      if (testerg==1) {
657        "(s) continue the development with a squarefree factor (*)";}
658      "(q) or just quit the algorithm (default action)";
659      "";"Please enter the letter of your choice:";
660      string str=read("")[1];
661    }
662    else { string str="r"; }      // printlevel <= 0: non-interactive behaviour
663    setring guenstig;
664    map polyhinueber=zweitring,x,y;
665    if (str=="r") {
666      poly test_sqr=squarefree(f);
667      if (test_sqr != f) {
668       if (printlevel>0) { "The given polynomial is in fact not squarefree."; }
669       else              { "The given polynomial is not squarefree!"; }
670       "I'll continue with the radical.";
671       if (printlevel>0) { pause("Hit RETURN to continue:"); }
672       f=test_sqr;
673      }
674      else {
675       dbprint(printlevel,
676        "everything is ok -- the polynomial is squarefree in char(k)=0");
677      }
678    }
679    else {
680      if ((str=="s") and (testerg==1)) {
681       "(*) attention: it could be that the factor is only one in char 32003!";
682        f=polyhinueber(test_sqr);
683      }
684      else {
685        if (str<>"c") {
686          setring altring;kill guenstig;kill zweitring;
687          return(return_error);}
688        else { "if the algorithm doesn't terminate, you were wrong...";}
689    }}
690    kill zweitring;
691    nm = getnm(f);             // N,M haben sich evtl. veraendert
692    N = nm[1]; M = nm[2];      // Berechne Schnittpunkte Newtonpolynom mit Achsen
693    if (defined(HNDebugOn)) {"I continue with the polynomial",f; }
694   }
695   else {
696     setring guenstig;
697     kill zweitring;
698   }
699  }
700 // ------------------- Fall Charakteristik > 0 -------------------------------
701  else {
702   poly test_sqr=squarefree(f);
703   if (test_sqr == 1) {
704    "The given polynomial is of the form g^"+string(p)+", therefore",
705    "reducible.";"Please try again.";
706    setring altring;
707    kill guenstig;
708    return(return_error);}
709   if (test_sqr != f) {
710    "The given polynomial is not squarefree. I'll continue with the radical.";
711    if (p != 0)
712     {"But if the polynomial contains a factor of the form g^"+string(p)+",";
713      "this factor will be lost.";}
714    if (printlevel>0) { pause("Hit RETURN to continue:"); }
715    f=test_sqr;
716    nm = getnm(f);              // N,M haben sich veraendert
717    N = nm[1]; M = nm[2];       // Berechne Schnittpunkte Newtonpolynom mit Achsen
718    if (defined(HNDebugOn)) {"I continue with the polynomial",f; }
719   }
720
721  }                             // endelse(p==0)
722
723  if (N==0) {
724    " Sorry. The remaining polynomial is a unit in the power series ring...";
725    setring altring;kill guenstig;return(return_error);
726  }
727 //---------------------- gewaehrleiste, dass x transvers ist -----------------
728  if (M < N)
729  { f = xytausch(f);            // Variablentausch : x jetzt transvers
730    getauscht = 1;              // den Tausch merken
731    M = M+N; N = M-N; M = M-N;  // M, N auch vertauschen
732  }
733  if (defined(HNDebugOn)) {
734   if (getauscht) {"x<->y were exchanged; polynomial is now ",f;}
735   else           {"x , y were not exchanged";}
736   "M resp. N are now",M,N;
737  }
738 }                              // end(if Abbruch==0)
739
740 ideal a(0);
741 while (Abbruch==0) {
742
743 //================= Beginn der Schleife (eigentliche Entwicklung) ============
744
745 //------------------- ist das Newtonpolygon eine gerade Linie? ---------------
746  if (testreducible(f,N,M)) {
747    dbprint(printlevel+1," The given polynomial is not irreducible");
748    kill guenstig;
749    setring altring;
750    return(return_error);       // Abbruch der Prozedur!
751  }
752  R = M%N;
753  Q = M / N;
754
755 //-------------------- Fall Rest der Division R = 0 : ------------------------
756  if (R == 0) {
757    c = koeff(f,0,N);
758    if (c == 0) {"Something has gone wrong! I didn't get N correctly!"; exit;}
759    e = gcd(M,N);
760 //----------------- Test, ob leitf = c*(y^N - delta*x^(m/e))^e ist -----------
761    if (p==0) {
762      delt = koeff(f,M/ e,N - N/ e) / (-1*e*c);
763      if (defined(HNDebugOn)) {"quasihomogeneous leading form:",
764         leit(f,N,M)," = ",c,"* (y -",delt,"* x^"+string(M/ e)+")^",e," ?";}
765      if (leit(f,N,M) != c*(y^(N/ e) - delt*x^(M/ e))^e) {
766        dbprint(printlevel+1," The given polynomial is reducible !");
767        Abbruch=1; Ausgabe=1; }
768    }
769    else {                     // p!=0
770      if (e%p != 0) {
771        delt = koeff(f,M/ e,N - N/ e) / (-1*e*c);
772        if (defined(HNDebugOn)) {"quasihomogeneous leading form:",
773           leit(f,N,M)," = ",c,"* (y -",delt,"* x^"+string(M/ e)+")^",e," ?";}
774        if (leit(f,N,M) != c*(y^(N/ e) - delt*x^(M/ e))^e) {
775           dbprint(printlevel+1," The given polynomial is reducible !");
776           Abbruch=1; Ausgabe=1; }
777      }
778
779      else {                   // e%p == 0
780        eps = e;
781        for (l = 0; eps%p == 0; l=l+1) { eps=eps/ p;}
782        if (defined(HNDebugOn)) {e," -> ",eps,"*",p,"^",l;}
783        delt = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c);
784
785        if ((ringchar != string(p)) and (delt != 0)) {
786 //- coeff. field is not Z/pZ => we`ve to correct delta by taking (p^l)th root-
787          if (delt == generat) {exponent=1;}
788          else {
789           if (delt == 1) {exponent=0;}
790           else {
791            exponent=pardeg(delt);
792
793 //-- an dieser Stelle kann ein Fehler auftreten, wenn wir eine transzendente -
794 //-- Erweiterung von Z/pZ haben: dann ist das hinzuadjungierte Element kein  -
795 //-- Erzeuger der mult. Gruppe, d.h. in Z/pZ (a) gibt es i.allg. keinen      -
796 //-- Exponenten mit z.B. a2+a = a^exp                                        -
797 //----------------------------------------------------------------------------
798          }}
799          delt = generat^(extgcd(n_elements-1,p^l)[3]*exponent);
800        }
801
802        if (defined(HNDebugOn)) {"quasihomogeneous leading form:",
803          leit(f,N,M)," = ",c,"* (y^"+string(N/ e),"-",delt,"* x^"
804          +string(M/ e)+")^",e,"  ?";}
805        if (leit(f,N,M) != c*(y^(N/ e) - delt*x^(M/ e))^e) {
806          dbprint(printlevel+1," The given polynomial is reducible !");
807          Abbruch=1; Ausgabe=1; }
808      }
809    }
810    if (Abbruch == 0) {
811      f = T1_Transform(f,delt,M/ e);
812      dbprint(printlevel-voice+2,"a("+string(zeile)+","+string(Q)+") = "
813              +string(delt));
814      a(zeile)[Q]=delt;
815      if (defined(HNDebugOn)) {"transformed polynomial: ",f;}}
816
817      nm=getnm(f); N=nm[1]; M=nm[2];        // Neuberechnung des Newtonpolygons
818  }
819 //--------------------------- Fall R > 0 : -----------------------------------
820  else {
821    dbprint(printlevel-voice+2, "h("+string(zeile)+ ") ="+string(Q));
822    hqs[zeile+1]=Q;                  // denn zeile beginnt mit dem Wert 0
823    a(zeile)[Q+1]=x;                 // Markierung des Zeilenendes der HNE
824    maxspalte=maxspalte*((Q+1) < maxspalte) + (Q+1)*((Q+1) >= maxspalte);
825                                     // Anpassung der Sp.zahl der HNE-Matrix
826    f = T_Transform(f,Q,N);
827    if (defined(HNDebugOn)) {"transformed polynomial: ",f;}
828    zeile=zeile+1;
829 //------------ Bereitstellung von Speicherplatz fuer eine neue Zeile: --------
830    ideal a(zeile);
831    M=N;N=R;
832  }
833
834 //--------------- schneidet das Newtonpolygon beide Achsen? ------------------
835  if (M==-1) {
836     dbprint(printlevel-voice+2,"The HNE is finite!");
837     a(zeile)[Q+1]=x;   // Markiere das Ende der Zeile
838     hqs[zeile+1]=Q;
839     maxspalte=maxspalte*((Q+1) < maxspalte) + (Q+1)*((Q+1) >= maxspalte);
840     if (N <> 1) { einzweig=0; }
841     f=0;               // transformiertes Polynom wird nicht mehr gebraucht
842     Abbruch=1;
843  }
844  else {if (M<N) {"Something has gone wrong: M<N";}}
845  if(defined(HNDebugOn)) {"new M,N:",M,N;}
846
847 //----------------- Abbruchbedingungen fuer die Schleife: --------------------
848  if ((N==1) and (Abbruch!=1) and ((M > maxspalte) or (minimalHNE==1))
849      and (size(a(zeile))>0))
850 //----------------------------------------------------------------------------
851 // Abbruch, wenn die Matrix so voll ist, dass eine neue Spalte angefangen
852 // werden muesste und die letzte Zeile nicht nur Nullen enthaelt
853 // oder wenn die Matrix nicht voll gemacht werden soll (minimale Information)
854 //----------------------------------------------------------------------------
855   { Abbruch=1; hqs[zeile+1]=-1;
856     if (maxspalte < ncols(a(zeile))) { maxspalte=ncols(a(zeile));}
857     if ((minimalHNE==1) and (M <= maxspalte)) {
858 // teile param mit, dass Eintraege der letzten Zeile nur teilw. richtig sind:-
859       hqs[zeile+1]=-M;
860 //------------- markiere den Rest der Zeile als unbekannt: -------------------
861       for (R=M; R <= maxspalte; R++) { a(zeile)[R]=y;}
862     }                  // R wird nicht mehr gebraucht
863   }
864 //========================= Ende der Schleife ================================
865
866 }
867 setring altring;
868 if (Ausgabe == 0) {
869 //-------------------- Ergebnis in den alten Ring transferieren: -------------
870   map zurueck=guenstig,maxideal(1)[1],maxideal(1)[2];
871   matrix amat[zeile+1][maxspalte];
872   ideal uebergabe;
873   for (e=0; e<=zeile; e=e+1) {
874     uebergabe=zurueck(a(e));
875     if (ncols(uebergabe) > 1) {
876      amat[e+1,1..ncols(uebergabe)]=uebergabe;}
877     else {amat[e+1,1]=uebergabe[1];}
878   }
879   if (ringwechsel) {
880     if (nvars(altring)==2) { f=fetch(guenstig,f); }
881     else                   { f=zurueck(f); }
882   }
883 }
884
885 kill guenstig;
886 if ((einzweig==0) && (voice==2) && (printlevel > -1)) {
887    "// Note: The curve is reducible, but we were able to compute a HNE.";
888    "// This means the result is only one of several existing HNE's.";
889 }
890 if (Ausgabe == 0) { return(list(amat,hqs,getauscht,f,einzweig));}
891 if (Ausgabe == 1) { return(return_error);}             // error has occurred
892 if (Ausgabe == 2) { return(list(matrix(ideal(0,x)),intvec(1),getauscht,
893                                 poly(0),einzweig));}   // HNE is x=0 or y=0
894}
895example
896{ "EXAMPLE:"; echo = 2;
897  ring exring = 7,(x,y),ds;
898  list Hne=develop(4x98+2x49y7+x11y14+2y14);
899  print(Hne[1]);
900  // therefore the HNE is:
901  // z(-1)= 3*z(0)^7 + z(0)^7*z(1),
902  // z(0) = z(1)*z(2),       (there is 1 zero in the 2nd row before x)
903  // z(1) = z(2)^3*z(3),     (there are 3 zeroes in the 3rd row)
904  // z(2) = z(3)*z(4),
905  // z(3) = -z(4)^2 + 0*z(4)^3 +...+ 0*z(4)^8 + ?*z(4)^9 + ...
906  // (the missing x in the last line indicates that it is not complete.)
907  Hne[2];
908  param(Hne);
909  // parametrization:   x(t)= -t^14+O(t^21),  y(t)= -3t^98+O(t^105)
910  // (the term -t^109 in y may have a wrong coefficient)
911  displayHNE(Hne);
912}
913
914///////////////////////////////////////////////////////////////////////////////
915//               procedures to extract information out of HNE                //
916///////////////////////////////////////////////////////////////////////////////
917
918proc param (list L, list #)
919"USAGE:  param(L [,s]); L list, s any type (optional)
920ASSUME:  L is the output of @code{develop(f)}, or of
921        @code{extdevelop(develop(f),n)}, or (one entry in) the list of HN
922        data created by @code{hnexpansion(f[,\"ess\"])}.
923RETURN: If L are the HN data of an irreducible plane curve singularity f: a
924        parametrization for f in the following format: @*
925        - if only the list L is given, the result is an ideal of two
926        polynomials p[1],p[2]: if the HNE was finite then f(p[1],p[2])=0};
927        if not, the true parametrization will be given by two power series,
928        and p[1],p[2] are truncations of these series.@*
929        - if the optional parameter s is given, the result is a list l:
930        l[1]=param(L) (ideal) and l[2]=intvec with two entries indicating
931        the highest degree up to which the coefficients of the monomials in
932        l[1] are exact (entry -1 means that the corresponding parametrization
933        is exact).
934        If L collects the HN data of a reducible plane curve singularity f,
935        the return value is a list of parametrizations in the respective
936        format.
937NOTE:   If the basering has only 2 variables, the first variable is chosen
938        as indefinite. Otherwise, the 3rd variable is chosen.
939SEE ALSO: develop, extdevelop
940KEYWORDS: parametrization
941EXAMPLE: example param;     shows an example
942         example develop;   shows another example
943"
944{
945 //-------------------------- Initialisierungen -------------------------------
946 int return_list;
947 if (size(#)>0) { return_list=1; }
948
949 if (typeof(L[1])=="list") { // output of hnexpansion (> 1 branch)
950   list Ergebnis;
951   for (int i=1; i<=size(L); i++) {
952     dbprint(printlevel-voice+4,"// Parametrization of branch number "
953       +string(i)+" computed.");
954     printlevel=printlevel+1;
955     if (return_list==1) { Ergebnis[i]=param(L[i],1); }
956     else                { Ergebnis[i]=param(L[i]); }
957     printlevel=printlevel-1;
958   }
959   return(Ergebnis);
960 }
961 else {
962   matrix m=L[1];
963   intvec v=L[2];
964   int switch=L[3];
965 }
966 if (switch==-1) {
967   "An error has occurred in develop, so there is no HNE.";
968   return(ideal(0,0));
969 }
970 int fehler,fehlervor,untergrad,untervor,beginn,i,zeile,hilf;
971
972 if (nvars(basering) > 2) { poly z(size(v)+1)=var(3); }
973 else                     { poly z(size(v)+1)=var(1); }
974 poly z(size(v));
975 zeile=size(v);
976 //------------- Parametrisierung der untersten Zeile der HNE -----------------
977 if (v[zeile] > 0) {
978   fehler=0;           // die Parametrisierung wird exakt werden
979   for (i=1; i<=v[zeile]; i++) {
980     z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
981   }
982 }
983 else {
984   untervor=1;         // = Untergrad der vorhergehenden Zeile
985   if (v[zeile]==-1) {
986     fehler=ncols(m)+1;
987     for (i=1; i<=ncols(m); i++) {
988       z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
989       if ((untergrad==0) and (m[zeile,i]!=0)) {untergrad=i;}
990                       // = Untergrad der aktuellen Zeile
991     }
992   }
993   else {
994     fehler= -v[zeile];
995     for (i=1; i<-v[zeile]; i++) {
996       z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
997       if ((untergrad==0) and (m[zeile,i]!=0)) {untergrad=i;}
998     }
999   }
1000 }
1001 //------------- Parametrisierung der restlichen Zeilen der HNE ---------------
1002 for (zeile=size(v)-1; zeile>0; zeile--) {
1003   poly z(zeile);
1004   beginn=0;             // Beginn der aktuellen Zeile
1005   for (i=1; i<=v[zeile]; i++) {
1006     z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
1007     if ((beginn==0) and (m[zeile,i]!=0)) { beginn=i;}
1008   }
1009   z(zeile)=z(zeile) + z(zeile+1)^v[zeile] * z(zeile+2);
1010   if (beginn==0) {
1011     if (fehler>0) {     // damit fehler=0 bleibt bei exakter Param.
1012     fehlervor=fehler;   // Fehler der letzten Zeile
1013     fehler=fehler+untergrad*(v[zeile]-1)+untervor;   // Fehler dieser Zeile
1014     hilf=untergrad;
1015     untergrad=untergrad*v[zeile]+untervor;
1016     untervor=hilf;}     // untervor = altes untergrad
1017   }
1018   else {
1019     fehlervor=fehler;
1020     fehler=fehler+untergrad*(beginn-1);
1021     untervor=untergrad;
1022     untergrad=untergrad*beginn;
1023   }
1024 }
1025 //--------------------- Ausgabe der Fehlerabschaetzung -----------------------
1026 if (switch==0) {
1027   if (fehler>0) {
1028     if (fehlervor>0) {
1029       dbprint(printlevel-voice+4,""+
1030         "// ** Warning: result is exact up to order "+string(fehlervor-1)+
1031         " in "+ string(var(1))+" and "+string(fehler-1)+" in " +
1032         string(var(2))+" !");
1033     }
1034     else {
1035       dbprint(printlevel-voice+4,""+
1036         "// ** Warning: result is exact up to order "+ string(fehler-1)+
1037         " in "+string(var(2))+" !");
1038     }
1039   }
1040   if (return_list==0) { return(ideal(z(2),z(1))); }
1041   else   { return(list(ideal(z(2),z(1)),intvec(fehlervor-1,fehler-1))); }
1042 }
1043 else {
1044   if (fehler>0) {
1045     if (fehlervor>0) {
1046       dbprint(printlevel-voice+4,""+
1047         "// ** Warning: result is exact up to order "+string(fehler-1)+
1048         " in "+ string(var(1))+" and "+string(fehlervor-1)+" in " +
1049         string(var(2))+" !");
1050     }
1051     else {
1052       dbprint(printlevel-voice+4,""+
1053        "// ** Warning: result is exact up to order "+ string(fehler-1)+
1054         " in "+string(var(1))+" !");
1055     }
1056   }
1057   if (return_list==0) { return(ideal(z(1),z(2))); }
1058   else   { return(list(ideal(z(1),z(2)),intvec(fehler-1,fehlervor-1))); }
1059 }
1060}
1061example
1062{ "EXAMPLE:"; echo = 2;
1063 ring exring=0,(x,y,t),ds;
1064 poly f=x3+2xy2+y2;
1065 list Hne=develop(f);
1066 list hne_extended=extdevelop(Hne,10);
1067            //   compare the HNE matrices ...
1068 print(Hne[1]);
1069 print(hne_extended[1]);
1070            // ... and the resulting parametrizations:
1071 param(Hne);
1072 param(hne_extended);
1073 param(hne_extended,0);
1074
1075 // An example with more than one branch:
1076 list L=hnexpansion(f*(x2+y4));
1077 def HNring = L[1]; setring HNring;
1078 param(hne);
1079}
1080
1081///////////////////////////////////////////////////////////////////////////////
1082
1083proc invariants
1084"USAGE:   invariants(INPUT); INPUT list or poly
1085ASSUME:  @code{INPUT} is the output of @code{develop(f)}, or of
1086         @code{extdevelop(develop(f),n)}, or one entry of the list of HN data
1087         computed by @code{hnexpansion(f[,\"ess\"])}.
1088RETURN:  list @code{INV} of the following format:
1089@format
1090    INV[1]:  intvec    (characteristic exponents)
1091    INV[2]:  intvec    (generators of the semigroup)
1092    INV[3]:  intvec    (Puiseux pairs, 1st components)
1093    INV[4]:  intvec    (Puiseux pairs, 2nd components)
1094    INV[5]:  int       (degree of the conductor)
1095    INV[6]:  intvec    (sequence of multiplicities)
1096@end format
1097         If @code{INPUT} contains no valid HN expansion, the empty list is
1098         returned.
1099ASSUME:  @code{INPUT} is a bivariate polynomial f, or the output of
1100         @code{hnexpansion(f)}, or the list of HN data computed by
1101         @code{hnexpansion(f [,\"ess\"])}.
1102RETURN:  list @code{INV}, such that @code{INV[i]} coincides with the output of
1103         @code{invariants(develop(f[i]))}, where f[i] is the i-th branch of
1104         f, and the last entry of @code{INV} contains further invariants of
1105         f in the format:
1106@format
1107    INV[last][1] : intmat    (contact matrix of the branches)
1108    INV[last][2] : intmat    (intersection multiplicities of the branches)
1109    INV[last][3] : int       (delta invariant of f)
1110@end format
1111NOTE:    In case the Hamburger-Noether expansion of the curve f is needed
1112         for other purposes as well it is better to calculate this first
1113         with the aid of @code{hnexpansion} and use it as input instead of
1114         the polynomial itself.
1115SEE ALSO: hnexpansion, develop, displayInvariants, multsequence, intersection
1116KEYWORDS: characteristic exponents; semigroup of values; Puiseux pairs;
1117          conductor, degree; multiplicities, sequence of
1118EXAMPLE:  example invariants; shows an example
1119"
1120{
1121 //---- INPUT = poly, or HNEring, or hne of reducible curve  -----------------
1122 if (typeof(#[1])!="matrix") {
1123   if (typeof(#[1])=="poly") {
1124      list L=hnexpansion(#[1]);
1125      if (typeof(L[1])=="ring") {
1126        def altring = basering;
1127        def HNring = L[1]; setring HNring;
1128        list Ergebnis = invariants(hne);
1129        setring altring;
1130        kill HNring;
1131        return(Ergebnis);
1132      }
1133      else {
1134        return(invariants(L));
1135      }
1136   }
1137   if (typeof(#[1])=="ring") {
1138     def altring = basering;
1139     def HNring = #[1]; setring HNring;
1140     list Ergebnis = invariants(hne);
1141     setring altring;
1142     kill HNring;
1143     return(Ergebnis);
1144   }
1145   if (typeof(#[1])=="list") {
1146     list hne=#;
1147     list Ergebnis;
1148     for (int lauf=1;lauf<=size(hne);lauf++) {
1149       Ergebnis[lauf]=invariants(hne[lauf]);
1150     }
1151     // Calculate the intersection matrix and the intersection multiplicities.
1152     intmat contact[size(hne)][size(hne)];
1153     intmat intersectionmatrix[size(hne)][size(hne)];
1154     int Lauf;
1155     for (lauf=1;lauf<=size(hne);lauf++) {
1156       for (Lauf=lauf+1;Lauf<=size(hne);Lauf++) {
1157         contact[lauf,Lauf]=separateHNE(hne[lauf],hne[Lauf]);
1158         contact[Lauf,lauf]=contact[lauf,Lauf];
1159         intersectionmatrix[lauf,Lauf]=intersection(hne[lauf],hne[Lauf]);
1160         intersectionmatrix[Lauf,lauf]=intersectionmatrix[lauf,Lauf];
1161       }
1162     }
1163     // Calculate the delta invariant.
1164     int inters;
1165     int del=Ergebnis[size(hne)][5]/2;
1166     for(lauf=1;lauf<=size(hne)-1;lauf++) {
1167       del=del+Ergebnis[lauf][5]/2;
1168       for(Lauf=lauf+1;Lauf<=size(hne);Lauf++) {
1169         inters=inters+intersectionmatrix[lauf,Lauf];
1170       }
1171     }
1172     del=del+inters;
1173     list LAST=contact,intersectionmatrix,del;
1174     Ergebnis[size(hne)+1]=LAST;
1175     return(Ergebnis);
1176   }
1177 }
1178 //-------------------------- Initialisierungen -------------------------------
1179 matrix m=#[1];
1180 intvec v=#[2];
1181 int switch=#[3];
1182 list ergebnis;
1183 if (switch==-1) {
1184   "An error has occurred in develop, so there is no HNE.";
1185   return(ergebnis);
1186 }
1187 intvec beta,s,svorl,ordnung,multseq,mpuiseux,npuiseux,halbgr;
1188 int genus,zeile,i,j,k,summe,conductor,ggT;
1189 string Ausgabe;
1190 int nc=ncols(m); int nr=nrows(m);
1191 ordnung[nr]=1;
1192         // alle Indizes muessen (gegenueber [Ca]) um 1 erhoeht werden,
1193         // weil 0..r nicht als Wertebereich erlaubt ist (aber nrows(m)==r+1)
1194
1195 //---------------- Bestimme den Untergrad der einzelnen Zeilen ---------------
1196 for (zeile=nr; zeile>1; zeile--) {
1197   if ((size(ideal(m[zeile,1..nc])) > 1) or (zeile==nr)) { // keine Nullzeile
1198      k=1;
1199      while (m[zeile,k]==0) {k++;}
1200      ordnung[zeile-1]=k*ordnung[zeile]; // vgl. auch Def. von untergrad in
1201      genus++;                           // proc param
1202      svorl[genus]=zeile;} // werden gerade in umgekehrter Reihenfolge abgelegt
1203   else {
1204      ordnung[zeile-1]=v[zeile]*ordnung[zeile]+ordnung[zeile+1];
1205 }}
1206 //----------------- charakteristische Exponenten (beta) ----------------------
1207 s[1]=1;
1208 for (k=1; k <= genus; k++) { s[k+1]=svorl[genus-k+1];} // s[2]==s(1), u.s.w.
1209 beta[1]=ordnung[1]; //charakt. Exponenten: Index wieder verschoben
1210 for (k=1; k <= genus; k++) {
1211   summe=0;
1212   for (i=1; i <= s[k]; i++) {summe=summe+v[i]*ordnung[i];}
1213   beta[k+1]=summe+ordnung[s[k]]+ordnung[s[k]+1]-ordnung[1];
1214 }
1215 //--------------------------- Puiseuxpaare -----------------------------------
1216 int produkt=1;
1217 for (i=1; i<=genus; i++) {
1218   ggT=gcd(beta[1],beta[i+1]*produkt);
1219   mpuiseux[i]=beta[i+1]*produkt / ggT;
1220   npuiseux[i]=beta[1] / ggT;
1221   produkt=produkt*npuiseux[i];
1222 }
1223 //---------------------- Grad des Konduktors ---------------------------------
1224 summe=1-ordnung[1];
1225 if (genus > 0) {
1226   for (i=2; i <= genus+1; i++) {
1227     summe=summe + beta[i] * (ordnung[s[i-1]] - ordnung[s[i]]);
1228   }                              // n.b.: Indizierung wieder um 1 verschoben
1229 }
1230 conductor=summe;
1231 //------------------- Erzeuger der Halbgruppe: -------------------------------
1232 halbgr=puiseux2generators(mpuiseux,npuiseux);
1233
1234 //------------------- Multiplizitaetensequenz: -------------------------------
1235 k=1;
1236 for (i=1; i<size(v); i++) {
1237   for (j=1; j<=v[i]; j++) {
1238     multseq[k]=ordnung[i];
1239     k++;
1240 }}
1241 multseq[k]=1;
1242 //--- fuelle die Multipl.seq. mit den notwendigen Einsen auf -- T.Keilen ----
1243 int tester=k;
1244 while((multseq[tester]==1) and (tester>1))
1245 {
1246   tester=tester-1;
1247 }
1248 if ((multseq[tester]!=1) and (multseq[tester]!=k-tester))
1249 {
1250   for (i=k+1; i<=tester+multseq[tester]; i++)
1251   {
1252     multseq[i]=1;
1253   }
1254 }
1255 //--- Ende T.Keilen --- 06.05.02
1256 //------------------------- Rueckgabe ----------------------------------------
1257 ergebnis=beta,halbgr,mpuiseux,npuiseux,conductor,multseq;
1258 return(ergebnis);
1259}
1260example
1261{ "EXAMPLE:"; echo = 2;
1262 ring exring=0,(x,y),dp;
1263 list Hne=develop(y4+2x3y2+x6+x5y);
1264 list INV=invariants(Hne);
1265 INV[1];                   // the characteristic exponents
1266 INV[2];                   // the generators of the semigroup of values
1267 INV[3],INV[4];            // the Puiseux pairs in packed form
1268 INV[5] / 2;               // the delta-invariant
1269 INV[6];                   // the sequence of multiplicities
1270                           // To display the invariants more 'nicely':
1271 displayInvariants(Hne);
1272 /////////////////////////////
1273 INV=invariants((x2-y3)*(x3-y5));
1274 INV[1][1];                // the characteristic exponents of the first branch
1275 INV[2][6];                // the sequence of multiplicities of the second branch
1276 print(INV[size(INV)][1]);         // the contact matrix of the branches
1277 print(INV[size(INV)][2]);         // the intersection numbers of the branches
1278 INV[size(INV)][3];                // the delta invariant of the curve
1279}
1280
1281///////////////////////////////////////////////////////////////////////////////
1282
1283proc displayInvariants
1284"USAGE:  displayInvariants(INPUT); INPUT list or poly
1285ASSUME:  @code{INPUT} is a bivariate polynomial, or the output of
1286         @code{develop(f)}, resp. of @code{extdevelop(develop(f),n)}, or (one
1287         entry of) the list of HN data computed by
1288         @code{hnexpansion(f[,\"ess\"])}.
1289RETURN:  none
1290DISPLAY: invariants of the corresponding branch, resp. of all branches,
1291         in a better readable form.
1292NOTE:    If the Hamburger-Noether expansion of the curve f is needed
1293         for other purposes as well it is better to calculate this first
1294         with the aid of @code{hnexpansion} and use it as input instead of
1295         the polynomial itself.
1296SEE ALSO: invariants, intersection, develop, hnexpansion
1297EXAMPLE: example displayInvariants;  shows an example
1298"
1299{
1300 // INPUT = polynomial or ring
1301 if (typeof(#[1])=="poly") {
1302   list L=hnexpansion(#[1]);
1303   if (typeof(L[1])=="ring") {
1304     def HNring = L[1]; setring HNring;
1305     displayInvariants(hne);
1306     return();
1307   }
1308   else {
1309     displayInvariants(L);
1310     return();
1311   }
1312 }
1313 if (typeof(#[1])=="ring")
1314 {
1315   def HNring = #[1]; setring HNring;
1316   displayInvariants(hne);
1317   return();
1318 }
1319 // INPUT = hne of a plane curve
1320 int i,j,k,mul;
1321 string Ausgabe;
1322 list ergebnis;
1323 //-- entferne ueberfluessige Daten zur Erhoehung der Rechengeschwindigkeit: --
1324 #=stripHNE(#);
1325 //-------------------- Ausgabe eines Zweiges ---------------------------------
1326 if (typeof(#[1])=="matrix") {
1327   ergebnis=invariants(#);
1328   if (size(ergebnis)!=0) {
1329    " characteristic exponents  :",ergebnis[1];
1330    " generators of semigroup   :",ergebnis[2];
1331    if (size(ergebnis[1])>1) {
1332     for (i=1; i<=size(ergebnis[3]); i++) {
1333       Ausgabe=Ausgabe+"("+string(ergebnis[3][i])+","
1334       +string(ergebnis[4][i])+")";
1335    }}
1336    " Puiseux pairs             :",Ausgabe;
1337    " degree of the conductor   :",ergebnis[5];
1338    " delta invariant           :",ergebnis[5]/2;
1339    " sequence of multiplicities:",ergebnis[6];
1340 }}
1341 //-------------------- Ausgabe aller Zweige ----------------------------------
1342 else {
1343  ergebnis=invariants(#);
1344  intmat contact=ergebnis[size(#)+1][1];
1345  intmat intersectionmatrix=ergebnis[size(#)+1][2];
1346  for (j=1; j<=size(#); j++) {
1347    " --- invariants of branch number",j,": ---";
1348    " characteristic exponents  :",ergebnis[j][1];
1349    " generators of semigroup   :",ergebnis[j][2];
1350    Ausgabe="";
1351    if (size(ergebnis[j][1])>1) {
1352     for (i=1; i<=size(ergebnis[j][3]); i++) {
1353       Ausgabe=Ausgabe+"("+string(ergebnis[j][3][i])+","
1354       +string(ergebnis[j][4][i])+")";
1355    }}
1356    " Puiseux pairs             :",Ausgabe;
1357    " degree of the conductor   :",ergebnis[j][5];
1358    " delta invariant           :",ergebnis[j][5]/2;
1359    " sequence of multiplicities:",ergebnis[j][6];
1360    "";
1361  }
1362  if (size(#)>1)
1363  {
1364    " -------------- contact numbers : -------------- ";"";
1365    Ausgabe="branch |   ";
1366    for (j=size(#); j>1; j--)
1367    {
1368      if (size(string(j))==1) { Ausgabe=Ausgabe+" "+string(j)+"    "; }
1369      else                    { Ausgabe=Ausgabe+string(j)+"    "; }
1370    }
1371    Ausgabe;
1372    Ausgabe="-------+";
1373    for (j=2; j<size(#); j++) { Ausgabe=Ausgabe+"------"; }
1374    Ausgabe=Ausgabe+"-----";
1375    Ausgabe;
1376  }
1377  for (j=1; j<size(#); j++)
1378  {
1379    if (size(string(j))==1) { Ausgabe="    "+string(j)+"  |"; }
1380    else                    { Ausgabe="   " +string(j)+"  |"; }
1381    for (k=size(#); k>j; k--)
1382    {
1383      mul=contact[j,k];//separateHNE(#[j],#[k]);
1384      for (i=1; i<=5-size(string(mul)); i++) { Ausgabe=Ausgabe+" "; }
1385      Ausgabe=Ausgabe+string(mul);
1386      if (k>j+1) { Ausgabe=Ausgabe+","; }
1387    }
1388    Ausgabe;
1389  }
1390  "";
1391  if (size(#)>1)
1392  {
1393    " -------------- intersection multiplicities : -------------- ";"";
1394    Ausgabe="branch |   ";
1395    for (j=size(#); j>1; j--)
1396    {
1397      if (size(string(j))==1) { Ausgabe=Ausgabe+" "+string(j)+"    "; }
1398      else                    { Ausgabe=Ausgabe+string(j)+"    "; }
1399    }
1400    Ausgabe;
1401    Ausgabe="-------+";
1402    for (j=2; j<size(#); j++) { Ausgabe=Ausgabe+"------"; }
1403    Ausgabe=Ausgabe+"-----";
1404    Ausgabe;
1405  }
1406  for (j=1; j<size(#); j++)
1407  {
1408    if (size(string(j))==1) { Ausgabe="    "+string(j)+"  |"; }
1409    else                    { Ausgabe="   " +string(j)+"  |"; }
1410    for (k=size(#); k>j; k--)
1411    {
1412      mul=intersectionmatrix[j,k];//intersection(#[j],#[k]);
1413      for (i=1; i<=5-size(string(mul)); i++) { Ausgabe=Ausgabe+" "; }
1414      Ausgabe=Ausgabe+string(mul);
1415      if (k>j+1) { Ausgabe=Ausgabe+","; }
1416    }
1417    Ausgabe;
1418  }
1419  "";
1420  " -------------- delta invariant of the curve : ",ergebnis[size(#)+1][3];
1421
1422 }
1423 return();
1424}
1425example
1426{ "EXAMPLE:"; echo = 2;
1427 ring exring=0,(x,y),dp;
1428 list Hne=develop(y4+2x3y2+x6+x5y);
1429 displayInvariants(Hne);
1430}
1431///////////////////////////////////////////////////////////////////////////////
1432
1433proc multiplicities
1434"USAGE:   multiplicities(L); L list
1435ASSUME:  L is the output of @code{develop(f)}, or of
1436         @code{extdevelop(develop(f),n)}, or one entry in the list @code{hne}
1437         in the ring created by @code{hnexpansion(f[,\"ess\"])}.
1438RETURN:  intvec of the different multiplicities that occur when successively
1439         blowing-up the curve singularity corresponding to f.
1440SEE ALSO: multsequence, develop
1441EXAMPLE: example multiplicities;  shows an example
1442"
1443{
1444 matrix m=#[1];
1445 intvec v=#[2];
1446 int switch=#[3];
1447 list ergebnis;
1448 if (switch==-1) {
1449   "An error has occurred in develop, so there is no HNE.";
1450   return(intvec(0));
1451 }
1452 intvec ordnung;
1453 int zeile,k;
1454 int nc=ncols(m); int nr=nrows(m);
1455 ordnung[nr]=1;
1456 //---------------- Bestimme den Untergrad der einzelnen Zeilen ---------------
1457 for (zeile=nr; zeile>1; zeile--) {
1458   if ((size(ideal(m[zeile,1..nc])) > 1) or (zeile==nr)) { // keine Nullzeile
1459      k=1;
1460      while (m[zeile,k]==0) {k++;}
1461      ordnung[zeile-1]=k*ordnung[zeile];
1462   }
1463   else {
1464      ordnung[zeile-1]=v[zeile]*ordnung[zeile]+ordnung[zeile+1];
1465 }}
1466 return(ordnung);
1467}
1468example
1469{ "EXAMPLE:"; echo = 2;
1470  int p=printlevel; printlevel=-1;
1471  ring r=0,(x,y),dp;
1472  multiplicities(develop(x5+y7));
1473  // The first value is the multiplicity of the curve itself, here it's 5
1474  printlevel=p;
1475}
1476///////////////////////////////////////////////////////////////////////////////
1477
1478proc puiseux2generators (intvec m, intvec n)
1479"USAGE:   puiseux2generators(m,n); m,n intvec
1480ASSUME:  m, resp. n, represent the 1st, resp. 2nd, components of Puiseux pairs
1481         (e.g., @code{m=invariants(L)[3]}, @code{n=invariants(L)[4]}).
1482RETURN:  intvec of the generators of the semigroup of values.
1483SEE ALSO: invariants
1484EXAMPLE: example puiseux2generators;  shows an example
1485"
1486{
1487 intvec beta;
1488 int q=1;
1489 //------------ glatte Kurve (eigentl. waeren m,n leer): ----------------------
1490 if (m==0) { return(intvec(1)); }
1491 //------------------- singulaere Kurve: --------------------------------------
1492 for (int i=1; i<=size(n); i++) { q=q*n[i]; }
1493 beta[1]=q; // == q_0
1494 m=1,m; n=1,n; // m[1] ist damit m_0 usw., genau wie beta[1]==beta_0
1495 for (i=2; i<=size(n); i++) {
1496  beta[i]=m[i]*q/n[i] - m[i-1]*q + n[i-1]*beta[i-1];
1497  q=q/n[i]; // == q_i
1498 }
1499 return(beta);
1500}
1501example
1502{ "EXAMPLE:"; echo = 2;
1503  // take (3,2),(7,2),(15,2),(31,2),(63,2),(127,2) as Puiseux pairs:
1504  puiseux2generators(intvec(3,7,15,31,63,127),intvec(2,2,2,2,2,2));
1505}
1506///////////////////////////////////////////////////////////////////////////////
1507
1508proc intersection (list hn1, list hn2)
1509"USAGE:   intersection(hne1,hne2); hne1, hne2 lists
1510ASSUME: @code{hne1, hne2} represent an HN expansion of an irreducible plane
1511        curve singularity (that is, are the output of @code{develop(f)}, or of
1512        @code{extdevelop(develop(f),n)}, or one entry of the list of HN data
1513        computed by @code{hnexpansion(f[,\"ess\"])}).
1514RETURN:  int, the intersection multiplicity of the irreducible plane curve
1515         singularities corresponding to @code{hne1} and @code{hne2}.
1516SEE ALSO: hnexpansion, displayInvariants
1517KEYWORDS: intersection multiplicity
1518EXAMPLE: example intersection;  shows an example
1519"
1520{
1521 //------------------ `intersect' ist schon reserviert ... --------------------
1522 int i,j,s,sum,schnitt,unterschied;
1523 matrix a1=hn1[1];
1524 matrix a2=hn2[1];
1525 intvec h1=hn1[2];
1526 intvec h2=hn2[2];
1527 intvec n1=multiplicities(hn1);
1528 intvec n2=multiplicities(hn2);
1529 if (hn1[3]!=hn2[3]) {
1530 //-- die jeweils erste Zeile von hn1,hn2 gehoert zu verschiedenen Parametern -
1531 //---------------- d.h. beide Kurven schneiden sich transversal --------------
1532   schnitt=n1[1]*n2[1];        // = mult(hn1)*mult(hn2)
1533 }
1534 else {
1535 //--------- die jeweils erste Zeile gehoert zum gleichen Parameter -----------
1536   unterschied=0;
1537   for (s=1; (h1[s]==h2[s]) && (s<size(h1)) && (s<size(h2))
1538              && (unterschied==0); s++) {
1539     for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]); i++) {;}
1540     if (i<=h1[s]) {
1541       unterschied=1;
1542       s--;                    // um s++ am Schleifenende wieder auszugleichen
1543     }
1544   }
1545   if (unterschied==0) {
1546     if ((s<size(h1)) && (s<size(h2))) {
1547       for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]) && (i<=h2[s]); i++) {;}
1548     }
1549     else {
1550 //-------------- Sonderfall: Unterschied in letzter Zeile suchen -------------
1551 // Beachte: Es koennen undefinierte Stellen auftreten, bei abbrechender HNE
1552 // muss die Ende-Markierung weg, h_[r] ist unendlich, die Matrix muss mit
1553 // Nullen fortgesetzt gedacht werden
1554 //----------------------------------------------------------------------------
1555       if (ncols(a1)>ncols(a2)) { j=ncols(a1); }
1556       else                     { j=ncols(a2); }
1557       unterschied=0;
1558       if ((h1[s]>0) && (s==size(h1))) {
1559         a1[s,h1[s]+1]=0;
1560         if (ncols(a1)<=ncols(a2)) { unterschied=1; }
1561       }
1562       if ((h2[s]>0) && (s==size(h2))) {
1563         a2[s,h2[s]+1]=0;
1564         if (ncols(a2)<=ncols(a1)) { unterschied=1; }
1565       }
1566       if (unterschied==1) {                   // mind. eine HNE war endlich
1567         matrix ma1[1][j]=a1[s,1..ncols(a1)];  // und bedarf der Fortsetzung
1568         matrix ma2[1][j]=a2[s,1..ncols(a2)];  // mit Nullen
1569       }
1570       else {
1571         if (ncols(a1)>ncols(a2)) { j=ncols(a2); }
1572         else                     { j=ncols(a1); }
1573         matrix ma1[1][j]=a1[s,1..j];          // Beschr. auf vergleichbaren
1574         matrix ma2[1][j]=a2[s,1..j];          // Teil (der evtl. y's enth.)
1575       }
1576       for (i=1; (ma1[1,i]==ma2[1,i]) && (i<j) && (ma1[1,i]!=var(2)); i++) {;}
1577       if (ma1[1,i]==ma2[1,i]) {
1578         "//** The two HNE's are identical!";
1579         "//** You have either tried to intersect a branch with itself,";
1580         "//** or the two branches have been developed separately.";
1581         "//   In the latter case use `extdevelop' to extend the HNE's until",
1582         "they differ.";
1583         return(-1);
1584       }
1585       if ((ma1[1,i]==var(2)) || (ma2[1,i]==var(2))) {
1586         "//** The two HNE's are (so far) identical. This is because they",
1587         "have been";
1588         "//** computed separately. I need more data; use `extdevelop' to",
1589         "extend them,";
1590         if (ma1[1,i]==var(2)) {"//** at least the first one.";}
1591         else                  {"//** at least the second one.";}
1592         return(-1);
1593       }
1594     }
1595   }
1596   sum=0;
1597   h1[size(h1)]=ncols(a1)+42;        // Ersatz fuer h1[r]=infinity
1598   h2[size(h2)]=ncols(a2)+42;
1599   for (j=1; j<s; j++) {sum=sum+h1[j]*n1[j]*n2[j];}
1600   if ((i<=h1[s]) && (i<=h2[s]))    { schnitt=sum+i*n1[s]*n2[s]; }
1601   if (i==h2[s]+1) { schnitt=sum+h2[s]*n1[s]*n2[s]+n2[s+1]*n1[s]; }
1602   if (i==h1[s]+1) { schnitt=sum+h1[s]*n2[s]*n1[s]+n1[s+1]*n2[s]; }
1603 // "s:",s-1,"i:",i,"S:",sum;
1604 }
1605 return(schnitt);
1606}
1607example
1608{
1609  "EXAMPLE:"; echo = 2;
1610  ring r=0,(x,y),dp;
1611  list Hne=hnexpansion((x2-y3)*(x2+y3));
1612  intersection(Hne[1],Hne[2]);
1613}
1614///////////////////////////////////////////////////////////////////////////////
1615
1616proc multsequence
1617"USAGE:   multsequence(INPUT); INPUT list or poly
1618ASSUME:  @code{INPUT} is the output of @code{develop(f)}, or of
1619         @code{extdevelop(develop(f),n)}, or one entry of the list of HN data
1620         computed by @code{hnexpansion(f[,\"ess\"])}.
1621RETURN:  intvec corresponding to the multiplicity sequence of the irreducible
1622         plane curve singularity described by the HN data (return value
1623         coincides with @code{invariants(INPUT)[6]}).
1624
1625ASSUME:  @code{INPUT} is a bivariate polynomial f, or the output of
1626         @code{hnexpansion(f)}, or the list of HN data computed by
1627         @code{hnexpansion(f [,\"ess\"])}.
1628RETURN:  list of two integer matrices:
1629@texinfo
1630@table @asis
1631@item  @code{multsequence(INPUT)[1][i,*]}
1632   contains the multiplicities of the branches at their infinitely near point
1633   of 0 in its (i-1) order neighbourhood (i.e., i=1: multiplicity of the
1634   branches themselves, i=2: multiplicity of their 1st quadratic transform,
1635   etc., @*
1636   Hence, @code{multsequence(INPUT)[1][*,j]} is the multiplicity sequence
1637   of branch j.
1638@item  @code{multsequence(INPUT)[2][i,*]}:
1639   contains the information which of these infinitely near points coincide.
1640@end table
1641@end texinfo
1642NOTE:  The order of the elements of the list of HN data obtained from
1643       @code{hnexpansion(f [,\"ess\"])} must not be changed (because otherwise
1644       the coincident infinitely near points couldn't be grouped together,
1645       see the meaning of the 2nd intmat in the example).
1646       Hence, it is not wise to compute the HN expansion of polynomial factors
1647       separately, put them into a list INPUT and call
1648       @code{multsequence(INPUT)}. @*
1649       Use @code{displayMultsequence} to produce a better readable output for
1650       reducible curves on the screen. @*
1651       In case the Hamburger-Noether expansion of the curve f is needed
1652       for other purposes as well it is better to calculate this first
1653       with the aid of @code{hnexpansion} and use it as input instead of
1654       the polynomial itself.
1655SEE ALSO: displayMultsequence, develop, hnexpansion, separateHNE
1656KEYWORDS: multiplicity sequence
1657EXAMPLE: example multsequence;  shows an example
1658"
1659{
1660 //---- INPUT = poly, or HNEring --------------------
1661 if (typeof(#[1])=="poly") {
1662   list L=hnexpansion(#[1]);
1663   if (typeof(L[1])=="ring") {
1664     def altring = basering;
1665     def HNring = L[1]; setring HNring;
1666     list Ergebnis = multsequence(hne);
1667     setring altring;
1668     kill HNring;
1669     return(Ergebnis);
1670   }
1671   else {
1672     return(multsequence(L));
1673   }
1674 }
1675 if (typeof(#[1])=="ring") {
1676   def altring = basering;
1677   def HNring = #[1]; setring HNring;
1678   list Ergebnis = multsequence(hne);
1679   setring altring;
1680   kill HNring;
1681   return(Ergebnis);
1682 }
1683 //-- entferne ueberfluessige Daten zur Erhoehung der Rechengeschwindigkeit: --
1684 #=stripHNE(#);
1685 int k,i,j;
1686 //----------------- Multiplizitaetensequenz eines Zweiges --------------------
1687 if (typeof(#[1])=="matrix") {
1688  intvec v=#[2];
1689  list ergebnis;
1690  if (#[3]==-1) {
1691    "An error has occurred in develop, so there is no HNE.";
1692   return(intvec(0));
1693  }
1694  intvec multips,multseq;
1695  multips=multiplicities(#);
1696  k=1;
1697  for (i=1; i<size(v); i++) {
1698    for (j=1; j<=v[i]; j++) {
1699      multseq[k]=multips[i];
1700      k++;
1701  }}
1702  multseq[k]=1;
1703  //--- fuelle die Multipl.seq. mit den notwendigen Einsen auf -- T.Keilen ----
1704  int tester=k;
1705  while((multseq[tester]==1) and (tester>1))
1706  {
1707    tester=tester-1;
1708  }
1709  if((multseq[tester]!=1) and (multseq[tester]!=k-tester))
1710  {
1711    for (i=k+1; i<=tester+multseq[tester]; i++)
1712    {
1713      multseq[i]=1;
1714    }
1715  }
1716  //--- Ende T.Keilen --- 06.05.02
1717  return(multseq);
1718 }
1719 //---------------------------- mehrere Zweige --------------------------------
1720 else {
1721   list HNEs=#;
1722   int anzahl=size(HNEs);
1723   int maxlength=0;
1724   int bisher;
1725   intvec schnitt,ones;
1726   ones[anzahl]=0;
1727   ones=ones+1;                  // = 1,1,...,1
1728   for (i=1; i<anzahl; i++) {
1729     schnitt[i]=separateHNE(HNEs[i],HNEs[i+1]);
1730     j=size(multsequence(HNEs[i]));
1731     maxlength=maxlength*(j < maxlength) + j*(j >= maxlength);
1732     maxlength=maxlength*(schnitt[i]+1 < maxlength)
1733               + (schnitt[i]+1)*(schnitt[i]+1 >= maxlength);
1734   }
1735   j=size(multsequence(HNEs[anzahl]));
1736   maxlength=maxlength*(j < maxlength) + j*(j >= maxlength);
1737
1738//-------------- Konstruktion der ersten zu berechnenden Matrix ---------------
1739   intmat allmults[maxlength][anzahl];
1740   for (i=1; i<=maxlength; i++)  { allmults[i,1..anzahl]=ones[1..anzahl]; }
1741   for (i=1; i<=anzahl; i++) {
1742     ones=multsequence(HNEs[i]);
1743     allmults[1..size(ones),i]=ones[1..size(ones)];
1744   }
1745//---------------------- Konstruktion der zweiten Matrix ----------------------
1746   intmat separate[maxlength][anzahl];
1747   for (i=1; i<=maxlength; i++) {
1748     k=1;
1749     bisher=0;
1750     if (anzahl==1) { separate[i,1]=1; }
1751     for (j=1; j<anzahl; j++)   {
1752       if (schnitt[j]<i) {
1753         separate[i,k]=j-bisher;
1754         bisher=j;
1755         k++;
1756       }
1757       separate[i,k]=anzahl-bisher;
1758     }
1759   }
1760  return(list(allmults,separate));
1761 }
1762}
1763example
1764{
1765  "EXAMPLE:"; echo = 2;
1766  ring r=0,(x,y),dp;
1767  list Hne=hnexpansion((x6-y10)*(x+y2-y3)*(x+y2+y3));
1768  multsequence(Hne[1]),"  |  ",multsequence(Hne[2]),"  |  ",
1769  multsequence(Hne[3]),"  |  ",multsequence(Hne[4]);
1770  multsequence(Hne);
1771  // The meaning of the entries of the 2nd matrix is as follows:
1772  displayMultsequence(Hne);
1773}
1774///////////////////////////////////////////////////////////////////////////////
1775
1776proc displayMultsequence
1777"USAGE:   displayMultsequence(INPUT); INPUT list or poly
1778ASSUME:  @code{INPUT} is a bivariate polynomial, or the output of
1779         @code{develop(f)}, resp. of @code{extdevelop(develop(f),n)}, or (one
1780         entry of) the list of HN data computed by @code{hnexpansion(f[,\"ess\"])},
1781         or the output of @code{hnexpansion(f)}.
1782RETURN:  nothing
1783DISPLAY: the sequence of multiplicities:
1784@format
1785 - if @code{INPUT=develop(f)} or @code{INPUT=extdevelop(develop(f),n)} or @code{INPUT=hne[i]}:
1786                      @code{a , b , c , ....... , 1}
1787 - if @code{INPUT=f} or @code{INPUT=hnexpansion(f)} or @code{INPUT=hne}:
1788                      @code{[(a_1, .... , b_1 , .... , c_1)],}
1789                      @code{[(a_2, ... ), ... , (... , c_2)],}
1790                      @code{ ........................................ ,}
1791                      @code{[(a_n),(b_n), ....., (c_n)]}
1792     with:
1793       @code{a_1 , ... , a_n} the sequence of multiplicities of the 1st branch,
1794       @code{[...]} the multiplicities of the j-th transform of all branches,
1795       @code{(...)} indicating branches meeting in an infinitely near point.
1796@end format
1797NOTE:  The Same restrictions as in @code{multsequence} apply for the input.@*
1798       In case the Hamburger-Noether expansion of the curve f is needed
1799       for other purposes as well it is better to calculate this first
1800       with the aid of @code{hnexpansion} and use it as input instead of
1801       the polynomial itself.
1802SEE ALSO: multsequence, develop, hnexpansion, separateHNE
1803EXAMPLE: example displayMultsequence;  shows an example
1804"
1805{
1806 //---- INPUT = poly, or HNEring --------------------
1807 if (typeof(#[1])=="poly") {
1808   list L=hnexpansion(#[1]);
1809   if (typeof(L[1])=="ring") {
1810     def HNring = L[1]; setring HNring;
1811     displayMultsequence(hne);
1812     return();
1813   }
1814   else {
1815     displayMultsequence(L);
1816     return();
1817   }
1818 }
1819 if (typeof(#[1])=="ring") {
1820   def HNring = #[1]; setring HNring;
1821   displayMultsequence(hne);
1822   return();
1823 }
1824
1825 //-- entferne ueberfluessige Daten zur Erhoehung der Rechengeschwindigkeit: --
1826 #=stripHNE(#);
1827 //----------------- Multiplizitaetensequenz eines Zweiges --------------------
1828 if (typeof(#[1])=="matrix") {
1829   if (#[3]==-1) {
1830     "An error has occurred in develop, so there is no HNE.";
1831   }
1832   else {
1833     "The sequence of multiplicities is  ",multsequence(#);
1834 }}
1835 //---------------------------- mehrere Zweige --------------------------------
1836 else {
1837   list multips=multsequence(#);
1838   int i,j,k,l;
1839   string output;
1840   for (i=1; i<=nrows(multips[1]); i++) {
1841     output="[";
1842     k=1;
1843     for (l=1; k<=ncols(multips[1]); l++) {
1844       output=output+"(";
1845       for (j=1; j<=multips[2][i,l]; j++) {
1846         output=output+string(multips[1][i,k]);
1847         k++;
1848         if (j<multips[2][i,l]) { output=output+","; }
1849       }
1850       output=output+")";
1851       if ((k-1) < ncols(multips[1])) { output=output+","; }
1852      }
1853     output=output+"]";
1854     if (i<nrows(multips[1])) { output=output+","; }
1855     output;
1856   }
1857 }
1858}
1859example
1860{
1861  "EXAMPLE:"; echo = 2;
1862  ring r=0,(x,y),dp;
1863  // Example 1: Input = output of develop
1864  displayMultsequence(develop(x3-y5));
1865
1866  // Example 2: Input = bivariate polynomial
1867  displayMultsequence((x6-y10)*(x+y2-y3)*(x+y2+y3));
1868}
1869
1870///////////////////////////////////////////////////////////////////////////////
1871
1872proc separateHNE (list hn1,list hn2)
1873"USAGE:    separateHNE(hne1,hne2);  hne1, hne2 lists
1874ASSUME:   hne1, hne2 are HNEs (=output of
1875          @code{develop(f)}, @code{extdevelop(develop(f),n)}, or
1876          one entry in the list @code{hne} in the ring created by
1877          @code{hnexpansion(f[,\"ess\"])}.
1878RETURN:   number of quadratic transformations needed to separate both curves
1879          (branches).
1880SEE ALSO: develop, hnexpansion, multsequence, displayMultsequence
1881EXAMPLE:  example separateHNE;  shows an example
1882"
1883{
1884 int i,j,s,unterschied,separated;
1885 matrix a1=hn1[1];
1886 matrix a2=hn2[1];
1887 intvec h1=hn1[2];
1888 intvec h2=hn2[2];
1889 if (hn1[3]!=hn2[3]) {
1890 //-- die jeweils erste Zeile von hn1,hn2 gehoert zu verschiedenen Parametern -
1891 //---------------- d.h. beide Kurven schneiden sich transversal --------------
1892   separated=1;
1893 }
1894 else {
1895 //--------- die jeweils erste Zeile gehoert zum gleichen Parameter -----------
1896   unterschied=0;
1897   for (s=1; (h1[s]==h2[s]) && (s<size(h1)) && (s<size(h2))
1898              && (unterschied==0); s++) {
1899     for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]); i++) {;}
1900     if (i<=h1[s]) {
1901       unterschied=1;
1902       s--;                    // um s++ am Schleifenende wieder auszugleichen
1903     }
1904   }
1905   if (unterschied==0) {
1906     if ((s<size(h1)) && (s<size(h2))) {
1907       for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]) && (i<=h2[s]); i++) {;}
1908     }
1909     else {
1910 //-------------- Sonderfall: Unterschied in letzter Zeile suchen -------------
1911 // Beachte: Es koennen undefinierte Stellen auftreten, bei abbrechender HNE
1912 // muss die Ende-Markierung weg, h_[r] ist unendlich, die Matrix muss mit
1913 // Nullen fortgesetzt gedacht werden
1914 //----------------------------------------------------------------------------
1915       if (ncols(a1)>ncols(a2)) { j=ncols(a1); }
1916       else                     { j=ncols(a2); }
1917       unterschied=0;
1918       if ((h1[s]>0) && (s==size(h1))) {
1919         a1[s,h1[s]+1]=0;
1920         if (ncols(a1)<=ncols(a2)) { unterschied=1; }
1921       }
1922       if ((h2[s]>0) && (s==size(h2))) {
1923         a2[s,h2[s]+1]=0;
1924         if (ncols(a2)<=ncols(a1)) { unterschied=1; }
1925       }
1926       if (unterschied==1) {                   // mind. eine HNE war endlich
1927         matrix ma1[1][j]=a1[s,1..ncols(a1)];  // und bedarf der Fortsetzung
1928         matrix ma2[1][j]=a2[s,1..ncols(a2)];  // mit Nullen
1929       }
1930       else {
1931         if (ncols(a1)>ncols(a2)) { j=ncols(a2); }
1932         else                     { j=ncols(a1); }
1933         matrix ma1[1][j]=a1[s,1..j];          // Beschr. auf vergleichbaren
1934         matrix ma2[1][j]=a2[s,1..j];          // Teil (der evtl. y's enth.)
1935       }
1936       for (i=1; (ma1[1,i]==ma2[1,i]) && (i<j) && (ma1[1,i]!=var(2)); i++) {;}
1937       if (ma1[1,i]==ma2[1,i]) {
1938         "//** The two HNE's are identical!";
1939         "//** You have either tried to compare a branch with itself,";
1940         "//** or the two branches have been developed separately.";
1941         "//   In the latter case use `extdevelop' to extend the HNE's until",
1942         "they differ.";
1943         return(-1);
1944       }
1945       if ((ma1[1,i]==var(2)) || (ma2[1,i]==var(2))) {
1946         "//** The two HNE's are (so far) identical. This is because they",
1947         "have been";
1948         "//** computed separately. I need more data; use `extdevelop' to",
1949         "extend them,";
1950         if (ma1[1,i]==var(2)) {"//** at least the first one.";}
1951         else                  {"//** at least the second one.";}
1952         return(-1);
1953       }
1954     }
1955   }
1956   separated=i;
1957   for (j=1; j<s; j++) { separated=separated+h1[j]; }
1958 }
1959 return(separated);
1960}
1961example
1962{ "EXAMPLE:"; echo = 2;
1963  int p=printlevel; printlevel=-1;
1964  ring r=0,(x,y),dp;
1965  list hne1=develop(x);
1966  list hne2=develop(x+y);
1967  list hne3=develop(x+y2);
1968  separateHNE(hne1,hne2);  // two transversal lines
1969  separateHNE(hne1,hne3);  // one quadratic transform. gives 1st example
1970  printlevel=p;
1971}
1972///////////////////////////////////////////////////////////////////////////////
1973
1974proc displayHNE(list ldev,list #)
1975"USAGE:   displayHNE(L[,n]); L list, n int
1976ASSUME:  L is the output of @code{develop(f)}, or of @code{exdevelop(f,n)},
1977         or of @code{hnexpansion(f[,\"ess\"])}, or (one entry in) the list
1978         @code{hne} in the ring created by @code{hnexpansion(f[,\"ess\"])}.
1979RETURN:  - if only one argument is given and if the input are the HN data
1980         of an irreducible plane curve singularity, no return value, but
1981         display an ideal HNE of the following form:
1982     @example
1983       y = []*x^1+[]*x^2   +...+x^<>*z(1)
1984       x =        []*z(1)^2+...+z(1)^<>*z(2)
1985       z(1) =     []*z(2)^2+...+z(2)^<>*z(3)
1986       .......             ..........................
1987       z(r-1) =   []*z(r)^2+[]*z(r)^3+......
1988     @end example
1989        where @code{x},@code{y} are the first 2 variables of the basering.
1990        The values of @code{[]} are the coefficients of the Hamburger-Noether
1991        matrix, the values of @code{<>} are represented by @code{x} in the
1992        HN matrix.@*
1993        - if a second argument is given and if the input are the HN data
1994        of an irreducible plane curve singularity, return a ring containing
1995        an ideal @code{HNE} as described above.@*
1996        - if L corresponds to the output of @code{hnexpansion(f)}
1997        or to the list of HN data computed by @code{hnexpansion(f[,\"ess\"])},
1998        @code{displayHNE(L[,n])} shows the HNE's of all branches of f in the
1999        format described above. The optional parameter is then ignored.
2000NOTE:  The 1st line of the above ideal (i.e., @code{HNE[1]}) means that
2001     @code{y=[]*z(0)^1+...}, the 2nd line (@code{HNE[2]}) means that
2002     @code{x=[]*z(1)^2+...}, so you can see which indeterminate
2003     corresponds to which line (it's also possible that @code{x} corresponds
2004     to the 1st line and @code{y} to the 2nd).
2005
2006SEE ALSO: develop, hnexpansion
2007EXAMPLE: example displayHNE; shows an example
2008"
2009{
2010 if ((typeof(ldev[1])=="list") || (typeof(ldev[1])=="none")) {
2011   for (int i=1; i<=size(ldev); i++) {
2012     "// Hamburger-Noether development of branch nr."+string(i)+":";
2013     displayHNE(ldev[i]);"";
2014   }
2015   return();
2016 }
2017 //--------------------- Initialisierungen und Ringwechsel --------------------
2018 matrix m=ldev[1];
2019 intvec v=ldev[2];
2020 int switch=ldev[3];
2021 if (switch==-1) {
2022   "An error has occurred throughout the expansion, so there is no HNE.";
2023   return(ideal(0));
2024 }
2025 def altring=basering;
2026 /////////////////////////////////////////////////////////
2027 //  Change by T. Keilen 08.06.2002
2028 //  ring + ring does not work if one ring is an algebraic extension
2029/*
2030 if (parstr(basering)!="") {
2031   if (charstr(basering)!=string(char(basering))+","+parstr(basering)) {
2032     execute
2033      ("ring dazu=("+charstr(basering)+"),z(0.."+string(size(v)-1)+"),ls;");
2034   }
2035   else { ring dazu=char(altring),z(0..size(v)-1),ls; }
2036 }
2037 else   { ring dazu=char(altring),z(0..size(v)-1),ls; }
2038 def displayring=dazu+altring;
2039*/
2040 execute("ring displayring=("+charstr(basering)+"),(z(0.."+string(size(v)-1)+"),"+varstr(basering)+"),(ls("+string(size(v))+"),"+ordstr(basering)+");");
2041 // End change by T. Keilen
2042 //////////////////////////////////////////////////////////////
2043 setring displayring;
2044 map holematrix=altring,0;        // mappt nur die Monome vom Grad Null
2045 matrix m=holematrix(m);
2046 int i,j;
2047
2048 // lossen: check the last row for finiteness (06/2004)
2049 int rowM=nrows(m);
2050 int colM=ncols(m);
2051 int undef_bd=v[size(v)];
2052 if ( undef_bd<-1 ){
2053   for (j=-undef_bd; j<=colM; j++) { m[rowM,j]=0; }
2054 }
2055
2056 //--------------------- Erzeuge Matrix n mit n[i,j]=z(j-1)^i -----------------
2057 matrix n[colM][rowM];
2058 for (j=1; j<=rowM; j++) {
2059    for (i=1; i<=colM; i++) { n[i,j]=z(j-1)^i; }
2060 }
2061 matrix displaymat=m*n;
2062 ideal HNE;
2063 for (i=1; i<rowM; i++) { HNE[i]=displaymat[i,i]+z(i)*z(i-1)^v[i]; }
2064 HNE[rowM]=displaymat[rowM,rowM];
2065
2066 // lossen: output modified (06/2004)
2067 if (size(#) == 0)
2068 {
2069   if (switch==0) {
2070     HNE=subst(HNE,z(0),var(size(v)+1));
2071   }
2072   else  {
2073     HNE=subst(HNE,z(0),var(size(v)+2));
2074   }
2075
2076   for (j=1; j<=ncols(HNE); j++){
2077     string stHNE(j)=string(HNE[j]);
2078   }
2079   if (undef_bd<-1)
2080   {
2081     stHNE(size(v))=stHNE(size(v))+" + ..... (terms of degree >="
2082                                 +string(-undef_bd)+")";
2083   }
2084   if (undef_bd==-1)
2085   {
2086     stHNE(size(v))=stHNE(size(v))+" + ..... (terms of degree >="
2087                                 +string(colM+1)+")";
2088   }
2089
2090   if (switch==0) {
2091     stHNE(1) = "  "+string(var(size(v)+2))+" = "+stHNE(1);
2092   }
2093   else {
2094     stHNE(1) = "  "+string(var(size(v)+1))+" = "+stHNE(1);
2095   }
2096   stHNE(1);
2097   if (ncols(HNE)==1) {return();}
2098
2099   if (switch==0) {
2100     stHNE(2) = "  "+string(var(size(v)+1))+" = "+stHNE(2);
2101   }
2102   else {
2103     stHNE(2) = "  "+string(var(size(v)+2))+" = "+stHNE(2);
2104   }
2105   stHNE(2);
2106
2107   for (j=3; j<=ncols(HNE); j++){
2108     stHNE(j)= "  "+"z(" +string(j-2)+ ") = "+stHNE(j);
2109     stHNE(j);
2110   }
2111   return();
2112 }
2113
2114 if (rowM<2) { HNE[2]=z(0); }
2115
2116 if (switch==0) {
2117    HNE[1] = HNE[1]-var(size(v)+2);
2118    HNE[2] = HNE[2]-var(size(v)+1);
2119 }
2120 else {
2121    HNE[1] = HNE[1]-var(size(v)+1);
2122    HNE[2] = HNE[2]-var(size(v)+2);
2123 }
2124if (size(#) == 0) {
2125   HNE;
2126   return();
2127 }
2128if (size(#) != 0) {
2129   HNE;
2130   export(HNE);
2131   return(displayring);
2132 }
2133}
2134example
2135{ "EXAMPLE:"; echo = 2;
2136  ring r=0,(x,y),dp;
2137  poly f=x3+2xy2+y2;
2138  list hn=develop(f);
2139  displayHNE(hn);
2140}
2141///////////////////////////////////////////////////////////////////////////////
2142//                      procedures for reducible curves                      //
2143///////////////////////////////////////////////////////////////////////////////
2144
2145// proc newtonhoehne (poly f)
2146// USAGE:   newtonhoehne(f);   f poly
2147// ASSUME:  basering = ...,(x,y),ds  or ls
2148// RETURN:  list of intvec(x,y) of coordinates of the newtonpolygon of f
2149// NOTE:    This proc is only available in versions of Singular that know the
2150//         command  system("newton",f);  f poly
2151// {
2152// intvec nm = getnm(f);
2153//  if ((nm[1]>0) && (nm[2]>0)) { f=jet(f,nm[1]*nm[2],nm); }
2154//  list erg=system("newton",f);
2155//  int i; list Ausgabe;
2156//  for (i=1; i<=size(erg); i++) { Ausgabe[i]=leadexp(erg[i]); }
2157// return(Ausgabe);
2158// }
2159///////////////////////////////////////////////////////////////////////////////
2160
2161proc newtonpoly (poly f, int #)
2162"USAGE:   newtonpoly(f);   f poly
2163ASSUME:  basering has exactly two variables; @*
2164         f is convenient, that is, f(x,0) != 0 != f(0,y).
2165RETURN:  list of intvecs (= coordinates x,y of the Newton polygon of f).
2166NOTE:    Procedure uses @code{execute}; this can be avoided by calling
2167         @code{newtonpoly(f,1)} if the ordering of the basering is @code{ls}.
2168KEYWORDS: Newton polygon
2169EXAMPLE: example newtonpoly;  shows an example
2170"
2171{
2172  if (size(#)>=1)
2173  {
2174    if (typeof(#[1])=="int")
2175    {
2176      // this is done to avoid the "execute" command for procedures in
2177      //  hnoether.lib
2178      def is_ls=#[1];
2179    }
2180  }
2181  if (defined(is_ls)<=0)
2182  {
2183    def @Rold=basering;
2184    execute("ring @RR=("+charstr(basering)+"),("+varstr(basering)+"),ls;");
2185    poly f=imap(@Rold,f);
2186  }
2187  intvec A=(0,ord(subst(f,var(1),0)));
2188  intvec B=(ord(subst(f,var(2),0)),0);
2189  intvec C,H; list L;
2190  int abbruch,i;
2191  poly hilf;
2192  L[1]=A;
2193  f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]));
2194  if (defined(is_ls))
2195  {
2196    map xytausch=basering,var(2),var(1);
2197  }
2198  else
2199  {
2200    map xytausch=@RR,var(2),var(1);
2201  }
2202  for (i=2; f!=0; i++)
2203  {
2204     abbruch=0;
2205     while (abbruch==0)
2206     {
2207        C=leadexp(f);
2208        if(jet(f,A[2]*C[1]-A[1]*C[2]-1,intvec(A[2]-C[2],C[1]-A[1]))==0)
2209        {
2210           abbruch=1;
2211        }
2212        else
2213        {
2214           f=jet(f,-C[1]-1,intvec(-1,0));
2215        }
2216    }
2217    hilf=jet(f,A[2]*C[1]-A[1]*C[2],intvec(A[2]-C[2],C[1]-A[1]));
2218    H=leadexp(xytausch(hilf));
2219    A=H[2],H[1];
2220    L[i]=A;
2221    f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]-A[1]));
2222  }
2223  L[i]=B;
2224  if (defined(is_ls))
2225  {
2226    return(L);
2227  }
2228  else
2229  {
2230    setring @Rold;
2231    return(L);
2232  }
2233}
2234example
2235{
2236 "EXAMPLE:"; echo = 2;
2237  ring r=0,(x,y),ls;
2238  poly f=x5+2x3y-x2y2+3xy5+y6-y7;
2239  newtonpoly(f);
2240}
2241///////////////////////////////////////////////////////////////////////////////
2242
2243proc is_NND (poly f, list #)
2244"USAGE:   is_NND(f[,mu,NP]);   f poly, mu int, NP list of intvecs
2245ASSUME:  f is convenient, that is, f(x,0) != 0 != f(0,y);@*
2246         mu (optional) is Milnor number of f.@*
2247         NP (optional) is output of @code{newtonpoly(f)}.
2248RETURN:  int: 1 if f is Newton non-degenerate, 0 otherwise.
2249SEE ALSO: newtonpoly
2250KEYWORDS: Newton non-degenerate; Newton polygon
2251EXAMPLE: example is_NND;  shows examples
2252"
2253{
2254  int i;
2255  int i_print=printlevel-voice+2;
2256
2257  if (size(#)==0)
2258  {
2259    int mu=milnor(f);
2260    list NP=newtonpoly(f);
2261  }
2262  else
2263  {
2264    if (typeof(#[1])=="int")
2265    {
2266      def mu=#[1];
2267      def NP=#[2];
2268      for (i=1;i<=size(NP);i++)
2269      {
2270        if (typeof(NP[i])!="intvec")
2271        {
2272          print("third input cannot be Newton polygon ==> ignored ")
2273          NP=newtonpoly(f);
2274          i=size(NP)+1;
2275        }
2276      }
2277    }
2278    else
2279    {
2280      print("second input cannot be Milnor number ==> ignored ")
2281      int mu=milnor(f);
2282      NP=newtonpoly(f);
2283    }
2284  }
2285
2286  // computation of the Newton number:
2287  int s=size(NP);
2288  int nN=-NP[1][2]-NP[s][1]+1;
2289  intmat m[2][2];
2290  for(i=1;i<=s-1;i++)
2291  {
2292    m=NP[i+1],NP[i];
2293    nN=nN+det(m);
2294  }
2295
2296  if(mu==nN)
2297  { // the Newton-polygon is non-degenerate
2298    // REFERENCE? (tfuer mehr als 2 Variable gilt nicht, dass mu=nu impliziert,
2299    // dass NP nicht ausgeartet ist!, Siehe KOMMENTAR in equising.lib in esIdeal)
2300    return(1);
2301  }
2302  else
2303  {
2304    return(0);
2305  }
2306}
2307example
2308{
2309 "EXAMPLE:"; echo = 2;
2310  ring r=0,(x,y),ls;
2311  poly f=x5+y3;
2312  is_NND(f);
2313  poly g=(x-y)^5+3xy5+y6-y7;
2314  is_NND(g);
2315
2316  // if already computed, one should give the Minor number and Newton polygon
2317  // as second and third input:
2318  int mu=milnor(g);
2319  list NP=newtonpoly(g);
2320  is_NND(g,mu,NP);
2321}
2322
2323
2324///////////////////////////////////////////////////////////////////////////////
2325
2326proc charPoly(poly f, int M, int N)
2327"USAGE:  charPoly(f,M,N);  f bivariate poly,  M,N int: length and height
2328                          of Newton polygon of f, which has to be only one line
2329RETURN:  the characteristic polynomial of f
2330EXAMPLE: example charPoly;  shows an example
2331"
2332{
2333 poly charp;
2334 int Np=N/ gcd(M,N);
2335 f=subst(f,var(1),1);
2336 for(charp=0; f<>0; f=f-lead(f))
2337  { charp=charp+leadcoef(f)*var(2)^(leadexp(f)[2]/ Np);}
2338 return(charp);
2339}
2340example
2341{ "EXAMPLE:"; echo = 2;
2342  ring exring=0,(x,y),dp;
2343  charPoly(y4+2y3x2-yx6+x8,8,4);
2344  charPoly(y6+3y3x2-x4,4,6);
2345}
2346///////////////////////////////////////////////////////////////////////////////
2347
2348proc find_in_list(list L,int p)
2349"USAGE:   find_in_list(L,p); L: list of intvec(x,y)
2350         (sorted in y: L[1][2]>=L[2][2]), int p >= 0
2351RETURN:  int i: L[i][2]=p if existent; otherwise i with L[i][2]<p if existent;
2352         otherwise i = size(L)+1;
2353EXAMPLE: example find_in_list;  shows an example
2354"
2355{
2356 int i;
2357 L[size(L)+1]=intvec(0,-1);          // falls p nicht in L[.][2] vorkommt
2358 for (i=1; L[i][2]>p; i++) {;}
2359 return(i);
2360}
2361example
2362{ "EXAMPLE:"; echo = 2;
2363  list L = intvec(0,4), intvec(1,2), intvec(2,1), intvec(4,0);
2364  find_in_list(L,1);
2365  L[find_in_list(L,2)];
2366}
2367///////////////////////////////////////////////////////////////////////////////
2368
2369proc get_last_divisor(int M, int N)
2370"USAGE:   get_last_divisor(M,N); int M,N
2371RETURN:  int Q: M=q1*N+r1, N=q2*r1+r2, ..., ri=Q*r(i+1) (Euclidean alg.)
2372EXAMPLE: example get_last_divisor; shows an example
2373"
2374{
2375 int R=M%N; int Q=M / N;
2376 while (R!=0) {M=N; N=R; R=M%N; Q=M / N;}
2377 return(Q)
2378}
2379example
2380{ "EXAMPLE"; echo = 2;
2381  ring r=0,(x,y),dp;
2382  get_last_divisor(12,10);
2383}
2384///////////////////////////////////////////////////////////////////////////////
2385proc redleit (poly f,intvec S, intvec E)
2386"USAGE:   redleit(f,S,E);  f poly, S,E intvec(x,y)
2387         S,E are two different points on a line in the Newton diagram of f
2388RETURN:  poly g: all monomials of f which lie on or below that line
2389NOTE:    The main purpose is that if the line defined by S and E is part of the
2390         Newton polygon, the result is the quasihomogeneous leading form of f
2391         w.r.t. that line.
2392SEE ALSO: newtonpoly
2393EXAMPLE: example redleit;  shows an example
2394"
2395{
2396 if (E[1]<S[1]) { intvec H=E; E=S; S=H; } // S,E verkehrt herum eingegeben
2397 return(jet(f,E[1]*S[2]-E[2]*S[1],intvec(S[2]-E[2],E[1]-S[1])));
2398}
2399example
2400{ "EXAMPLE"; echo = 2;
2401  ring exring=0,(x,y),dp;
2402  redleit(y6+xy4-2x3y2+x4y+x6,intvec(3,2),intvec(4,1));
2403}
2404///////////////////////////////////////////////////////////////////////////////
2405
2406
2407proc extdevelop (list l, int Exaktheit)
2408"USAGE:   extdevelop(L,N); list L, int N
2409ASSUME:  L is the output of @code{develop(f)}, or of @code{extdevelop(l,n)},
2410         or one entry in the list @code{hne} in the ring created by
2411          @code{hnexpansion(f[,\"ess\"])}.
2412RETURN:  an extension of the Hamburger-Noether development of f as a list
2413         in the same format as L has (up to the last entry in the output
2414         of @code{develop(f)}).@*
2415         Type @code{help develop;}, resp. @code{help hnexpansion;} for more
2416         details.
2417NOTE:    The new HN-matrix will have at least N columns (if the HNE is not
2418         finite). In particular, if f is irreducible then (in most cases)
2419         @code{extdevelop(develop(f),N)} will produce the same result as
2420         @code{develop(f,N)}.@*
2421         If the matrix M of L has n columns then, compared with
2422         @code{parametrization(L)}, @code{paramametrize(extdevelop(L,N))} will increase the
2423         exactness by at least (N-n) more significant monomials.
2424SEE ALSO: develop, hnexpansion, param
2425EXAMPLE: example extdevelop;  shows an example
2426"
2427{
2428 //------------ Initialisierungen und Abfangen unzulaessiger Aufrufe ----------
2429 matrix m=l[1];
2430 intvec v=l[2];
2431 int switch=l[3];
2432 if (nvars(basering) < 2) {
2433   " Sorry. I need two variables in the ring.";
2434   return(list(matrix(maxideal(1)[1]),intvec(0),-1,poly(0)));}
2435 if (switch==-1) {
2436   "An error has occurred in develop, so there is no HNE and no extension.";
2437   return(l);
2438 }
2439 poly f=l[4];
2440 if (f==0) {
2441   " No extension is possible";
2442   return(l);
2443 }
2444 int Q=v[size(v)];
2445 if (Q>0) {
2446   " The HNE was already exact";
2447   return(l);
2448 }
2449 else {
2450   if (Q==-1) { Q=ncols(m); }
2451   else { Q=-Q-1; }
2452 }
2453 int zeile=nrows(m);
2454 int spalten,i,M;
2455 ideal lastrow=m[zeile,1..Q];
2456 int ringwechsel=(varstr(basering)!="x,y") or (ordstr(basering)!="ls(2),C");
2457
2458 //------------------------- Ringwechsel, falls noetig ------------------------
2459 if (ringwechsel) {
2460  def altring = basering;
2461  int p = char(basering);
2462  if (charstr(basering)!=string(p)) {
2463     string tststr=charstr(basering);
2464     tststr=tststr[1..find(tststr,",")-1];     //-> "p^k" bzw. "p"
2465     if (tststr==string(p)) {
2466       if (size(parstr(basering))>1) {         // ring (p,a,..),...
2467        execute("ring extdguenstig=("+charstr(basering)+"),(x,y),ls;");
2468       }
2469       else {                                  // ring (p,a),...
2470        string mipl=string(minpoly);
2471        ring extdguenstig=(p,`parstr(basering)`),(x,y),ls;
2472        if (mipl!="0") { execute("minpoly="+mipl+";"); }
2473       }
2474     }
2475     else {
2476       execute("ring extdguenstig=("+charstr(basering)+"),(x,y),ls;");
2477     }
2478  }
2479  else {                               // charstr(basering)== p : no parameter
2480     ring extdguenstig=p,(x,y),ls;
2481  }
2482  export extdguenstig;
2483  map hole=altring,x,y;
2484 //----- map kann sehr zeitaufwendig sein, daher Vermeidung, wo moeglich: -----
2485  if (nvars(altring)==2) { poly f=fetch(altring,f); }
2486  else                   { poly f=hole(f);          }
2487  ideal a=hole(lastrow);
2488 }
2489 else { ideal a=lastrow; }
2490 list Newton=newtonpoly(f,1);
2491 int M1=Newton[size(Newton)-1][1];     // konstant
2492 number delt;
2493 if (Newton[size(Newton)-1][2]!=1) {
2494    " *** The transformed polynomial was not valid!!";}
2495 else {
2496 //--------------------- Fortsetzung der HNE ----------------------------------
2497  while (Q<Exaktheit) {
2498    M=ord(subst(f,y,0));
2499    Q=M-M1;
2500 //------ quasihomogene Leitform ist c*x^M1*y+d*x^(M1+Q) => delta=-d/c: -------
2501    delt=-koeff(f,M,0)/koeff(f,M1,1);
2502    a[Q]=delt;
2503    dbprint(printlevel-voice+2,"a("+string(zeile-1)+","+string(Q)+") = "+string(delt));
2504    if (Q<Exaktheit) {
2505     f=T1_Transform(f,delt,Q);
2506     if (defined(HNDebugOn)) { "transformed polynomial:",f; }
2507     if (subst(f,y,0)==0) {
2508       dbprint(printlevel-voice+2,"The HNE is finite!");
2509       a[Q+1]=x; Exaktheit=Q;
2510       f=0;                        // Speicherersparnis: f nicht mehr gebraucht
2511     }
2512    }
2513  }
2514 }
2515 //------- Wechsel in alten Ring, Zusammensetzung alte HNE + Erweiterung ------
2516 if (ringwechsel) {
2517  setring altring;
2518  map zurueck=extdguenstig,var(1),var(2);
2519  if (nvars(altring)==2) { f=fetch(extdguenstig,f); }
2520  else                   { f=zurueck(f);            }
2521  lastrow=zurueck(a);
2522 }
2523 else { lastrow=a; }
2524 if (ncols(lastrow)>ncols(m)) { spalten=ncols(lastrow); }
2525 else { spalten=ncols(m); }
2526 matrix mneu[zeile][spalten];
2527 for (i=1; i<nrows(m); i++) {
2528  mneu[i,1..ncols(m)]=m[i,1..ncols(m)];
2529 }
2530 mneu[zeile,1..ncols(lastrow)]=lastrow;
2531 if (lastrow[ncols(lastrow)]!=var(1)) {
2532  if (ncols(lastrow)==spalten) { v[zeile]=-1; }  // keine undefinierten Stellen
2533  else {
2534   v[zeile]=-Q-1;
2535   for (i=ncols(lastrow)+1; i<=spalten; i++) {
2536    mneu[zeile,i]=var(2);           // fuelle nicht def. Stellen der Matrix auf
2537 }}}
2538 else { v[zeile]=Q; }               // HNE war exakt
2539 if (ringwechsel)
2540 {
2541   kill extdguenstig;
2542 }
2543
2544 return(list(mneu,v,switch,f));
2545}
2546example
2547{
2548  "EXAMPLE:"; echo = 2;
2549  ring exring=0,(x,y),dp;
2550  list Hne=hnexpansion(x14-3y2x11-y3x10-y2x9+3y4x8+y5x7+3y4x6+x5*(-y6+y5)
2551                      -3y6x3-y7x2+y8);
2552  displayHNE(Hne);    // HNE of 1st,3rd branch is finite
2553  print(extdevelop(Hne[1],5)[1]);
2554  list ehne=extdevelop(Hne[2],5);
2555  displayHNE(ehne);
2556  param(Hne[2]);
2557  param(ehne);
2558
2559}
2560///////////////////////////////////////////////////////////////////////////////
2561
2562proc stripHNE (list l)
2563"USAGE:   stripHNE(L);  L list
2564ASSUME:  L is the output of @code{develop(f)}, or of
2565         @code{extdevelop(develop(f),n)}, or (one entry of) the list
2566         @code{hne} in the ring created by @code{hnexpansion(f[,\"ess\"])}.
2567RETURN:  list in the same format as L, but all polynomials L[4], resp.
2568         L[i][4], are set to zero.
2569NOTE:    The purpose of this procedure is to remove huge amounts of data
2570         no longer needed. It is useful, if one or more of the polynomials
2571         in L consume much memory. It is still possible to compute invariants,
2572         parametrizations etc. with the stripped HNE(s), but it is not possible
2573         to use @code{extdevelop} with them.
2574SEE ALSO: develop, hnexpansion, extdevelop
2575EXAMPLE: example stripHNE;  shows an example
2576"
2577{
2578 list h;
2579 if (typeof(l[1])=="matrix") { l[4]=poly(0); }
2580 else {
2581  for (int i=1; i<=size(l); i++) {
2582    h=l[i];
2583    h[4]=poly(0);
2584    l[i]=h;
2585  }
2586 }
2587 return(l);
2588}
2589example
2590{
2591  "EXAMPLE:"; echo = 2;
2592  ring r=0,(x,y),dp;
2593  list Hne=develop(x2+y3+y4);
2594  Hne;
2595  stripHNE(Hne);
2596}
2597///////////////////////////////////////////////////////////////////////////////
2598static proc extractHNEs(list HNEs, int transvers)
2599"USAGE:  extractHNEs(HNEs,transvers);  list HNEs (output from HN),
2600        int transvers: 1 if x,y were exchanged, 0 else
2601RETURN: list of Hamburger-Noether-Extensions in the form of hne in hnexpansion
2602NOTE:   This procedure is only for internal purpose; examples don't make sense
2603"
2604{
2605 int i,maxspalte,hspalte,hnezaehler;
2606 list HNEaktu,Ergebnis;
2607 for (hnezaehler=1; hnezaehler<=size(HNEs); hnezaehler++) {
2608  maxspalte=0;
2609  HNEaktu=HNEs[hnezaehler];
2610  if (defined(HNDebugOn)) {"To process:";HNEaktu;}
2611  if (size(HNEaktu)!=size(HNEaktu[1])+2) {
2612     "The ideals and the hqs in HNEs[",hnezaehler,"] don't match!!";
2613     HNEs[hnezaehler];
2614  }
2615 //------------ ermittle  maximale Anzahl benoetigter Spalten: ----------------
2616  for (i=2; i<size(HNEaktu); i++) {
2617    hspalte=ncols(HNEaktu[i]);
2618    maxspalte=maxspalte*(hspalte < maxspalte)+hspalte*(hspalte >= maxspalte);
2619  }
2620 //------------- schreibe Ausgabe fuer hnezaehler-ten Zweig: ------------------
2621  matrix ma[size(HNEaktu)-2][maxspalte];
2622  for (i=1; i<=(size(HNEaktu)-2); i++) {
2623    if (ncols(HNEaktu[i+1]) > 1) {
2624      ma[i,1..ncols(HNEaktu[i+1])]=HNEaktu[i+1]; }
2625    else { ma[i,1]=HNEaktu[i+1][1];}
2626  }
2627  Ergebnis[hnezaehler]=list(ma,HNEaktu[1],transvers,HNEaktu[size(HNEaktu)]);
2628  kill ma;
2629 }
2630 return(Ergebnis);
2631}
2632///////////////////////////////////////////////////////////////////////////////
2633
2634proc factorfirst(poly f, int M, int N)
2635"USAGE : factorfirst(f,M,N); f poly, M,N int
2636RETURN: number d such that f=const*(y^(N/e) - d*x^(M/e))^e, where e=gcd(M,N),
2637        0 if such a d does not exist
2638EXAMPLE: example factorfirst;  shows an example
2639"
2640{
2641 number c = koeff(f,0,N);
2642 number delt;
2643 int eps,l;
2644 int p=char(basering);
2645 string ringchar=charstr(basering);
2646
2647 if (c == 0) {"Something has gone wrong! I didn't get N correctly!"; exit;}
2648 int e = gcd(M,N);
2649
2650 if (p==0) { delt = koeff(f,M/ e,N - N/ e) / (-1*e*c); }
2651 else {
2652   if (e%p != 0) { delt = koeff(f,M/ e,N - N/ e) / (-1*e*c); }
2653   else {
2654     eps = e;
2655     for (l = 0; eps%p == 0; l=l+1) { eps=eps/ p;}
2656     if (defined(HNDebugOn)) {e," -> ",eps,"*",p,"^",l;}
2657     delt = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c);
2658
2659     if ((charstr(basering) != string(p)) and (delt != 0)) {
2660 //------ coefficient field is not Z/pZ => (p^l)th root is not identity -------
2661       delt=0;
2662       if (defined(HNDebugOn)) {
2663         "trivial factorization not implemented for",
2664         "parameters---I've to use 'factorize'";
2665       }
2666     }
2667   }
2668 }
2669 if (defined(HNDebugOn)) {"quasihomogeneous leading form:",f," = ",c,
2670        "* (y^"+string(N/ e),"-",delt,"* x^"+string(M/ e)+")^",e," ?";}
2671 if (f != c*(var(2)^(N/ e) - delt*var(1)^(M/ e))^e) {return(0);}
2672 else {return(delt);}
2673}
2674example
2675{ "EXAMPLE:"; echo = 2;
2676  ring exring=7,(x,y),dp;
2677  factorfirst(2*(y3-3x4)^5,20,15);
2678  factorfirst(x14+y7,14,7);
2679  factorfirst(x14+x8y3+y7,14,7);
2680}
2681
2682///////////////////////////////////////////////////////////////////////////
2683
2684proc hnexpansion(poly f,list #)
2685"USAGE:   hnexpansion(f[,\"ess\"]);   f poly
2686ASSUME:  f is a bivariate polynomial (in the first 2 ring variables)
2687RETURN:  list @code{L}, containing Hamburger-Noether data of @code{f}:
2688         If the computation of the HNE required no field extension, @code{L}
2689         is a list of lists @code{L[i]} (corresponding to the output of
2690         @code{develop}, applied to a branch of @code{f}, but the last entry
2691         being omitted):
2692@texinfo
2693@table @asis
2694@item @code{L[i][1]}; matrix:
2695         Each row contains the coefficients of the corresponding line of the
2696         Hamburger-Noether expansion (HNE) for the i-th branch. The end of
2697         the line is marked in the matrix by the first ring variable
2698         (usually x).
2699@item @code{L[i][2]}; intvec:
2700         indicating the length of lines of the HNE
2701@item @code{L[i][3]}; int:
2702         0  if the 1st ring variable was transversal (with respect to the
2703            i-th branch), @*
2704         1  if the variables were changed at the beginning of the
2705            computation, @*
2706        -1  if an error has occurred.
2707@item @code{L[i][4]}; poly:
2708         the transformed equation of the i-th branch to make it possible
2709         to extend the Hamburger-Noether data a posteriori without having
2710         to do all the previous calculation once again (0 if not needed).
2711@end table
2712@end texinfo
2713         If the computation of the HNE required a field extension, the first
2714         entry @code{L[1]} of the list is a ring, in which a list @code{hne}
2715         of lists (the HN data, as above) and a polynomial @code{f} (image of
2716         @code{f} over the new field) are stored.
2717         @*
2718         If called with an additional input parameter, @code{hnexpansion}
2719         computes only one representative for each class of conjugate
2720         branches (over the ground field active when calling the procedure).
2721         In this case, the returned list @code{L} always has only two
2722         entries: @code{L[1]} is either a list of lists (the HN data) or a
2723         ring (as above), and @code{L[2]} is an integer vector (the number
2724         of branches in the respective conjugacy classes).
2725
2726NOTE:    If f is known to be irreducible as a power series, @code{develop(f)}
2727         could be chosen instead to avoid a change of basering during the
2728         computations. @*
2729         Increasing  @code{printlevel} leads to more and more comments. @*
2730         Having defined a variable @code{HNDebugOn} leads to a maximum
2731         number of comments.
2732
2733SEE ALSO: develop, extdevelop, param, displayHNE
2734EXAMPLE: example hnexpansion;  shows an example
2735"
2736{
2737 int essential;
2738 if (size(#)==1) { essential=1; }
2739 int field_ext;
2740 def altring=basering;
2741
2742 //--------- Falls Ring (p^k,a),...: Wechsel in (p,a),... + minpoly -----------
2743 if ( (find(charstr(basering),string(char(basering)))!=1) &&
2744      (charstr(basering)<>"real")&&
2745      (charstr(basering)<>"complex") ) {
2746   string strmip=string(minpoly);
2747   string strf=string(f);
2748   execute("ring tempr=("+string(char(basering))+","+parstr(basering)+"),("
2749           +varstr(basering)+"),dp;");
2750   execute("minpoly="+strmip+";");
2751   execute("poly f="+strf+";");
2752   field_ext=1;
2753   def L=pre_HN(f,essential);
2754   if (size(L)==0) { return(list()); }
2755   def HNEring=L[1];
2756   setring HNEring;
2757   if ((typeof(hne[1])=="ideal")) { return(list()); }
2758   if ((voice==2) && (printlevel > -1)) {
2759     "// Attention: The parameter",par(1),"may have changed its meaning!";
2760     "// It needs no longer be a generator of the cyclic group of unities!";
2761   }
2762   dbprint(printlevel-voice+2,
2763     "// result: "+string(size(hne))+" branch(es) successfully computed.");
2764 }
2765 else {
2766   def L=pre_HN(f,essential);
2767   if (size(L)==0) { return(list()); }
2768   if (L[2]==1) { field_ext=1; }
2769   intvec hne_conj=L[3];
2770   def HNEring=L[1];
2771   setring HNEring;
2772   if ((typeof(hne[1])=="ideal")) { return(list()); }
2773   dbprint(printlevel-voice+2,
2774      "// result: "+string(size(hne))+" branch(es) successfully computed.");
2775 }
2776
2777 // ----- Lossen 10/02 : the branches have to be resorted to be able to
2778 // -----                display the multsequence in a nice way
2779 if (size(hne)>2)
2780 {
2781   int i,j,k,m;
2782   list dummy;
2783   int nbsave;
2784   int no_br = size(hne);
2785   intmat nbhd[no_br][no_br];
2786   for (i=1;i<no_br;i++)
2787   {
2788     for (j=i+1;j<=no_br;j++)
2789     {
2790       nbhd[i,j]=separateHNE(hne[i],hne[j]);
2791       k=i+1;
2792       while ( (nbhd[i,k] >= nbhd[i,j]) and (k<j) )
2793       {
2794         k++;
2795       }
2796       if (k<j)  // branches have to be resorted
2797       {
2798         dummy=hne[j];
2799         nbsave=nbhd[i,j];
2800         for (m=k; m<j; m++)
2801         {
2802           hne[m+1]=hne[m];
2803           nbhd[i,m+1]=nbhd[i,m];
2804         }
2805         hne[k]=dummy;
2806         nbhd[i,k]=nbsave;
2807       }
2808     }
2809   }
2810 }
2811 // -----
2812
2813 if (field_ext==1) {
2814   dbprint(printlevel-voice+3,"
2815// 'hnexpansion' created a list of one ring.
2816// To see the ring and the data stored in the ring, type (if you assigned
2817// the name L to the list):
2818     show(L);
2819// To display the computed HN expansion, type
2820     def HNring = L[1]; setring HNring;  displayHNE(hne); ");
2821   if (essential==1) {
2822     dbprint(printlevel-voice+3,""+
2823"// As second entry of the returned list L, you obtain an integer vector,
2824// indicating the number of conjugates for each of the computed branches.");
2825     return(list(HNEring,hne_conj));
2826   }
2827   return(list(HNEring));
2828 }
2829 else { // no change of basering necessary --> map data to original ring
2830   setring altring;
2831   if ((npars(altring)==1) and (minpoly!=0)) {
2832     ring HNhelpring=char(altring),(a,x,y),ls;
2833     list hne=imap(HNEring,hne);
2834     setring altring;
2835     map mmm=HNhelpring,par(1),var(1),var(2);
2836     list hne=mmm(hne);
2837     kill mmm,HNhelpring;
2838   }
2839   else {
2840     list hne=fetch(HNEring,hne);
2841   }
2842   kill HNEring;
2843   if (essential==1) {
2844     dbprint(printlevel-voice+3,""+
2845"// No change of ring necessary, return value is a list:
2846//   first entry  =  list :  HN expansion of essential branches.
2847//   second entry =  intvec: numbers of conjugated branches ");
2848     return(list(hne,hne_conj));
2849   }
2850   else {
2851     dbprint(printlevel-voice+3,""+
2852"// No change of ring necessary, return value is HN expansion.");
2853     return(hne);
2854   }
2855 }
2856}
2857example
2858{
2859  "EXAMPLE:"; echo = 2;
2860  ring r=0,(x,y),dp;
2861  // First, an example which requires no field extension:
2862  list Hne=hnexpansion(x4-y6);
2863  size(Hne);           // number of branches
2864  displayHNE(Hne);     // HN expansion of branches
2865  param(Hne[1]);       // parametrization of 1st branch
2866  param(Hne[2]);       // parametrization of 2nd branch
2867
2868  // An example which requires a field extension:
2869  list L=hnexpansion((x4-y6)*(y2+x4));
2870  def R=L[1]; setring R; displayHNE(hne);
2871  basering;
2872  setring r; kill R;
2873
2874  // Computing only one representative per conjugacy class:
2875  L=hnexpansion((x4-y6)*(y2+x4),"ess");
2876  def R=L[1]; setring R; displayHNE(hne);
2877  L[2];     // number of branches in respective conjugacy classes
2878}
2879
2880///////////////////////////////////////////////////////////////////////////////
2881
2882static proc pre_HN (poly f, int essential)
2883"NOTE: This procedure is only for internal use, it is called via
2884       hnexpansion
2885RETURN: list:  first entry = HNEring  (containing list hne, poly f)
2886               second entry = 0  if no change of base ring necessary
2887                              1  if change of base ring necessary
2888               third entry = numbers of conjugates ( if essential = 1 )
2889        if some error has occured, the empty list is returned
2890"
2891{
2892 def altring = basering;
2893 int p = char(basering);
2894 int field_ext;
2895 intvec hne_conj;
2896
2897 //-------------------- Tests auf Zulaessigkeit von basering ------------------
2898 if (charstr(basering)=="real") {
2899   " Singular cannot factorize over 'real' as ground field";
2900   return(list());
2901 }
2902 if (size(maxideal(1))<2) {
2903   " A univariate polynomial ring makes no sense !";
2904   return(list());
2905 }
2906 if (size(maxideal(1))>2) {
2907   dbprint(printlevel-voice+2,
2908   " Warning: all variables except the first two will be ignored!");
2909 }
2910 if (find(charstr(basering),string(char(basering)))!=1) {
2911   " ring of type (p^k,a) not implemented";
2912 //----------------------------------------------------------------------------
2913 // weder primitives Element noch factorize noch map "char p^k" -> "char -p"
2914 // [(p^k,a)->(p,a) ground field] noch fetch
2915 //----------------------------------------------------------------------------
2916   return(list());
2917 }
2918 //----------------- Definition eines neuen Ringes: HNEring -------------------
2919 string namex=varstr(1); string namey=varstr(2);
2920 if (string(char(altring))==charstr(altring)) { // kein Parameter, nicht 'real'
2921   ring HNEring = char(altring),(x,y),ls;
2922   map m=altring,x,y;
2923   poly f=m(f);
2924   export f;
2925   kill m;
2926 }
2927 else {
2928   string mipl=string(minpoly);
2929   if (mipl=="0") {
2930     "// ** WARNING: Algebraic extension of given ground field not possible!";
2931     "// ** We try to develop this polynomial, but if the need for a field";
2932     "// ** extension occurs during the calculation, we cannot proceed with";
2933     "// ** the corresponding branches.";
2934     execute("ring HNEring=("+charstr(basering)+"),(x,y),ls;");
2935   }
2936   else {
2937    string pa=parstr(altring);
2938    ring HNhelpring=p,`pa`,dp;
2939    execute("poly mipo="+mipl+";");  // Minimalpolynom in Polynom umgewandelt
2940    ring HNEring=(p,a),(x,y),ls;
2941    map getminpol=HNhelpring,a;
2942    mipl=string(getminpol(mipo));    // String umgewandelt mit 'a' als Param.
2943    execute("minpoly="+mipl+";");    // "minpoly=poly is not supported"
2944    kill HNhelpring, getminpol;
2945   }
2946   if (nvars(altring)==2) {
2947     poly f=fetch(altring,f);
2948     export f;
2949   }
2950   else {
2951     if (defined(pa)) { // Parameter hatte vorher anderen Namen als 'a'
2952       ring HNhelpring=p,(`pa`,x,y),ls;
2953       poly f=imap(altring,f);
2954       setring HNEring;
2955       map m=HNhelpring,a,x,y;
2956       poly f=m(f);
2957       kill HNhelpring;
2958     }
2959     else {
2960       map m=altring,x,y;
2961       poly f=m(f);
2962     }
2963     export f;
2964     kill m;
2965   }
2966 }
2967
2968 if (defined(HNDebugOn))
2969 {"received polynomial: ",f,", with x =",namex,", y =",namey;}
2970
2971 //----------------------- Variablendefinitionen ------------------------------
2972 int Abbruch,i,NullHNEx,NullHNEy;
2973 string str;
2974 list Newton,hne;
2975
2976 // --- changed for SINGULAR 3: ---
2977 hne=ideal(0);
2978 export hne;
2979
2980 //====================== Tests auf Zulaessigkeit des Polynoms ================
2981
2982 //-------------------------- Test, ob Einheit oder Null ----------------------
2983 if (subst(subst(f,x,0),y,0)!=0) {
2984   dbprint(printlevel+1,
2985           "The given polynomial is a unit in the power series ring!");
2986   setring altring; kill HNEring;
2987   return(list());                   // there are no HNEs
2988 }
2989 if (f==0) {
2990   dbprint(printlevel+1,"The given polynomial is zero!");
2991   setring altring; kill HNEring;
2992   return(list());                   // there are no HNEs
2993 }
2994
2995 //-----------------------  Test auf Quadratfreiheit --------------------------
2996
2997 if ((p==0) and (size(charstr(basering))==1)) {
2998
2999 //-------- Fall basering==0,... : Wechsel in Ring mit char >0 ----------------
3000 // weil squarefree eine Standardbasis berechnen muss (verwendet Syzygien)
3001 // -- wenn f in diesem Ring quadratfrei ist, dann erst recht im Ring HNEring
3002 //----------------------------------------------------------------------------
3003  int testerg=(polytest(f)==0);
3004  ring zweitring = 32003,(x,y),dp;
3005
3006  map polyhinueber=HNEring,x,y;         // fetch geht nicht
3007  poly f=polyhinueber(f);
3008  poly test_sqr=squarefree(f);
3009  if (test_sqr != f) {
3010   if (printlevel>0) {
3011     "Most probably the given polynomial is not squarefree. But the test was";
3012     "made in characteristic 32003 and not 0 to improve speed. You can";
3013     "(r) redo the test in char 0 (but this may take some time)";
3014     "(c) continue the development, if you're sure that the polynomial IS",
3015     "squarefree";
3016     if (testerg==1) {
3017       "(s) continue the development with a squarefree factor (*)";}
3018     "(q) or just quit the algorithm (default action)";
3019     "";"Please enter the letter of your choice:";
3020     str=read("")[1];             // reads one character
3021   }
3022   else { str="r"; }              // printlevel <= 0: non-interactive behaviour
3023   setring HNEring;
3024   map polyhinueber=zweitring,x,y;
3025   if (str=="r") {
3026     poly test_sqr=squarefree(f);
3027     if (test_sqr != f) {
3028      if (printlevel>0) { "The given polynomial is in fact not squarefree."; }
3029      else              { "The given polynomial is not squarefree!"; }
3030      "I'll continue with the radical.";
3031      f=test_sqr;
3032     }
3033     else {
3034      dbprint(printlevel,
3035       "everything is ok -- the polynomial is squarefree in characteristic 0");
3036     }
3037   }
3038   else {
3039     if ((str=="s") and (testerg==1)) {
3040       "(*)attention: it could be that the factor is only one in char 32003!";
3041       f=polyhinueber(test_sqr);
3042     }
3043     else {
3044       if (str<>"c") {
3045         setring altring;
3046         kill HNEring;kill zweitring;
3047         return(list());}
3048       else { "if the algorithm doesn't terminate, you were wrong...";}
3049   }}
3050   kill zweitring;
3051   if (defined(HNDebugOn)) {"I continue with the polynomial",f; }
3052  }
3053  else {
3054    setring HNEring;
3055    kill zweitring;
3056  }
3057 }
3058 //------------------ Fall Char > 0 oder Ring hat Parameter -------------------
3059 else {
3060  poly test_sqr=squarefree(f);
3061  if (test_sqr != f) {
3062   if (printlevel>0) {
3063    if (test_sqr == 1) {
3064     "The given polynomial is of the form g^"+string(p)+",";
3065     "therefore not squarefree.  You can:";
3066     " (q) quit the algorithm (recommended) or";
3067     " (f) continue with the full radical (using a factorization of the";
3068     "     pure power part; this could take much time)";
3069     "";"Please enter the letter of your choice:";
3070     str=read("")[1];
3071     if (str<>"f") { str="q"; }
3072    }
3073    else {
3074     "The given polynomial is not squarefree.";
3075     if (p != 0)
3076      {
3077       " You can:";
3078       " (c) continue with a squarefree divisor (but factors of the form g^"
3079       +string(p);
3080       "     are lost; this is recommended, takes no extra time)";
3081       " (f) continue with the full radical (using a factorization of the";
3082       "     pure power part; this could take some time)";
3083       " (q) quit the algorithm";
3084       "";"Please enter the letter of your choice:";
3085       str=read("")[1];
3086       if ((str<>"f") && (str<>"q")) { str="c"; }
3087      }
3088     else { "I'll continue with the radical."; str="c"; }
3089    }                                // endelse (test_sqr!=1)
3090   }
3091   else {
3092     "//** Error: The given polynomial is not squarefree!";
3093     "//** Since the global variable `printlevel' has the value",printlevel,
3094       "we stop here.";
3095     "//   Either call me again with a squarefree polynomial f or assign";
3096     "            printlevel=1;";
3097     "//   before calling me with a non-squarefree f.";
3098     "//   If printlevel > 0, I present some possibilities how to proceed.";
3099     str="q";
3100   }
3101   if (str=="q") {
3102    setring altring;kill HNEring;
3103    return(list());
3104   }
3105   if (str=="c") { f=test_sqr; }
3106   if (str=="f") { f=allsquarefree(f,test_sqr); }
3107  }
3108  if (defined(HNDebugOn)) {"I continue with the polynomial",f; }
3109
3110 }
3111 //====================== Ende Test auf Quadratfreiheit =======================
3112 if (subst(subst(f,x,0),y,0)!=0) {
3113   "The polynomial is a unit in the power series ring. No HNE computed.";
3114   setring altring;kill HNEring;
3115   return(list());
3116 }
3117 //---------------------- Test, ob f teilbar durch x oder y -------------------
3118 if (subst(f,y,0)==0) {
3119   f=f/y; NullHNEy=1; }             // y=0 is a solution
3120 if (subst(f,x,0)==0) {
3121   f=f/x; NullHNEx=1; }             // x=0 is a solution
3122
3123 Newton=newtonpoly(f,1);
3124 i=1; Abbruch=0;
3125 //----------------------------------------------------------------------------
3126 // finde Eckpkt. des Newtonpolys, der den Teil abgrenzt, fuer den x transvers:
3127 // Annahme: Newton ist sortiert, s.d. Newton[1]=Punkt auf der y-Achse,
3128 // Newton[letzt]=Punkt auf der x-Achse
3129 //----------------------------------------------------------------------------
3130 while ((i<size(Newton)) and (Abbruch==0)) {
3131  if ((Newton[i+1][1]-Newton[i][1])>=(Newton[i][2]-Newton[i+1][2]))
3132   {Abbruch=1;}
3133  else {i=i+1;}
3134 }
3135 int grenze1=Newton[i][2];
3136 int grenze2=Newton[i][1];
3137 //----------------------------------------------------------------------------
3138 // Stelle Ring bereit zur Uebertragung der Daten im Fall einer Koerperer-
3139 // weiterung. Definiere Objekte, die spaeter uebertragen werden.
3140 // Binde die Listen (azeilen,...) an den Ring (um sie nicht zu ueberschreiben
3141 // bei Def. in einem anderen Ring).
3142 //----------------------------------------------------------------------------
3143 ring HNE_noparam = char(altring),(a,x,y),ls;
3144 poly f;
3145 list azeilen=ideal(0);
3146 list HNEs=ideal(0);
3147 list aneu=ideal(0);
3148 list faktoren=ideal(0);
3149
3150 ideal deltais;
3151 poly delt;
3152
3153 //----- hier steht die Anzahl bisher benoetigter Ringerweiterungen drin: -----
3154 int EXTHNEnumber=0;
3155
3156 list EXTHNEring;
3157 list HNE_RingDATA;
3158 int number_of_letztring;
3159 setring HNEring;
3160 number_of_letztring=0;
3161
3162 // ================= Die eigentliche Berechnung der HNE: =====================
3163
3164 // ------- Berechne HNE von allen Zweigen, fuer die x transversal ist: -------
3165 if (defined(HNDebugOn))
3166   {"1st step: Treat Newton polygon until height",grenze1;}
3167 if (grenze1>0) {
3168  if (EXTHNEnumber>0){ EXTHNEring = EXTHNEring(1..EXTHNEnumber); }
3169  HNE_RingDATA = list(HNEring, HNE_noparam, EXTHNEnumber, EXTHNEring,
3170                      number_of_letztring);
3171
3172  list hilflist=HN(HNE_RingDATA,f,grenze1,1,essential,0,hne_conj,1);
3173  kill HNEring, HNE_noparam;
3174  if (EXTHNEnumber>0) { kill EXTHNEring(1..EXTHNEnumber);}
3175  def HNEring = hilflist[1][1];
3176  def HNE_noparam = hilflist[1][2];
3177  EXTHNEnumber = hilflist[1][3];
3178  for (i=1; i<=EXTHNEnumber; i++) { def EXTHNEring(i)=hilflist[1][4][i]; }
3179  if (hilflist[2]==0) { setring HNEring; number_of_letztring=0; }
3180  else                { setring EXTHNEring(hilflist[2]);}
3181  if (hilflist[3]==1){field_ext=1;}
3182  hne_conj=hilflist[5];
3183
3184  if (number_of_letztring != hilflist[2])
3185  {  // Ringwechsel in Prozedur HN
3186     map hole=HNE_noparam,transfproc,x,y;
3187     setring HNE_noparam;
3188     if (not(defined(f))) {poly f;}
3189     f=imap(HNEring,f);
3190     setring EXTHNEring(EXTHNEnumber);
3191     if (not(defined(f))) {poly f; f=hole(f); export f;}
3192     else                 {f=hole(f);}
3193  }
3194  number_of_letztring = hilflist[2];
3195  kill hilflist;
3196 }
3197
3198 if (NullHNEy==1) {
3199  if ((typeof(hne[1])=="ideal")) { hne=list(); }
3200  hne=hne+list(list(matrix(ideal(0,x)),intvec(1),int(0),poly(0)));
3201  if (hne_conj==0) { hne_conj=1; }
3202  else { hne_conj = hne_conj, 1; }
3203 }
3204 // --------------- Berechne HNE von allen verbliebenen Zweigen: --------------
3205 if (defined(HNDebugOn))
3206    {"2nd step: Treat Newton polygon until height",grenze2;}
3207 if (grenze2>0) {
3208
3209  if (EXTHNEnumber>0){ EXTHNEring = EXTHNEring(1..EXTHNEnumber); }
3210
3211  if (essential==1) { number_of_letztring=0; }
3212  if (number_of_letztring==0) { setring HNEring; }
3213  else                        { setring EXTHNEring(number_of_letztring); }
3214  map xytausch=basering,y,x;
3215
3216  HNE_RingDATA = list(HNEring, HNE_noparam, EXTHNEnumber, EXTHNEring,
3217                      number_of_letztring);
3218  list hilflist=HN(HNE_RingDATA,xytausch(f),grenze2,1,essential,1,hne_conj,1);
3219  kill HNEring, HNE_noparam;
3220  if (EXTHNEnumber>0){ kill EXTHNEring(1..EXTHNEnumber); }
3221  def HNEring = hilflist[1][1];
3222  def HNE_noparam = hilflist[1][2];
3223  EXTHNEnumber = hilflist[1][3];
3224  for (i=1; i<=EXTHNEnumber; i++) { def EXTHNEring(i)=hilflist[1][4][i]; }
3225  if (hilflist[2]==0) { setring HNEring; number_of_letztring=0; }
3226  else                { setring EXTHNEring(hilflist[2]);
3227                        number_of_letztring=hilflist[2]; }
3228  if (hilflist[3]==1){field_ext=1;}
3229  hne_conj=hilflist[5];
3230  kill hilflist;
3231 }
3232 if (NullHNEx==1) {
3233  if ((typeof(hne[1])=="ideal")) { hne=list(); }
3234  hne=hne+list(list(matrix(ideal(0,x)),intvec(1),int(1),poly(0)));
3235  if (hne_conj==0) { hne_conj=1; }
3236  else { hne_conj = hne_conj, 1; }
3237 }
3238
3239
3240 // --- aufraeumen ---
3241 if (defined(HNEakut)){
3242   kill HNEakut,faktoren,deltais,transformiert,teiler,leitf;
3243 }
3244 if (defined(hilflist)) {kill hilflist;}
3245 if (defined(erg)) {kill erg;}
3246 if (defined(delt)) {kill delt;}
3247 if (defined(azeilen)) { kill azeilen;}
3248 if (defined(aneu)) { kill aneu;}
3249 if (defined(transfproc)) { kill transfproc;}
3250 if (defined(transf)) { kill transf;}
3251 if (not(defined(f))) { poly f = imap(HNEring,f); export f; }
3252
3253 return(list(basering,field_ext,hne_conj));
3254}
3255
3256//////////////////////////////////////////////////////////////////////////////
3257proc essdevelop (poly f)
3258"USAGE:   essdevelop(f); f poly
3259NOTE:     command is obsolete, use hnexpansion(f,\"ess\") instead.
3260SEE ALSO: hnexpansion, develop, extdevelop, param
3261"
3262{
3263 printlevel=printlevel+1;
3264 list Ergebnis=hnexpansion(f,1);
3265 printlevel=printlevel-1;
3266 return(Ergebnis);
3267}
3268
3269///////////////////////////////////////////////////////////////////////////////
3270static proc HN (list HNE_RingDATA,poly fneu,int grenze,Aufruf_Ebene,
3271                     essential,getauscht,intvec hne_conj,int conj_factor)
3272"NOTE: This procedure is only for internal use, it is called via pre_HN
3273RETURN: list: first entry = list of HNErings,
3274              second entry = number of new base ring (0 for HNEring,
3275                                                      -1 for HNE_noparam,
3276                                                      i for EXTHNEring(i))
3277              third entry = 0 if no field extension necessary
3278                            1 if field extension necessary
3279              forth entry = HNEs (only if no change of basering)
3280"
3281{
3282 //---------- Variablendefinitionen fuer den unverzweigten Teil: --------------
3283 if (defined(HNDebugOn)) {"procedure HN",Aufruf_Ebene;}
3284 int Abbruch,ende,i,j,k,e,M,N,Q,R,zeiger,zeile,zeilevorher,dd,ii;
3285 intvec hqs;
3286 int field_ext;
3287 int ring_changed, hneshift;
3288 intvec conjugates,conj2,conj1;
3289
3290 list EXTHNEring;
3291 def HNEring = HNE_RingDATA[1];
3292 def HNE_noparam = HNE_RingDATA[2];
3293 int EXTHNEnumber = HNE_RingDATA[3];
3294 for (i=1; i<=EXTHNEnumber; i++) { def EXTHNEring(i)=HNE_RingDATA[4][i]; }
3295 int number_of_letztring = HNE_RingDATA[5];
3296 if (defined(basering))
3297 {
3298   if (number_of_letztring==0) { kill HNEring; def HNEring=basering; }
3299   else                 { kill EXTHNEring(number_of_letztring);
3300                          def EXTHNEring(number_of_letztring)=basering; }
3301 }
3302 else
3303 {
3304   if ( number_of_letztring==0) { setring HNEring; }
3305   else                         { setring EXTHNEring(number_of_letztring); }
3306 }
3307 if (not(defined(hne))) {list hne;}
3308 poly fvorher;
3309 list erg=ideal(0); list HNEs=ideal(0); // um die Listen an den Ring zu binden
3310
3311 //-------------------- Bedeutung von Abbruch: --------------------------------
3312 //------- 0:keine Verzweigung | 1:Verzweigung,nicht fertig | 2:fertig --------
3313 //
3314 // Struktur von HNEs : Liste von Listen L (fuer jeden Zweig) der Form
3315 // L[1]=intvec (hqs), L[2],L[3],... ideal (die Zeilen (0,1,...) der HNE)
3316 // L[letztes]=poly (transformiertes f)
3317 //----------------------------------------------------------------------------
3318 list Newton;
3319 number delt;
3320 int p = char(basering);                // Ringcharakteristik
3321 list azeilen=ideal(0);
3322
3323 ideal hilfid; intvec hilfvec;
3324
3325 // ======================= der unverzweigte Teil: ============================
3326 while (Abbruch==0) {
3327  Newton=newtonpoly(fneu,1);
3328  zeiger=find_in_list(Newton,grenze);
3329  if (Newton[zeiger][2] != grenze)
3330    {"Didn't find an edge in the Newton polygon!";}
3331  if (zeiger==size(Newton)-1) {
3332    if (defined(HNDebugOn)) {"only one relevant side in Newton polygon";}
3333    M=Newton[zeiger+1][1]-Newton[zeiger][1];
3334    N=Newton[zeiger][2]-Newton[zeiger+1][2];
3335    R = M%N;
3336    Q = M / N;
3337
3338    //-------- 1. Versuch: ist der quasihomogene Leitterm reine Potenz ? ------
3339    //              (dann geht alles wie im irreduziblen Fall)
3340    //-------------------------------------------------------------------------
3341    e = gcd(M,N);
3342    delt=factorfirst(redleit(fneu,Newton[zeiger],Newton[zeiger+1])
3343                      /x^Newton[zeiger][1],M,N);
3344    if (delt==0) {
3345      if (defined(HNDebugOn)) {" The given polynomial is reducible !";}
3346      Abbruch=1;
3347    }
3348    if (Abbruch==0) {
3349      //----------- fneu,zeile retten fuer den Spezialfall (###): -------------
3350      fvorher=fneu;zeilevorher=zeile;
3351      if (R==0) {
3352        //-------- transformiere fneu mit T1, wenn kein Abbruch nachher: ------
3353        if (N>1) { fneu = T1_Transform(fneu,delt,M/ e); }
3354        else     { ende=1; }
3355        if (defined(HNDebugOn)) {"a("+string(zeile)+","+string(Q)+") =",delt;}
3356        azeilen[zeile+1][Q]=delt;
3357      }
3358      else {
3359        //------------- R > 0 : transformiere fneu mit T2 ---------------------
3360        erg=T2_Transform(fneu,delt,M,N,referencepoly(Newton));
3361        fneu=erg[1];delt=erg[2];
3362        //----- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: ----
3363        while (R!=0) {
3364         if (defined(HNDebugOn)) { "h("+string(zeile)+") =",Q; }
3365         hqs[zeile+1]=Q;         // denn zeile beginnt mit dem Wert 0
3366         //--------------- markiere das Zeilenende der HNE: -------------------
3367         azeilen[zeile+1][Q+1]=x;
3368         zeile=zeile+1;
3369         //-------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ----
3370         azeilen[zeile+1]=ideal(0);
3371         M=N; N=R; R=M%N; Q=M / N;
3372        }
3373        if (defined(HNDebugOn)) {"a("+string(zeile)+","+string(Q)+") =",delt;}
3374        azeilen[zeile+1][Q]=delt;
3375      }
3376      if (defined(HNDebugOn)) {"transformed polynomial: ",fneu;}
3377      grenze=e;
3378      //----------------------- teste Abbruchbedingungen: ---------------------
3379      if (subst(fneu,y,0)==0) {              // <==> y|fneu
3380        dbprint(printlevel-voice+3,"finite HNE of one branch found");
3381           // voice abzufragen macht bei rekursiven procs keinen Sinn
3382        azeilen[zeile+1][Q+1]=x;
3383        //----- Q wird nur in hqs eingetragen, wenn der Spezialfall nicht
3384        //      eintritt (siehe unten) -----
3385        Abbruch=2;
3386        if (grenze>1) {
3387         if (jet(fneu,1,intvec(0,1))==0) {
3388           //- jet(...)=alle Monome von fneu, die nicht durch y2 teilbar sind -
3389           "THE TEST FOR SQUAREFREENESS WAS BAD!!";
3390           " The polynomial was NOT squarefree!!!";}
3391         else {
3392           //----------------------- Spezialfall (###): -----------------------
3393           // Wir haben das Problem, dass die HNE eines Zweiges hier abbricht,
3394           // aber ein anderer Zweig bis hierher genau die gleiche HNE hat, die
3395           // noch weiter geht
3396           // Loesung: mache Transform. rueckgaengig und behandle fneu im
3397           // Verzweigungsteil
3398           //------------------------------------------------------------------
3399          Abbruch=1;
3400          fneu=fvorher;zeile=zeilevorher;grenze=Newton[zeiger][2];
3401        }}
3402        else {fneu=0;}     // fneu nicht mehr gebraucht - spare Speicher
3403        if (Abbruch==2) { hqs[zeile+1]=Q; }
3404      }                 // Spezialfall nicht eingetreten
3405      else {
3406        if (ende==1) {
3407          dbprint(printlevel-voice+2,"HNE of one branch found");
3408          Abbruch=2; hqs[zeile+1]=-Q-1;}
3409      }
3410    }                   // end(if Abbruch==0)
3411  }                     // end(if zeiger...)
3412  else { Abbruch=1;}
3413 }                      // end(while Abbruch==0)
3414
3415 // ===================== der Teil bei Verzweigung: ===========================
3416 if (Abbruch==1) {
3417  //---------- Variablendefinitionen fuer den verzweigten Teil: ---------------
3418  poly leitf,teiler,transformiert;
3419  list aneu=ideal(0);
3420  list faktoren;
3421  ideal deltais;
3422  list HNEakut=ideal(0);
3423  intvec eis;
3424  int zaehler,hnezaehler,zl,zl1,M1,N1,R1,Q1,needext;
3425  int numberofRingchanges,lastRingnumber,ringischanged,flag;
3426  string letztringname;
3427
3428  zeiger=find_in_list(Newton,grenze);
3429  if (defined(HNDebugOn)) {
3430    "Branching part reached---Newton polygon :",Newton;
3431    "relevant part until height",grenze,", from",Newton[zeiger],"on";
3432  }
3433  azeilen=list(hqs)+azeilen; // hat jetzt Struktur von HNEs: hqs in der 1.Zeile
3434
3435  //======= Schleife fuer jede zu betrachtende Seite des Newtonpolygons: ======
3436  for(i=zeiger; i<size(Newton); i++) {
3437   if ((essential==1) and (EXTHNEnumber>number_of_letztring)) {
3438     // ----- setze ring zurueck fuer neue Kante  -----
3439     // ---- (damit konjugierte Zweige erkennbar) -----
3440     hneshift=hneshift+hnezaehler;
3441     hnezaehler=0;
3442     ring_changed=0;
3443     def SaveRing = EXTHNEring(EXTHNEnumber);
3444     setring SaveRing;
3445     if (not(defined(HNEs))) { // HN wurde zum 2.Mal von pre_HN aufgerufen
3446       list HNEs=ideal(0);
3447     }
3448     for (k=number_of_letztring+1; k<=EXTHNEnumber; k++) { kill EXTHNEring(k);}
3449     EXTHNEnumber=number_of_letztring;
3450     if (EXTHNEnumber==0) { setring HNEring; }
3451     else                 { setring EXTHNEring(EXTHNEnumber); }
3452     if (not(defined(HNEs))) { list HNEs; }
3453     HNEs=ideal(0);
3454     deltais=0;
3455     delt=0;
3456     if (defined(zerlege)) { kill zerlege; }
3457   }
3458
3459   if (defined(HNDebugOn)) { "we consider side",Newton[i],Newton[i+1]; }
3460   M=Newton[i+1][1]-Newton[i][1];
3461   N=Newton[i][2]-Newton[i+1][2];
3462   R = M%N;
3463   Q = M / N;
3464   e=gcd(M,N);
3465   needext=1;
3466   letztringname=nameof(basering);
3467   lastRingnumber=EXTHNEnumber;
3468   faktoren=list(ideal(charPoly(redleit(fneu,Newton[i],Newton[i+1])
3469                       /(x^Newton[i][1]*y^Newton[i+1][2]),M,N)  ),
3470                 intvec(1));                  // = (zu faktoriserendes Poly, 1)
3471   conjugates=conj_factor;
3472
3473   //-- wechsle so lange in Ringerweiterungen, bis Leitform vollstaendig
3474   //   in Linearfaktoren zerfaellt -----
3475   for (numberofRingchanges=1; needext==1; numberofRingchanges++) {
3476    leitf=redleit(fneu,Newton[i],Newton[i+1])/
3477                     (x^Newton[i][1]*y^Newton[i+1][2]);
3478    delt=factorfirst(leitf,M,N);
3479    needext=0;
3480    if (delt==0) {
3481     //---------- Sonderbehandlung: faktorisiere einige Polynome ueber Q(a): --
3482     if ((charstr(basering)=="0,a") and (essential==0)) {
3483        // ====================================================
3484        // neu CL:  06.10.05
3485        poly CHPOLY=charPoly(leitf,M,N);
3486        poly tstpoly;
3487        if (defined(faktoren)!=0) {
3488          // Test, damit kein Fehler eingebaut (vermutlich nicht notwendig)
3489          tstpoly = faktoren[1][1]^faktoren[2][1];
3490          for (k=2; k<=size(faktoren[1]); k++) {
3491            tstpoly = tstpoly * faktoren[1][k]^faktoren[2][k];
3492          }
3493          tstpoly = CHPOLY-tstpoly*(CHPOLY/tstpoly);
3494          kill CHPOLY;
3495        }
3496        if ((numberofRingchanges>1) and (defined(faktoren)!=0) and (tstpoly==0)) {
3497          def L_help=factorlist(faktoren,conjugates);
3498          faktoren=L_help[1];
3499          conjugates=L_help[2];
3500          kill L_help;
3501        }
3502        else {
3503          faktoren=factorize(charPoly(leitf,M,N));
3504          conjugates=conj_factor;
3505          for (k=2;k<=size(faktoren[2]);k++) {conjugates=conjugates,conj_factor;}
3506        }
3507        kill tstpoly;
3508        // Ende neu (vorher nur else Fall)
3509        // ====================================================
3510     }
3511     else {
3512       //------------------ faktorisiere das charakt. Polynom: ----------------
3513       if ((numberofRingchanges==1) or (essential==0)) {
3514         def L_help=factorlist(faktoren,conjugates);
3515         faktoren=L_help[1];
3516         conjugates=L_help[2];
3517         kill L_help;
3518       }
3519       else {     // eliminiere alle konjugierten Nullstellen bis auf eine:
3520         ideal hilf_id;
3521         for (zaehler=1; zaehler<=size(faktoren[1]); zaehler++) {
3522           hilf_id=factorize(faktoren[1][zaehler],1);
3523           if (size(hilf_id)>1) {
3524             poly fac=hilf_id[1];
3525             dd=deg(fac);
3526             // Zur Sicherheit:
3527             if (deg(fac)==0) { fac=hilf_id[2]; }
3528             for (k=2;k<=size(hilf_id);k++) {
3529               dd=dd+deg(hilf_id[k]);
3530               if (deg(hilf_id[k])<deg(fac)) { fac=hilf_id[k]; }
3531             }
3532             faktoren[1][zaehler]=fac;
3533             kill fac;
3534             if (conjugates[zaehler]==conj_factor) {
3535               // number of conjugates not yet determined for this factor
3536               conjugates[zaehler]=conjugates[zaehler]*dd;
3537             }
3538           }
3539           else {
3540             faktoren[1][zaehler]=hilf_id[1];
3541           }
3542         }
3543       }
3544     }
3545
3546     zaehler=1; eis=0;
3547     for (j=1; j<=size(faktoren[2]); j++) {
3548      teiler=faktoren[1][j];
3549      if (teiler/y != 0) {         // sonst war's eine Konstante --> wegwerfen!
3550        if (defined(HNDebugOn)) {"factor of leading form found:",teiler;}
3551        if (teiler/y2 == 0) {      // --> Faktor hat die Form cy+d
3552          deltais[zaehler]=-subst(teiler,y,0)/koeff(teiler,0,1); //=-d/c
3553          eis[zaehler]=faktoren[2][j];
3554          conj2[zaehler]=conjugates[j];
3555          zaehler++;
3556        }
3557        else {
3558          dbprint(printlevel-voice+2,
3559             " Change of basering (field extension) necessary!");
3560          if (defined(HNDebugOn)) { teiler,"is not yet properly factorized!"; }
3561          if (needext==0) { poly zerlege=teiler; }
3562          needext=1;
3563          field_ext=1;
3564        }
3565      }
3566     }  // end(for j)
3567    }
3568    else { deltais=ideal(delt); eis=e; conj2=conj_factor; }
3569    if (defined(HNDebugOn)) {"roots of char. poly:";deltais;
3570                             "with multiplicities:",eis;}
3571    if (needext==1) {
3572      //--------------------- fuehre den Ringwechsel aus: ---------------------
3573      ringischanged=1;
3574      if ((size(parstr(basering))>0) && string(minpoly)=="0") {
3575        " ** We've had bad luck! The HNE cannot be calculated completely!";
3576        // HNE in transzendenter Erweiterung fehlgeschlagen
3577        kill zerlege;
3578        ringischanged=0; break;    // weiter mit gefundenen Faktoren
3579      }
3580      if (parstr(basering)=="") {
3581        EXTHNEnumber++;
3582        def EXTHNEring(EXTHNEnumber) = splitring(zerlege);
3583        setring EXTHNEring(EXTHNEnumber);
3584
3585        poly transf=0;
3586        poly transfproc=0;
3587        ring_changed=1;
3588        export transfproc;
3589      }
3590      else {
3591        if (numberofRingchanges>1) {  // ein Ringwechsel hat nicht gereicht
3592         def helpring = splitring(zerlege,list(transf,transfproc,faktoren));
3593         kill EXTHNEring(EXTHNEnumber);
3594         def EXTHNEring(EXTHNEnumber)=helpring;
3595         setring EXTHNEring(EXTHNEnumber);
3596         kill helpring;
3597
3598         poly transf=erg[1];
3599         poly transfproc=erg[2];
3600         ring_changed=1;
3601         list faktoren=erg[3];
3602         export transfproc;
3603         kill erg;
3604        }
3605        else {
3606         if (ring_changed==1) { // in dieser proc geschah schon Ringwechsel
3607          EXTHNEnumber++;
3608          def EXTHNEring(EXTHNEnumber) = splitring(zerlege,list(a,transfproc));
3609          setring EXTHNEring(EXTHNEnumber);
3610          poly transf=erg[1];
3611          poly transfproc=erg[2];
3612          export transfproc;
3613          kill erg;
3614         }
3615         else { // parameter war vorher da
3616          EXTHNEnumber++;
3617          def EXTHNEring(EXTHNEnumber) = splitring(zerlege,a);
3618          setring EXTHNEring(EXTHNEnumber);
3619          poly transf=erg[1];
3620          poly transfproc=transf;
3621          ring_changed=1;
3622          export transfproc;
3623          kill erg;
3624        }}
3625      }
3626      //-----------------------------------------------------------------------
3627      // transf enthaelt jetzt den alten Parameter des Ringes, der aktiv war
3628      // vor Beginn der Schleife (evtl. also ueber mehrere Ringwechsel
3629      // weitergereicht),
3630      // transfproc enthaelt den alten Parameter des Ringes, der aktiv war zu
3631      // Beginn der Prozedur, und der an die aufrufende Prozedur zurueckgegeben
3632      // werden muss
3633      // transf ist Null, falls der alte Ring keinen Parameter hatte,
3634      // das gleiche gilt fuer transfproc
3635      //-----------------------------------------------------------------------
3636
3637      //---- Neudef. von Variablen, Uebertragung bisher errechneter Daten: ----
3638      poly leitf,teiler,transformiert;
3639      list aneu=ideal(0);
3640      ideal deltais;
3641      number delt;
3642      setring HNE_noparam;
3643      if (defined(letztring)) { kill letztring; }
3644      if (EXTHNEnumber>1) { def letztring=EXTHNEring(EXTHNEnumber-1); }
3645      else                { def letztring=HNEring; }
3646      if (not defined(fneu)) {poly fneu;}
3647      fneu=imap(letztring,fneu);
3648      if (not defined(f)) {poly f;}
3649      f=imap(letztring,f);
3650      if (not defined(hne)) {list hne;}
3651      hne=imap(letztring,hne);
3652
3653      if (not defined(faktoren)) {list faktoren; }
3654      faktoren=imap(letztring,faktoren);
3655
3656      if (not(defined(azeilen))){list azeilen,HNEs;}
3657      azeilen=imap(letztring,azeilen);
3658      HNEs=imap(letztring,HNEs);
3659
3660      setring EXTHNEring(EXTHNEnumber);
3661      if (not(defined(hole))) { map hole; }
3662      hole=HNE_noparam,transf,x,y;
3663      poly fneu=hole(fneu);
3664      if (not defined(faktoren)) {
3665        list faktoren;
3666        faktoren=hole(faktoren);
3667      }
3668      if (not(defined(f)))
3669      {
3670        poly f=hole(f);
3671        list hne=hole(hne);
3672        export f,hne;
3673      }
3674    }
3675   }    // end (while needext==1) bzw. for (numberofRingchanges)
3676
3677   if (eis==0) { i++; continue; }
3678   if (ringischanged==1) {
3679    list erg,HNEakut;            // dienen nur zum Sp. von Zwi.erg.
3680
3681    ideal hilfid;
3682    erg=ideal(0); HNEakut=erg;
3683
3684    hole=HNE_noparam,transf,x,y;
3685    setring HNE_noparam;
3686    if (not(defined(azeilen))){list azeilen,HNEs;}
3687    azeilen=imap(letztring,azeilen);
3688    HNEs=imap(letztring,HNEs);
3689
3690    setring EXTHNEring(EXTHNEnumber);
3691    list azeilen=hole(azeilen);
3692    list HNEs=hole(HNEs);
3693    kill letztring;
3694    ringischanged=0;
3695   }
3696
3697   //============ Schleife fuer jeden gefundenen Faktor der Leitform: =========
3698   for (j=1; j<=size(eis); j++) {
3699     //---- Mache Transformation T1 oder T2, trage Daten in HNEs ein,
3700     //     falls HNE abbricht: ----
3701
3702    //------------------------ Fall R==0: -------------------------------------
3703    if (R==0) {
3704      transformiert = T1_Transform(fneu,number(deltais[j]),M/ e);
3705      if (defined(HNDebugOn)) {
3706        "a("+string(zeile)+","+string(Q)+") =",deltais[j];
3707        "transformed polynomial: ",transformiert;
3708      }
3709      if (subst(transformiert,y,0)==0) {
3710       dbprint(printlevel-voice+3,"finite HNE found");
3711       hnezaehler++;
3712       //-------- trage deltais[j],x ein in letzte Zeile, fertig: -------------
3713       HNEakut=azeilen+list(poly(0));        // =HNEs[hnezaehler];
3714       hilfid=HNEakut[zeile+2]; hilfid[Q]=deltais[j]; hilfid[Q+1]=x;
3715       HNEakut[zeile+2]=hilfid;
3716       HNEakut[1][zeile+1]=Q;                // aktualisiere Vektor mit den hqs
3717       HNEs[hnezaehler]=HNEakut;
3718       conj1[hneshift+hnezaehler]=conj2[j];
3719       if (eis[j]>1) {
3720        transformiert=transformiert/y;
3721        if (subst(transformiert,y,0)==0){"THE TEST FOR SQUAREFREENESS WAS BAD!"
3722                                  +"! The polynomial was NOT squarefree!!!";}
3723        else {
3724          //--- Spezialfall (###) eingetreten: Noch weitere Zweige vorhanden --
3725          eis[j]=eis[j]-1;
3726        }
3727       }
3728      }
3729    }
3730    else {
3731      //------------------------ Fall R <> 0: ---------------------------------
3732      erg=T2_Transform(fneu,number(deltais[j]),M,N,referencepoly(Newton));
3733      transformiert=erg[1];delt=erg[2];
3734      if (defined(HNDebugOn)) {"transformed polynomial: ",transformiert;}
3735      if (subst(transformiert,y,0)==0) {
3736       dbprint(printlevel-voice+3,"finite HNE found");
3737       hnezaehler++;
3738       //---------------- trage endliche HNE in HNEs ein: ---------------------
3739       HNEakut=azeilen;           // dupliziere den gemeins. Anfang der HNE's
3740       zl=2;                      // (kommt schliesslich nach HNEs[hnezaehler])
3741       //----------------------------------------------------------------------
3742       // Werte von:  zeile: aktuelle Zeilennummer der HNE (gemeinsamer Teil)
3743       //             zl : die HNE spaltet auf; zeile+zl ist der Index fuer die
3744       //                  Zeile des aktuellen Zweigs; (zeile+zl-2) ist die
3745       //                  tatsaechl. Zeilennr. (bei 0 angefangen) der HNE
3746       //                  ([1] <- intvec(hqs), [2] <- 0. Zeile usw.)
3747       //----------------------------------------------------------------------
3748
3749       //----- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: -----
3750       M1=M;N1=N;R1=R;Q1=M1/ N1;
3751       while (R1!=0) {
3752        if (defined(HNDebugOn)) { "h("+string(zeile+zl-2)+") =",Q1; }
3753        HNEakut[1][zeile+zl-1]=Q1;
3754        HNEakut[zeile+zl][Q1+1]=x;
3755                                  // markiere das Zeilenende der HNE
3756        zl=zl+1;
3757        //----- Bereitstellung von Speicherplatz fuer eine neue Zeile: --------
3758        HNEakut[zeile+zl]=ideal(0);
3759
3760        M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1;
3761       }
3762       if (defined(HNDebugOn)) {
3763         "a("+string(zeile+zl-2)+","+string(Q1)+") =",delt;
3764       }
3765       HNEakut[zeile+zl][Q1]  =delt;
3766       HNEakut[zeile+zl][Q1+1]=x;
3767       HNEakut[1][zeile+zl-1] =Q1;     // aktualisiere Vektor mit hqs
3768       HNEakut[zeile+zl+1]=poly(0);
3769       HNEs[hnezaehler]=HNEakut;
3770       conj1[hneshift+hnezaehler]=conj2[j];
3771
3772       //-------------------- Ende der Eintragungen in HNEs -------------------
3773
3774       if (eis[j]>1) {
3775        transformiert=transformiert/y;
3776        if (subst(transformiert,y,0)==0){"THE TEST FOR SQUAREFREENESS WAS BAD!"
3777                               +" The polynomial was NOT squarefree!!!";}
3778         else {
3779          //--- Spezialfall (###) eingetreten: Noch weitere Zweige vorhanden --
3780          eis[j]=eis[j]-1;
3781       }}
3782      }                           // endif (subst()==0)
3783    }                             // endelse (R<>0)
3784
3785    //========== Falls HNE nicht abbricht: Rekursiver Aufruf von HN: ==========
3786    //------------------- Berechne HNE von transformiert ----------------------
3787    if (subst(transformiert,y,0)!=0) {
3788     lastRingnumber=EXTHNEnumber;
3789
3790     if (EXTHNEnumber>0){ EXTHNEring = EXTHNEring(1..EXTHNEnumber); }
3791     HNE_RingDATA = list( HNEring, HNE_noparam, EXTHNEnumber, EXTHNEring,
3792                          lastRingnumber);
3793     if (defined(HNerg)) {kill HNerg;}
3794     list HNerg=HN(HNE_RingDATA,transformiert,eis[j],Aufruf_Ebene+1,
3795                                essential,getauscht,hne_conj,conj2[j]);
3796     HNE_RingDATA = HNerg[1];
3797     if (conj1==0) { conj1=HNerg[5]; }
3798     else  { conj1 = conj1,HNerg[5]; }
3799
3800     if (HNerg[3]==1) { field_ext=1; }
3801     if (HNerg[2]==lastRingnumber) { // kein Ringwechsel in HN aufgetreten
3802       if (not(defined(aneu))) { list aneu; }
3803       aneu = HNerg[4];
3804     }
3805     else { // Ringwechsel aufgetreten
3806       if (defined(HNDebugOn))
3807          {" ring change in HN(",Aufruf_Ebene+1,") detected";}
3808       EXTHNEnumber = HNerg[1][3];
3809       for (ii=lastRingnumber+1; ii<=EXTHNEnumber; ii++) {
3810         def EXTHNEring(ii)=HNerg[1][4][ii];
3811       }
3812       if (HNerg[2]==0) { setring HNEring; }
3813       else             { setring EXTHNEring(HNerg[2]); }
3814       def tempRing=HNerg[4];
3815       def aneu=imap(tempRing,HNEs);
3816       kill tempRing;
3817
3818       //--- stelle lokale Variablen im neuen Ring wieder her, und rette
3819       //    gegebenenfalls ihren Inhalt: ----
3820       list erg,faktoren,HNEakut;
3821       ideal hilfid;
3822       erg=ideal(0); faktoren=erg; HNEakut=erg;
3823       poly leitf,teiler,transformiert;
3824       map hole=HNE_noparam,transfproc,x,y;
3825       setring HNE_noparam;
3826       if (lastRingnumber>0) { def letztring=EXTHNEring(lastRingnumber); }
3827       else                  { def letztring=HNEring; }
3828       if (not defined(HNEs)) {list HNEs;}
3829       HNEs=imap(letztring,HNEs);
3830       if (not defined(azeilen)) {list azeilen;}
3831       azeilen=imap(letztring,azeilen);
3832       if (not defined(deltais)) {ideal deltais;}
3833       deltais=imap(letztring,deltais);
3834       if (not defined(delt)) {poly delt;}
3835       delt=imap(letztring,delt);
3836       if (not defined(fneu)) {poly fneu;}
3837       fneu=imap(letztring,fneu);
3838       if (not defined(f)) {poly f;}
3839       f=imap(letztring,f);
3840       if (not defined(hne)) {list hne;}
3841       hne=imap(letztring,hne);
3842
3843       setring EXTHNEring(EXTHNEnumber);
3844       list HNEs=hole(HNEs);
3845       list azeilen=hole(azeilen);
3846       ideal deltais=hole(deltais);
3847       number delt=number(hole(delt));
3848       poly fneu=hole(fneu);
3849       if (not(defined(f)))
3850       {
3851         poly f=hole(f);
3852         list hne=hole(hne);
3853         export f,hne;
3854       }
3855     }
3856
3857     //========== Verknuepfe bisherige HNE mit von HN gelieferten HNEs: ======
3858     if (R==0) {
3859       HNEs,hnezaehler=constructHNEs(HNEs,hnezaehler,aneu,azeilen,zeile,
3860                       deltais,Q,j);
3861       kill aneu;
3862     }
3863     else {
3864      for (zaehler=1; zaehler<=size(aneu); zaehler++) {
3865       hnezaehler++;
3866       HNEakut=azeilen;          // dupliziere den gemeinsamen Anfang der HNE's
3867       zl=2;                     // (kommt schliesslich nach HNEs[hnezaehler])
3868       //------------ Trage Beitrag dieser Transformation T2 ein: -------------
3869       //------ Zur Bedeutung von zeile, zl: siehe Kommentar weiter oben ------
3870
3871       //----- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: -----
3872       M1=M;N1=N;R1=R;Q1=M1/ N1;
3873       while (R1!=0) {
3874        if (defined(HNDebugOn)) { "h("+string(zeile+zl-2)+") =",Q1; }
3875        HNEakut[1][zeile+zl-1]=Q1;
3876        HNEakut[zeile+zl][Q1+1]=x;    // Markierung des Zeilenendes der HNE
3877        zl=zl+1;
3878        //----- Bereitstellung von Speicherplatz fuer eine neue Zeile: --------
3879        HNEakut[zeile+zl]=ideal(0);
3880        M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1;
3881       }
3882       if (defined(HNDebugOn)) {
3883         "a("+string(zeile+zl-2)+","+string(Q1)+") =",delt;
3884       }
3885       HNEakut[zeile+zl][Q1]=delt;
3886
3887       //-- Daten aus T2_Transform sind eingetragen; haenge Daten von HN an: --
3888       hilfid=HNEakut[zeile+zl];
3889       for (zl1=Q1+1; zl1<=ncols(aneu[zaehler][2]); zl1++) {
3890        hilfid[zl1]=aneu[zaehler][2][zl1];
3891       }
3892       HNEakut[zeile+zl]=hilfid;
3893       // ------ vorher HNEs[.][zeile+zl]<-aneu[.][2],
3894       //        jetzt [zeile+zl+1] <- [3] usw.: --------
3895       for (zl1=3; zl1<=size(aneu[zaehler]); zl1++) {
3896         HNEakut[zeile+zl+zl1-2]=aneu[zaehler][zl1];
3897       }
3898       //--- setze hqs zusammen: HNEs[hnezaehler][1]=HNEs[..][1],aneu[..][1] --
3899       hilfvec=HNEakut[1],aneu[zaehler][1];
3900       HNEakut[1]=hilfvec;
3901       //-------- weil nicht geht: liste[1]=liste[1],aneu[zaehler][1] ---------
3902       HNEs[hnezaehler]=HNEakut;
3903      }                     // end(for zaehler)
3904     kill aneu;
3905     }                      // endelse (R<>0)
3906    }                       // endif (subst()!=0)  (weiteres Aufblasen mit HN)
3907
3908   }                        // end(for j) (Behandlung der einzelnen delta_i)
3909
3910
3911   // ------------------------- new for essdevelop ----------------------------
3912   if ((essential==1) and (defined(SaveRing))) {
3913     // ----- uebertrage Daten in gemeinsame Koerpererweiterung ---------------
3914     if (EXTHNEnumber>number_of_letztring) {
3915       // ----- fuer aktuelle Kante war Koerpererweiterung erforderlich -------
3916       EXTHNEnumber++;
3917       if (not(defined(minPol))) { poly miniPol; }
3918       miniPol=minpoly;
3919       setring SaveRing;
3920       if (not(defined(miniPol))) { poly miniPol; }
3921       miniPol=minpoly;
3922       setring HNE_noparam;
3923       if (not(defined(a_x))){ map a_x,a_y; poly mp_save, mp_new; }
3924       mp_save=imap(SaveRing,miniPol);
3925       mp_new=imap(EXTHNEring(EXTHNEnumber-1),miniPol);
3926       if (mp_save==mp_new) { // Sonderfall: wieder gleicher Ring
3927         def EXTHNEring(EXTHNEnumber)=SaveRing;
3928         setring EXTHNEring(EXTHNEnumber);
3929         if (not(defined(f))) {poly f; f=hole(f); export f;}
3930         list dummyL=imap(EXTHNEring(EXTHNEnumber-1),HNEs);
3931         if (not(defined(HNEs))) { def HNEs=list(); }
3932         if ((size(HNEs)==1) and (typeof(HNEs[1])=="ideal")) {HNEs=list();}
3933         HNEs[size(HNEs)+1..size(HNEs)+size(dummyL)]=dummyL[1..size(dummyL)];
3934         kill dummyL,SaveRing;
3935       }
3936       else { // verschiedene Ringe
3937         a_x=HNE_noparam,x,0,0;
3938         a_y=HNE_noparam,y,0,0;
3939         mp_save=a_x(mp_save); // minpoly aus SaveRing mit a --> x
3940         mp_new=a_y(mp_new);   // minpoly aus SaveRing mit a --> y
3941         setring SaveRing;
3942         poly mp_new=imap(HNE_noparam,mp_new);
3943         list Lfac=factorize(mp_new,1);
3944         poly fac=Lfac[1][1];
3945         for (k=2;k<=size(Lfac[1]);k++) {
3946           if (deg(Lfac[1][k])<deg(fac)) { fac=Lfac[1][k]; }
3947         }
3948
3949         if (deg(fac)==1) { // keine Erweiterung noetig
3950           def EXTHNEring(EXTHNEnumber)=SaveRing;
3951           setring HNE_noparam;
3952           HNEs=imap(EXTHNEring(EXTHNEnumber-1),HNEs);
3953           setring EXTHNEring(EXTHNEnumber);
3954           if (not(defined(f))) {poly f; f=hole(f); export f;}
3955           map phi=HNE_noparam,-subst(fac,y,0)/koeff(fac,0,1),x,y;
3956           list dummyL=phi(HNEs);
3957           if (not(defined(HNEs))) { def HNEs=list(); }
3958           if ((size(HNEs)==1) and (typeof(HNEs[1])=="ideal")) {HNEs=list();}
3959           HNEs[size(HNEs)+1..size(HNEs)+size(dummyL)]=dummyL[1..size(dummyL)];
3960           kill dummyL,phi,SaveRing;
3961         }
3962         else { // Koerpererweiterung noetig
3963           def EXTHNEring(EXTHNEnumber) = splitring(fac,list(a,transfproc));
3964           setring EXTHNEring(EXTHNEnumber);
3965           poly transf=erg[1];  // image of parameter from SaveRing
3966           poly transfproc=erg[2];
3967           poly transb=erg[3];  // image of parameter from EXTHNEring(..)
3968           export transfproc;
3969           kill erg;
3970           setring HNE_noparam;
3971           if (not(defined(HNEs1))) { list HNEs1=ideal(0); }
3972           HNEs1=imap(EXTHNEring(EXTHNEnumber-1),HNEs);
3973           if (not(defined(hne))) { list hne=ideal(0); }
3974           hne=imap(SaveRing,hne);
3975           HNEs=imap(SaveRing,HNEs);
3976           setring EXTHNEring(EXTHNEnumber);
3977           map hole=HNE_noparam,transf,x,y;
3978           poly fneu=hole(fneu);
3979           poly f=hole(f);
3980           map phi=HNE_noparam,transb,x,y;
3981           list HNEs=hole(HNEs);
3982           list hne=hole(hne);
3983           export f,hne;
3984           if ((size(HNEs)==1) and (typeof(HNEs[1])=="ideal")) {HNEs=list();}
3985           list dummyL=phi(HNEs1);
3986           HNEs[size(HNEs)+1..size(HNEs)+size(dummyL)]=dummyL[1..size(dummyL)];
3987           kill dummyL,phi,SaveRing;
3988         }
3989       }
3990     }
3991     else { // nur bei letzter Kante muss was getan werden
3992       if (i==size(Newton)-1) {
3993         EXTHNEnumber++;
3994         if (number_of_letztring==0) { def letztring=HNEring; }
3995         else       { def letztring=EXTHNEring(EXTHNEnumber); }
3996         if (minpoly==0) {
3997           def EXTHNEring(EXTHNEnumber)=SaveRing;
3998           setring EXTHNEring(EXTHNEnumber);
3999           if (not(defined(f))) {poly f; f=hole(f); export f;}
4000           if ((size(HNEs)==1) and (typeof(HNEs[1])=="ideal")) {HNEs=list();}
4001           list dummyL=imap(letztring,HNEs);
4002           HNEs[size(HNEs)+1..size(HNEs)+size(dummyL)]=dummyL[1..size(dummyL)];
4003           kill dummyL,letztring,SaveRing;
4004         }
4005         else { // muessen Daten nach SaveRing uebertragen;
4006           setring HNE_noparam;
4007           if (not(defined(HNEs))) { list HNEs; }
4008           HNEs=imap(letztring,HNEs);
4009           def EXTHNEring(EXTHNEnumber)=SaveRing;
4010           setring EXTHNEring(EXTHNEnumber);
4011           if (not(defined(hole))) { map hole; }
4012           hole=HNE_noparam,transfproc,x,y;
4013           list dummyL=hole(HNEs);
4014           if (not(defined(HNEs))) { def HNEs=list(); }
4015           if ((size(HNEs)==1) and (typeof(HNEs[1])=="ideal")) {HNEs=list();}
4016           HNEs[size(HNEs)+1..size(HNEs)+size(dummyL)]=dummyL[1..size(dummyL)];
4017           kill dummyL, letztring,SaveRing;
4018         }
4019       }
4020     }
4021   }
4022   // -----------------end of new part (loop for essential=1) ----------------
4023  } // end (Loop uber Kanten)
4024  if (defined(SaveRing)) { kill SaveRing; }
4025 }
4026 else {
4027  HNEs[1]=list(hqs)+azeilen+list(fneu); // fneu ist transform. Polynom oder Null
4028  conj1[1]=conj_factor;
4029 }
4030
4031 if (Aufruf_Ebene == 1)
4032 {
4033   if ((number_of_letztring!=EXTHNEnumber) and (not(defined(hne))))
4034   {
4035     //----- falls Zweige in transz. Erw. berechnet werden konnten ---------
4036     if (defined(transfproc))
4037     { // --- Ringwechsel hat stattgefunden ---
4038       if (defined(HNDebugOn)) {" ring change in HN(",1,") detected";}
4039       if (not(defined(hole))) { map hole; }
4040       hole=HNE_noparam,transfproc,x,y;
4041       setring HNE_noparam;
4042       f=imap(HNEring,f);
4043       if (number_of_letztring==0) { def letztring=HNEring; }
4044       else                        { def letztring=EXTHNEring(EXTHNEnumber); }
4045       if (not(defined(hne))) { list hne; }
4046       hne=imap(letztring,hne);
4047       setring EXTHNEring(EXTHNEnumber);
4048       if (not(defined(f))) { poly f=hole(f); export f; }
4049       list hne=hole(hne);
4050       export hne;
4051     }
4052   }
4053   if (size(HNEs)>0) {
4054     if ((size(HNEs)>1) or (typeof(HNEs[1])!="ideal") or (size(HNEs[1])>0)) {
4055       if ((typeof(hne[1])=="ideal")) { hne=list(); }
4056       hne=hne+extractHNEs(HNEs,getauscht);
4057       if (hne_conj==0) { hne_conj=conj1; }
4058       else { hne_conj = hne_conj, conj1; }
4059     }
4060   }
4061 }
4062 else
4063 { // HN wurde rekursiv aufgerufen
4064   if (number_of_letztring!=EXTHNEnumber)
4065   { // Ringwechsel hatte stattgefunden
4066     string mipl_alt = string(minpoly);
4067     execute("ring tempRing = ("+charstr(basering)+"),("+varstr(basering)+
4068                              "),("+ordstr(basering)+");");
4069     execute("minpoly="+ mipl_alt +";");
4070     list HNEs=imap(EXTHNEring(EXTHNEnumber),HNEs);
4071     export HNEs;
4072     if (defined(HNDebugOn)) {" ! tempRing defined ! ";}
4073   }
4074   if (conj1!=0) { hne_conj=conj1; }
4075   else          { hne_conj=conj_factor; }
4076 }
4077 if (EXTHNEnumber>0){ EXTHNEring = EXTHNEring(1..EXTHNEnumber); }
4078 HNE_RingDATA = list(HNEring, HNE_noparam, EXTHNEnumber, EXTHNEring);
4079 if (number_of_letztring==EXTHNEnumber) {
4080   return(list(HNE_RingDATA,EXTHNEnumber,field_ext,HNEs,hne_conj));
4081 }
4082 else {
4083   if (defined(tempRing)) {
4084     return(list(HNE_RingDATA,EXTHNEnumber,field_ext,tempRing,hne_conj));
4085   }
4086   return(list(HNE_RingDATA,EXTHNEnumber,field_ext,0,hne_conj));
4087 }
4088}
4089
4090///////////////////////////////////////////////////////////////////////////////
4091
4092static proc constructHNEs (list HNEs,int hnezaehler,list aneu,list azeilen,
4093                    int zeile,ideal deltais,int Q,int j)
4094"NOTE: This procedure is only for internal use, it is called via HN"
4095{
4096  int zaehler,zl;
4097  ideal hilfid;
4098  list hilflist;
4099  intvec hilfvec;
4100  for (zaehler=1; zaehler<=size(aneu); zaehler++) {
4101     hnezaehler++;
4102     // HNEs[hnezaehler]=azeilen;            // dupliziere gemeins. Anfang
4103 //----------------------- trage neu berechnete Daten ein ---------------------
4104     hilfid=azeilen[zeile+2];
4105     hilfid[Q]=deltais[j];
4106     for (zl=Q+1; zl<=ncols(aneu[zaehler][2]); zl++) {
4107      hilfid[zl]=aneu[zaehler][2][zl];
4108     }
4109     hilflist=azeilen; hilflist[zeile+2]=hilfid;
4110 //------------------ haenge uebrige Zeilen von aneu[] an: --------------------
4111     for (zl=3; zl<=size(aneu[zaehler]); zl++) {
4112      hilflist[zeile+zl]=aneu[zaehler][zl];
4113     }
4114 //--- setze die hqs zusammen: HNEs[hnezaehler][1]=HNEs[..][1],aneu[..][1] ----
4115     if (hilflist[1]==0) { hilflist[1]=aneu[zaehler][1]; }
4116     else { hilfvec=hilflist[1],aneu[zaehler][1]; hilflist[1]=hilfvec; }
4117     HNEs[hnezaehler]=hilflist;
4118  }
4119  return(HNEs,hnezaehler);
4120}
4121///////////////////////////////////////////////////////////////////////////////
4122
4123proc referencepoly (list newton)
4124"USAGE:   referencepoly(newton);
4125         newton is list of intvec(x,y) which represents points in the Newton
4126         diagram (e.g. output of the proc newtonpoly)
4127RETURN:  a polynomial which has newton as Newton diagram
4128SEE ALSO: newtonpoly
4129EXAMPLE: example referencepoly;  shows an example
4130"
4131{
4132 poly f;
4133 for (int i=1; i<=size(newton); i++) {
4134   f=f+var(1)^newton[i][1]*var(2)^newton[i][2];
4135 }
4136 return(f);
4137}
4138example
4139{ "EXAMPLE:"; echo = 2;
4140 ring exring=0,(x,y),ds;
4141 referencepoly(list(intvec(0,4),intvec(2,3),intvec(5,1),intvec(7,0)));
4142}
4143///////////////////////////////////////////////////////////////////////////////
4144
4145proc factorlist (list L, list #)
4146"USAGE:   factorlist(L);   L a list in the format of `factorize'
4147RETURN:  the nonconstant irreducible factors of
4148         L[1][1]^L[2][1] * L[1][2]^L[2][2] *...* L[1][size(L[1])]^...
4149         with multiplicities (same format as factorize)
4150SEE ALSO: factorize
4151EXAMPLE: example factorlist;  shows an example
4152"
4153{
4154 int k;
4155 if ((size(#)>=1) and ((typeof(#[1])=="intvec") or (typeof(#[1])=="int"))) {
4156   int with_conj = 1; intvec C = #[1];
4157 }
4158 else {
4159   int with_conj = 0; intvec C = L[2];
4160 }
4161 // eine Sortierung der Faktoren eruebrigt sich, weil keine zwei versch.
4162 // red.Fakt. einen gleichen irred. Fakt. haben koennen (I.3.27 Diplarb.)
4163 int i,gross;
4164 list faktoren,hilf;
4165 intvec conjugates;
4166 ideal hil1,hil2;
4167 intvec v,w,hilf_conj;
4168 for (i=1; (L[1][i] == jet(L[1][i],0)) && (i<size(L[1])); i++) {;}
4169 if (L[1][i] != jet(L[1][i],0)) {
4170   hilf=factorize(L[1][i]);
4171 // Achtung!!! factorize(..,2) wirft entgegen der Beschreibung nicht nur
4172 // konstante Faktoren raus, sondern alle Einheiten in der LOKALISIERUNG nach
4173 // der Monomordnung!!! Im Beispiel unten verschwindet der Faktor x+y+1, wenn
4174 // man ds statt dp als Ordnung nimmt!
4175   hilf_conj=C[i];
4176   for (k=2;k<=size(hilf[2]);k++){ hilf_conj=hilf_conj,C[i]; }
4177   hilf[2]=hilf[2]*L[2][i];
4178   hil1=hilf[1];
4179   gross=size(hil1);
4180   if (gross>1) {
4181     v=hilf[2];
4182     faktoren=list(ideal(hil1[2..gross]),intvec(v[2..gross]));
4183     conjugates=intvec(hilf_conj[2..gross]);
4184   }
4185   else         { faktoren=hilf; conjugates=hilf_conj; }
4186 }
4187 else {
4188   faktoren=L;
4189   conjugates=C;
4190 }
4191
4192 for (i++; i<=size(L[2]); i++) {
4193 //------------------------- linearer Term -- irreduzibel ---------------------
4194   if (L[1][i] == jet(L[1][i],1)) {
4195     if (L[1][i] != jet(L[1][i],0)) {           // konst. Faktoren eliminieren
4196       hil1=faktoren[1];
4197       hil1[size(hil1)+1]=L[1][i];
4198       faktoren[1]=hil1;
4199       v=faktoren[2],L[2][i];
4200       conjugates=conjugates,C[i];
4201       faktoren[2]=v;
4202     }
4203   }
4204 //----------------------- nichtlinearer Term -- faktorisiere -----------------
4205   else {
4206     hilf=factorize(L[1][i]);
4207     hilf_conj=C[i];
4208     for (k=2;k<=size(hilf[2]);k++){ hilf_conj=hilf_conj,C[i]; }
4209     hilf[2]=hilf[2]*L[2][i];
4210     hil1=faktoren[1];
4211     hil2=hilf[1];
4212     gross=size(hil2);
4213       // hil2[1] ist konstant, wird weggelassen:
4214     hil1[(size(hil1)+1)..(size(hil1)+gross-1)]=hil2[2..gross];
4215       // ideal+ideal does not work due to simplification;
4216       // ideal,ideal not allowed
4217     faktoren[1]=hil1;
4218     w=hilf[2];
4219     v=faktoren[2],w[2..gross];
4220     conjugates=conjugates,hilf_conj[2..gross];
4221     faktoren[2]=v;
4222   }
4223 }
4224 if (with_conj==0) { return(faktoren); }
4225 else { return(list(faktoren,conjugates)); }  // for essential development
4226}
4227example
4228{ "EXAMPLE:"; echo = 2;
4229 ring exring=0,(x,y),ds;
4230 list L=list(ideal(x,(x-y)^2*(x+y+1),x+y),intvec(2,2,1));
4231 L;
4232 factorlist(L);
4233}
4234
4235///////////////////////////////////////////////////////////////////////////////
4236
4237proc delta
4238"USAGE:  delta(INPUT);  INPUT a polynomial defining an isolated plane curve
4239         singularity at 0, or the Hamburger-Noether expansion thereof, i.e.
4240         the output of @code{develop(f)}, or the output of @code{hnexpansion(f)},
4241         or the list of HN data computed by @code{hnexpansion(f)}.
4242RETURN:  int, the delta invariant of the singularity at 0, that is, the vector
4243         space dimension of R~/R, (R~ the normalization of the local ring of
4244         the singularity).
4245NOTE:    In case the Hamburger-Noether expansion of the curve f is needed
4246         for other purposes as well it is better to calculate this first
4247         with the aid of @code{hnexpansion} and use it as input instead of
4248         the polynomial itself.
4249SEE ALSO: deltaLoc, invariants
4250KEYWORDS: delta invariant
4251EXAMPLE: example delta;  shows an example
4252"
4253{
4254  if (typeof(#[1])=="poly") { // INPUT = polynomial defining the singularity
4255    list L=hnexpansion(#[1]);
4256    if (typeof(L[1])=="ring") {
4257      def altring = basering;
4258      def HNring = L[1]; setring HNring;
4259      return(delta(hne));
4260    }
4261    else {
4262      return(delta(L));
4263    }
4264  }
4265  if (typeof(#[1])=="ring") { // INPUT = HNEring of curve
4266    def HNring = #[1]; setring HNring;
4267    return(delta(hne));
4268  }
4269  if (typeof(#[1])=="matrix")
4270  { // INPUT = hne of an irreducible curve
4271    return(invariants(#)[5]/2);
4272  }
4273  else
4274  { // INPUT = hne of a reducible curve
4275    list INV=invariants(#);
4276    return(INV[size(INV)][3]);
4277  }
4278}
4279example
4280{ "EXAMPLE:"; echo = 2;
4281  ring r = 32003,(x,y),ds;
4282  poly f = x25+x24-4x23-1x22y+4x22+8x21y-2x21-12x20y-4x19y2+4x20+10x19y
4283           +12x18y2-24x18y-20x17y2-4x16y3+x18+60x16y2+20x15y3-9x16y
4284           -80x14y3-10x13y4+36x14y2+60x12y4+2x11y5-84x12y3-24x10y5
4285           +126x10y4+4x8y6-126x8y5+84x6y6-36x4y7+9x2y8-1y9;
4286  delta(f);
4287}
4288
4289///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.