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

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