source: git/Singular/LIB/hnoether.lib @ 2e7859

spielwiese
Last change on this file since 2e7859 was 2e7859, checked in by Olaf Bachmann <obachman@…>, 24 years ago
* SUB_PROCEDURES to AUXILIARY PROCEDURES for correct parsing git-svn-id: file:///usr/local/Singular/svn/trunk@4356 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 128.0 KB
Line 
1// $Id: hnoether.lib,v 1.21 2000-05-19 11:41:05 obachman Exp $
2// last change:           27.04.00
3///////////////////////////////////////////////////////////////////////////////
4// This library is for Singular 1-3-7 or newer
5
6version="$Id: hnoether.lib,v 1.21 2000-05-19 11:41:05 obachman Exp $";
7info="
8LIBRARY:  hnoether.lib   Procedures for the Hamburger-Noether (Puiseux) development
9
10AUTHOR:   Martin Lamm,   lamm@mathematik.uni-kl.de
11
12OVERVIEW:
13           A library for computing the Hamburger-Noether resp. Puiseux
14           development of a plane curve singularity following A. Campillo:
15           Algebroid curves in positive characteristic, Lecture Notes in Math.
16           813, Springer (1980).
17           The library contains also procedures for computing the (topological)
18           numerical invariants of plane curve singularities.
19
20MAIN PROCEDURES:
21 HNdevelop(f);       Hamburger-Noether development of f
22 develop(f [,n]);    Hamburger-Noether development of irreducible curves
23 reddevelop(f);      Hamburger-Noether development of reducible curves
24 essdevelop(f);      Hamburger-Noether development of essential branches
25 extdevelop(hne,n);  extension of the Hamburger-Noether development hne of f
26 param(hne [,x]);    returns a parametrization of f (input=output(develop))
27 displayHNE(hne);    display Hamburger-Noether development as an ideal
28 invariants(hne);    invariants of f, e.g. the characteristic exponents
29 displayInvariants(hne);  display invariants of f
30 intersection(hne1,hne2); intersection multiplicity of two curves
31 multsequence(hne);  sequence of multiplicities
32 displayMultsequence(hne);  display sequence of multiplicities
33 stripHNE(hne);      reduce amount of memory consumed by hne
34 is_irred(f);        test for irreducibility
35
36AUXILIARY PROCEDURES:
37 puiseux2generators(m,n); convert Puiseux pairs to generators of semigroup
38 multiplicities(hne);     multiplicities of blowed up curves
39 separateHNE(hne1,hne2);  number of quadratic transf. needed for separation
40 squarefree(f);      returns a squarefree divisor of the poly f
41 allsquarefree(f,l); returns the maximal squarefree divisor of the poly f
42 further_hn_proc();  show further procedures useful for interactive use
43
44KEYWORDS: Hamburger-Noether expansion; curve singularities
45";
46///////////////////////////////////////////////////////////////////////////////
47LIB "primitiv.lib";
48LIB "inout.lib";
49
50///////////////////////////////////////////////////////////////////////////////
51
52proc further_hn_proc()
53"USAGE: further_hn_proc();
54   The library `hnoether.lib' contains some more procedures which are not
55   shown when typing @code{help hnoether.lib}. They may be useful for
56   interactive use (e.g. if you want to do the calculation of a HNE
57   \"by hand\" to see the intermediate results), and they can be enumerated by
58   calling further_hn_proc. Use @code{help <procedure>;} for detailed
59   information about each of them.
60"
61{
62 "
63 The following procedures are also part of the library `hnoether.lib':
64
65 newtonpoly(f);      Newton polygon of polynom f
66 getnm(f);           intersection pts. of Newton polygon with the axes
67 T_Transform(f,Q,N); returns f(y,xy^Q)/y^NQ (f: poly, Q,N: int)
68 T1_Transform(f,d,M); returns f(x,y+d*x^M)  (f: poly,d:number,M:int)
69 T2_Transform(f,d,M,N,ref);   a composition of T1 & T
70 koeff(f,I,J);       gets coefficient of indicated monomial of poly f
71 redleit(f,S,E);     restriction of monomials of f to line (S-E)
72 leit(f,n,m);        special case of redleit (for irred. polynomials)
73 testreducible(f,n,m); tests whether f is reducible
74 charPoly(f,M,N);    characteristic polynomial of f
75 find_in_list(L,p);  find int p in list L
76 get_last_divisor(M,N); last divisor in Euclid's algorithm
77 factorfirst(f,M,N); try to factor f in a trivial way without `factorize'
78 factorlist(L);      factorize a list L of polynomials
79 referencepoly(D);   a polynomial f s.t. D is the Newton diagram of f";
80
81// extrafactor(f,M,N); try to factor charPoly(f,M,N) where 'factorize' cannot
82// (will hopefully become obsolete soon)
83//
84//       static procedures not useful for interactive use:
85// polytest(f);        tests coefficients and exponents of poly f
86// extractHNEs(H,t);   extracts output H of HN to output of reddevelop
87// HN(f,grenze);       recursive subroutine for reddevelop
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, 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/ggt; N=N/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 / 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   "delta =",d,"Q =",Q,"tau,sigma =",tau,sigma;
218  }
219 //------------------- Durchfuehrung der Transformation T2 --------------------
220  f=Tstar(f);
221  refpoly=Tstar(refpoly);
222 //--- dividiere f so lange durch x & y, wie die Div. aufgeht, benutze ein  ---
223 //--- Referenzpolynom mit gleichem Newtonpolynom wie f zur Beschleunigung: ---
224  for (hilf=refpoly/xp; hilf*xp==refpoly; hilf=refpoly/xp) {refpoly=hilf; s++;}
225  for (hilf=refpoly/yp; hilf*yp==refpoly; hilf=refpoly/yp) {refpoly=hilf; t++;}
226  f=f/(xp^s*yp^t);
227  return(list(T1(f),d));
228}
229example
230{ "EXAMPLE:"; echo = 2;
231  ring exrg=0,(x,y),ds;
232  export exrg;
233  poly f=y2-2x2y+x6-x5y+x4y2;
234  T2_Transform(f,1/2,4,1,f);
235  T2_Transform(f,1/2,4,1,referencepoly(newtonpoly(f)));
236  // if  size(referencepoly) << size(f)  the 2nd example would be faster
237  referencepoly(newtonpoly(f));
238  kill exrg;
239}
240///////////////////////////////////////////////////////////////////////////////
241
242proc koeff (poly f, int I, int J)
243"USAGE:   koeff(f,I,J); f bivariate polynomial, I,J integers
244RETURN:  if f = sum(a(i,j)*x^i*y^j), then koeff(f,I,J)= a(I,J) (of type number)
245NOTE:    J must be in the range of the exponents of the 2nd ring variable
246EXAMPLE: example koeff;  shows an example
247"
248{
249  matrix mat = coeffs(coeffs(f,var(2))[J+1,1],var(1));
250  if (size(mat) <= I) { return(0);}
251  else { return(leadcoef(mat[I+1,1]));}
252}
253example
254{ "EXAMPLE:"; echo = 2;
255  ring r=0,(x,y),dp;
256  koeff(x2+2xy+3xy2-x2y-2y3,1,2);
257}
258///////////////////////////////////////////////////////////////////////////////
259
260proc squarefree (poly f)
261"USAGE:  squarefree(f);  f a bivariate polynomial
262RETURN: a squarefree divisor of f.
263        Normally the return value is a greatest squarefree divisor, but there
264        is an exception: factors with a p-th root, p the ring characteristic,
265        are lost.
266SEE ALSO: allsquarefree
267EXAMPLE: example squarefree; shows some examples.
268"
269{
270 //----------------- Wechsel in geeigneten Ring & Variablendefinition ---------
271  def altring = basering;
272  int e;
273  int gcd_ok=1;
274  string mipl="0";
275  if (size(parstr(altring))==1) { mipl=string(minpoly); }
276 //---- test: char = (p^k,a) (-> gcd not implemented) or (p,a) (gcd works) ----
277  if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) {
278    string tststr=charstr(basering);
279    tststr=tststr[1..find(tststr,",")-1];           //-> "p^k" bzw. "p"
280    gcd_ok=(tststr==string(char(basering)));
281  }
282  execute("ring rsqrf = ("+charstr(altring)+"),(x,y),dp;");
283  if ((gcd_ok!=0) && (mipl!="0")) { execute("minpoly="+mipl+";"); }
284  poly f=fetch(altring,f);
285  poly dif,g,l;
286  if ((char(basering)==0) and (charstr(basering)!=string(char(basering)))
287      and (mipl!="0")) {
288    gcd_ok=0;                   // since Singular 1.2 gcd no longer implemented
289  }
290  if (gcd_ok!=0) {
291 //--------------------- Berechne f/ggT(f,df/dx,df/dy) ------------------------
292    dif=diff(f,x);
293    if (dif==0) { g=f; }        // zur Beschleunigung
294    else { g=gcd(f,dif); }
295    if (g!=1) {                 // sonst schon sicher, dass f quadratfrei
296     dif=diff(f,y);
297     if (dif!=0) { g=gcd(g,dif); }
298    }
299    if (g!=1) {
300     e=0;
301     if (g==f) { l=1; }         // zur Beschleunigung
302     else {
303       module m=syz(ideal(g,f));
304       if (deg(m[2,1])>0) {
305         "!! The Singular command 'syz' has returned a wrong result !!";
306         l=1;                   // Division f/g muss aufgehen
307       }
308       else { l=m[1,1]; }
309     }
310    }
311    else { e=1; }
312  }
313  else {
314 //------------------- Berechne syz(f,df/dx) oder syz(f,df/dy) ----------------
315 //-- Achtung: Ist f reduzibel, koennen Faktoren mit Ableitung Null verloren --
316 //-- gehen! Ist aber nicht weiter schlimm, weil char (p^k,a) nur im irred.  --
317 //-- Fall vorkommen kann. Wenn f nicht g^p ist, wird auf jeden Fall         --
318 //------------------------ ein Faktor gefunden. ------------------------------
319    dif=diff(f,x);
320    if (dif == 0) {
321     dif=diff(f,y);
322     if (dif==0) { e=2; l=1; } // f is of power divisible by char of basefield
323     else { l=syz(ideal(dif,f))[1,1];  // x^p+y^(p-1) abgedeckt
324            if (subst(f,x,0)==0) { l=l*x; }
325            if (deg(l)==deg(f))  { e=1;}
326            else {e=0;}
327     }
328    }
329    else { l=syz(ideal(dif,f))[1,1];
330           if (subst(f,y,0)==0) { l=l*y; }
331           if (deg(l)==deg(f))  { e=1;}
332           else {e=0;}
333    }
334  }
335 //--------------- Wechsel in alten Ring und Rueckgabe des Ergebnisses --------
336  setring altring;
337  if (e==1) { return(f); }    // zur Beschleunigung
338  else {
339   poly l=fetch(rsqrf,l);
340   return(l);
341  }
342}
343example
344{ "EXAMPLE:"; echo = 2;
345 ring exring=3,(x,y),dp;
346 squarefree(x3+y);
347 squarefree(x3);
348 squarefree(x2);
349 squarefree(xy+x);
350 squarefree((x+y)^3*(x-y)^2); // (x+y)^3 is lost
351 squarefree((x+y)^4*(x-y)^2); // result is (x+y)*(x-y)
352}
353///////////////////////////////////////////////////////////////////////////////
354
355proc allsquarefree (poly f, poly l)
356"USAGE : allsquarefree(f,l);  poly f,l
357        l is the squarefree divisor of f you get by @code{l=squarefree(f);}
358RETURN: the squarefree divisor of f consisting of all irreducible factors of f
359NOTE  : This proc uses factorize to get the missing factors of f not in l and
360        therefore may be slow.
361SEE ALSO: squarefree
362EXAMPLE: example allsquarefree;  shows an example
363"
364{
365 //------------------------ Wechsel in geeigneten Ring ------------------------
366 def altring = basering;
367 string mipl="0";
368 if (size(parstr(altring))==1) { mipl=string(minpoly); }
369 if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) {
370   string tststr=charstr(basering);
371   tststr=tststr[1..find(tststr,",")-1];           //-> "p^k" bzw. "p"
372   if (tststr!=string(char(basering))) {
373     " Sorry -- not implemented for this ring (gcd doesn't work)";
374     return(l);
375   }
376 }
377 execute("ring rsqrf = ("+charstr(altring)+"),(x,y),dp;");
378 if (mipl!="0") { execute("minpoly="+mipl+";"); }
379 poly f=fetch(altring,f);
380 poly l=fetch(altring,l);
381 //---------- eliminiere bereits mit squarefree gefundene Faktoren ------------
382 poly g=l;
383 while (deg(g)!=0) {
384   f=syz(ideal(g,f))[1,1];                         // f=f/g;
385   g=gcd(f,l);
386 }                                                 // jetzt f=h^p
387 //--------------- Berechne uebrige Faktoren mit factorize --------------------
388 if (deg(f)>0) {
389  g=1;
390  ideal factf=factorize(f,1);
391  for (int i=1; i<=size(factf); i++) { g=g*factf[i]; }
392  poly testp=squarefree(g);
393  if (deg(testp)<deg(g)) {
394    "!! factorize has not worked correctly !!";
395    if (testp==1) {" We cannot proceed ..."; g=1;}
396    else {" But we could recover some factors..."; g=testp;}
397  }
398  l=l*g;
399 }
400 //--------------- Wechsel in alten Ring und Rueckgabe des Ergebnisses --------
401 setring altring;
402 l=fetch(rsqrf,l);
403 return(l);
404}
405example
406{ "EXAMPLE:"; echo = 2;
407  ring exring=7,(x,y),dp;
408  poly f=(x+y)^7*(x-y)^8;
409  poly l=squarefree(f);
410  l; // x-y found because 7 does not divide 8
411  allsquarefree(f,l); // all factors (x+y)*(x-y) found
412}
413///////////////////////////////////////////////////////////////////////////////
414
415proc is_irred (poly f)
416"USAGE:   is_irred(f); f a squarefree bivariate polynomial
417RETURN:  an integer:
418 @*      is_irred(f) @math{= 1} if f is irreducible as a formal power series
419         over the algebraic closure of its coefficient field (f defines an
420         analytically irreducible curve at zero),
421 @*      is_irred(f) @math{= 0} else.
422NOTE:    @math{f = 0} and units in the ring of formal power series are not
423         considered irreducible.
424KEYWORDS: irreducible power series
425EXAMPLE: example is_irred;  shows an example
426"
427{
428  int pl=printlevel;
429  printlevel=-1;
430  list hnl=develop(f,-1);
431  printlevel=pl;
432  return(hnl[5]);
433}
434example
435{ "EXAMPLE:"; echo = 2;
436  ring exring=0,(x,y),ls;
437  is_irred(x2+y3);
438  is_irred(x2+y2);
439  is_irred(x2+y3+1);
440}
441///////////////////////////////////////////////////////////////////////////////
442
443static proc polytest(poly f)
444"USAGE : polytest(f); f poly in x and y
445RETURN: a monomial of f with |coefficient| > 16001
446          or exponent divisible by 32003, if there is one
447        0 else (in this case computing a squarefree divisor
448                in characteristic 32003 could make sense)
449NOTE:   this procedure is only useful in characteristic zero, because otherwise
450        there is no appropriate ordering of the leading coefficients
451"
452{
453 poly verbrecher=0;
454 intvec leitexp;
455 for (; (f<>0) and (verbrecher==0); f=f-lead(f)) {
456  if ((leadcoef(f)<-16001) or (leadcoef(f)>16001)) {verbrecher=lead(f);}
457  leitexp=leadexp(f);
458  if (( ((leitexp[1] % 32003) == 0)   and (leitexp[1]<>0))
459     or ( ((leitexp[2] % 32003) == 0) and (leitexp[2]<>0)) )
460       {verbrecher=lead(f);}
461 }
462 return(verbrecher);
463}
464
465//////////////////////////////////////////////////////////////////////////////
466
467
468proc develop
469"USAGE:   develop(f [,n]); f polynomial in two variables, n integer
470RETURN:  Hamburger-Noether development of f, in form of a list:
471  @*     [1]: Hamburger-Noether matrix:
472            Each row contains the coefficients of the corresponding line of the
473            Hamburger-Noether expansion (HNE). The end of the line is marked in
474            the matrix as the first ring variable (usually x).
475  @*     [2]: intvec indicating the length of lines of the HNE
476  @*     [3]: int:
477           0  if the 1st ring variable was transversal (with respect to f),
478  @*       1  if the variables were changed at the beginning of the
479              computation,
480           -1 if an error has occurred.
481  @*     [4]: the transformed polynomial of f to make it possible to extend the
482              Hamburger-Noether development a posteriori without having to do
483              all the previous calculation once again (0 if not needed)
484  @*     [5]: int:
485           1  if the curve has exactly one branch (i.e. is irreducible),
486  @*       0  else (i.e. the curve has more than one HNE or f is not valid).
487
488ASSUME:  f irreducible as a power series (for reducible f use
489         @code{reddevelop})
490DISPLAY: The (non zero) elements of the HNE (if not called by another proc).
491
492NOTE: If the optional parameter @math{n} is given, the HN-matrix will have
493      at least @math{n} columns. Otherwise the number of columns will be chosen
494      minimal s.t. the matrix contains all necessary information (i.e. all
495      lines of the HNE but the last (which is in general infinite) have place).
496      If @math{n} is negative, the algorithm is stopped as soon as possible,
497      i.e. the information computed is enough for `@code{invariants}', but the
498      HN-matrix may contain undetermined elements, which are marked with the
499      2nd variable (of the basering).
500      In any case, the optional parameter only affects the calculation of
501      the LAST line of the HNE; @code{develop}(f) gives already all necessary
502      information for the procedure `@code{invariants}'. A negative value of
503      @math{n} will result in a very poor parametrization, but it can make
504      `@code{develop}' faster; a positive value will improve the exactness of
505      the parametrization.
506
507      For time critical computations it is recommended to use
508      \"@code{ring ...,(x,y),ls}\" as basering - it increases the algorithm's
509      speed.
510
511SEE ALSO: reddevelop, extdevelop, displayHNE
512EXAMPLES: example develop; shows an example
513          example param;   shows an example for using the 2nd parameter
514"
515{
516 //--------- Abfangen unzulaessiger Ringe: 1) nur eine Unbestimmte ------------
517 poly f=#[1];
518 if (size(#) > 1) {int maxspalte=#[2];}
519 else             {int maxspalte= 1 ; }
520 if (nvars(basering) < 2) {
521   " Sorry. I need two variables in the ring.";
522   return(list(matrix(maxideal(1)[1]),intvec(0),-1,poly(0),0));}
523 if (nvars(basering) > 2) {
524   " Warning! You have defined too many variables!";
525   " All variables except the first two will be ignored!";}
526
527 string namex=varstr(1); string namey=varstr(2);
528 list return_error=matrix(maxideal(1)[2]),intvec(0),int(-1),poly(0),int(0);
529
530 //------------- 2) mehrere Unbestimmte, weitere unzulaessige Ringe -----------
531 // Wir koennen einheitlichen Rueckgabewert benutzen, aus dem ersichtlich ist,
532 // dass ein Fehler aufgetreten ist: return_error.
533 //----------------------------------------------------------------------------
534
535 if (charstr(basering)=="real") {
536  " The algorithm doesn't work with 'real' as coefficient field.";
537                     // denn : map from characteristic -1 to -1 not implemented
538  return(return_error);
539 }
540 if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) {
541 //-- teste, ob char = (p^k,a) (-> a primitiv; erlaubt) oder (p,a[,b,...]) ----
542    string tststr=charstr(basering);
543    tststr=tststr[1..find(tststr,",")-1];           //-> "p^k" bzw. "p"
544    int primit=(tststr==string(char(basering)));
545    if (primit!=0) {
546      " Such extensions of Z/p are not implemented.";
547      " Please try (p^k,a) as ground field or use `reddevelop'.";
548      return(return_error);
549    }
550 }
551 //---- Ende der unzulaessigen Ringe; Ringwechsel in einen guenstigen Ring: ---
552
553 int ringwechsel=(varstr(basering)!="x,y") or (ordstr(basering)!="ls(2),C");
554
555 def altring = basering;
556 if (ringwechsel) {
557   string mipl=string(minpoly);
558   execute("ring guenstig = ("+charstr(altring)+"),(x,y),ls;");
559   if ((char(basering)==0) && (mipl!="0")) {
560     execute("minpoly="+mipl+";");
561   }}
562 else { def guenstig=basering; }
563 export guenstig;
564
565 //-------------------------- Initialisierungen -------------------------------
566 map m=altring,x,y;
567 if (ringwechsel) { poly f=m(f); }
568 if (defined(HNDebugOn))
569 {"received polynomial: ",f,", where x =",namex,", y =",namey;}
570
571 int M,N,Q,R,l,e,hilf,eps,getauscht,Abbruch,zeile,exponent,Ausgabe;
572
573 // Werte von Ausgabe: 0 : normale HNE-Matrix,
574 // 1 : Fehler aufgetreten - Matrix (namey) zurueck
575 // 2 : Die HNE ist eine Nullzeile - Matrix (0) zurueck
576 // int maxspalte=1; geaendert: wird jetzt am Anfang gesetzt
577
578 int minimalHNE=0;          // Flag fuer minimale HNE-Berechnung
579 int einzweig=1;            // Flag fuer Irreduzibilit"at
580 intvec hqs;                // erhaelt die Werte von h(zeile)=Q;
581
582 if (maxspalte<0) {
583   minimalHNE=1;
584   maxspalte=1;
585 }
586
587 number c,delta;
588 int p = char(basering);
589 string ringchar=charstr(basering);
590 map xytausch = basering,y,x;
591 if ((p!=0) and (ringchar != string(p))) {
592                            // coefficient field is extension of Z/pZ
593   execute("int n_elements="+ringchar[1,size(ringchar)-2]+";");
594                            // number of elements of actual ring
595   number generat=par(1);   // generator of the coefficient field of the ring
596 }
597
598
599 //========= Abfangen von unzulaessigen oder trivialen Eingaben ===============
600 //------------ Nullpolynom oder Einheit im Potenzreihenring: -----------------
601 if (f == 0) {
602   dbprint(printlevel+1,"You have given me the zero-polynomial !");
603   Abbruch=1; Ausgabe=1;
604 }
605 else {
606   intvec nm = getnm(f);
607   N = nm[1]; M = nm[2]; // Berechne Schnittpunkte Newtonpolygon mit Achsen
608   if (N == 0) {
609     dbprint(printlevel+1,"The given polynomial is a unit as power series !");
610     Abbruch=1; Ausgabe=1;
611   }
612   else {
613    if (N == -1) {
614      if ((voice==2) && (printlevel > -1)) { "The HNE is x = 0"; }
615      Abbruch=1; Ausgabe=2; getauscht=1;
616      if (M <> 1) { einzweig=0; }
617    }
618    else {
619     if (M == -1) {
620       if ((voice==2) && (printlevel > -1)) { "The HNE is y = 0"; }
621       Abbruch=1; Ausgabe=2;
622       if (N <> 1) { einzweig=0; }
623   }}}
624 }
625 //--------------------- Test auf Quadratfreiheit -----------------------------
626 if (Abbruch==0) {
627
628 //-------- Fall basering==0,... : Wechsel in Ring mit char >0 ----------------
629 // weil squarefree eine Standardbasis berechnen muss (verwendet Syzygien)
630 // -- wenn f in diesem Ring quadratfrei ist, dann erst recht im Ring guenstig
631 //----------------------------------------------------------------------------
632
633  if ((p==0) and (size(charstr(basering))==1)) {
634   int testerg=(polytest(f)==0);
635   ring zweitring = 32003,(x,y),dp;
636   map polyhinueber=guenstig,x,y;     // fetch geht nicht
637   poly f=polyhinueber(f);
638   poly test_sqr=squarefree(f);
639   if (test_sqr != f) {
640    if (printlevel>0) {
641      "Most probably the given polynomial is not squarefree. But the test was";
642      "made in characteristic 32003 and not 0 to improve speed. You can";
643      "(r) redo the test in char 0 (but this may take some time)";
644      "(c) continue the development, if you're sure that the polynomial",
645      "IS squarefree";
646      if (testerg==1) {
647        "(s) continue the development with a squarefree factor (*)";}
648      "(q) or just quit the algorithm (default action)";
649      "";"Please enter the letter of your choice:";
650      string str=read("")[1];
651    }
652    else { string str="r"; }      // printlevel <= 0: non-interactive behaviour
653    setring guenstig;
654    map polyhinueber=zweitring,x,y;
655    if (str=="r") {
656      poly test_sqr=squarefree(f);
657      if (test_sqr != f) {
658       if (printlevel>0) { "The given polynomial is in fact not squarefree."; }
659       else              { "The given polynomial is not squarefree!"; }
660       "I'll continue with the radical.";
661       if (printlevel>0) { pause("Hit RETURN to continue:"); }
662       f=test_sqr;
663      }
664      else {
665       dbprint(printlevel,
666        "everything is ok -- the polynomial is squarefree in char(k)=0");
667      }
668    }
669    else {
670      if ((str=="s") and (testerg==1)) {
671       "(*) attention: it could be that the factor is only one in char 32003!";
672        f=polyhinueber(test_sqr);
673      }
674      else {
675        if (str<>"c") {
676          setring altring;kill guenstig;kill zweitring;
677          return(return_error);}
678        else { "if the algorithm doesn't terminate, you were wrong...";}
679    }}
680    kill zweitring;
681    nm = getnm(f);             // N,M haben sich evtl. veraendert
682    N = nm[1]; M = nm[2];      // Berechne Schnittpunkte Newtonpoly mit Achsen
683    if (defined(HNDebugOn)) {"I continue with the polynomial",f; }
684   }
685   else {
686     setring guenstig;
687     kill zweitring;
688   }
689  }
690 // ------------------- Fall Charakteristik > 0 -------------------------------
691  else {
692   poly test_sqr=squarefree(f);
693   if (test_sqr == 1) {
694    "The given polynomial is of the form g^"+string(p)+", therefore",
695    "reducible.";"Please try again.";
696    setring altring;
697    kill guenstig;
698    return(return_error);}
699   if (test_sqr != f) {
700    "The given polynomial is not squarefree. I'll continue with the radical.";
701    if (p != 0)
702     {"But if the polynomial contains a factor of the form g^"+string(p)+",";
703      "this factor will be lost.";}
704    if (printlevel>0) { pause("Hit RETURN to continue:"); }
705    f=test_sqr;
706    nm = getnm(f);              // N,M haben sich veraendert
707    N = nm[1]; M = nm[2];       // Berechne Schnittpunkte Newtonpoly mit Achsen
708    if (defined(HNDebugOn)) {"I continue with the polynomial",f; }
709   }
710
711  }                             // endelse(p==0)
712
713  if (N==0) {
714    " Sorry. The remaining polynomial is a unit in the power series ring...";
715    setring altring;kill guenstig;return(return_error);
716  }
717 //---------------------- gewaehrleiste, dass x transvers ist -----------------
718  if (M < N)
719  { f = xytausch(f);            // Variablentausch : x jetzt transvers
720    getauscht = 1;              // den Tausch merken
721    M = M+N; N = M-N; M = M-N;  // M, N auch vertauschen
722  }
723  if (defined(HNDebugOn)) {
724   if (getauscht) {"x<->y were exchanged; poly is now ",f;}
725   else           {"x , y were not exchanged";}
726   "M resp. N are now",M,N;
727  }
728 }                              // end(if Abbruch==0)
729
730 ideal a(0);
731 while (Abbruch==0) {
732
733 //================= Beginn der Schleife (eigentliche Entwicklung) ============
734
735 //------------------- ist das Newtonpolygon eine gerade Linie? ---------------
736  if (testreducible(f,N,M)) {
737    dbprint(printlevel+1," The given polynomial is not irreducible");
738    kill guenstig;
739    setring altring;
740    return(return_error);       // Abbruch der Prozedur!
741  }
742  R = M%N;
743  Q = M / N;
744
745 //-------------------- Fall Rest der Division R = 0 : ------------------------
746  if (R == 0) {
747    c = koeff(f,0,N);
748    if (c == 0) {"Something has gone wrong! I didn't get N correctly!"; exit;}
749    e = gcd(M,N);
750 //----------------- Test, ob leitf = c*(y^N - delta*x^(m/e))^e ist -----------
751    if (p==0) {
752      delta = koeff(f,M/ e,N - N/ e) / (-1*e*c);
753      if (defined(HNDebugOn)) {"quasihomogeneous leading form:",
754         leit(f,N,M)," = ",c,"* (y -",delta,"* x^"+string(M/ e)+")^",e," ?";}
755      if (leit(f,N,M) != c*(y^(N/ e) - delta*x^(M/ e))^e) {
756        dbprint(printlevel+1," The given polynomial is reducible !");
757        Abbruch=1; Ausgabe=1; }
758    }
759    else {                     // p!=0
760      if (e%p != 0) {
761        delta = koeff(f,M/ e,N - N/ e) / (-1*e*c);
762        if (defined(HNDebugOn)) {"quasihomogeneous leading form:",
763           leit(f,N,M)," = ",c,"* (y -",delta,"* x^"+string(M/ e)+")^",e," ?";}
764        if (leit(f,N,M) != c*(y^(N/ e) - delta*x^(M/ e))^e) {
765           dbprint(printlevel+1," The given polynomial is reducible !");
766           Abbruch=1; Ausgabe=1; }
767      }
768
769      else {                   // e%p == 0
770        eps = e;
771        for (l = 0; eps%p == 0; l=l+1) { eps=eps/ p;}
772        if (defined(HNDebugOn)) {e," -> ",eps,"*",p,"^",l;}
773        delta = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c);
774
775        if ((ringchar != string(p)) and (delta != 0)) {
776 //- coeff. field is not Z/pZ => we`ve to correct delta by taking (p^l)th root-
777          if (delta == generat) {exponent=1;}
778          else {
779           if (delta == 1) {exponent=0;}
780           else {
781            exponent=pardeg(delta);
782
783 //-- an dieser Stelle kann ein Fehler auftreten, wenn wir eine transzendente -
784 //-- Erweiterung von Z/pZ haben: dann ist das hinzuadjungierte Element kein  -
785 //-- Erzeuger der mult. Gruppe, d.h. in Z/pZ (a) gibt es i.allg. keinen      -
786 //-- Exponenten mit z.B. a2+a = a^exp                                        -
787 //----------------------------------------------------------------------------
788          }}
789          delta = generat^(extgcd(n_elements-1,p^l)[3]*exponent);
790        }
791
792        if (defined(HNDebugOn)) {"quasihomogeneous leading form:",
793          leit(f,N,M)," = ",c,"* (y^"+string(N/ e),"-",delta,"* x^"
794          +string(M/ e)+")^",e,"  ?";}
795        if (leit(f,N,M) != c*(y^(N/ e) - delta*x^(M/ e))^e) {
796          dbprint(printlevel+1," The given polynomial is reducible !");
797          Abbruch=1; Ausgabe=1; }
798      }
799    }
800    if (Abbruch == 0) {
801      f = T1_Transform(f,delta,M/ e);
802      if ((voice==2) && (printlevel > -1))
803       { "a("+string(zeile)+","+string(Q)+") =",delta; }
804      a(zeile)[Q]=delta;
805      if (defined(HNDebugOn)) {"transformed polynomial: ",f;}}
806
807      nm=getnm(f); N=nm[1]; M=nm[2];        // Neuberechnung des Newtonpolygons
808  }
809 //--------------------------- Fall R > 0 : -----------------------------------
810  else {
811    if ((voice==2) && (printlevel > -1)) { "h("+string(zeile)+") =",Q; }
812    hqs[zeile+1]=Q;                  // denn zeile beginnt mit dem Wert 0
813    a(zeile)[Q+1]=x;                 // Markierung des Zeilenendes der HNE
814    maxspalte=maxspalte*((Q+1) < maxspalte) + (Q+1)*((Q+1) >= maxspalte);
815                                     // Anpassung der Sp.zahl der HNE-Matrix
816    f = T_Transform(f,Q,N);
817    if (defined(HNDebugOn)) {"transformed polynomial: ",f;}
818    zeile=zeile+1;
819 //------------ Bereitstellung von Speicherplatz fuer eine neue Zeile: --------
820    ideal a(zeile);
821    M=N;N=R;
822  }
823
824 //--------------- schneidet das Newtonpolygon beide Achsen? ------------------
825  if (M==-1) {
826     if ((voice==2) && (printlevel > -1)) { "The HNE is finite!"; }
827     a(zeile)[Q+1]=x;   // Markiere das Ende der Zeile
828     hqs[zeile+1]=Q;
829     maxspalte=maxspalte*((Q+1) < maxspalte) + (Q+1)*((Q+1) >= maxspalte);
830     if (N <> 1) { einzweig=0; }   
831     f=0;               // transformiertes Polynom wird nicht mehr gebraucht
832     Abbruch=1;
833  }
834  else {if (M<N) {"Something has gone wrong: M<N";}}
835  if(defined(HNDebugOn)) {"new M,N:",M,N;}
836
837 //----------------- Abbruchbedingungen fuer die Schleife: --------------------
838  if ((N==1) and (Abbruch!=1) and ((M > maxspalte) or (minimalHNE==1))
839      and (size(a(zeile))>0))
840 //----------------------------------------------------------------------------
841 // Abbruch, wenn die Matrix so voll ist, dass eine neue Spalte angefangen
842 // werden muesste und die letzte Zeile nicht nur Nullen enthaelt
843 // oder wenn die Matrix nicht voll gemacht werden soll (minimale Information)
844 //----------------------------------------------------------------------------
845   { Abbruch=1; hqs[zeile+1]=-1;
846     if (maxspalte < ncols(a(zeile))) { maxspalte=ncols(a(zeile));}
847     if ((minimalHNE==1) and (M <= maxspalte)) {
848 // teile param mit, dass Eintraege der letzten Zeile nur teilw. richtig sind:-
849       hqs[zeile+1]=-M;
850 //------------- markiere den Rest der Zeile als unbekannt: -------------------
851       for (R=M; R <= maxspalte; R++) { a(zeile)[R]=y;}
852     }                  // R wird nicht mehr gebraucht
853   }
854 //========================= Ende der Schleife ================================
855
856 }
857 setring altring;
858 if (Ausgabe == 0) {
859 //-------------------- Ergebnis in den alten Ring transferieren: -------------
860   map zurueck=guenstig,maxideal(1)[1],maxideal(1)[2];
861   matrix amat[zeile+1][maxspalte];
862   ideal uebergabe;
863   for (e=0; e<=zeile; e=e+1) {
864     uebergabe=zurueck(a(e));
865     if (ncols(uebergabe) > 1) {
866      amat[e+1,1..ncols(uebergabe)]=uebergabe;}
867     else {amat[e+1,1]=uebergabe[1];}
868   }
869   if (ringwechsel) {
870     if (nvars(altring)==2) { f=fetch(guenstig,f); }
871     else                   { f=zurueck(f); }
872   }
873 }
874
875 kill guenstig;
876 if ((einzweig==0) && (voice==2) && (printlevel > -1)) {
877    "// Note: The curve is reducible, but we were able to compute a HNE.";
878    "// This means the result is only one of several existing HNE's.";
879 }
880 if (Ausgabe == 0) { return(list(amat,hqs,getauscht,f,einzweig));}
881 if (Ausgabe == 1) { return(return_error);}             // error has occurred
882 if (Ausgabe == 2) { return(list(matrix(ideal(0,x)),intvec(1),getauscht,
883                                 poly(0),einzweig));}   // HNE is x=0 or y=0
884}
885example
886{ "EXAMPLE:"; echo = 2;
887  ring exring = 7,(x,y),ds;
888  list hne=develop(4x98+2x49y7+x11y14+2y14);
889  print(hne[1]);
890  // therefore the HNE is:
891  // z(-1)= 3*z(0)^7 + z(0)^7*z(1),
892  // z(0) = z(1)*z(2),
893  //        (note that there is 1 zero in the 2nd row before x)
894  // z(1) = z(2)^3*z(3),    (there are 3 zeroes in the 3rd row)
895  // z(2) = z(3)*z(4),
896  // z(3) = -z(4)^2 + 0*z(4)^3 +...+ 0*z(4)^8 + ?*z(4)^9 + ...
897  //   (the missing x in the matrix indicates that this line
898  //    is not complete. It can only occur in the last line of the
899  //    HNE, and normally does.)
900  hne[2];
901  pause("Hit RETURN to continue:");
902  param(hne);
903  // returns the parametrization
904  // x(t)= -t^14+O(t^21),  y(t)= -3t^98+O(t^105)
905  // (the term -t^109 in y may have a wrong coefficient)
906  displayHNE(hne);
907}
908
909///////////////////////////////////////////////////////////////////////////////
910//               procedures to extract information out of HNE                //
911///////////////////////////////////////////////////////////////////////////////
912
913proc param
914"USAGE:  param(l [,x]) takes the output of develop(f)
915        (list l (matrix m, intvec v, int s[,poly g,...]))
916        and gives a parametrization for f; the first variable of the active
917        ring is chosen as indefinite. If the ring contains more than
918        two variables, the 3rd variable is chosen (remember that develop takes
919        the first two variables and therefore the other vars should be unused).
920 @*     x is an optional parameter of any type.
921RETURN: If only the list l is given:
922        ideal of two polynomials: if the HNE was finite, f(param[1],param[2])
923        will be  zero. If not, the real parametrization will be
924        two power series; then param will return a truncation of these series.
925
926        If the optional parameter x is given:
927        list L with L[1]=param(l) (the ideal containing the parametrization)
928        and L[2]=intvec with two entries indicating the highest degree unto
929        which the coefficients of the monomials in L[1] are exact (a value of
930        -1 is for infinity, i.e. the corresponding polynomial is exact).
931SEE ALSO: develop, extdevelop
932KEYWORDS: parametrization
933EXAMPLE: example param;     shows an example
934         example develop;   shows another example
935"
936{
937 //-------------------------- Initialisierungen -------------------------------
938 if (typeof(#[1])=="list") {
939   list Li=#[1];
940   matrix m=Li[1];
941   intvec v=Li[2];
942   int switch=Li[3];
943   int return_list=1;
944 }
945 else {
946   matrix m=#[1];
947   intvec v=#[2];
948   int switch=#[3];
949   int return_list=0;
950 }
951 if (switch==-1) {
952   "An error has occurred in develop, so there is no HNE.";
953   return(ideal(0,0));
954 }
955 int fehler,fehlervor,untergrad,untervor,beginn,i,zeile,hilf;
956
957 if (nvars(basering) > 2) { poly z(size(v)+1)=var(3); }
958 else                     { poly z(size(v)+1)=var(1); }
959 poly z(size(v));
960 zeile=size(v);
961 //------------- Parametrisierung der untersten Zeile der HNE -----------------
962 if (v[zeile] > 0) {
963   fehler=0;           // die Parametrisierung wird exakt werden
964   for (i=1; i<=v[zeile]; i++) {
965     z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
966 }}
967 else {
968   untervor=1;         // = Untergrad der vorhergehenden Zeile
969   if (v[zeile]==-1) {
970     fehler=ncols(m)+1;
971     for (i=1; i<=ncols(m); i++) {
972       z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
973       if ((untergrad==0) and (m[zeile,i]!=0)) {untergrad=i;}
974                       // = Untergrad der aktuellen Zeile
975   }}
976   else {
977     fehler= -v[zeile];
978     for (i=1; i<-v[zeile]; i++) {
979       z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
980       if ((untergrad==0) and (m[zeile,i]!=0)) {untergrad=i;}
981   }}
982 }
983 //------------- Parametrisierung der restlichen Zeilen der HNE ---------------
984 for (zeile=size(v)-1; zeile>0; zeile--) {
985   poly z(zeile);
986   beginn=0;             // Beginn der aktuellen Zeile
987   for (i=1; i<=v[zeile]; i++) {
988     z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
989     if ((beginn==0) and (m[zeile,i]!=0)) { beginn=i;}
990   }
991   z(zeile)=z(zeile) + z(zeile+1)^v[zeile] * z(zeile+2);
992   if (beginn==0) {
993     if (fehler>0) {     // damit fehler=0 bleibt bei exakter Param.
994     fehlervor=fehler;   // Fehler der letzten Zeile
995     fehler=fehler+untergrad*(v[zeile]-1)+untervor;   // Fehler dieser Zeile
996     hilf=untergrad;
997     untergrad=untergrad*v[zeile]+untervor;
998     untervor=hilf;}}    // untervor = altes untergrad
999   else {
1000     fehlervor=fehler;
1001     fehler=fehler+untergrad*(beginn-1);
1002     untervor=untergrad;
1003     untergrad=untergrad*beginn;
1004 }}
1005 //--------------------- Ausgabe der Fehlerabschaetzung -----------------------
1006 if (switch==0) {
1007   if (fehler>0) {
1008     if (fehlervor>0) {
1009       if ((voice==2) && (printlevel > -1)) {
1010         "// ** Warning: result is exact up to order",fehlervor-1,"in",
1011         maxideal(1)[1],"and",fehler-1,"in",maxideal(1)[2],"!";
1012     }}
1013     else {
1014       if ((voice==2) && (printlevel > -1)) {
1015         "// ** Warning: result is exact up to order",fehler-1,"in",
1016         maxideal(1)[2],"!";
1017     }}
1018   }
1019   if (return_list==0) { return(ideal(z(2),z(1))); }
1020   else   { return(list(ideal(z(2),z(1)),intvec(fehlervor-1,fehler-1))); }
1021 }
1022 else {
1023   if (fehler>0) {
1024     if (fehlervor>0) {
1025      if ((voice==2) && (printlevel > -1)) {
1026        "// ** Warning: result is exact up to order",fehler-1,"in",
1027        maxideal(1)[1],"and",fehlervor-1,"in",maxideal(1)[2],"!";
1028     }}
1029     else {
1030      if ((voice==2) && (printlevel > -1)) {
1031        "// ** Warning: result is exact up to order",fehler-1,"in",
1032        maxideal(1)[1],"!";
1033     }}
1034   }
1035   if (return_list==0) { return(ideal(z(1),z(2))); }
1036   else   { return(list(ideal(z(1),z(2)),intvec(fehler-1,fehlervor-1))); }
1037 }
1038}
1039example
1040{ "EXAMPLE:"; echo = 2;
1041 ring exring=0,(x,y,t),ds;
1042 poly f=x3+2xy2+y2;
1043 list hne=develop(f);
1044 list hne_extended1=develop(f,6);
1045 list hne_extended2=develop(f,10);
1046 pause("Hit RETURN to continue:");
1047                                   //   compare the different matrices ...
1048 print(hne[1]);
1049 print(hne_extended1[1]);
1050 print(hne_extended2[1]);
1051                                   // ... and the resulting parametrizations:
1052 param(hne);
1053 param(hne_extended1);
1054 param(hne_extended2);
1055 param(hne_extended2,0);
1056}
1057///////////////////////////////////////////////////////////////////////////////
1058
1059proc invariants
1060"USAGE:   invariants(l) takes the output of develop(f)
1061         (list l (matrix m, intvec v, int s[,poly g,...]))
1062         and computes the following irreducible curve invariants:
1063RETURN:  a list, if l contains a valid HNE:
1064      @format
1065    - invariants(l)[1]: intvec         (characteristic exponents)
1066    - invariants(l)[2]: intvec         (generators of the semigroup)
1067    - invariants(l)[3],[4]: 2 x intvec (Puiseux pairs, 1st and 2nd components)
1068    - invariants(l)[5]: int            (degree of the conductor)
1069    - invariants(l)[6]: intvec         (sequence of multiplicities)
1070      @end format
1071         an empty list, if an error occurred in the procedure develop.
1072SEE ALSO: develop, displayInvariants, multsequence, intersection
1073KEYWORDS: characteristic exponents; semigroup of values; Puiseux pairs;
1074          conductor, degree; multiplicities, sequence of
1075EXAMPLE:  example invariants; shows an example
1076"
1077{
1078 //-------------------------- Initialisierungen -------------------------------
1079 matrix m=#[1];
1080 intvec v=#[2];
1081 int switch=#[3];
1082 list ergebnis;
1083 if (switch==-1) {
1084   "An error has occurred in develop, so there is no HNE.";
1085   return(ergebnis);
1086 }
1087 intvec beta,s,svorl,ordnung,multseq,mpuiseux,npuiseux,halbgr;
1088 int genus,zeile,i,j,k,summe,conductor,ggT;
1089 string Ausgabe;
1090 int nc=ncols(m); int nr=nrows(m);
1091 ordnung[nr]=1;
1092         // alle Indizes muessen (gegenueber [Ca]) um 1 erhoeht werden,
1093         // weil 0..r nicht als Wertebereich erlaubt ist (aber nrows(m)==r+1)
1094
1095 //---------------- Bestimme den Untergrad der einzelnen Zeilen ---------------
1096 for (zeile=nr; zeile>1; zeile--) {
1097   if ((size(ideal(m[zeile,1..nc])) > 1) or (zeile==nr)) { // keine Nullzeile
1098      k=1;
1099      while (m[zeile,k]==0) {k++;}
1100      ordnung[zeile-1]=k*ordnung[zeile]; // vgl. auch Def. von untergrad in
1101      genus++;                           // proc param
1102      svorl[genus]=zeile;} // werden gerade in umgekehrter Reihenfolge abgelegt
1103   else {
1104      ordnung[zeile-1]=v[zeile]*ordnung[zeile]+ordnung[zeile+1];
1105 }}
1106 //----------------- charakteristische Exponenten (beta) ----------------------
1107 s[1]=1;
1108 for (k=1; k <= genus; k++) { s[k+1]=svorl[genus-k+1];} // s[2]==s(1), u.s.w.
1109 beta[1]=ordnung[1]; //charakt. Exponenten: Index wieder verschoben
1110 for (k=1; k <= genus; k++) {
1111   summe=0;
1112   for (i=1; i <= s[k]; i++) {summe=summe+v[i]*ordnung[i];}
1113   beta[k+1]=summe+ordnung[s[k]]+ordnung[s[k]+1]-ordnung[1];
1114 }
1115 //--------------------------- Puiseuxpaare -----------------------------------
1116 int produkt=1;
1117 for (i=1; i<=genus; i++) {
1118   ggT=gcd(beta[1],beta[i+1]*produkt);
1119   mpuiseux[i]=beta[i+1]*produkt / ggT;
1120   npuiseux[i]=beta[1] / ggT;
1121   produkt=produkt*npuiseux[i];
1122 }
1123 //---------------------- Grad des Konduktors ---------------------------------
1124 summe=1-ordnung[1];
1125 if (genus > 0) {
1126   for (i=2; i <= genus+1; i++) {
1127     summe=summe + beta[i] * (ordnung[s[i-1]] - ordnung[s[i]]);
1128   }                              // n.b.: Indizierung wieder um 1 verschoben
1129 }
1130 conductor=summe;
1131 //------------------- Erzeuger der Halbgruppe: -------------------------------
1132 halbgr=puiseux2generators(mpuiseux,npuiseux);
1133
1134 //------------------- Multiplizitaetensequenz: -------------------------------
1135 k=1;
1136 for (i=1; i<size(v); i++) {
1137   for (j=1; j<=v[i]; j++) {
1138     multseq[k]=ordnung[i];
1139     k++;
1140 }}
1141 multseq[k]=1;
1142 //------------------------- Rueckgabe ----------------------------------------
1143 ergebnis=beta,halbgr,mpuiseux,npuiseux,conductor,multseq;
1144 return(ergebnis);
1145}
1146example
1147{ "EXAMPLE:"; echo = 2;
1148 ring exring=0,(x,y),dp;
1149 list hne=develop(y4+2x3y2+x6+x5y);
1150 list erg=invariants(hne);
1151 erg[1];                   // the characteristic exponents
1152 erg[2];                   // the generators of the semigroup of values
1153 erg[3],erg[4];            // the Puiseux pairs in packed form
1154 erg[5] / 2;               // the delta-invariant
1155 erg[6];                   // the sequence of multiplicities
1156                           // To display the invariants more 'nicely':
1157 displayInvariants(hne);
1158}
1159///////////////////////////////////////////////////////////////////////////////
1160
1161proc displayInvariants
1162"USAGE:   displayInvariants(l); takes the output both of
1163         develop(f) and reddevelop(f)
1164          ( list l (matrix m, intvec v, int s[,poly g,...])
1165            or list of lists in the form l )
1166RETURN:  nothing
1167DISPLAY: invariants of the corresponding branch, resp. of all branches,
1168         in a better readable form.
1169SEE ALSO: invariants, intersection, develop, reddevelop
1170EXAMPLE: example displayInvariants;  shows an example
1171"
1172{
1173 int i,j,k,mul;
1174 string Ausgabe;
1175 list ergebnis;
1176 //-- entferne ueberfluessige Daten zur Erhoehung der Rechengeschwindigkeit: --
1177 #=stripHNE(#);
1178 //-------------------- Ausgabe eines Zweiges ---------------------------------
1179 if (typeof(#[1])=="matrix") {
1180   ergebnis=invariants(#);
1181   if (size(ergebnis)!=0) {
1182    " characteristic exponents  :",ergebnis[1];
1183    " generators of semigroup   :",ergebnis[2];
1184    if (size(ergebnis[1])>1) {
1185     for (i=1; i<=size(ergebnis[3]); i++) {
1186       Ausgabe=Ausgabe+"("+string(ergebnis[3][i])+","
1187       +string(ergebnis[4][i])+")";
1188    }}
1189    " Puiseux pairs             :",Ausgabe;
1190    " degree of the conductor   :",ergebnis[5];
1191    " delta invariant           :",ergebnis[5]/2;
1192    " sequence of multiplicities:",ergebnis[6];
1193 }}
1194 //-------------------- Ausgabe aller Zweige ----------------------------------
1195 else {
1196  for (j=1; j<=size(#); j++) {
1197    ergebnis=invariants(#[j]);
1198    " --- invariants of branch number",j,": ---";
1199    " characteristic exponents  :",ergebnis[1];
1200    " generators of semigroup   :",ergebnis[2];
1201    Ausgabe="";
1202    if (size(ergebnis[1])>1) {
1203     for (i=1; i<=size(ergebnis[3]); i++) {
1204       Ausgabe=Ausgabe+"("+string(ergebnis[3][i])+","
1205       +string(ergebnis[4][i])+")";
1206    }}
1207    " Puiseux pairs             :",Ausgabe;
1208    " degree of the conductor   :",ergebnis[5];
1209    " delta invariant           :",ergebnis[5]/2;
1210    " sequence of multiplicities:",ergebnis[6];
1211    "";
1212  }
1213  if (size(#)>1) {
1214    " -------------- intersection multiplicities : -------------- ";"";
1215    Ausgabe="branch |   ";
1216    for (j=size(#); j>1; j--) {
1217      if (size(string(j))==1) { Ausgabe=Ausgabe+" "+string(j)+"    "; }
1218      else                    { Ausgabe=Ausgabe+string(j)+"    "; }
1219    }
1220    Ausgabe;
1221    Ausgabe="-------+";
1222    for (j=2; j<size(#); j++) { Ausgabe=Ausgabe+"------"; }
1223    Ausgabe=Ausgabe+"-----";
1224    Ausgabe;
1225  }
1226  for (j=1; j<size(#); j++) {
1227    if (size(string(j))==1) { Ausgabe="    "+string(j)+"  |"; }
1228    else                    { Ausgabe="   " +string(j)+"  |"; }
1229    for (k=size(#); k>j; k--) {
1230      mul=intersection(#[j],#[k]);
1231      for (i=1; i<=5-size(string(mul)); i++) { Ausgabe=Ausgabe+" "; }
1232      Ausgabe=Ausgabe+string(mul);
1233      if (k>j+1) { Ausgabe=Ausgabe+","; }
1234    }
1235    Ausgabe;
1236  }
1237 }
1238 return();
1239}
1240example
1241{ "EXAMPLE:"; echo = 2;
1242 ring exring=0,(x,y),dp;
1243 list hne=develop(y4+2x3y2+x6+x5y);
1244 displayInvariants(hne);
1245}
1246///////////////////////////////////////////////////////////////////////////////
1247
1248proc multiplicities
1249"USAGE:   multiplicities(l) takes the output of develop(f)
1250         (list l (matrix m, intvec v, int s[,poly g,...]))
1251RETURN:  intvec of the different multiplicities that occur during the
1252         successive blowing up of the curve corresponding to f
1253SEE ALSO: develop
1254EXAMPLE: example multiplicities;  shows an example
1255"
1256{
1257 matrix m=#[1];
1258 intvec v=#[2];
1259 int switch=#[3];
1260 list ergebnis;
1261 if (switch==-1) {
1262   "An error has occurred in develop, so there is no HNE.";
1263   return(intvec(0));
1264 }
1265 intvec ordnung;
1266 int zeile,k;
1267 int nc=ncols(m); int nr=nrows(m);
1268 ordnung[nr]=1;
1269 //---------------- Bestimme den Untergrad der einzelnen Zeilen ---------------
1270 for (zeile=nr; zeile>1; zeile--) {
1271   if ((size(ideal(m[zeile,1..nc])) > 1) or (zeile==nr)) { // keine Nullzeile
1272      k=1;
1273      while (m[zeile,k]==0) {k++;}
1274      ordnung[zeile-1]=k*ordnung[zeile];
1275   }
1276   else {
1277      ordnung[zeile-1]=v[zeile]*ordnung[zeile]+ordnung[zeile+1];
1278 }}
1279 return(ordnung);
1280}
1281example
1282{ "EXAMPLE:"; echo = 2;
1283  int p=printlevel; printlevel=-1;
1284  ring r=0,(x,y),dp;
1285  multiplicities(develop(x5+y7));
1286  // The first value is the multiplicity of the curve itself, here it's 5
1287  printlevel=p;
1288}
1289///////////////////////////////////////////////////////////////////////////////
1290
1291proc puiseux2generators (intvec m, intvec n)
1292"USAGE:   puiseux2generators(m,n);
1293         m,n intvecs with 1st resp. 2nd components of Puiseux pairs
1294RETURN:  intvec of the generators of the semigroup of values
1295SEE ALSO: invariants
1296EXAMPLE: example puiseux2generators;  shows an example
1297"
1298{
1299 intvec beta;
1300 int q=1;
1301 //------------ glatte Kurve (eigentl. waeren m,n leer): ----------------------
1302 if (m==0) { return(intvec(1)); }
1303 //------------------- singulaere Kurve: --------------------------------------
1304 for (int i=1; i<=size(n); i++) { q=q*n[i]; }
1305 beta[1]=q; // == q_0
1306 m=1,m; n=1,n; // m[1] ist damit m_0 usw., genau wie beta[1]==beta_0
1307 for (i=2; i<=size(n); i++) {
1308  beta[i]=m[i]*q/n[i] - m[i-1]*q + n[i-1]*beta[i-1];
1309  q=q/n[i]; // == q_i
1310 }
1311 return(beta);
1312}
1313example
1314{ "EXAMPLE:"; echo = 2;
1315  // take (3,2),(7,2),(15,2),(31,2),(63,2),(127,2) as Puiseux pairs:
1316  puiseux2generators(intvec(3,7,15,31,63,127),intvec(2,2,2,2,2,2));
1317}
1318///////////////////////////////////////////////////////////////////////////////
1319
1320proc intersection (list hn1, list hn2)
1321"USAGE:   intersection(hne1,hne2);
1322         hne1,hne2: two lists representing a HNE (normally two entries out of
1323         the output of reddevelop), i.e. list(matrix,intvec,int[,poly]).
1324RETURN:  The intersection multiplicity of the branches corresponding to
1325         hne1 & hne2 (of type int).
1326SEE ALSO: reddevelop, displayInvariants
1327KEYWORDS: intersection multiplicity
1328EXAMPLE: example intersection;  shows an example
1329"
1330{
1331 //------------------ `intersect' ist schon reserviert ... --------------------
1332 int i,j,s,sum,schnitt,unterschied;
1333 matrix a1=hn1[1];
1334 matrix a2=hn2[1];
1335 intvec h1=hn1[2];
1336 intvec h2=hn2[2];
1337 intvec n1=multiplicities(hn1);
1338 intvec n2=multiplicities(hn2);
1339 if (hn1[3]!=hn2[3]) {
1340 //-- die jeweils erste Zeile von hn1,hn2 gehoert zu verschiedenen Parametern -
1341 //---------------- d.h. beide Kurven schneiden sich transversal --------------
1342   schnitt=n1[1]*n2[1];        // = mult(hn1)*mult(hn2)
1343 }
1344 else {
1345 //--------- die jeweils erste Zeile gehoert zum gleichen Parameter -----------
1346   unterschied=0;
1347   for (s=1; (h1[s]==h2[s]) && (s<size(h1)) && (s<size(h2))
1348              && (unterschied==0); s++) {
1349     for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]); i++) {;}
1350     if (i<=h1[s]) {
1351       unterschied=1;
1352       s--;                    // um s++ am Schleifenende wieder auszugleichen
1353     }
1354   }
1355   if (unterschied==0) {
1356     if ((s<size(h1)) && (s<size(h2))) {
1357       for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]) && (i<=h2[s]); i++) {;}
1358     }
1359     else {
1360 //-------------- Sonderfall: Unterschied in letzter Zeile suchen -------------
1361 // Beachte: Es koennen undefinierte Stellen auftreten, bei abbrechender HNE
1362 // muss die Ende-Markierung weg, h_[r] ist unendlich, die Matrix muss mit
1363 // Nullen fortgesetzt gedacht werden
1364 //----------------------------------------------------------------------------
1365       if (ncols(a1)>ncols(a2)) { j=ncols(a1); }
1366       else                     { j=ncols(a2); }
1367       unterschied=0;
1368       if ((h1[s]>0) && (s==size(h1))) {
1369         a1[s,h1[s]+1]=0;
1370         if (ncols(a1)<=ncols(a2)) { unterschied=1; }
1371       }
1372       if ((h2[s]>0) && (s==size(h2))) {
1373         a2[s,h2[s]+1]=0;
1374         if (ncols(a2)<=ncols(a1)) { unterschied=1; }
1375       }
1376       if (unterschied==1) {                   // mind. eine HNE war endlich
1377         matrix ma1[1][j]=a1[s,1..ncols(a1)];  // und bedarf der Fortsetzung
1378         matrix ma2[1][j]=a2[s,1..ncols(a2)];  // mit Nullen
1379       }
1380       else {
1381         if (ncols(a1)>ncols(a2)) { j=ncols(a2); }
1382         else                     { j=ncols(a1); }
1383         matrix ma1[1][j]=a1[s,1..j];          // Beschr. auf vergleichbaren
1384         matrix ma2[1][j]=a2[s,1..j];          // Teil (der evtl. y's enth.)
1385       }
1386       for (i=1; (ma1[1,i]==ma2[1,i]) && (i<j) && (ma1[1,i]!=var(2)); i++) {;}
1387       if (ma1[1,i]==ma2[1,i]) {
1388         "//** The two HNE's are identical!";
1389         "//** You have either tried to intersect a branch with itself,";
1390         "//** or the two branches have been developed separately.";
1391         "//   In the latter case use `extdevelop' to extend the HNE's until",
1392         "they differ.";
1393         return(-1);
1394       }
1395       if ((ma1[1,i]==var(2)) || (ma2[1,i]==var(2))) {
1396         "//** The two HNE's are (so far) identical. This is because they",
1397         "have been";
1398         "//** computed separately. I need more data; use `extdevelop' to",
1399         "extend them,";
1400         if (ma1[1,i]==var(2)) {"//** at least the first one.";}
1401         else                  {"//** at least the second one.";}
1402         return(-1);
1403       }
1404     }
1405   }
1406   sum=0;
1407   h1[size(h1)]=ncols(a1)+42;        // Ersatz fuer h1[r]=infinity
1408   h2[size(h2)]=ncols(a2)+42;
1409   for (j=1; j<s; j++) {sum=sum+h1[j]*n1[j]*n2[j];}
1410   if ((i<=h1[s]) && (i<=h2[s]))    { schnitt=sum+i*n1[s]*n2[s]; }
1411   if (i==h2[s]+1) { schnitt=sum+h2[s]*n1[s]*n2[s]+n2[s+1]*n1[s]; }
1412   if (i==h1[s]+1) { schnitt=sum+h1[s]*n2[s]*n1[s]+n1[s+1]*n2[s]; }
1413 // "s:",s-1,"i:",i,"S:",sum;
1414 }
1415 return(schnitt);
1416}
1417example
1418{
1419  // -------- prepare for example ---------
1420  if (nameof(basering)=="HNEring") {
1421   def rettering=HNEring;
1422   kill HNEring;
1423  }
1424  // ------ the example starts here -------
1425  "EXAMPLE:"; echo = 2;
1426  ring r=0,(x,y),dp;
1427  list hne=reddevelop((x2-y3)*(x2+y3));
1428  intersection(hne[1],hne[2]);
1429  kill HNEring,r;
1430  echo = 0;
1431  // --- restore HNEring if previously defined ---
1432  if (defined(rettering)) {
1433   setring rettering;
1434   def HNEring=rettering;
1435   export HNEring;
1436  }
1437}
1438///////////////////////////////////////////////////////////////////////////////
1439
1440proc multsequence
1441"USAGE:   multsequence(l);
1442         l is the output of either develop(f) or reddevelop(f)
1443RETURN:  - if l=@code{develop}(f) or l=@code{reddevelop}(f)[i]:
1444           intvec of the sequence of multiplicities of the curve
1445           resp. of the i-th branch (the same as @code{invariants}(l)[6]; )
1446
1447         - if l=@code{reddevelop}(f) : list of two integer matrices:
1448  @* @code{multsequence}(l)[1][i,*] contains the multiplicities of the branches
1449           at their infinitely near point of 0 in its (i-1) order
1450           neighbourhood (i.e. i=1: multiplicity of the branches themselves,
1451                 i=2: multiplicity of their 1st quadratic transformed etc.,
1452       @*  or to say it in another way: @code{multsequence}(l)[1][*,j] is the
1453                 sequence of multiplicities of the j-th branch).
1454  @* @code{multsequence}(l)[2][i,*] contains the information which of these
1455           infinitely near points coincide.
1456NOTE:  The order of elements of the list obtained from @code{reddevelop} must
1457       not be changed (because then the coincident infinitely near points
1458       couldn't be grouped together, cf. meaning of 2nd intmat in example).
1459       Hence it is not wise to compute the HNE of several polynomials
1460       separately, put them into a list l and call @code{multsequence}(l).
1461
1462       Use @code{displayMultsequence} to produce a better readable output for
1463       reducible curves on the screen.
1464SEE ALSO: displayMultsequence, develop, reddevelop, separateHNE
1465EXAMPLE: example multsequence;  shows an example
1466"
1467{
1468 //-- entferne ueberfluessige Daten zur Erhoehung der Rechengeschwindigkeit: --
1469 #=stripHNE(#);
1470 int k,i,j;
1471 //----------------- Multiplizitaetensequenz eines Zweiges --------------------
1472 if (typeof(#[1])=="matrix") {
1473  intvec v=#[2];
1474  list ergebnis;
1475  if (#[3]==-1) {
1476    "An error has occurred in develop, so there is no HNE.";
1477   return(intvec(0));
1478  }
1479  intvec multips,multseq;
1480  multips=multiplicities(#);
1481  k=1;
1482  for (i=1; i<size(v); i++) {
1483    for (j=1; j<=v[i]; j++) {
1484      multseq[k]=multips[i];
1485      k++;
1486  }}
1487  multseq[k]=1;
1488  return(multseq);
1489 }
1490 //---------------------------- mehrere Zweige --------------------------------
1491 else {
1492   list HNEs=#;
1493   int anzahl=size(HNEs);
1494   int maxlength=0;
1495   int bisher;
1496   intvec schnitt,ones;
1497   ones[anzahl]=0;
1498   ones=ones+1;                  // = 1,1,...,1
1499   for (i=1; i<anzahl; i++) {
1500     schnitt[i]=separateHNE(HNEs[i],HNEs[i+1]);
1501     j=size(multsequence(HNEs[i]));
1502     maxlength=maxlength*(j < maxlength) + j*(j >= maxlength);
1503     maxlength=maxlength*(schnitt[i]+1 < maxlength)
1504               + (schnitt[i]+1)*(schnitt[i]+1 >= maxlength);
1505   }
1506   j=size(multsequence(HNEs[anzahl]));
1507   maxlength=maxlength*(j < maxlength) + j*(j >= maxlength);
1508
1509//-------------- Konstruktion der ersten zu berechnenden Matrix ---------------
1510   intmat allmults[maxlength][anzahl];
1511   for (i=1; i<=maxlength; i++)  { allmults[i,1..anzahl]=ones[1..anzahl]; }
1512   for (i=1; i<=anzahl; i++) {
1513     ones=multsequence(HNEs[i]);
1514     allmults[1..size(ones),i]=ones[1..size(ones)];
1515   }
1516//---------------------- Konstruktion der zweiten Matrix ----------------------
1517   intmat separate[maxlength][anzahl];
1518   for (i=1; i<=maxlength; i++) {
1519     k=1;
1520     bisher=0;
1521     if (anzahl==1) { separate[i,1]=1; }
1522     for (j=1; j<anzahl; j++)   {
1523       if (schnitt[j]<i) {
1524         separate[i,k]=j-bisher;
1525         bisher=j;
1526         k++;
1527       }
1528       separate[i,k]=anzahl-bisher;
1529     }
1530   }
1531  return(list(allmults,separate));
1532 }
1533}
1534example
1535{
1536  // -------- prepare for example ---------
1537  if (nameof(basering)=="HNEring") {
1538   def rettering=HNEring;
1539   kill HNEring;
1540  }
1541  // ------ the example starts here -------
1542  "EXAMPLE:"; echo = 2;
1543  ring r=0,(x,y),dp;
1544  list hne=reddevelop((x6-y10)*(x+y2-y3)*(x+y2+y3));
1545
1546  multsequence(hne[1]),"  |  ",multsequence(hne[2]),"  |  ",
1547  multsequence(hne[3]),"  |  ",multsequence(hne[4]);
1548  multsequence(hne);
1549  // The meaning of the entries of the 2nd matrix is as follows:
1550  displayMultsequence(hne);
1551  kill HNEring,r;
1552  echo = 0;
1553  // --- restore HNEring if previously defined ---
1554  if (defined(rettering)) {
1555   setring rettering;
1556   def HNEring=rettering;
1557   export HNEring;
1558  }
1559}
1560///////////////////////////////////////////////////////////////////////////////
1561
1562proc displayMultsequence
1563"USAGE:   displayMultsequence(l);
1564         l is the output of either develop(f) or reddevelop(f)
1565RETURN:  nothing
1566DISPLAY: the sequence of multiplicities:
1567       if l=develop(f) or l=reddevelop(f)[i] : just the sequence a,b,c,...,1
1568 @*    if l=reddevelop(f) : a sequence of the following form:
1569
1570 @*    @math{[(a_1,....,b_1,....,c_1)]},
1571 @*    @math{[(a_2,...),...,(...,c_2)]},
1572 @*    @math{   ...................}
1573 @*    @math{[(a_n),(b_n),.....,(c_n)]}
1574
1575 @*    with @math{a_1,...,a_n} the sequence of multiplicities of the 1st branch
1576 @*    [...] the multiplicities of the j-th transformed of all branches
1577 @*    (...) the multiplicities of those branches meeting in the same
1578             infinitely near point.
1579NOTE:    The same restrictions for l as in `@code{multsequence}' apply.
1580SEE ALSO: multsequence, develop, reddevelop, separateHNE
1581EXAMPLE: example displayMultsequence;  shows an example
1582"
1583{
1584 //-- entferne ueberfluessige Daten zur Erhoehung der Rechengeschwindigkeit: --
1585 #=stripHNE(#);
1586 //----------------- Multiplizitaetensequenz eines Zweiges --------------------
1587 if (typeof(#[1])=="matrix") {
1588   if (#[3]==-1) {
1589     "An error has occurred in develop, so there is no HNE.";
1590   }
1591   else {
1592     "The sequence of multiplicities is  ",multiplicities(#);
1593 }}
1594 //---------------------------- mehrere Zweige --------------------------------
1595 else {
1596   list multips=multsequence(#);
1597   int i,j,k,l;
1598   string output;
1599   for (i=1; i<=nrows(multips[1]); i++) {
1600     output="[";
1601     k=1;
1602     for (l=1; k<=ncols(multips[1]); l++) {
1603       output=output+"(";
1604       for (j=1; j<=multips[2][i,l]; j++) {
1605         output=output+string(multips[1][i,k]);
1606         k++;
1607         if (j<multips[2][i,l]) { output=output+","; }
1608       }
1609       output=output+")";
1610       if ((k-1) < ncols(multips[1])) { output=output+","; }
1611      }
1612     output=output+"]";
1613     if (i<nrows(multips[1])) { output=output+","; }
1614     output;
1615   }
1616 }
1617} // example multsequence; geht wegen echo nicht (muesste auf 3 gesetzt werden)
1618example
1619{
1620  // -------- prepare for example ---------
1621  if (nameof(basering)=="HNEring") {
1622   def rettering=HNEring;
1623   kill HNEring;
1624  }
1625  // ------ the example starts here -------
1626  "EXAMPLE:"; echo = 2;
1627  ring r=0,(x,y),dp;
1628  list hne=reddevelop((x6-y10)*(x+y2-y3)*(x+y2+y3));
1629
1630  displayMultsequence(hne[1]);
1631  displayMultsequence(hne);
1632  kill HNEring,r;
1633  echo = 0;
1634  // --- restore HNEring if previously defined ---
1635  if (defined(rettering)) {
1636   setring rettering;
1637   def HNEring=rettering;
1638   export HNEring;
1639  }
1640}
1641///////////////////////////////////////////////////////////////////////////////
1642
1643proc separateHNE (list hn1,list hn2)
1644"USAGE:    separateHNE(hne1,hne2);  list hne1,hne2: output of develop
1645RETURN:   number of quadratic transformations needed to separate both curves
1646SEE ALSO: develop, reddevelop, multsequence, displayMultsequence
1647EXAMPLE:  example separateHNE;  shows an example
1648"
1649{
1650 int i,j,s,unterschied,separated;
1651 matrix a1=hn1[1];
1652 matrix a2=hn2[1];
1653 intvec h1=hn1[2];
1654 intvec h2=hn2[2];
1655 if (hn1[3]!=hn2[3]) {
1656 //-- die jeweils erste Zeile von hn1,hn2 gehoert zu verschiedenen Parametern -
1657 //---------------- d.h. beide Kurven schneiden sich transversal --------------
1658   separated=1;
1659 }
1660 else {
1661 //--------- die jeweils erste Zeile gehoert zum gleichen Parameter -----------
1662   unterschied=0;
1663   for (s=1; (h1[s]==h2[s]) && (s<size(h1)) && (s<size(h2))
1664              && (unterschied==0); s++) {
1665     for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]); i++) {;}
1666     if (i<=h1[s]) {
1667       unterschied=1;
1668       s--;                    // um s++ am Schleifenende wieder auszugleichen
1669     }
1670   }
1671   if (unterschied==0) {
1672     if ((s<size(h1)) && (s<size(h2))) {
1673       for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]) && (i<=h2[s]); i++) {;}
1674     }
1675     else {
1676 //-------------- Sonderfall: Unterschied in letzter Zeile suchen -------------
1677 // Beachte: Es koennen undefinierte Stellen auftreten, bei abbrechender HNE
1678 // muss die Ende-Markierung weg, h_[r] ist unendlich, die Matrix muss mit
1679 // Nullen fortgesetzt gedacht werden
1680 //----------------------------------------------------------------------------
1681       if (ncols(a1)>ncols(a2)) { j=ncols(a1); }
1682       else                     { j=ncols(a2); }
1683       unterschied=0;
1684       if ((h1[s]>0) && (s==size(h1))) {
1685         a1[s,h1[s]+1]=0;
1686         if (ncols(a1)<=ncols(a2)) { unterschied=1; }
1687       }
1688       if ((h2[s]>0) && (s==size(h2))) {
1689         a2[s,h2[s]+1]=0;
1690         if (ncols(a2)<=ncols(a1)) { unterschied=1; }
1691       }
1692       if (unterschied==1) {                   // mind. eine HNE war endlich
1693         matrix ma1[1][j]=a1[s,1..ncols(a1)];  // und bedarf der Fortsetzung
1694         matrix ma2[1][j]=a2[s,1..ncols(a2)];  // mit Nullen
1695       }
1696       else {
1697         if (ncols(a1)>ncols(a2)) { j=ncols(a2); }
1698         else                     { j=ncols(a1); }
1699         matrix ma1[1][j]=a1[s,1..j];          // Beschr. auf vergleichbaren
1700         matrix ma2[1][j]=a2[s,1..j];          // Teil (der evtl. y's enth.)
1701       }
1702       for (i=1; (ma1[1,i]==ma2[1,i]) && (i<j) && (ma1[1,i]!=var(2)); i++) {;}
1703       if (ma1[1,i]==ma2[1,i]) {
1704         "//** The two HNE's are identical!";
1705         "//** You have either tried to compare a branch with itself,";
1706         "//** or the two branches have been developed separately.";
1707         "//   In the latter case use `extdevelop' to extend the HNE's until",
1708         "they differ.";
1709         return(-1);
1710       }
1711       if ((ma1[1,i]==var(2)) || (ma2[1,i]==var(2))) {
1712         "//** The two HNE's are (so far) identical. This is because they",
1713         "have been";
1714         "//** computed separately. I need more data; use `extdevelop' to",
1715         "extend them,";
1716         if (ma1[1,i]==var(2)) {"//** at least the first one.";}
1717         else                  {"//** at least the second one.";}
1718         return(-1);
1719       }
1720     }
1721   }
1722   separated=i;
1723   for (j=1; j<s; j++) { separated=separated+h1[j]; }
1724 }
1725 return(separated);
1726}
1727example
1728{ "EXAMPLE:"; echo = 2;
1729  int p=printlevel; printlevel=-1;
1730  ring r=0,(x,y),dp;
1731  list hne1=develop(x);
1732  list hne2=develop(x+y);
1733  list hne3=develop(x+y2);
1734  separateHNE(hne1,hne2);  // two transversal lines
1735  separateHNE(hne1,hne3);  // one quadratic transform. gives 1st example
1736  printlevel=p;
1737}
1738///////////////////////////////////////////////////////////////////////////////
1739
1740proc displayHNE(list ldev,list #)
1741"USAGE:   displayHNE(ldev,[,n]); ldev=list (the output of develop(f) or
1742          reddevelop(f)), n=int
1743RETURN:  - if only one argument is given, no return value, but
1744         display an ideal HNE of the following form:
1745     @example
1746     HNE[1]=-y+[]*z(0)^1+[]*z(0)^2+...+z(0)^<>*z(1)
1747     HNE[2]=-x+          []*z(1)^2+...+z(1)^<>*z(2)
1748     HNE[3]=             []*z(2)^2+...+z(2)^<>*z(3)
1749     .......             ..........................
1750     HNE[r+1]=           []*z(r)^2+[]*z(r)^3+......
1751     @end example
1752     where x,y are the indeterminates of the basering. The values of @code{[]}
1753     are the coefficients of the Hamburger-Noether matrix, the values of
1754     @code{<>} are represented in the HN-matrix as `x'.
1755     @*  The 1st line (@code{HNE[1]}) means that
1756     @math{y=}@code{[]}@math{*z(0)^1+...},
1757       the 2nd line (@code{HNE[2]}) means that
1758     @math{x=}@code{[]}@math{*z(1)^2+...},
1759     so you can see which indeterminate corresponds to which line
1760     (it's also possible that x corresponds to the 1st line and y to the 2nd).
1761     @*  - if a second argument is given, create and export a new ring with
1762         name `displayring' containing an ideal `HNE' as described above.
1763
1764     If ldev contains the output of @code{reddevelop}(f), @code{displayHNE}
1765     shows the HNE's of all branches of f in the form described above.
1766     The optional parameter is then ignored.
1767
1768SEE ALSO: develop, reddevelop
1769EXAMPLE: example displayHNE; shows an example
1770"
1771{
1772 if ((typeof(ldev[1])=="list") || (typeof(ldev[1])=="none")) {
1773   for (int i=1; i<=size(ldev); i++) {
1774     "// Hamburger-Noether development of branch nr."+string(i)+":";
1775     displayHNE(ldev[i]);"";
1776   }
1777   return();
1778 }
1779 //--------------------- Initialisierungen und Ringwechsel --------------------
1780 matrix m=ldev[1];
1781 intvec v=ldev[2];
1782 int switch=ldev[3];
1783 if (switch==-1) {
1784   "An error has occurred in develop, so there is no HNE.";
1785   return(ideal(0));
1786 }
1787 def altring=basering;
1788 if (parstr(basering)!="") {
1789   if (charstr(basering)!=string(char(basering))+","+parstr(basering)) {
1790     execute
1791      ("ring dazu=("+charstr(basering)+"),z(0.."+string(size(v)-1)+"),ls;");
1792   }
1793   else { ring dazu=char(altring),z(0..size(v)-1),ls; }
1794 }
1795 else   { ring dazu=char(altring),z(0..size(v)-1),ls; }
1796 def displayring=dazu+altring;
1797 setring displayring;
1798 if (size(#) != 0) {
1799    export displayring;
1800  }
1801 map holematrix=altring,0;        // mappt nur die Monome vom Grad Null
1802 matrix m=holematrix(m);
1803 //--------------------- Erzeuge Matrix n mit n[i,j]=z(j-1)^i -----------------
1804 int i;
1805 matrix n[ncols(m)][nrows(m)];
1806 for (int j=1; j<=nrows(m); j++) {
1807    for (i=1; i<=ncols(m); i++) { n[i,j]=z(j-1)^i; }
1808 }
1809 matrix displaymat=m*n;
1810 ideal HNE;
1811 for (i=1; i<nrows(m); i++) { HNE[i]=displaymat[i,i]+z(i)*z(i-1)^v[i]; }
1812 HNE[nrows(m)]=displaymat[nrows(m),nrows(m)];
1813 if (nrows(m)<2) { HNE[2]=z(0); }
1814 if (switch==0) {
1815    HNE[1] = HNE[1]-var(size(v)+2);
1816    HNE[2] = HNE[2]-var(size(v)+1);
1817 }
1818 else {
1819    HNE[1] = HNE[1]-var(size(v)+1);
1820    HNE[2] = HNE[2]-var(size(v)+2);
1821 }
1822if (size(#) == 0) {
1823   HNE;
1824   return();
1825 }
1826if (size(#) != 0) {
1827   "// basering is now 'displayring' containing ideal 'HNE'";
1828   keepring(displayring);
1829   export(HNE);
1830   return(HNE);
1831 }
1832}
1833example
1834{ "EXAMPLE:"; echo = 2;
1835  ring r=0,(x,y),dp;
1836  poly f=x3+2xy2+y2;
1837  list hn=develop(f);
1838  displayHNE(hn);
1839}
1840///////////////////////////////////////////////////////////////////////////////
1841//                      procedures for reducible curves                      //
1842///////////////////////////////////////////////////////////////////////////////
1843
1844// proc newtonhoehne (poly f)
1845// USAGE:   newtonhoehne(f);   f poly
1846// ASSUME:  basering = ...,(x,y),ds  or ls
1847// RETURN:  list of intvec(x,y) of coordinates of the newtonpolygon of f
1848// NOTE:    This proc is only available in versions of Singular that know the
1849//         command  system("newton",f);  f poly
1850// {
1851// intvec nm = getnm(f);
1852//  if ((nm[1]>0) && (nm[2]>0)) { f=jet(f,nm[1]*nm[2],nm); }
1853//  list erg=system("newton",f);
1854//  int i; list Ausgabe;
1855//  for (i=1; i<=size(erg); i++) { Ausgabe[i]=leadexp(erg[i]); }
1856// return(Ausgabe);
1857// }
1858///////////////////////////////////////////////////////////////////////////////
1859
1860proc newtonpoly (poly f)
1861"USAGE:   newtonpoly(f);   f poly
1862RETURN:  list of intvec(x,y) of coordinates of the Newton polygon of f
1863ASSUME: For performance reasons, newtonpoly assumes that
1864      - the basering has ordering `ls'
1865      - f(x,0) != 0 != f(0,y), f(0,0) = 0
1866EXAMPLE: example newtonpoly;  shows an example
1867"
1868{
1869 intvec A=(0,ord(subst(f,var(1),0)));
1870 intvec B=(ord(subst(f,var(2),0)),0);
1871 intvec C,H; list L;
1872 int abbruch,i;
1873 poly hilf;
1874
1875 L[1]=A;
1876 //-------- wirf alle Monome auf oder oberhalb der Geraden AB raus: -----------
1877 f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]));
1878 map xytausch=basering,var(2),var(1);
1879 for (i=2; f!=0; i++) {
1880   abbruch=0;
1881   while (abbruch==0) {
1882 // finde den Punkt aus {verbliebene Pkte (a,b) mit a minimal} mit b minimal: -
1883
1884     C=leadexp(f);           // Ordnung ls ist wesentlich!
1885
1886     if (jet(f,A[2]*C[1]-A[1]*C[2]-1,intvec(A[2]-C[2],C[1]-A[1]))==0)
1887       { abbruch=1; }        // keine Monome unterhalb der Geraden AC
1888
1889 // ----- alle Monome auf der Parallelen zur y-Achse durch C wegwerfen: -------
1890 // ------------------ (links von C gibt es sowieso keine mehr) ---------------
1891     else { f=jet(f,-C[1]-1,intvec(-1,0)); }
1892   }
1893 //- finde alle Monome auf der Geraden durch A und C (unterhalb gibt's keine) -
1894   hilf=jet(f,A[2]*C[1]-A[1]*C[2],intvec(A[2]-C[2],C[1]-A[1]));
1895
1896   H=leadexp(xytausch(hilf));
1897   A=H[2],H[1];
1898
1899 // die Alternative waere ein Ringwechsel nach ..,(y,x),ds gewesen
1900 // A_neu ist der naechste Eckpunkt (unterster Punkt auf Geraden durch A,C)
1901
1902   L[i]=A;
1903 //----------------- alle Monome auf oder unterhalb AB raus -------------------
1904   f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]-A[1]));
1905 }
1906 L[i]=B;
1907 return(L);
1908}
1909example
1910{ "EXAMPLE:";
1911  ring @exring_Newt=0,(x,y),ls;
1912  export @exring_Newt;
1913  "  ring exring=0,(x,y),ls;";
1914  echo = 2;
1915  poly f=x5+2x3y-x2y2+3xy5+y6-y7;
1916  newtonpoly(f);
1917  echo = 0;
1918  kill @exring_Newt;
1919}
1920///////////////////////////////////////////////////////////////////////////////
1921
1922proc charPoly(poly f, int M, int N)
1923"USAGE:  charPoly(f,M,N);  f bivariate poly,  M,N int: length and height
1924                          of Newton polygon of f, which has to be only one line
1925RETURN:  the characteristic polynomial of f
1926EXAMPLE: example charPoly;  shows an example
1927"
1928{
1929 poly charp;
1930 int Np=N/ gcd(M,N);
1931 f=subst(f,var(1),1);
1932 for(charp=0; f<>0; f=f-lead(f))
1933  { charp=charp+leadcoef(f)*var(2)^(leadexp(f)[2]/ Np);}
1934 return(charp);
1935}
1936example
1937{ "EXAMPLE:"; echo = 2;
1938  ring exring=0,(x,y),dp;
1939  charPoly(y4+2y3x2-yx6+x8,8,4);
1940  charPoly(y6+3y3x2-x4,4,6);
1941}
1942///////////////////////////////////////////////////////////////////////////////
1943
1944proc find_in_list(list L,int p)
1945"USAGE:   find_in_list(L,p); L: list of intvec(x,y)
1946         (sorted in y: L[1][2]>=L[2][2]), int p >= 0
1947RETURN:  int i: L[i][2]=p if existent; otherwise i with L[i][2]<p if existent;
1948         otherwise i = size(L)+1;
1949EXAMPLE: example find_in_list;  shows an example
1950"
1951{
1952 int i;
1953 L[size(L)+1]=intvec(0,-1);          // falls p nicht in L[.][2] vorkommt
1954 for (i=1; L[i][2]>p; i++) {;}
1955 return(i);
1956}
1957example
1958{ "EXAMPLE:"; echo = 2;
1959  list L = intvec(0,4), intvec(1,2), intvec(2,1), intvec(4,0);
1960  find_in_list(L,1);
1961  L[find_in_list(L,2)];
1962}
1963///////////////////////////////////////////////////////////////////////////////
1964
1965proc get_last_divisor(int M, int N)
1966"USAGE:   get_last_divisor(M,N); int M,N
1967RETURN:  int Q: M=q1*N+r1, N=q2*r1+r2, ..., ri=Q*r(i+1) (Euclidean alg.)
1968EXAMPLE: example get_last_divisor; shows an example
1969"
1970{
1971 int R=M%N; int Q=M / N;
1972 while (R!=0) {M=N; N=R; R=M%N; Q=M / N;}
1973 return(Q)
1974}
1975example
1976{ "EXAMPLE"; echo = 2;
1977  ring r=0,(x,y),dp;
1978  get_last_divisor(12,10);
1979}
1980///////////////////////////////////////////////////////////////////////////////
1981proc redleit (poly f,intvec S, intvec E)
1982"USAGE:   redleit(f,S,E);  f poly, S,E intvec(x,y)
1983         S,E are two different points on a line in the Newton diagram of f
1984RETURN:  poly g: all monomials of f which lie on or below that line
1985NOTE:    The main purpose is that if the line defined by S and E is part of the
1986         Newton polygon, the result is the quasihomogeneous leading form of f
1987         wrt. that line.
1988SEE ALSO: newtonpoly
1989EXAMPLE: example redleit;  shows an example
1990"
1991{
1992 if (E[1]<S[1]) { intvec H=E; E=S; S=H; } // S,E verkehrt herum eingegeben
1993 return(jet(f,E[1]*S[2]-E[2]*S[1],intvec(S[2]-E[2],E[1]-S[1])));
1994}
1995example
1996{ "EXAMPLE"; echo = 2;
1997  ring exring=0,(x,y),dp;
1998  redleit(y6+xy4-2x3y2+x4y+x6,intvec(3,2),intvec(4,1));
1999}
2000///////////////////////////////////////////////////////////////////////////////
2001
2002
2003proc extdevelop (list l, int Exaktheit)
2004"USAGE:   extdevelop(l,n);
2005         list l (matrix m, intvec v, int s, poly g[,int b]), int n
2006 @*      takes the output l of develop(f) ( or extdevelop(L,N),
2007         or one entry of the output of reddevelop(f))  and
2008RETURN:  an extension of the Hamburger-Noether development of f in the same
2009         form as l (i.e. list(matrix, intvec, int, poly) ).
2010         The new HN-matrix will have at least n columns
2011         (if the HNE isn't finite).
2012         Thus if f is irreducible, @code{extdevelop}(@code{develop}(f),n);
2013         (in most cases) will produce the same result as @code{develop}(f,n).
2014         Type `help develop;'  for more details.
2015NOTE:    If the matrix m of l has N columns, the exactness of
2016         @code{param}(@code{extdevelop}(l,n))  will be increased by at least
2017         (n-N) more significant monomials compared with @code{param}(l).
2018SEE ALSO: develop, reddevelop, param
2019EXAMPLE: example extdevelop;  shows an example
2020"
2021{
2022 //------------ Initialisierungen und Abfangen unzulaessiger Aufrufe ----------
2023 matrix m=l[1];
2024 intvec v=l[2];
2025 int switch=l[3];
2026 if (nvars(basering) < 2) {
2027   " Sorry. I need two variables in the ring.";
2028   return(list(matrix(maxideal(1)[1]),intvec(0),-1,poly(0)));}
2029 if (switch==-1) {
2030   "An error has occurred in develop, so there is no HNE and no extension.";
2031   return(l);
2032 }
2033 poly f=l[4];
2034 if (f==0) {
2035   " No extension is possible";
2036   return(l);
2037 }
2038 int Q=v[size(v)];
2039 if (Q>0) {
2040   " The HNE was already exact";
2041   return(l);
2042 }
2043 else {
2044   if (Q==-1) { Q=ncols(m); }
2045   else { Q=-Q-1; }
2046 }
2047 int zeile=nrows(m);
2048 int spalten,i,M;
2049 ideal lastrow=m[zeile,1..Q];
2050 int ringwechsel=(varstr(basering)!="x,y") or (ordstr(basering)!="ls(2),C");
2051
2052 //------------------------- Ringwechsel, falls noetig ------------------------
2053 if (ringwechsel) {
2054  def altring = basering;
2055  int p = char(basering);
2056  if (charstr(basering)!=string(p)) {
2057     string tststr=charstr(basering);
2058     tststr=tststr[1..find(tststr,",")-1];     //-> "p^k" bzw. "p"
2059     if (tststr==string(p)) {
2060       if (size(parstr(basering))>1) {         // ring (p,a,..),...
2061        execute("ring extdguenstig=("+charstr(basering)+"),(x,y),ls;");
2062       }
2063       else {                                  // ring (p,a),...
2064        string mipl=string(minpoly);
2065        ring extdguenstig=(p,`parstr(basering)`),(x,y),ls;
2066        if (mipl!="0") { execute("minpoly="+mipl+";"); }
2067       }
2068     }
2069     else {
2070       execute("ring extdguenstig=("+charstr(basering)+"),(x,y),ls;");
2071     }
2072  }
2073  else {                               // charstr(basering)== p : no parameter
2074     ring extdguenstig=p,(x,y),ls;
2075  }
2076  export extdguenstig;
2077  map hole=altring,x,y;
2078 //----- map kann sehr zeitaufwendig sein, daher Vermeidung, wo moeglich: -----
2079  if (nvars(altring)==2) { poly f=fetch(altring,f); }
2080  else                   { poly f=hole(f);          }
2081  ideal a=hole(lastrow);
2082 }
2083 else { ideal a=lastrow; }
2084 list Newton=newtonpoly(f);
2085 int M1=Newton[size(Newton)-1][1];     // konstant
2086 number delta;
2087 if (Newton[size(Newton)-1][2]!=1) {
2088    " *** The transformed polynomial was not valid!!";}
2089 else {
2090 //--------------------- Fortsetzung der HNE ----------------------------------
2091  while (Q<Exaktheit) {
2092    M=ord(subst(f,y,0));
2093    Q=M-M1;
2094 //------ quasihomogene Leitform ist c*x^M1*y+d*x^(M1+Q) => delta=-d/c: -------
2095    delta=-koeff(f,M,0)/koeff(f,M1,1);
2096    a[Q]=delta;
2097    if ((voice==2) && (printlevel > -1))
2098     { "a("+string(zeile-1)+","+string(Q)+") =",delta; }
2099    if (Q<Exaktheit) {
2100     f=T1_Transform(f,delta,Q);
2101     if (defined(HNDebugOn)) { "transformed polynomial:",f; }
2102     if (subst(f,y,0)==0) {
2103       if ((voice==2) && (printlevel > -1)) { "The HNE is finite!"; }
2104       a[Q+1]=x; Exaktheit=Q;
2105       f=0;                        // Speicherersparnis: f nicht mehr gebraucht
2106     }
2107    }
2108  }
2109 }
2110 //------- Wechsel in alten Ring, Zusammensetzung alte HNE + Erweiterung ------
2111 if (ringwechsel) {
2112  setring altring;
2113  map zurueck=extdguenstig,var(1),var(2);
2114  if (nvars(altring)==2) { f=fetch(extdguenstig,f); }
2115  else                   { f=zurueck(f);            }
2116  lastrow=zurueck(a);
2117 }
2118 else { lastrow=a; }
2119 if (ncols(lastrow)>ncols(m)) { spalten=ncols(lastrow); }
2120 else { spalten=ncols(m); }
2121 matrix mneu[zeile][spalten];
2122 for (i=1; i<nrows(m); i++) {
2123  mneu[i,1..ncols(m)]=m[i,1..ncols(m)];
2124 }
2125 mneu[zeile,1..ncols(lastrow)]=lastrow;
2126 if (lastrow[ncols(lastrow)]!=var(1)) {
2127  if (ncols(lastrow)==spalten) { v[zeile]=-1; }  // keine undefinierten Stellen
2128  else {
2129   v[zeile]=-Q-1;
2130   for (i=ncols(lastrow)+1; i<=spalten; i++) {
2131    mneu[zeile,i]=var(2);           // fuelle nicht def. Stellen der Matrix auf
2132 }}}
2133 else { v[zeile]=Q; }               // HNE war exakt
2134 if (ringwechsel)
2135 {
2136   if(system("with","Namespaces")) { kill Top::extdguenstig; }
2137   kill extdguenstig;
2138 }
2139
2140 return(list(mneu,v,switch,f));
2141}
2142example
2143{
2144  // -------- prepare for example ---------
2145  if (nameof(basering)=="HNEring") {
2146   def rettering=HNEring;
2147   kill HNEring;
2148  }
2149  // ------ the example starts here -------
2150  "EXAMPLE:"; echo = 2;
2151  ring exring=0,(x,y),dp;
2152  list hne=reddevelop(x14-3y2x11-y3x10-y2x9+3y4x8+y5x7+3y4x6+x5*(-y6+y5)
2153                      -3y6x3-y7x2+y8);
2154  print(hne[1][1]); // finite HNE
2155  print(extdevelop(hne[1],5)[1]);
2156  // ------------------------------------------
2157  print(hne[2][1]); // HNE that can be extended
2158  list ehne=extdevelop(hne[2],5);
2159  print(ehne[1]); // new HN-matrix has 5 columns
2160  param(hne[2]);
2161  param(ehne);
2162  kill HNEring,exring;
2163  echo = 0;
2164  // --- restore HNEring if previously defined ---
2165  if (defined(rettering)) {
2166   setring rettering;
2167   def HNEring=rettering;
2168   export HNEring;
2169  }
2170}
2171///////////////////////////////////////////////////////////////////////////////
2172
2173proc stripHNE (list l)
2174"USAGE:   stripHNE(l);    takes the output both of develop(f) and reddevelop(f)
2175 @*        ( list l (matrix m, intvec v, int s[,poly g,...])
2176             or list of lists in the form l )
2177RETURN:  a list in the same format as l, but all polynomials g are set to zero
2178NOTE:    The purpose of this procedure is to remove huge amounts of data
2179         no longer needed. It is useful, if one or more of the polynomials
2180         in l consumes much memory. It is still possible to compute invariants,
2181         parametrizations etc. with the stripped HNE(s), but it is not possible
2182         to use extdevelop with them.
2183SEE ALSO: develop, reddevelop, extdevelop
2184EXAMPLE: example stripHNE;  shows an example
2185"
2186{
2187 list h;
2188 if (typeof(l[1])=="matrix") { l[4]=poly(0); }
2189 else {
2190  for (int i=1; i<=size(l); i++) {
2191    h=l[i];
2192    h[4]=poly(0);
2193    l[i]=h;
2194  }
2195 }
2196 return(l);
2197}
2198example
2199{ "EXAMPLE:"; echo = 2;
2200  ring r=0,(x,y),dp;
2201  list hne=develop(x2+y3+y4);
2202  hne;
2203  stripHNE(hne);
2204}
2205///////////////////////////////////////////////////////////////////////////////
2206static
2207proc extractHNEs(list HNEs, int transvers)
2208"USAGE:  extractHNEs(HNEs,transvers);  list HNEs (output from HN),
2209        int transvers: 1 if x,y were exchanged, 0 else
2210RETURN: list of Hamburger-Noether-Extensions in the form of reddevelop
2211NOTE:   This procedure is only for internal purpose; examples don't make sense
2212"
2213{
2214 int i,maxspalte,hspalte,hnezaehler;
2215 list HNEaktu,Ergebnis;
2216 for (hnezaehler=1; hnezaehler<=size(HNEs); hnezaehler++) {
2217  maxspalte=0;
2218  HNEaktu=HNEs[hnezaehler];
2219  if (defined(HNDebugOn)) {"To process:";HNEaktu;}
2220  if (size(HNEaktu)!=size(HNEaktu[1])+2) {
2221     "The ideals and the hqs in HNEs[",hnezaehler,"] don't match!!";
2222     HNEs[hnezaehler];
2223  }
2224 //------------ ermittle  maximale Anzahl benoetigter Spalten: ----------------
2225  for (i=2; i<size(HNEaktu); i++) {
2226    hspalte=ncols(HNEaktu[i]);
2227    maxspalte=maxspalte*(hspalte < maxspalte)+hspalte*(hspalte >= maxspalte);
2228  }
2229 //------------- schreibe Ausgabe fuer hnezaehler-ten Zweig: ------------------
2230  matrix ma[size(HNEaktu)-2][maxspalte];
2231  for (i=1; i<=(size(HNEaktu)-2); i++) {
2232    if (ncols(HNEaktu[i+1]) > 1) {
2233      ma[i,1..ncols(HNEaktu[i+1])]=HNEaktu[i+1]; }
2234    else { ma[i,1]=HNEaktu[i+1][1];}
2235  }
2236  Ergebnis[hnezaehler]=list(ma,HNEaktu[1],transvers,HNEaktu[size(HNEaktu)]);
2237  kill ma;
2238 }
2239 return(Ergebnis);
2240}
2241///////////////////////////////////////////////////////////////////////////////
2242
2243proc factorfirst(poly f, int M, int N)
2244"USAGE : factorfirst(f,M,N); f poly, M,N int
2245RETURN: number d: f=c*(y^(N/e) - d*x^(M/e))^e with e=gcd(M,N), number c fitting
2246        0 if d does not exist
2247EXAMPLE: example factorfirst;  shows an example
2248"
2249{
2250 number c = koeff(f,0,N);
2251 number delta;
2252 int eps,l;
2253 int p=char(basering);
2254 string ringchar=charstr(basering);
2255
2256 if (c == 0) {"Something has gone wrong! I didn't get N correctly!"; exit;}
2257 int e = gcd(M,N);
2258
2259 if (p==0) { delta = koeff(f,M/ e,N - N/ e) / (-1*e*c); }
2260 else {
2261   if (e%p != 0) { delta = koeff(f,M/ e,N - N/ e) / (-1*e*c); }
2262   else {
2263     eps = e;
2264     for (l = 0; eps%p == 0; l=l+1) { eps=eps/ p;}
2265     if (defined(HNDebugOn)) {e," -> ",eps,"*",p,"^",l;}
2266     delta = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c);
2267
2268     if ((charstr(basering) != string(p)) and (delta != 0)) {
2269 //------ coefficient field is not Z/pZ => (p^l)th root is not identity -------
2270       delta=0;
2271       if (defined(HNDebugOn)) {
2272         "trivial factorization not implemented for",
2273         "parameters---I've to use 'factorize'";
2274       }
2275     }
2276   }
2277 }
2278 if (defined(HNDebugOn)) {"quasihomogeneous leading form:",f," = ",c,
2279        "* (y^"+string(N/ e),"-",delta,"* x^"+string(M/ e)+")^",e," ?";}
2280 if (f != c*(var(2)^(N/ e) - delta*var(1)^(M/ e))^e) {return(0);}
2281 else {return(delta);}
2282}
2283example
2284{ "EXAMPLE:"; echo = 2;
2285  ring exring=7,(x,y),dp;
2286  factorfirst(2*(y3-3x4)^5,20,15);
2287  factorfirst(x14+y7,14,7);
2288  factorfirst(x14+x8y3+y7,14,7);
2289}
2290///////////////////////////////////////////////////////////////////////////////
2291
2292proc HNdevelop (poly f)
2293"USAGE:   HNdevelop(f); f poly
2294RETURN:  Hamburger-Noether development of f:
2295         A list of lists in the form of develop(f)
2296         (matrix,intvec,int,poly -- the last int indicating irreducibility is
2297         omitted); each entry contains the data for one of the branches of f.
2298         For more details type `help develop;'
2299CREATE:  a ring with name `HNEring', variables `x,y' and ordering `ls' over
2300         a field extension of the current basering's ground field.
2301         As the Hamburger-Noether development normally does not exist
2302         in the originally given basering, @code{HNdevelop} always defines
2303         HNEring and changes to it. The field extension is chosen minimal.
2304NOTE:    @code{HNdevelop} decides which procedure (@code{develop} or
2305         @code{reddevelop}) applies best to the given problem and calls it.
2306SEE ALSO: develop, reddevelop, extdevelop
2307EXAMPLE: example HNdevelop;  shows an example
2308"
2309{
2310 int irred=0;
2311 //--------- Falls Ring (p^k,a),...: Wechsel in (p,a),... + minpoly -----------
2312 if ((find(charstr(basering),string(char(basering)))!=1) &&
2313     (charstr(basering)<>"real")) {
2314   string strmip=string(minpoly);
2315   string strf=string(f);
2316   execute("ring tempr=("+string(char(basering))+","+parstr(basering)+"),("
2317           +varstr(basering)+"),dp;");
2318   execute("minpoly="+strmip+";");
2319   execute("poly f="+strf+";");
2320   list hne=reddevelop(f);
2321   if ((voice==2) && (printlevel > -1)) {
2322     "// Attention: The parameter",par(1),"has changed its meaning!";
2323     "// It need no longer be a generator of the cyclic group of unities!";
2324   }
2325 }
2326 else {
2327 //--- Falls Ring (0,a),... + minpoly : solange factorize nicht in Singular ---
2328 //------- implementiert ist, develop aufrufen (kann spaeter entfallen) -------
2329   if ((char(basering)==0) && (npars(basering)==1)) {
2330     if (string(minpoly)<>"0") { irred=1; }
2331   }
2332 //------------------ Aufruf der geeigneten Prozedur --------------------------
2333   if (irred==0) {
2334     list hne=pre_HN(f,0);       // = reddevelop(f);
2335     dbprint(printlevel+1,
2336        "// result: "+string(size(hne))+" branch(es) successfully computed,",
2337        "//         basering has changed to HNEring");
2338   }
2339   else {
2340     def altring=basering;
2341     string strmip=string(minpoly);
2342     ring HNEring=(char(altring),`parstr(altring)`),(x,y),ls;
2343     execute("minpoly="+strmip+";");
2344     export HNEring;
2345     poly f=fetch(altring,f);
2346     list hn=develop(f,-1);
2347     list hne;
2348     if (hn[3] <> -1) {
2349       hne[1]=list(hn[1],hn[2],hn[3],hn[4]);
2350       if (hn[5] <> 1) {
2351   " ** WARNING : The curve is reducible, but only one branch could be found!";
2352       }
2353     }
2354     else { " ** Sorry -- could not find a HNE."; }
2355     dbprint(printlevel+1,"// note: basering has changed to HNEring");
2356   }
2357 }
2358 keepring basering;
2359 return(hne);
2360}
2361example
2362{
2363  // -------- prepare for example ---------
2364  if (nameof(basering)=="HNEring") {
2365   def rettering=HNEring;
2366   kill HNEring;
2367  }
2368  // ------ the example starts here -------
2369  "EXAMPLE:"; echo = 2;
2370  ring r=0,(x,y),dp;
2371  list hne=HNdevelop(x4-y6);
2372  nameof(basering);
2373  size(hne);           // i.e. x4-y6 has two branches
2374  print(hne[1][1]);    // HN-matrix of 1st branch
2375  param(hne[1]);
2376  param(hne[2]);
2377  displayInvariants(hne);
2378  kill HNEring,r;
2379  echo = 0;
2380  // --- restore HNEring if previously defined ---
2381  if (defined(rettering)) {
2382   setring rettering;
2383   def HNEring=rettering;
2384   export HNEring;
2385  }
2386}
2387///////////////////////////////////////////////////////////////////////////////
2388
2389proc reddevelop (poly f)
2390"USAGE:   reddevelop(f); f poly
2391RETURN:  Hamburger-Noether development of f:
2392         A list of lists in the form of develop(f)
2393         (matrix,intvec,int,poly -- the last int indicating irreducibility is
2394         omitted); each entry contains the data for one of the branches of f.
2395         For more details type `help develop;'
2396CREATE:  a ring with name `HNEring', variables `x,y' and ordering `ls' over
2397         a field extension of the current basering's ground field.
2398         As the Hamburger-Noether development of a reducible curve normally
2399         does not exist in the given basering, @code{reddevelop} always defines
2400         HNEring and changes to it. The field extension is chosen minimal.
2401SEE ALSO: develop, extdevelop, displayHNE
2402EXAMPLE: example reddevelop;  shows an example
2403"
2404{
2405 list Ergebnis=pre_HN(f,0);
2406 if (size(Ergebnis)>0) {     // otherwise an error may have occurred
2407  dbprint(printlevel+1,
2408   "// result: "+string(size(Ergebnis))+" branch(es) successfully computed,",
2409   "//         basering has changed to HNEring");
2410 }
2411 keepring basering;
2412 return(Ergebnis);
2413}
2414example
2415{
2416  // -------- prepare for example ---------
2417  if (nameof(basering)=="HNEring") {
2418   def rettering=HNEring;
2419   kill HNEring;
2420  }
2421  // ------ the example starts here -------
2422  "EXAMPLE:"; echo = 2;
2423  ring r=0,(x,y),dp;
2424  list hne=reddevelop(x4-y6);
2425  size(hne);           // i.e. x4-y6 has two branches
2426  print(hne[1][1]);    // HN-matrix of 1st branch
2427  param(hne[1]);
2428  param(hne[2]);
2429  displayInvariants(hne);
2430  kill HNEring,r;
2431  pause("Hit RETURN to continue:");
2432  // ----------------- a more interesting example: --------------------
2433  ring r = 32003,(x,y),dp;
2434  poly f = x25+x24-4x23-1x22y+4x22+8x21y-2x21-12x20y-4x19y2+4x20+10x19y
2435          +12x18y2-24x18y-20x17y2-4x16y3+x18+60x16y2+20x15y3-9x16y
2436          -80x14y3-10x13y4+36x14y2+60x12y4+2x11y5-84x12y3-24x10y5
2437          +126x10y4+4x8y6-126x8y5+84x6y6-36x4y7+9x2y8-1y9;
2438  list hne=reddevelop(f);
2439  size(hne);
2440  print(hne[1][1]);
2441  print(hne[4][1]);
2442  // a ring change was necessary, a is a parameter
2443  HNEring;
2444  kill HNEring,r;
2445  echo = 0;
2446  // --- restore HNEring if previously defined ---
2447  if (defined(rettering)) {
2448   setring rettering;
2449   def HNEring=rettering;
2450   export HNEring;
2451  }
2452}
2453///////////////////////////////////////////////////////////////////////////////
2454
2455static proc pre_HN (poly f, int essential)
2456"NOTE: This procedure is only for internal use, it is called via
2457      reddevelop or essdevelop"
2458{
2459 def altring = basering;
2460 int p = char(basering);                 // Ringcharakteristik
2461
2462 //-------------------- Tests auf Zulaessigkeit von basering ------------------
2463 if (charstr(basering)=="real") {
2464   " Singular cannot factorize over 'real' as ground field";
2465   return(list());
2466 }
2467 if (size(maxideal(1))<2) {
2468   " A univariate polynomial ring makes no sense !";
2469   return(list());
2470 }
2471 if ((size(maxideal(1))>2) && (printlevel > -1)) {
2472   " Warning: all but the first two variables are ignored!";
2473 }
2474 if (find(charstr(basering),string(char(basering)))!=1) {
2475   " ring of type (p^k,a) not implemented";
2476 //----------------------------------------------------------------------------
2477 // weder primitives Element noch factorize noch map "char p^k" -> "char -p"
2478 // [(p^k,a)->(p,a) ground field] noch fetch
2479 //----------------------------------------------------------------------------
2480   return(list());
2481 }
2482 //----------------- Definition eines neuen Ringes: HNEring -------------------
2483 string namex=varstr(1); string namey=varstr(2);
2484 if (string(char(altring))==charstr(altring)) {       // kein Parameter
2485   ring HNEring = char(altring),(x,y),ls;
2486   map m=altring,x,y;
2487   poly f=m(f);
2488   kill m;
2489 }
2490 else {
2491   string mipl=string(minpoly);
2492   if (mipl=="0") {
2493     " ** WARNING: No algebraic extension of this ground field is possible!";
2494     " ** We try to develop this polynomial, but if the need for an extension";
2495     " ** occurs during the calculation, we cannot proceed with the";
2496     " ** corresponding branches ...";
2497     execute("ring HNEring=("+charstr(basering)+"),(x,y),ls;");
2498 //--- ring ...=(char(.),`parstr()`),... geht nicht, wenn mehr als 1 Param. ---
2499   }
2500   else {
2501    string pa=parstr(altring);
2502    ring HNhelpring=p,`pa`,dp;
2503    execute("poly mipo="+mipl+";");  // Minimalpolynom in Polynom umgewandelt
2504    ring HNEring=(p,a),(x,y),ls;
2505    map getminpol=HNhelpring,a;
2506    mipl=string(getminpol(mipo));    // String umgewandelt mit 'a' als Param.
2507    execute("minpoly="+mipl+";");     // "minpoly=poly is not supported"
2508    kill HNhelpring, getminpol;
2509   }
2510   if (nvars(altring)==2) { poly f=fetch(altring,f); }
2511   else {
2512     map m=altring,x,y;
2513     poly f=m(f);
2514     kill m;
2515   }
2516 }
2517 export HNEring;
2518
2519 if (defined(HNDebugOn))
2520 {"received polynomial: ",f,", with x =",namex,", y =",namey;}
2521
2522 //----------------------- Variablendefinitionen ------------------------------
2523 int Abbruch,i,NullHNEx,NullHNEy;
2524 string str;
2525 list Newton,Ergebnis,hilflist;
2526
2527 //====================== Tests auf Zulaessigkeit des Polynoms ================
2528
2529 //-------------------------- Test, ob Einheit oder Null ----------------------
2530 if (subst(subst(f,x,0),y,0)!=0) {
2531   dbprint(printlevel+1,
2532           "The given polynomial is a unit in the power series ring!");
2533   keepring HNEring;
2534   return(list());                   // there are no HNEs
2535 }
2536 if (f==0) {
2537   dbprint(printlevel+1,"The given polynomial is zero!");
2538   keepring HNEring;
2539   return(list());                   // there are no HNEs
2540 }
2541
2542 //-----------------------  Test auf Quadratfreiheit --------------------------
2543
2544 if ((p==0) and (size(charstr(basering))==1)) {
2545
2546 //-------- Fall basering==0,... : Wechsel in Ring mit char >0 ----------------
2547 // weil squarefree eine Standardbasis berechnen muss (verwendet Syzygien)
2548 // -- wenn f in diesem Ring quadratfrei ist, dann erst recht im Ring HNEring
2549 //----------------------------------------------------------------------------
2550  int testerg=(polytest(f)==0);
2551  ring zweitring = 32003,(x,y),dp;
2552
2553  map polyhinueber=HNEring,x,y;         // fetch geht nicht
2554  poly f=polyhinueber(f);
2555  poly test_sqr=squarefree(f);
2556  if (test_sqr != f) {
2557   if (printlevel>0) {
2558     "Most probably the given polynomial is not squarefree. But the test was";
2559     "made in characteristic 32003 and not 0 to improve speed. You can";
2560     "(r) redo the test in char 0 (but this may take some time)";
2561     "(c) continue the development, if you're sure that the polynomial IS",
2562     "squarefree";
2563     if (testerg==1) {
2564       "(s) continue the development with a squarefree factor (*)";}
2565     "(q) or just quit the algorithm (default action)";
2566     "";"Please enter the letter of your choice:";
2567     str=read("")[1];             // reads one character
2568   }
2569   else { str="r"; }              // printlevel <= 0: non-interactive behaviour
2570   setring HNEring;
2571   map polyhinueber=zweitring,x,y;
2572   if (str=="r") {
2573     poly test_sqr=squarefree(f);
2574     if (test_sqr != f) {
2575      if (printlevel>0) { "The given polynomial is in fact not squarefree."; }
2576      else              { "The given polynomial is not squarefree!"; }
2577      "I'll continue with the radical.";
2578      f=test_sqr;
2579     }
2580     else {
2581      dbprint(printlevel,
2582       "everything is ok -- the polynomial is squarefree in characteristic 0");
2583     }
2584   }
2585   else {
2586     if ((str=="s") and (testerg==1)) {
2587       "(*)attention: it could be that the factor is only one in char 32003!";
2588       f=polyhinueber(test_sqr);
2589     }
2590     else {
2591       if (str<>"c") {
2592         setring altring;
2593         if(system("with","Namespaces")) { kill Top::HNEring; }
2594         kill HNEring;kill zweitring;
2595         return(list());}
2596       else { "if the algorithm doesn't terminate, you were wrong...";}
2597   }}
2598   kill zweitring;
2599   if (defined(HNDebugOn)) {"I continue with the polynomial",f; }
2600  }
2601  else {
2602    setring HNEring;
2603    kill zweitring;
2604  }
2605 }
2606 //------------------ Fall Char > 0 oder Ring hat Parameter -------------------
2607 else {
2608  poly test_sqr=squarefree(f);
2609  if (test_sqr != f) {
2610   if (printlevel>0) {
2611    if (test_sqr == 1) {
2612     "The given polynomial is of the form g^"+string(p)+",";
2613     "therefore not squarefree.  You can:";
2614     " (q) quit the algorithm (recommended) or";
2615     " (f) continue with the full radical (using a factorization of the";
2616     "     pure power part; this could take much time)";
2617     "";"Please enter the letter of your choice:";
2618     str=read("")[1];
2619     if (str<>"f") { str="q"; }
2620    }
2621    else {
2622     "The given polynomial is not squarefree.";
2623     if (p != 0)
2624      {
2625       " You can:";
2626       " (c) continue with a squarefree divisor (but factors of the form g^"
2627       +string(p);
2628       "     are lost; this is recommended, takes no more time)";
2629       " (f) continue with the full radical (using a factorization of the";
2630       "     pure power part; this could take much time)";
2631       " (q) quit the algorithm";
2632       "";"Please enter the letter of your choice:";
2633       str=read("")[1];
2634       if ((str<>"f") && (str<>"q")) { str="c"; }
2635      }
2636     else { "I'll continue with the radical."; str="c"; }
2637    }                                // endelse (test_sqr!=1)
2638   }
2639   else {
2640     "//** Error: The given polynomial is not squarefree!";
2641     "//** Since the global variable `printlevel' has the value",printlevel,
2642       "we stop here.";
2643     "//   Either call me again with a squarefree polynomial f or assign";
2644     "            printlevel=1;";
2645     "//   before calling me with a non-squarefree f.";
2646     "//   If printlevel > 0, I will present to you some possibilities how to",
2647       "proceed.";
2648     str="q";
2649   }
2650   if (str=="q") {
2651    if(system("with","Namespaces")) { kill Top::HNEring; }
2652    setring altring;kill HNEring;
2653    return(list());
2654   }
2655   if (str=="c") { f=test_sqr; }
2656   if (str=="f") { f=allsquarefree(f,test_sqr); }
2657  }
2658  if (defined(HNDebugOn)) {"I continue with the polynomial",f; }
2659
2660 }
2661 //====================== Ende Test auf Quadratfreiheit =======================
2662 if (subst(subst(f,x,0),y,0)!=0) {
2663   "Sorry. The remaining polynomial is a unit in the power series ring...";
2664   keepring HNEring;
2665   return(list());
2666 }
2667 //---------------------- Test, ob f teilbar durch x oder y -------------------
2668 if (subst(f,y,0)==0) {
2669   f=f/y; NullHNEy=1; }             // y=0 is a solution
2670 if (subst(f,x,0)==0) {
2671   f=f/x; NullHNEx=1; }             // x=0 is a solution
2672
2673 Newton=newtonpoly(f);
2674 i=1; Abbruch=0;
2675 //----------------------------------------------------------------------------
2676 // finde Eckpkt. des Newtonpolys, der den Teil abgrenzt, fuer den x transvers:
2677 // Annahme: Newton ist sortiert, s.d. Newton[1]=Punkt auf der y-Achse,
2678 // Newton[letzt]=Punkt auf der x-Achse
2679 //----------------------------------------------------------------------------
2680 while ((i<size(Newton)) and (Abbruch==0)) {
2681  if ((Newton[i+1][1]-Newton[i][1])>=(Newton[i][2]-Newton[i+1][2]))
2682   {Abbruch=1;}
2683  else {i=i+1;}
2684 }
2685 int grenze1=Newton[i][2];
2686 int grenze2=Newton[i][1];
2687 //----------------------------------------------------------------------------
2688 // Stelle Ring bereit zur Uebertragung der Daten im Fall einer Koerperer-
2689 // weiterung. Definiere Objekte, die spaeter uebertragen werden.
2690 // Binde die Listen (azeilen,...) an den Ring (um sie nicht zu ueberschreiben
2691 // bei Def. in einem anderen Ring).
2692 // Exportiere Objekte, damit sie auch in der proc HN noch da sind
2693 //----------------------------------------------------------------------------
2694 ring HNE_noparam = char(altring),(a,x,y),ls;
2695 export HNE_noparam;
2696 poly f;
2697 list azeilen=ideal(0);
2698 list HNEs=ideal(0);
2699 list aneu=ideal(0);
2700 list faktoren=ideal(0);
2701 ideal deltais;
2702 poly delta;                   // nicht number, weil delta von a abhaengen kann
2703 export f,azeilen,HNEs,aneu,faktoren,deltais,delta;
2704 //----- hier steht die Anzahl bisher benoetigter Ringerweiterungen drin: -----
2705 int EXTHNEnumber=0; export EXTHNEnumber;
2706 setring HNEring;
2707
2708 // ================= Die eigentliche Berechnung der HNE: =====================
2709
2710 // ------- Berechne HNE von allen Zweigen, fuer die x transversal ist: -------
2711 if (defined(HNDebugOn))
2712   {"1st step: Treat Newton polygon until height",grenze1;}
2713 if (grenze1>0) {
2714  hilflist=HN(f,grenze1,1,essential);
2715  if (typeof(hilflist[1][1])=="ideal") { hilflist[1]=list(); }
2716 //- fuer den Fall, dass keine Zweige in transz. Erw. berechnet werden konnten-
2717  Ergebnis=extractHNEs(hilflist[1],0);
2718  if (hilflist[2]!=-1) {
2719    if (defined(HNDebugOn)) {" ring change in HN(",1,") detected";}
2720    poly transfproc=hilflist[2];
2721    map hole=HNE_noparam,transfproc,x,y;
2722    setring HNE_noparam;
2723    f=imap(HNEring,f);
2724    setring EXTHNEring(EXTHNEnumber);
2725    poly f=hole(f);
2726  }
2727 }
2728 if (NullHNEy==1) {
2729  Ergebnis=Ergebnis+list(list(matrix(ideal(0,x)),intvec(1),int(0),poly(0)));
2730 }
2731 // --------------- Berechne HNE von allen verbliebenen Zweigen: --------------
2732 if (defined(HNDebugOn))
2733    {"2nd step: Treat Newton polygon until height",grenze2;}
2734 if (grenze2>0) {
2735  map xytausch=basering,y,x;
2736  kill hilflist;
2737  def letztring=basering;
2738  if (EXTHNEnumber==0) { setring HNEring; }
2739  else                 { setring EXTHNEring(EXTHNEnumber); }
2740  list hilflist=HN(xytausch(f),grenze2,1,essential);
2741  if (typeof(hilflist[1][1])=="ideal") { hilflist[1]=list(); }
2742  if (not defined(Ergebnis)) {
2743 //-- HN wurde schon mal ausgefuehrt; Ringwechsel beim zweiten Aufruf von HN --
2744    if (defined(HNDebugOn)) {" ring change in HN(",1,") detected";}
2745    poly transfproc=hilflist[2];
2746    map hole=HNE_noparam,transfproc,x,y;
2747    setring HNE_noparam;
2748    list Ergebnis=imap(letztring,Ergebnis);
2749    setring EXTHNEring(EXTHNEnumber);
2750    list Ergebnis=hole(Ergebnis);
2751  }
2752  Ergebnis=Ergebnis+extractHNEs(hilflist[1],1);
2753 }
2754 if (NullHNEx==1) {
2755  Ergebnis=Ergebnis+list(list(matrix(ideal(0,x)),intvec(1),int(1),poly(0)));
2756 }
2757 //------------------- Loesche globale, nicht mehr benoetigte Objekte: --------
2758 if (EXTHNEnumber>0) {
2759  if(system("with","Namespaces")) { kill Top::HNEring; }
2760  kill HNEring;
2761  def HNEring=EXTHNEring(EXTHNEnumber);
2762  setring HNEring;
2763  export HNEring;
2764  kill EXTHNEring(1..EXTHNEnumber);
2765 }
2766 kill HNE_noparam;
2767 kill EXTHNEnumber;
2768 keepring basering;
2769 return(Ergebnis);
2770}
2771
2772///////////////////////////////////////////////////////////////////////////////
2773proc essdevelop (poly f)
2774"USAGE:   essdevelop(f); f poly
2775RETURN:  Hamburger-Noether development of essential branches of f:
2776         A list of lists in the form of develop(f)
2777         (matrix,intvec,int,poly -- the last int indicating irreducibility is
2778         omitted); each entry contains the data for one of the branches of f.
2779         For more details type `help develop;'
2780CREATE:  a ring with name `HNEring', variables `x,y' and ordering `ls' over
2781         a field extension of the current basering's ground field.
2782         As the Hamburger-Noether development of a reducible curve normally
2783         does not exist in the given basering, @code{essdevelop} always defines
2784         HNEring and changes to it. The field extension is chosen minimal.
2785NOTE:    If the HNE needs a field extension, some of the branches will be
2786         conjugate. In this case @code{essdevelop} reduces the computation to
2787         one representative for each group of conjugate branches.
2788         Note that the degree of each branch is in general less than the degree
2789         of the field extension in which all HNEs can be put.
2790
2791         Use @code{reddevelop} or @code{HNdevelop} to compute a complete HNE,
2792         i.e. a HNE for all branches.
2793SEE ALSO: develop, reddevelop, HNdevelop, extdevelop
2794EXAMPLE: example essdevelop;  shows an example
2795"
2796{
2797 list Ergebnis=pre_HN(f,1);
2798 dbprint(printlevel+1,
2799    "// result: "+string(size(Ergebnis))+" branch(es) successfully computed;");
2800 if (string(minpoly) <> "0") {
2801   dbprint(printlevel+1,
2802    "// note that conjugate branches are omitted and that the number",
2803    "// of branches represented by each remaining one may vary!");
2804 }
2805 dbprint(printlevel+1,
2806    "// basering has changed to HNEring");
2807 keepring basering;
2808 return(Ergebnis);
2809}
2810example
2811{
2812  // -------- prepare for example ---------
2813  if (nameof(basering)=="HNEring") {
2814   def rettering=HNEring;
2815   kill HNEring;
2816  }
2817  // ------ the example starts here -------
2818  "EXAMPLE:"; echo = 2;
2819  ring r=2,(x,y),dp;
2820  poly f=(x4+x2y+y2)*(x3+xy2+y3);
2821  // --------- compute all branches: ---------
2822  list hne=reddevelop(f);
2823  displayHNE(hne);
2824  setring r;
2825  // --- compute only one of conjugate branches: ---
2826  list hne=essdevelop(f);
2827  displayHNE(hne);
2828  // nr. 1 of essdevelop represents nr. 1 - 3 of reddevelop and
2829  // nr. 2 of essdevelop represents nr. 4 + 5 of reddevelop
2830  echo = 0;
2831  // --- restore HNEring if previously defined ---
2832  if (defined(rettering)) {
2833   setring rettering;
2834   def HNEring=rettering;
2835   export HNEring;
2836  }
2837}
2838
2839///////////////////////////////////////////////////////////////////////////////
2840static
2841proc HN (poly f,int grenze, int Aufruf_Ebene, int essential)
2842"NOTE: This procedure is only for internal use, it is called via pre_HN"
2843{
2844 //---------- Variablendefinitionen fuer den unverzweigten Teil: --------------
2845 if (defined(HNDebugOn)) {"procedure HN",Aufruf_Ebene;}
2846 int Abbruch,ende,i,j,e,M,N,Q,R,zeiger,zeile,zeilevorher;
2847 intvec hqs;
2848 poly fvorher;
2849 list erg=ideal(0); list HNEs=ideal(0); // um die Listen an den Ring zu binden
2850
2851 //-------------------- Bedeutung von Abbruch: --------------------------------
2852 //------- 0:keine Verzweigung | 1:Verzweigung,nicht fertig | 2:fertig --------
2853 //
2854 // Struktur von HNEs : Liste von Listen L (fuer jeden Zweig) der Form
2855 // L[1]=intvec (hqs), L[2],L[3],... ideal (die Zeilen (0,1,...) der HNE)
2856 // L[letztes]=poly (transformiertes f)
2857 //----------------------------------------------------------------------------
2858 list Newton;
2859 number delta;
2860 int p = char(basering);                // Ringcharakteristik
2861 list azeilen=ideal(0);
2862 ideal hilfid; list hilflist=ideal(0); intvec hilfvec;
2863
2864 // ======================= der unverzweigte Teil: ============================
2865 while (Abbruch==0) {
2866  Newton=newtonpoly(f);
2867  zeiger=find_in_list(Newton,grenze);
2868  if (Newton[zeiger][2] != grenze)
2869    {"Didn't find an edge in the Newton polygon!";}
2870  if (zeiger==size(Newton)-1) {
2871    if (defined(HNDebugOn)) {"only one relevant side in Newton polygon";}
2872    M=Newton[zeiger+1][1]-Newton[zeiger][1];
2873    N=Newton[zeiger][2]-Newton[zeiger+1][2];
2874    R = M%N;
2875    Q = M / N;
2876
2877 //-------- 1. Versuch: ist der quasihomogene Leitterm reine Potenz ? ---------
2878 //              (dann geht alles wie im irreduziblen Fall)
2879 //----------------------------------------------------------------------------
2880    e = gcd(M,N);
2881    delta=factorfirst(redleit(f,Newton[zeiger],Newton[zeiger+1])
2882                      /x^Newton[zeiger][1],M,N);
2883    if (delta==0) {
2884      if (defined(HNDebugOn)) {" The given polynomial is reducible !";}
2885      Abbruch=1;
2886    }
2887    if (Abbruch==0) {
2888 //-------------- f,zeile retten fuer den Spezialfall (###): ------------------
2889      fvorher=f;zeilevorher=zeile;
2890      if (R==0) {
2891 //------------- transformiere f mit T1, wenn kein Abbruch nachher: -----------
2892        if (N>1) { f = T1_Transform(f,delta,M/ e); }
2893        else     { ende=1; }
2894        if (defined(HNDebugOn)) {"a("+string(zeile)+","+string(Q)+") =",delta;}
2895        azeilen[zeile+1][Q]=delta;
2896      }
2897      else {
2898 //------------- R > 0 : transformiere f mit T2 -------------------------------
2899        erg=T2_Transform(f,delta,M,N,referencepoly(Newton));
2900        f=erg[1];delta=erg[2];
2901 //------- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: ---------
2902        while (R!=0) {
2903         if (defined(HNDebugOn)) { "h("+string(zeile)+") =",Q; }
2904         hqs[zeile+1]=Q;         // denn zeile beginnt mit dem Wert 0
2905 //------------------ markiere das Zeilenende der HNE: ------------------------
2906         azeilen[zeile+1][Q+1]=x;
2907         zeile=zeile+1;
2908 //----------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ---------
2909         azeilen[zeile+1]=ideal(0);
2910         M=N; N=R; R=M%N; Q=M / N;
2911        }
2912        if (defined(HNDebugOn)) {"a("+string(zeile)+","+string(Q)+") =",delta;}
2913        azeilen[zeile+1][Q]=delta;
2914      }
2915      if (defined(HNDebugOn)) {"transformed polynomial: ",f;}
2916      grenze=e;
2917 //----------------------- teste Abbruchbedingungen: --------------------------
2918      if (subst(f,y,0)==0) {              // <==> y|f
2919        dbprint(printlevel+1,"finite HNE of one branch found");
2920           // voice abzufragen macht bei rekursiven procs keinen Sinn
2921        azeilen[zeile+1][Q+1]=x;
2922 //- Q wird nur in hqs eingetragen, wenn der Spezialfall nicht eintritt (s.u.)-
2923        Abbruch=2;
2924        if (grenze>1) {
2925         if (jet(f,1,intvec(0,1))==0) {
2926 //------ jet(...)=alle Monome von f, die nicht durch y2 teilbar sind ---------
2927 "THE TEST FOR SQUAREFREENESS WAS BAD!! The polynomial was NOT squarefree!!!";}
2928         else {
2929 //-------------------------- Spezialfall (###): ------------------------------
2930 // Wir haben das Problem, dass die HNE eines Zweiges hier abbricht, aber ein
2931 // anderer Zweig bis hierher genau die gleiche HNE hat, die noch weiter geht
2932 // Loesung: mache Transform. rueckgaengig und behandle f im Verzweigungsteil
2933 //----------------------------------------------------------------------------
2934          Abbruch=1;
2935          f=fvorher;zeile=zeilevorher;grenze=Newton[zeiger][2];
2936        }}
2937        else {f=0;}     // f nicht mehr gebraucht - spare Speicher
2938        if (Abbruch==2) { hqs[zeile+1]=Q; }
2939      }                 // Spezialfall nicht eingetreten
2940      else {
2941        if (ende==1) {
2942          dbprint(printlevel+1,"HNE of one branch found");
2943          Abbruch=2; hqs[zeile+1]=-Q-1;}
2944      }
2945    }                   // end(if Abbruch==0)
2946  }                     // end(if zeiger...)
2947  else { Abbruch=1;}
2948 }                      // end(while Abbruch==0)
2949
2950 // ===================== der Teil bei Verzweigung: ===========================
2951
2952 if (Abbruch==1) {
2953 //---------- Variablendefinitionen fuer den verzweigten Teil: ----------------
2954  poly leitf,teiler,transformiert;
2955  list aneu=ideal(0);
2956  list faktoren;
2957  list HNEakut=ideal(0);
2958  ideal deltais;
2959  intvec eis;
2960  int zaehler,hnezaehler,zl,zl1,M1,N1,R1,Q1,needext;
2961  int numberofRingchanges,lastRingnumber,ringischanged,flag;
2962  string letztringname;
2963
2964  zeiger=find_in_list(Newton,grenze);
2965  if (defined(HNDebugOn)) {
2966    "Branching part reached---Newton polygon :",Newton;
2967    "relevant part until height",grenze,", from",Newton[zeiger],"on";
2968  }
2969  azeilen=list(hqs)+azeilen; // hat jetzt Struktur von HNEs: hqs in der 1.Zeile
2970
2971 //======= Schleife fuer jede zu betrachtende Seite des Newtonpolygons: =======
2972  for(i=zeiger; i<size(Newton); i++) {
2973   if (defined(HNDebugOn)) { "we consider side",Newton[i],Newton[i+1]; }
2974   M=Newton[i+1][1]-Newton[i][1];
2975   N=Newton[i][2]-Newton[i+1][2];
2976   R = M%N;
2977   Q = M / N;
2978   e=gcd(M,N);
2979   needext=1;
2980   letztringname=nameof(basering);
2981   lastRingnumber=EXTHNEnumber;
2982   faktoren=list(ideal(charPoly(redleit(f,Newton[i],Newton[i+1])
2983                       /(x^Newton[i][1]*y^Newton[i+1][2]),M,N)  ),
2984                 intvec(1));                  // = (zu faktoriserendes Poly, 1)
2985
2986 //-- wechsle so lange in Ringerw., bis Leitform vollst. in Linearfakt. zerf.:-
2987   for (numberofRingchanges=1; needext==1; numberofRingchanges++) {
2988    leitf=redleit(f,Newton[i],Newton[i+1])/(x^Newton[i][1]*y^Newton[i+1][2]);
2989    delta=factorfirst(leitf,M,N);
2990    needext=0;
2991    if (delta==0) {
2992
2993 //---------- Sonderbehandlung: faktorisere einige Polynome ueber Q(a): -------
2994     if (charstr(basering)=="0,a") {
2995       faktoren,flag=extrafactor(leitf,M,N);  // damit funktion. Bsp. Baladi 5
2996       if (flag==0)
2997         { ERROR("Could not factorize in field of type (0,a)!"); }
2998     }
2999     else {
3000 //------------------ faktorisiere das charakt. Polynom: ----------------------
3001       if ((numberofRingchanges==1) or (essential==0)) {
3002         faktoren=factorlist(faktoren);
3003       }
3004       else {     // eliminiere alle konjugierten Nullstellen bis auf eine:
3005         ideal hilf_id;
3006         for (zaehler=1; zaehler<=size(faktoren[1]); zaehler++) {
3007           hilf_id=factorize(faktoren[1][zaehler],0)[1];
3008           if (size(hilf_id)>1) { faktoren[1][zaehler]=hilf_id[2]; }
3009           else                 { faktoren[1][zaehler]=hilf_id[1]; }
3010         }
3011       }
3012     }
3013
3014     zaehler=1; eis=0;
3015     for (j=1; j<=size(faktoren[2]); j++) {
3016      teiler=faktoren[1][j];
3017      if (teiler/y != 0) {         // sonst war's eine Einheit --> wegwerfen!
3018        if (defined(HNDebugOn)) {"factor of leading form found:",teiler;}
3019        if (teiler/y2 == 0) {      // --> Faktor hat die Form cy+d
3020          deltais[zaehler]=-subst(teiler,y,0)/koeff(teiler,0,1); //=-d/c
3021          eis[zaehler]=faktoren[2][j];
3022          zaehler++;
3023        }
3024        else {
3025          dbprint(printlevel+1,
3026             " Change of basering (field extension) necessary!");
3027          if (defined(HNDebugOn)) { teiler,"is not properly factored!"; }
3028          if (needext==0) { poly zerlege=teiler; }
3029          needext=1;
3030        }
3031      }
3032     }                             // end(for j)
3033    }
3034    else { deltais=ideal(delta); eis=e;}
3035    if (defined(HNDebugOn)) {"roots of char. poly:";deltais;
3036                             "with multiplicities:",eis;}
3037    if (needext==1) {
3038 //--------------------- fuehre den Ringwechsel aus: --------------------------
3039      ringischanged=1;
3040      if ((size(parstr(basering))>0) && string(minpoly)=="0") {
3041        " ** We've had bad luck! The HNE cannot completely be calculated!";
3042                                   // HNE in transzendenter Erw. fehlgeschlagen
3043        kill zerlege;
3044        ringischanged=0; break;    // weiter mit gefundenen Faktoren
3045      }
3046      if (parstr(basering)=="") {
3047        EXTHNEnumber++;
3048        splitring(zerlege,"EXTHNEring("+string(EXTHNEnumber)+")");
3049        poly transf=0;
3050        poly transfproc=0;
3051      }
3052      else {
3053        if (defined(translist)) { kill translist; } // Vermeidung einer Warnung
3054        if (numberofRingchanges>1) {  // ein Ringwechsel hat nicht gereicht
3055         list translist=splitring(zerlege,"",list(transf,transfproc,faktoren));
3056         poly transf=translist[1];
3057         poly transfproc=translist[2];
3058         list faktoren=translist[3];
3059        }
3060        else {
3061         if (defined(transfproc)) { // in dieser proc geschah schon Ringwechsel
3062          EXTHNEnumber++;
3063          list translist=splitring(zerlege,"EXTHNEring("
3064               +string(EXTHNEnumber)+")",list(a,transfproc));
3065          poly transf=translist[1];
3066          poly transfproc=translist[2];
3067         }
3068         else {
3069          EXTHNEnumber++;
3070          list translist=splitring(zerlege,"EXTHNEring("
3071               +string(EXTHNEnumber)+")",a);
3072          poly transf=translist[1];
3073          poly transfproc=transf;
3074        }}
3075      }
3076 //----------------------------------------------------------------------------
3077 // transf enthaelt jetzt den alten Parameter des Ringes, der aktiv war vor
3078 // Beginn der Schleife (evtl. also ueber mehrere Ringwechsel weitergereicht),
3079 // transfproc enthaelt den alten Parm. des R., der aktiv war zu Beginn der
3080 // Prozedur, und der an die aufrufende Prozedur zurueckgegeben werden muss
3081 // transf ist Null, falls der alte Ring keinen Parameter hatte,
3082 // das gleiche gilt fuer transfproc
3083 //----------------------------------------------------------------------------
3084
3085 //------ Neudef. von Variablen, Uebertragung bisher errechneter Daten: -------
3086      poly leitf,teiler,transformiert;
3087      list aneu=ideal(0);
3088      ideal deltais;
3089      number delta;
3090      setring HNE_noparam;
3091      if (defined(letztring)) { kill letztring; }
3092      if (lastRingnumber>0) { def letztring=EXTHNEring(lastRingnumber); }
3093      else                  { def letztring=HNEring; }
3094      f=imap(letztring,f);
3095      faktoren=imap(letztring,faktoren);
3096      setring EXTHNEring(EXTHNEnumber);
3097      map hole=HNE_noparam,transf,x,y;
3098      poly f=hole(f);
3099      if (not defined(faktoren)) {
3100        list faktoren=hole(faktoren);
3101      }
3102    }
3103   }    // end (while needext==1) bzw. for (numberofRingchanges)
3104
3105   if (eis==0) { i++; continue; }
3106   if (ringischanged==1) {
3107    list erg,hilflist,HNEakut;            // dienen nur zum Sp. von Zwi.erg.
3108    ideal hilfid;
3109    erg=ideal(0); hilflist=erg; HNEakut=erg;
3110
3111    hole=HNE_noparam,transf,x,y;
3112    setring HNE_noparam;
3113    azeilen=imap(letztring,azeilen);
3114    HNEs=imap(letztring,HNEs);
3115
3116    setring EXTHNEring(EXTHNEnumber);
3117    list azeilen=hole(azeilen);
3118    list HNEs=hole(HNEs);
3119    kill letztring;
3120    ringischanged=0;
3121   }
3122
3123 //============ Schleife fuer jeden gefundenen Faktor der Leitform: ===========
3124   for (j=1; j<=size(eis); j++) {
3125 //-- Mache Transf. T1 oder T2, trage Daten in HNEs ein, falls HNE abbricht: --
3126
3127 //------------------------ Fall R==0: ----------------------------------------
3128    if (R==0) {
3129      transformiert = T1_Transform(f,number(deltais[j]),M/ e);
3130      if (defined(HNDebugOn)) {
3131        "a("+string(zeile)+","+string(Q)+") =",deltais[j];
3132        "transformed polynomial: ",transformiert;
3133      }
3134      if (subst(transformiert,y,0)==0) {
3135       dbprint(printlevel+1,"finite HNE found");
3136       hnezaehler++;
3137 //------------ trage deltais[j],x ein in letzte Zeile, fertig: ---------------
3138       HNEakut=azeilen+list(poly(0));        // =HNEs[hnezaehler];
3139       hilfid=HNEakut[zeile+2]; hilfid[Q]=deltais[j]; hilfid[Q+1]=x;
3140       HNEakut[zeile+2]=hilfid;
3141       HNEakut[1][zeile+1]=Q;                // aktualisiere Vektor mit den hqs
3142       HNEs[hnezaehler]=HNEakut;
3143       if (eis[j]>1) {
3144        transformiert=transformiert/y;
3145        if (subst(transformiert,y,0)==0) {
3146 "THE TEST FOR SQUAREFREENESS WAS BAD!! The polynomial was NOT squarefree!!!";}
3147        else {
3148 //------ Spezialfall (###) eingetreten: Noch weitere Zweige vorhanden --------
3149          eis[j]=eis[j]-1;
3150        }
3151       }
3152      }
3153    }
3154    else {
3155 //------------------------ Fall R <> 0: --------------------------------------
3156      erg=T2_Transform(f,number(deltais[j]),M,N,referencepoly(Newton));
3157      transformiert=erg[1];delta=erg[2];
3158      if (defined(HNDebugOn)) {"transformed polynomial: ",transformiert;}
3159      if (subst(transformiert,y,0)==0) {
3160       dbprint(printlevel+1,"finite HNE found");
3161       hnezaehler++;
3162 //---------------- trage endliche HNE in HNEs ein: ---------------------------
3163       HNEakut=azeilen;           // dupliziere den gemeins. Anfang der HNE's
3164       zl=2;                      // (kommt schliesslich nach HNEs[hnezaehler])
3165 //----------------------------------------------------------------------------
3166 // Werte von:  zeile: aktuelle Zeilennummer der HNE (gemeinsamer Teil)
3167 //             zl   : die HNE spaltet auf; zeile+zl ist der Index fuer die
3168 // Zeile des aktuellen Zweigs; (zeile+zl-2) ist die tatsaechl. Zeilennr.
3169 // (bei 0 angefangen) der HNE  ([1] <- intvec(hqs), [2] <- 0. Zeile usw.)
3170 //----------------------------------------------------------------------------
3171
3172 //---------- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: ------
3173       M1=M;N1=N;R1=R;Q1=M1/ N1;
3174       while (R1!=0) {
3175        if (defined(HNDebugOn)) { "h("+string(zeile+zl-2)+") =",Q1; }
3176        HNEakut[1][zeile+zl-1]=Q1;
3177        HNEakut[zeile+zl][Q1+1]=x;
3178                                  // markiere das Zeilenende der HNE
3179        zl=zl+1;
3180 //-------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ------------
3181        HNEakut[zeile+zl]=ideal(0);
3182
3183        M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1;
3184       }
3185       if (defined(HNDebugOn)) {
3186         "a("+string(zeile+zl-2)+","+string(Q1)+") =",delta;
3187       }
3188       HNEakut[zeile+zl][Q1]  =delta;
3189       HNEakut[zeile+zl][Q1+1]=x;
3190       HNEakut[1][zeile+zl-1] =Q1;     // aktualisiere Vektor mit hqs
3191       HNEakut[zeile+zl+1]=poly(0);
3192       HNEs[hnezaehler]=HNEakut;
3193 //-------------------- Ende der Eintragungen in HNEs -------------------------
3194
3195       if (eis[j]>1) {
3196        transformiert=transformiert/y;
3197        if (subst(transformiert,y,0)==0) {
3198 "THE TEST FOR SQUAREFREENESS WAS BAD!! The polynomial was NOT squarefree!!!";}
3199         else {
3200 //--------- Spezialfall (###) eingetreten: Noch weitere Zweige vorhanden -----
3201          eis[j]=eis[j]-1;
3202       }}
3203      }                           // endif (subst()==0)
3204    }                             // endelse (R<>0)
3205
3206 //========== Falls HNE nicht abbricht: Rekursiver Aufruf von HN: =============
3207 //------------------- Berechne HNE von transformiert -------------------------
3208    if (subst(transformiert,y,0)!=0) {
3209     lastRingnumber=EXTHNEnumber;
3210     list HNerg=HN(transformiert,eis[j],Aufruf_Ebene+1,essential);
3211     if (HNerg[2]==-1) {          // kein Ringwechsel in HN aufgetreten
3212       aneu=HNerg[1];  }
3213     else {
3214       if (defined(HNDebugOn))
3215          {" ring change in HN(",Aufruf_Ebene+1,") detected";}
3216       list aneu=HNerg[1];
3217       poly transfproc=HNerg[2];
3218
3219 //- stelle lokale Var. im neuen Ring wieder her und rette ggf. ihren Inhalt: -
3220       list erg,hilflist,faktoren,HNEakut;
3221       ideal hilfid;
3222       erg=ideal(0); hilflist=erg; faktoren=erg; HNEakut=erg;
3223       poly leitf,teiler,transformiert;
3224
3225       map hole=HNE_noparam,transfproc,x,y;
3226       setring HNE_noparam;
3227       if (lastRingnumber>0) { def letztring=EXTHNEring(lastRingnumber); }
3228       else                  { def letztring=HNEring; }
3229       HNEs=imap(letztring,HNEs);
3230       azeilen=imap(letztring,azeilen);
3231       deltais=imap(letztring,deltais);
3232       delta=imap(letztring,delta);
3233       f=imap(letztring,f);
3234
3235       setring EXTHNEring(EXTHNEnumber);
3236       list HNEs=hole(HNEs);
3237       list azeilen=hole(azeilen);
3238       ideal deltais=hole(deltais);
3239       number delta=number(hole(delta));
3240       poly f=hole(f);
3241     }
3242     kill HNerg;
3243 //----------------------------------------------------------------------------
3244 // HNerg muss jedesmal mit "list" neu definiert werden, weil vorher noch nicht
3245 // ------- klar ist, ob der Ring nach Aufruf von HN noch derselbe ist --------
3246
3247 //============= Verknuepfe bisherige HNE mit von HN gelieferten HNEs: ========
3248     if (R==0) {
3249       HNEs,hnezaehler=constructHNEs(HNEs,hnezaehler,aneu,azeilen,zeile,
3250                       deltais,Q,j);
3251     }
3252     else {
3253      for (zaehler=1; zaehler<=size(aneu); zaehler++) {
3254       hnezaehler++;
3255       HNEakut=azeilen;          // dupliziere den gemeinsamen Anfang der HNE's
3256       zl=2;                     // (kommt schliesslich nach HNEs[hnezaehler])
3257 //---------------- Trage Beitrag dieser Transformation T2 ein: ---------------
3258 //--------- Zur Bedeutung von zeile, zl: siehe Kommentar weiter oben ---------
3259
3260 //--------- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: -------
3261       M1=M;N1=N;R1=R;Q1=M1/ N1;
3262       while (R1!=0) {
3263        if (defined(HNDebugOn)) { "h("+string(zeile+zl-2)+") =",Q1; }
3264        HNEakut[1][zeile+zl-1]=Q1;
3265        HNEakut[zeile+zl][Q1+1]=x;    // Markierung des Zeilenendes der HNE
3266        zl=zl+1;
3267 //-------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ------------
3268        HNEakut[zeile+zl]=ideal(0);
3269        M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1;
3270       }
3271       if (defined(HNDebugOn)) {
3272         "a("+string(zeile+zl-2)+","+string(Q1)+") =",delta;
3273       }
3274       HNEakut[zeile+zl][Q1]=delta;
3275
3276 //--- Daten aus T2_Transform sind eingetragen; haenge Daten von HN an: -------
3277       hilfid=HNEakut[zeile+zl];
3278       for (zl1=Q1+1; zl1<=ncols(aneu[zaehler][2]); zl1++) {
3279        hilfid[zl1]=aneu[zaehler][2][zl1];
3280       }
3281       HNEakut[zeile+zl]=hilfid;
3282 //--- vorher HNEs[.][zeile+zl]<-aneu[.][2], jetzt [zeile+zl+1] <- [3] usw.: --
3283       for (zl1=3; zl1<=size(aneu[zaehler]); zl1++) {
3284         HNEakut[zeile+zl+zl1-2]=aneu[zaehler][zl1];
3285       }
3286 //--- setze die hqs zusammen: HNEs[hnezaehler][1]=HNEs[..][1],aneu[..][1] ----
3287       hilfvec=HNEakut[1],aneu[zaehler][1];
3288       HNEakut[1]=hilfvec;
3289 //----------- weil nicht geht: liste[1]=liste[1],aneu[zaehler][1] ------------
3290       HNEs[hnezaehler]=HNEakut;
3291      }                     // end(for zaehler)
3292     }                      // endelse (R<>0)
3293    }                       // endif (subst()!=0)  (weiteres Aufblasen mit HN)
3294
3295   }                        // end(for j) (Behandlung der einzelnen delta_i)
3296
3297  }
3298  keepring basering;
3299  if (defined(transfproc)) { return(list(HNEs,transfproc)); }
3300  else                     { return(list(HNEs,poly(-1))); }
3301 // -1 als 2. Rueckgabewert zeigt an, dass kein Ringwechsel stattgefunden hat -
3302 }
3303 else {
3304  HNEs[1]=list(hqs)+azeilen+list(f); // f ist das transform. Poly oder Null
3305  keepring basering;
3306  return(list(HNEs,poly(-1)));
3307 //-- in dieser proc trat keine Verzweigung auf, also auch kein Ringwechsel ---
3308 }
3309}
3310///////////////////////////////////////////////////////////////////////////////
3311
3312static
3313proc constructHNEs (list HNEs,int hnezaehler,list aneu,list azeilen,
3314                    int zeile,ideal deltais,int Q,int j)
3315"NOTE: This procedure is only for internal use, it is called via HN"
3316{
3317  int zaehler,zl;
3318  ideal hilfid;
3319  list hilflist;
3320  intvec hilfvec;
3321  for (zaehler=1; zaehler<=size(aneu); zaehler++) {
3322     hnezaehler++;
3323     // HNEs[hnezaehler]=azeilen;            // dupliziere gemeins. Anfang
3324 //----------------------- trage neu berechnete Daten ein ---------------------
3325     hilfid=azeilen[zeile+2];
3326     hilfid[Q]=deltais[j];
3327     for (zl=Q+1; zl<=ncols(aneu[zaehler][2]); zl++) {
3328      hilfid[zl]=aneu[zaehler][2][zl];
3329     }
3330     hilflist=azeilen; hilflist[zeile+2]=hilfid;
3331 //------------------ haenge uebrige Zeilen von aneu[] an: --------------------
3332     for (zl=3; zl<=size(aneu[zaehler]); zl++) {
3333      hilflist[zeile+zl]=aneu[zaehler][zl];
3334     }
3335 //--- setze die hqs zusammen: HNEs[hnezaehler][1]=HNEs[..][1],aneu[..][1] ----
3336     if (hilflist[1]==0) { hilflist[1]=aneu[zaehler][1]; }
3337     else { hilfvec=hilflist[1],aneu[zaehler][1]; hilflist[1]=hilfvec; }
3338     HNEs[hnezaehler]=hilflist;
3339  }
3340  return(HNEs,hnezaehler);
3341}
3342///////////////////////////////////////////////////////////////////////////////
3343
3344proc referencepoly (list newton)
3345"USAGE:   referencepoly(newton);
3346         newton is list of intvec(x,y) which represents points in the Newton
3347         diagram (e.g. output of the proc newtonpoly)
3348RETURN:  a polynomial which has newton as Newton diagram
3349SEE ALSO: newtonpoly
3350EXAMPLE: example referencepoly;  shows an example
3351"
3352{
3353 poly f;
3354 for (int i=1; i<=size(newton); i++) {
3355   f=f+var(1)^newton[i][1]*var(2)^newton[i][2];
3356 }
3357 return(f);
3358}
3359example
3360{ "EXAMPLE:"; echo = 2;
3361 ring exring=0,(x,y),ds;
3362 referencepoly(list(intvec(0,4),intvec(2,3),intvec(5,1),intvec(7,0)));
3363}
3364///////////////////////////////////////////////////////////////////////////////
3365
3366proc factorlist (list L)
3367"USAGE:   factorlist(L);   L a list in the format of `factorize'
3368RETURN:  the nonconstant irreducible factors of
3369         L[1][1]^L[2][1] * L[1][2]^L[2][2] *...* L[1][size(L[1])]^...
3370         with multiplicities (same format as factorize)
3371SEE ALSO: factorize
3372EXAMPLE: example factorlist;  shows an example
3373"
3374{
3375 // eine Sortierung der Faktoren eruebrigt sich, weil keine zwei versch.
3376 // red.Fakt. einen gleichen irred. Fakt. haben koennen (I.3.27 Diplarb.)
3377 int i,gross;
3378 list faktoren,hilf;
3379 ideal hil1,hil2;
3380 intvec v,w;
3381 for (i=1; (L[1][i] == jet(L[1][i],0)) && (i<size(L[1])); i++) {;}
3382 if (L[1][i] != jet(L[1][i],0)) {
3383   hilf=factorize(L[1][i]);
3384 // Achtung!!! factorize(..,2) wirft entgegen der Beschreibung nicht nur
3385 // konstante Faktoren raus, sondern alle Einheiten in der LOKALISIERUNG nach
3386 // der Monomordnung!!! Im Beispiel unten verschwindet der Faktor x+y+1, wenn
3387 // man ds statt dp als Ordnung nimmt!
3388   hilf[2]=hilf[2]*L[2][i];
3389   hil1=hilf[1];
3390   gross=size(hil1);
3391   if (gross>1) {
3392       // faktoren=list(hilf[1][2..gross],hilf[2][2..gross]);
3393       // -->  `? indexed object must have a name'
3394     v=hilf[2];
3395     faktoren=list(ideal(hil1[2..gross]),intvec(v[2..gross]));
3396   }
3397   else         { faktoren=hilf; }
3398 }
3399 else {
3400   faktoren=L;
3401 }
3402
3403 for (i++; i<=size(L[2]); i++) {
3404 //------------------------- linearer Term -- irreduzibel ---------------------
3405   if (L[1][i] == jet(L[1][i],1)) {
3406     if (L[1][i] != jet(L[1][i],0)) {           // konst. Faktoren eliminieren
3407       hil1=faktoren[1];
3408       hil1[size(hil1)+1]=L[1][i];
3409       faktoren[1]=hil1;
3410       v=faktoren[2],L[2][i];
3411       faktoren[2]=v;
3412     }
3413   }
3414 //----------------------- nichtlinearer Term -- faktorisiere -----------------
3415   else {
3416     hilf=factorize(L[1][i]);
3417     hilf[2]=hilf[2]*L[2][i];
3418     hil1=faktoren[1];
3419     hil2=hilf[1];
3420     gross=size(hil2);
3421       // hil2[1] ist konstant, wird weggelassen:
3422     hil1[(size(hil1)+1)..(size(hil1)+gross-1)]=hil2[2..gross];
3423       // ideal+ideal does not work due to simplification;
3424       // ideal,ideal not allowed
3425     faktoren[1]=hil1;
3426     w=hilf[2];
3427     v=faktoren[2],w[2..gross];
3428     faktoren[2]=v;
3429   }
3430 }
3431 return(faktoren);
3432}
3433example
3434{ "EXAMPLE:"; echo = 2;
3435 ring exring=0,(x,y),ds;
3436 list L=ideal(x,(x-y)^2*(x+y+1),x+y),intvec(2,2,1);
3437 L;
3438 factorlist(L);
3439}
3440///////////////////////////////////////////////////////////////////////////////
3441
3442proc extrafactor (poly leitf, int M, int N)
3443"USAGE:   factors,flag=extrafactor(f,M,N);
3444         list factors, int flag,M,N, poly f in x,y
3445RETURN:  for some special f, factors is a list of the factors of
3446         charPoly(f,M,N), same format as factorize
3447         (see help charPoly; help factorize)
3448         In this case, flag!=0 (TRUE). If extrafactor cannot factorize f,
3449         flag will be 0 (FALSE), factors will be the empty list.
3450ASSUME:  basering is (0,a),(x,y),ds or ls
3451         Newton polygon of f is one side with height N, length M
3452NOTE:    This procedure is designed to make reddevelop able to compute some
3453         special cases that need a ring extension in char 0.
3454         It becomes obsolete as soon as `factorize' works also in such rings.
3455EXAMPLE: example extrafactor;  shows an example
3456"
3457{
3458 list faktoren;
3459
3460      if (a2==-1) {
3461       poly testpol=charPoly(leitf,M,N);
3462       if (testpol==1+2y2+y4) {
3463         faktoren=list(ideal(y+a,y-a),intvec(2,2));
3464         testpol=0;
3465       }
3466 //------------ ist Poly von der Form q*i+r*y^n, n in N, q,r in Q ?: ----------
3467       if ((size(testpol)==2) && (find(string(lead(testpol)),"a")>0)
3468           && (find(string(testpol-lead(testpol)),"a")==0)) {
3469        faktoren=list(ideal(testpol),intvec(1));
3470        testpol=0;
3471       }
3472      }
3473      if (a4==-625) {
3474        poly testpol=charPoly(leitf,M,N);
3475        if (testpol==4a2-4y2)
3476         { faktoren=list(ideal(-4,y+a,y-a),intvec(1,1,1)); testpol=0;}
3477        if (testpol==-4a2-4y2)
3478         { faktoren=list(ideal(-4,y+1/25*a3,y-1/25*a3),intvec(1,1,1));
3479           testpol=0;}
3480      }
3481      if (a2+a==-1) {
3482        poly testpol=charPoly(leitf,M,N);
3483        if (testpol==1+2y+3y2+2y3+y4)
3484         { faktoren=list(ideal(1,y-a,y+a+1),intvec(1,2,2)); testpol=0;}
3485      }
3486      if (defined(testpol)==0) { poly testpol=1; }
3487      if (testpol!=0) {
3488        "factorize not implemented in char (0,a)!";
3489        "could not factorize:",charPoly(leitf,M,N);
3490        if (printlevel>0) { pause("Hit RETURN to continue:"); }
3491      }
3492 return(faktoren,(testpol==0)); // Test: faktoren==list() geht leider nicht
3493}
3494example
3495{ "EXAMPLE:"; echo=2;
3496  ring r=(0,a),(x,y),ds;
3497  minpoly=a2+1;
3498  poly f=x4+2x2y2+y4;
3499  charPoly(f,4,4);
3500  list factors;
3501  int flag;
3502  factors,flag=extrafactor(f,4,4);
3503  if (flag!=0) {factors;}
3504}
Note: See TracBrowser for help on using the repository browser.