source: git/Singular/LIB/hnoether.lib @ 63be42

spielwiese
Last change on this file since 63be42 was fff61a7, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: modified info-string git-svn-id: file:///usr/local/Singular/svn/trunk@1999 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 99.5 KB
Line 
1// $Id: hnoether.lib,v 1.8 1998-05-27 10:54:24 Singular Exp $
2// author:  Martin Lamm,  email: lamm@mathematik.uni-kl.de
3// last change:           26.03.98
4///////////////////////////////////////////////////////////////////////////////
5
6version="$Id: hnoether.lib,v 1.8 1998-05-27 10:54:24 Singular Exp $";
7info="
8LIBRARY:  hnoether.lib   PROCEDURES FOR THE HAMBURGER-NOETHER-DEVELOPMENT
9
10           Important procedures:
11 develop(f [,n]);    Hamburger-Noether development from the irred. polynomial f
12 reddevelop(f);      Hamburger-Noether development from the (red.) polynomial f
13 extdevelop(hne,n);  extension of Hamburger-Noether development hne from f
14 param(hne);         returns a parametrization of f (input=output(develop))
15 displayHNE(hne);    display Hamburger-Noether development as an ideal
16 invariants(hne);    invariants of f, e.g. the characteristic exponents
17 displayInvariants(hne);  display invariants of f
18 generators(hne);    computes the generators of the semigroup of values
19 intersection(hne1,hne2); intersection multiplicity of two curves
20 stripHNE(hne);      reduce amount of memory consumed by hne
21";
22////////////////////////////////////////////////////////////////////////////
23//           perhaps useful procedures:
24// puiseux2generators(m,n); convert puiseux pairs to generators of semigroup
25// multiplicities(hne);     multiplicities of blowed up curves
26// newtonpoly(f);      newtonpolygon of polynom f
27// newtonhoehne(f);    same as newtonpoly, but uses internal procedure
28// getnm(f);           intersection points of the Newton-polygon with the axes
29// T_Transform(f,Q,N); returns f(y,xy^Q)/y^NQ (type f(x,y): poly, Q,N: int)
30// T1_Transform(f,d,M); returns f(x,y+d*x^M)  (type f(x,y): poly,d:number,M:int)
31// T2_Transform(f,d,M,N);   a composition of T1 & T
32// koeff(f,I,J);       gets coefficient of indicated monomial of poly f (I,J:int)
33// redleit(f,S,E);     restriction of monomials of f to line (S-E)
34// squarefree(f);      returns a squarefree divisor of the poly f
35// allsquarefree(f,l); returns the maximal squarefree divisor of the poly f
36// set_list(L,v,l);    managing interconnected lists
37//
38//           procedures (more or less) for internal use:
39// leit(f,n,m);        special case of redleit (for irred. polynomials)
40// testreducible(f,n,m); tests whether f is reducible
41// polytest(f);        tests coefficients and exponents of poly f
42// charPoly(f,M,N);    characteristic polynomial of f
43// find_in_list(L,p);  find int p in list L
44// get_last_divisor(M,N); last divisor in Euclid's algorithm
45// factorfirst(f,M,N); try to factor f in a trivial way before using 'factorize'
46// extrafactor(f,M,N); try to factor charPoly(f,M,N) where 'factorize' cannot
47// extractHNEs(H,t);   extracts output H of HN to output of reddevelop
48// HN(f,grenze);       recursive subroutine for reddevelop
49// constructHNEs(...); subroutine for HN
50///////////////////////////////////////////////////////////////////////////////
51LIB "primitiv.lib";
52///////////////////////////////////////////////////////////////////////////////
53
54proc getnm (poly f)
55"USAGE:   getnm(f); f polynomial in x,y
56ASSUME:  basering = ...,(x,y),ls;  or ds
57RETURN:  intvec(n,m) : (0,n) is the intersection point of the Newton
58         polygon of f with the y-axis, n=-1 if it doesn't exist
59         (m,0) is its intersection point with the x-axis,
60         m=-1 if this point doesn't exist
61EXAMPLE: example getnm; shows an example
62"
63{
64 // wir gehen davon aus, dass die Prozedur von develop aufgerufen wurde, also
65 // die Ringordnung ist ls (ds ginge auch)
66 return(ord(subst(f,x,0)),ord(subst(f,y,0)));
67}
68example
69{ "EXAMPLE:"; echo = 2;
70   ring r = 0,(x,y),ds;
71   poly f = x5+x4y3-y2+y4;
72   getnm(f);
73}
74///////////////////////////////////////////////////////////////////////////////
75
76proc leit (poly f, int n, int m)
77"// leit(poly f, int n, int m) returns all monomials on the line
78// from (0,n) to (m,0) in the Newton diagram
79"
80{
81 return(jet(f,m*n,intvec(n,m))-jet(f,m*n-1,intvec(n,m)))
82}
83///////////////////////////////////////////////////////////////////////////////
84proc testreducible (poly f, int n, int m)
85"// testreducible(poly f,int n,int m) returns true if there are points in the
86// Newton diagram below the line (0,n)-(m,0)
87"
88{
89 return(size(jet(f,m*n-1,intvec(n,m))) != 0)
90}
91///////////////////////////////////////////////////////////////////////////////
92proc T_Transform (poly f, int Q, int N)
93"// returns f(y,xy^Q)/y^NQ
94"
95{
96 map T = basering,y,x*y^Q;
97 return(T(f)/y^(N*Q));
98}
99///////////////////////////////////////////////////////////////////////////////
100proc T1_Transform (poly f, number d, int M)
101"// returns f(x,y+d*x^M)
102"
103{
104 map T1 = basering,x,y+d*x^M;
105 return(T1(f));
106}
107///////////////////////////////////////////////////////////////////////////////
108
109proc T2_Transform (poly f, number d, int M, int N)
110"USAGE: T2_Transform(f,d,M,N); f poly, d number; M,N int
111RETURN: list: poly T2(f,d',M,N), number d' in \{ d, 1/d \}
112"
113{
114 //---------------------- compute gcd and extgcd of N,M -----------------------
115  int ggt=gcd(M,N);
116  M=M/ggt; N=N/ggt;
117  list ts=extgcd(M,N);
118  int tau,sigma=ts[2],-ts[3];
119  if (sigma<0) { tau=-tau; sigma=-sigma;}
120 // es gilt: 0<=tau<=N, 0<=sigma<=M, |N*sigma-M*tau| = 1 = ggT(M,N)
121  if (N*sigma < M*tau) { d = 1/d; }
122 //--------------------------- euklid. Algorithmus ----------------------------
123  int R;
124  int M1,N1=M,N;
125  for ( R=M1%N1; R!=0; ) { M1=N1; N1=R; R=M1%N1;}
126  int Q=M1 / N1;
127  map T1 = basering,x,y+d*x^Q;
128  map Tstar=basering,x^(N-Q*tau)*y^tau,x^(M-sigma*Q)*y^sigma;
129  if (defined(Protokoll)) {
130   "Trafo. T2: x->x^"+string(N-Q*tau)+"*y^"+string(tau)+", y->x^"
131    +string(M-sigma*Q)+"*y^"+string(sigma);
132   "delta =",d,"Q =",Q,"tau,sigma =",tau,sigma;
133  }
134 //------------------- Durchfuehrung der Transformation T2 --------------------
135 // man koennte an T2_Transform noch eine Liste mit allen Eckpunkten des
136 // Newtonpolys uebergeben und nur die Transformierten der Eckpunkte betrachten
137 // um max. Exponenten e,f von x,y zu finden--noch einfacher: ein Referenzpo-
138 // lynom, das das gleiche Newtonpolygon hat (z.B. Einschraenkung auf die Mono-
139 // me in den Eckpunkten, oder alle Monome haben Koeff. 1, das ist einfacher)
140 //----------------------------------------------------------------------------
141
142  f=Tstar(f); // hier dann refpoly=Tstar(refpoly)
143  poly hilf;
144 // dividiere f so lange durch x, wie die Div. aufgeht:
145  for (hilf=f/x; hilf*x==f; hilf=f/x) {f=hilf;}
146  for (hilf=f/y; hilf*y==f; hilf=f/y) {f=hilf;} // gleiches fuer y
147  return(list(T1(f),d));
148}
149///////////////////////////////////////////////////////////////////////////////
150
151proc koeff (poly f, int I, int J)
152"USAGE:  koeff(f,I,J); f polynomial in x and y, I,J integers
153RETURN: if f = sum(a(i,j)*x^i*y^j), then koeff(f,I,J)= a(I,J) (of type number)
154EXAMPLE: example koeff;  shows an example
155"{
156  matrix mat = coeffs(coeffs(f,y)[J+1,1],x);
157  if (size(mat) <= I) { return(0);}
158  else { return(leadcoef(mat[I+1,1]));}
159}
160example
161{ "EXAMPLE:"; echo = 2;
162  ring r=0,(x,y),dp;
163  koeff(x2+2xy+3xy2-x2y-2y3,1,2);
164}
165///////////////////////////////////////////////////////////////////////////////
166
167proc squarefree (poly f)
168"USAGE:  squarefree(f); f poly in x and y
169RETURN: a squarefree divisor of f
170        Normally the return value is a greatest squarefree divisor, but there
171        is an exception: factors with a p-th root, p the ring characteristic,
172        are lost.
173EXAMPLE: example squarefree; shows some examples.
174"{
175 // es gibt noch ein Problem: syz gibt manchmal in Koerpererweiterung (p,a)
176 // nicht f/g zurueck, obwohl die Division aufgeht ... In dem Fall ist die
177 // Rueckgabe evtl. nicht quadratfrei! Bsp.: squarefree((1/a*x3+y7)^3);
178 // in char 7
179 // In Singular 1.1.3 funktioniert syz, dafuer geht das gleiche Bsp. schief,
180 // weil gcd auf einmal falsch wird! (ring (7,a),(x,y),dp; minpoly=a2+1; )
181
182 //----------------- Wechsel in geeigneten Ring & Variablendefinition ---------
183  def altring = basering;
184  int e;
185  int gcd_ok=1;
186  string mipl="0";
187  if (size(parstr(altring))==1) {mipl=string(minpoly);}
188 //---- test: char = (p^k,a) (-> gcd not implemented) or (p,a) (gcd works) ----
189  if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) {
190    string tststr=charstr(basering);
191    tststr=tststr[1..find(tststr,",")-1];           //-> "p^k" bzw. "p"
192    gcd_ok=(tststr==string(char(basering)));
193  }
194  execute("ring rsqrf = ("+charstr(altring)+"),(x,y),dp;");
195  if ((gcd_ok!=0) && (mipl!="0")) { execute "minpoly="+mipl+";"; }
196  poly f=fetch(altring,f);
197  poly dif,g,l;
198  if (gcd_ok!=0) {
199 //-------------------- Berechne f/ggT(f,df/dx,df/dy) ------------------------
200    dif=diff(f,x);
201    if (dif==0) { g=f; }        // zur Beschleunigung
202    else { g=gcd(f,dif); }
203    if (g!=1) {                 // sonst schon sicher, dass f quadratfrei
204     dif=diff(f,y);
205     if (dif!=0) { g=gcd(g,dif); }
206    }
207    if (g!=1) {
208     e=0;
209     if (g==f) { l=1; }         // zur Beschleunigung
210     else {
211       module m=syz(ideal(g,f));
212       if (deg(m[2,1])>0) {
213         "!! The Singular command 'syz' has returned a wrong result !!";
214         l=1;                   // Division f/g muss aufgehen
215       }
216       else { l=m[1,1]; }
217     }
218 //----------------------------------------------------------------------------
219 // Polynomdivision geht, wenn factory eingebunden ist
220 // dann also l=f/g; aber: ist ebenfalls fehlerhaft! (s.o. als Bsp.)
221 // wenn / zur Polynomdivision einsetzbar ist, kann der Ringwechsel gespart
222 // werden, sollte dann also zum Einsatz gebracht werden (Zeitersparnis)
223 //----------------------------------------------------------------------------
224    }
225    else { e=1; }
226  }
227  else {
228 //------------------- Berechne syz(f,df/dx) oder syz(f,df/dy) ----------------
229 //-- Achtung: Ist f reduzibel, koennen Faktoren mit Ableitung Null verloren --
230 //-- gehen! Ist aber nicht weiter schlimm, weil char (p^k,a) nur im irred.  --
231 //-- Fall vorkommen kann. Wenn f nicht g^p ist, wird auf jeden Fall         --
232 //------------------------ ein Faktor gefunden. ------------------------------
233    dif=diff(f,x);
234    if (dif == 0) {
235     dif=diff(f,y);
236     if (dif==0) { e=2; l=1; } // f is of power divisible by char of basefield
237     else { l=syz(ideal(dif,f))[1,1];  // x^p+y^(p-1) abgedeckt
238            if (deg(l)==deg(f)) { e=1;}
239            else {e=0;}
240     }
241    }
242    else { l=syz(ideal(dif,f))[1,1];
243           if (deg(l)==deg(f)) { e=1;}
244           else {e=0;}
245    }
246  }
247 //--------------- Wechsel in alten Ring und Rueckgabe des Ergebnisses --------
248  setring altring;
249  if (e==1) { return(f); }    // zur Beschleunigung
250  else {
251   poly l=fetch(rsqrf,l);
252   return(l);
253  }
254}
255example
256{ "EXAMPLE:"; echo = 2;
257 ring exring=3,(x,y),dp;
258 squarefree(x3+y);
259 squarefree(x3);
260 squarefree(x2);
261 squarefree(xy+x);
262 squarefree((x+y)^3*(x-y)^2); // (x+y)^3 is lost
263 squarefree((x+y)^4*(x-y)^2); // result is (x+y)*(x-y)
264}
265///////////////////////////////////////////////////////////////////////////////
266
267proc allsquarefree (poly f, poly l)
268"USAGE : allsquarefree(f,l);  poly f,l
269        l is the squarefree divisor of f you get by l=squarefree(f);
270RETURN: the squarefree divisor of f consisting of all irreducible factors of f
271NOTE  : * this proc uses factorize to get the missing factors of f not in l and
272          therefore may be slow. Furthermore factorize can still return
273          a factor more than once (what should not occur)
274        * this proc uses polynomial division of factory (which still has a bug
275         at the moment: it can return nonsense, if the basering has parameters)
276EXAMPLE: example allsquarefree;  shows an example
277"{
278 //---------- eliminiere bereits mit squarefree gefundene Faktoren ------------
279 poly g=l;                            // g = gcd(f,l);
280 while (deg(g)!=0) {
281  f=f/g;
282  g=gcd(f,l);
283 }                                    // jetzt f=h^p
284 //--------------- Berechne uebrige Faktoren mit factorize --------------------
285 if (deg(f)>0) {
286  g=1;
287  ideal factf=factorize(f,1);
288  for (int i=1; i<=size(factf); i++) { g=g*factf[i]; }
289  poly testp=squarefree(g);
290  if (deg(testp)<deg(g)) {
291    "!! factorize has not worked correctly !!";
292    if (testp==1) {" We cannot proceed ..."; g=1;}
293    else {" But we could recover the correct result..."; g=testp;}
294  }
295  l=l*g;
296 }
297 return(l);
298}
299example
300{ "EXAMPLE:"; echo = 2;
301  ring exring=7,(x,y),dp;
302  poly f=(x+y)^7*(x-y)^8;
303  poly l=squarefree(f);
304  l; // x-y found because 7 does not divide 8
305  allsquarefree(f,l); // all factors (x+y)*(x-y) found
306}
307///////////////////////////////////////////////////////////////////////////////
308
309proc polytest(poly f)
310"USAGE : polytest(f); f poly in x and y
311RETURN: a monomial of f with |coefficient| > 16001
312          or exponent divisible by 32003, if there is one
313        0 else (in this case computing a squarefree divisor
314                in characteristic 32003 could make sense)
315NOTE:   this procedure is only useful in characteristic zero, because otherwise
316        there is no appropriate ordering of the leading coefficients
317"{
318 poly verbrecher=0;
319 intvec leitexp;
320 for (; (f<>0) and (verbrecher==0); f=f-lead(f)) {
321  if ((leadcoef(f)<-16001) or (leadcoef(f)>16001)) {verbrecher=lead(f);}
322  leitexp=leadexp(f);
323  if (( ((leitexp[1] % 32003) == 0)   and (leitexp[1]<>0))
324     or ( ((leitexp[2] % 32003) == 0) and (leitexp[2]<>0)) )
325       {verbrecher=lead(f);}
326 }
327 return(verbrecher);
328}
329
330//////////////////////////////////////////////////////////////////////////////
331
332
333proc develop
334"USAGE:   develop(f [,n]); f polynomial in two variables, n integer
335ASSUME:  f irreducible as a power series (for reducible f use reddevelop)
336RETURN:  Hamburger-Noether development of f, in form of a list:
337         [1]: Hamburger-Noether matrix
338            each row contains the coefficients of the corresponding line of the
339            Hamburger-Noether extension (HNE). The end of the line is marked in
340            the matrix as the first ring variable (usually x)
341         [2]: intvec indicating the length of lines of the HNE
342         [3]: int:
343           0  if the 1st ring variable was transverse (with respect to f),
344           1  if the variables were changed at the beginning of the computation
345           -1 if an error has occurred
346         [4]: the transformed polynomial of f to make it possible to extend the
347              Hamburger-Noether development a posteriori without having to do
348              all the previous calculation once again (0 if not needed)
349DISPLAY: The (non zero) elements of the HNE
350
351NOTE:    If the optional parameter n is given, the HN-matrix will have at least
352         n columns. Otherwise the number of columns will be chosen minimal s.t.
353         the matrix contains all necessary information (i.e. all lines of the
354         HNE but the last (which is in general infinite) appear).
355         If n is negative, the algorithm is stopped as soon as possible, i.e.
356         the information computed is enough for 'invariants', but the HN-matrix
357         may contain undetermined elements, which are marked with the
358         2nd variable (of the basering).
359         In any case, the optional parameter only affects the calculation of
360         the LAST line of the HNE; develop(f) gives already all necessary
361         information for the procedure 'invariants'. A negative value of n will
362         result in a very poor parametrization, but it can make 'develop'
363         faster; a positive value will improve the exactness of the
364         parametrization.
365
366         For time critical computations it is recommended to use
367         \"ring ...,(x,y),ls\" as basering - it increases the algorithm's speed.
368
369EXAMPLES: example develop; shows an example
370          example param;   shows an example for using the 2nd parameter
371"{
372 //--------- Abfangen unzulaessiger Ringe: 1) nur eine Unbestimmte ------------
373 poly f=#[1];
374 if (size(#) > 1) {int maxspalte=#[2];}
375 else             {int maxspalte= 1 ; }
376 if (nvars(basering) < 2) {
377   " Sorry. I need two variables in the ring.";
378   return(list(matrix(maxideal(1)[1]),intvec(0),-1,poly(0)));}
379 if (nvars(basering) > 2) {
380   " Warning! You have defined too many variables!";
381   " All variables except the first two will be ignored!";}
382
383 string namex=varstr(1); string namey=varstr(2);
384 list return_error=matrix(maxideal(1)[2]),intvec(0),int(-1),poly(0);
385
386 //------------- 2) mehrere Unbestimmte, weitere unzulaessige Ringe -----------
387 // Wir koennen einheitlichen Rueckgabewert benutzen, aus dem ersichtlich ist,
388 // dass ein Fehler aufgetreten ist: return_error.
389 //----------------------------------------------------------------------------
390
391 if (charstr(basering)=="real") {
392  "The algorithm doesn't work with 'real' as coefficient field.";
393                     // denn : map from characteristic -1 to -1 not implemented
394  return(return_error);
395 }
396 if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) {
397 //-- teste, ob char = (p^k,a) (-> a primitiv; erlaubt) oder (p,a[,b,...]) ----
398    string tststr=charstr(basering);
399    tststr=tststr[1..find(tststr,",")-1];           //-> "p^k" bzw. "p"
400    int primit=(tststr==string(char(basering)));
401    if (primit!=0) {
402      "such extensions of Z/p are not implemented.",
403      "Please try (p^k,a) as ground field.";
404      return(return_error);
405    }
406 }
407 //---- Ende der unzulaessigen Ringe; Ringwechsel in einen guenstigen Ring: ---
408
409 int ringwechsel=(varstr(basering)!="x,y") or (ordstr(basering)!="ls(2),C");
410
411 def altring = basering;
412 if (ringwechsel) {
413   string mipl=string(minpoly);
414   execute("ring guenstig = ("+charstr(altring)+"),(x,y),ls;");
415   if ((char(basering)==0) && (mipl!="0")) {
416     execute "minpoly="+mipl+";";
417   }}
418 else { def guenstig=basering; }
419 export guenstig;
420
421 //-------------------------- Initialisierungen -------------------------------
422 map m=altring,x,y;
423 if (ringwechsel) { poly f=m(f); }
424 if (defined(Protokoll))
425 {"received polynomial: ",f,", where x =",namex,", y =",namey;}
426
427 int M,N,Q,R,l,e,hilf,eps,getauscht,Abbruch,zeile,exponent,Ausgabe;
428
429 // Werte von Ausgabe: 0 : normale HNE-Matrix,
430 // 1 : Fehler aufgetreten - Matrix (namey) zurueck
431 // 2 : Die HNE ist eine Nullzeile - Matrix (0) zurueck
432 // int maxspalte=1; geaendert: wird jetzt am Anfang gesetzt
433
434 int minimalHNE=0;          // Flag fuer minimale HNE-Berechnung
435 intvec hqs;                // erhaelt die Werte von h(zeile)=Q;
436
437 if (maxspalte<0) {
438   minimalHNE=1;
439   maxspalte=1;
440 }
441
442 number c,delta;
443 int p = char(basering);
444 string ringchar=charstr(basering);
445 map xytausch = basering,y,x;
446 if ((p!=0) and (ringchar != string(p))) {
447                            // coefficient field is extension of Z/pZ
448   execute "int n_elements="+ringchar[1,size(ringchar)-2]+";";
449                            // number of elements of actual ring
450   number generat=par(1);   // generator of the coefficient field of the ring
451 }
452
453
454 //========= Abfangen von unzulaessigen oder trivialen Eingaben ===============
455 //------------ Nullpolynom oder Einheit im Potenzreihenring: -----------------
456 if (f == 0) {"You have given me the zero-polynomial !"; Abbruch=1; Ausgabe=1;}
457 else {
458   intvec nm = getnm(f);
459   N = nm[1]; M = nm[2]; //Berechne Schnittpunkte Newtonpolygon mit Achsen
460   if (N == 0)
461   {" The given polynomial is a unit as power series!!"; Abbruch=1; Ausgabe=1;}
462   else {
463    if (N == -1) {"The HNE is x = 0"; Abbruch=1; Ausgabe=2; getauscht=1;}
464    else {
465     if (M == -1) {"The HNE is y = 0"; Abbruch=1; Ausgabe=2;}
466   }}
467   if (N+M==-2) {"Remark: The polynomial is divisible both by x and y,",
468                 "so it is not irreducible.";
469                 "The result of develop is only one of the existing HNE's.";
470   }
471 }
472 //--------------------- Test auf Quadratfreiheit -----------------------------
473 if (Abbruch==0) {
474
475 //-------- Fall basering==0,... : Wechsel in Ring mit char >0 ----------------
476 // weil squarefree eine Standardbasis berechnen muss (verwendet Syzygien)
477 // -- wenn f in diesem Ring quadratfrei ist, dann erst recht im Ring guenstig
478 //----------------------------------------------------------------------------
479
480  if ((p==0) and (size(charstr(basering))==1)) {
481   int testerg=(polytest(f)==0);
482   ring zweitring = 32003,(x,y),dp;
483   map polyhinueber=guenstig,x,y;     // fetch geht nicht
484   poly f=polyhinueber(f);
485   poly test_sqr=squarefree(f);
486   if (test_sqr != f) {
487    "Most probably the given polynomial is not squarefree. But the test was";
488    "made in characteristic 32003 and not 0 to improve speed. You can";
489    "(r) redo the test in char 0 (but this may take some time)";
490    "(c) continue the development, if you're sure that the polynomial",
491    "IS squarefree";
492    if (testerg==1) {
493      "(s) continue the development with a squarefree factor (*)";}
494    "(q) or just quit the algorithm (default action)";
495    "";"Please enter the letter of your choice:";
496    string str=read("")[1];
497    setring guenstig;
498    map polyhinueber=zweitring,x,y;
499    if (str=="r") {
500      poly test_sqr=squarefree(f);
501      if (test_sqr != f) {
502       "The given polynomial is in fact not squarefree.";
503       "I'll continue with the radical.";
504       pause;
505       f=test_sqr;
506      }
507      else {"everything is ok -- the polynomial is squarefree in char(k)=0";}
508    }
509    else {
510      if ((str=="s") and (testerg==1)) {
511       "(*) attention: it could be that the factor is only one in char 32003!";
512        f=polyhinueber(test_sqr);
513      }
514      else {
515        if (str<>"c") {
516          setring altring;kill guenstig;kill zweitring;
517          return(return_error);}
518        else { "if the algorithm doesn't terminate, you were wrong...";}
519    }}
520    kill zweitring;
521    nm = getnm(f);             // N,M haben sich evtl. veraendert
522    N = nm[1]; M = nm[2];      // Berechne Schnittpunkte Newtonpoly mit Achsen
523    if (defined(Protokoll)) {"I continue with the polynomial",f; }
524   }
525   else {
526     setring guenstig;
527     kill zweitring;
528   }
529  }
530 // ------------------- Fall Charakteristik > 0 -------------------------------
531  else {
532   poly test_sqr=squarefree(f);
533   if (test_sqr == 1) {
534    "The given polynomial is of the form g^"+string(p)+", therefore",
535    "reducible.";"Please try again.";
536    setring altring;
537    kill guenstig;
538    return(return_error);}
539   if (test_sqr != f) {
540    "The given polynomial is not squarefree. I'll continue with the radical.";
541    if (p != 0)
542     {"But if the polynomial contains a factor of the form g^"+string(p)+",";
543      "this factor will be lost.";}
544    pause;
545    f=test_sqr;
546    nm = getnm(f);              // N,M haben sich veraendert
547    N = nm[1]; M = nm[2];       // Berechne Schnittpunkte Newtonpoly mit Achsen
548    if (defined(Protokoll)) {"I continue with the polynomial",f; }
549   }
550
551  }                             // endelse(p==0)
552
553  if (N==0) {
554    " Sorry. The remaining polynomial is a unit in the power series ring...";
555    setring altring;kill guenstig;return(return_error);
556  }
557 //---------------------- gewaehrleiste, dass x transvers ist -----------------
558  if (M < N)
559  { f = xytausch(f);            // Variablentausch : x jetzt transvers
560    getauscht = 1;              // den Tausch merken
561    M = M+N; N = M-N; M = M-N;  // M, N auch vertauschen
562  }
563  if (defined(Protokoll)) {
564   if (getauscht) {"x<->y were exchanged; poly is now ",f;}
565   else           {"x , y were not exchanged";}
566   "M resp. N are now",M,N;
567  }
568 }                              // end(if Abbruch==0)
569
570 ideal a(0);
571 while (Abbruch==0) {
572
573 //================= Beginn der Schleife (eigentliche Entwicklung) ============
574
575 //------------------- ist das Newtonpolygon eine gerade Linie? ---------------
576  if (testreducible(f,N,M)) {
577    " The given polynomial is not irreducible";
578    kill guenstig;
579    setring altring;
580    return(return_error);       // Abbruch der Prozedur!
581  }
582  R = M%N;
583  Q = M / N;
584
585 //-------------------- Fall Rest der Division R = 0 : ------------------------
586  if (R == 0) {
587    c = koeff(f,0,N);
588    if (c == 0) {"Something has gone wrong! I didn't get N correctly!"; exit;}
589    e = gcd(M,N);
590 //----------------- Test, ob leitf = c*(y^N - delta*x^(m/e))^e ist -----------
591    if (p==0) {
592      delta = koeff(f,M/ e,N - N/ e) / (-1*e*c);
593      if (defined(Protokoll)) {"quasihomogeneous leading form:",
594         leit(f,N,M)," = ",c,"* (y -",delta,"* x^"+string(M/ e)+")^",e," ?";}
595      if (leit(f,N,M) != c*(y^(N/ e) - delta*x^(M/ e))^e)
596        { " The given Polynomial is reducible !"; Abbruch=1; Ausgabe=1;}
597    }
598    else {                     // p!=0
599      if (e%p != 0) {
600        delta = koeff(f,M/ e,N - N/ e) / (-1*e*c);
601        if (defined(Protokoll)) {"quasihomogeneous leading form:",
602           leit(f,N,M)," = ",c,"* (y -",delta,"* x^"+string(M/ e)+")^",e," ?";}
603        if (leit(f,N,M) != c*(y^(N/ e) - delta*x^(M/ e))^e)
604         { " The given Polynomial is reducible !"; Abbruch=1; Ausgabe=1;}
605      }
606
607      else {                   // e%p == 0
608        eps = e;
609        for (l = 0; eps%p == 0; l=l+1) { eps=eps/ p;}
610        if (defined(Protokoll)) {e," -> ",eps,"*",p,"^",l;}
611        delta = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c);
612
613        if ((ringchar != string(p)) and (delta != 0)) {
614 //- coeff. field is not Z/pZ => we`ve to correct delta by taking (p^l)th root-
615          if (delta == generat) {exponent=1;}
616          else {
617           if (delta == 1) {exponent=0;}
618           else {
619            exponent=pardeg(delta);
620
621 //-- an dieser Stelle kann ein Fehler auftreten, wenn wir eine transzendente -
622 //-- Erweiterung von Z/pZ haben: dann ist das hinzuadjungierte Element kein  -
623 //-- Erzeuger der mult. Gruppe, d.h. in Z/pZ (a) gibt es i.allg. keinen      -
624 //-- Exponenten mit z.B. a2+a = a^exp                                        -
625 //----------------------------------------------------------------------------
626          }}
627          delta = generat^(extgcd(n_elements-1,p^l)[3]*exponent);
628        }
629
630        if (defined(Protokoll)) {"quasihomogeneous leading form:",
631          leit(f,N,M)," = ",c,"* (y^"+string(N/ e),"-",delta,"* x^"
632          +string(M/ e)+")^",e,"  ?";}
633        if (leit(f,N,M) != c*(y^(N/ e) - delta*x^(M/ e))^e)
634          { " The given Polynomial is reducible !"; Abbruch=1; Ausgabe=1;}
635      }
636    }
637    if (Abbruch == 0) {
638      f = T1_Transform(f,delta,M/ e);
639      "a("+string(zeile)+","+string(Q)+") =",delta;
640      a(zeile)[Q]=delta;
641      if (defined(Protokoll)) {"transformed polynomial: ",f;}}
642
643      nm=getnm(f); N=nm[1]; M=nm[2];        // Neuberechnung des Newtonpolygons
644  }
645 //--------------------------- Fall R > 0 : -----------------------------------
646  else {
647    "h("+string(zeile)+") =",Q;
648    hqs[zeile+1]=Q;                  // denn zeile beginnt mit dem Wert 0
649    a(zeile)[Q+1]=x;                 // Markierung des Zeilenendes der HNE
650    maxspalte=maxspalte*((Q+1) < maxspalte) + (Q+1)*((Q+1) >= maxspalte);
651                                     // Anpassung der Sp.zahl der HNE-Matrix
652    f = T_Transform(f,Q,N);
653    if (defined(Protokoll)) {"transformed polynomial: ",f;}
654    zeile=zeile+1;
655 //------------ Bereitstellung von Speicherplatz fuer eine neue Zeile: --------
656    ideal a(zeile);
657    M=N;N=R;
658  }
659
660 //--------------- schneidet das Newtonpolygon beide Achsen? ------------------
661  if (M==-1) {
662     "The HNE is finite!";
663     a(zeile)[Q+1]=x;   // Markiere das Ende der Zeile
664     hqs[zeile+1]=Q;
665     maxspalte=maxspalte*((Q+1) < maxspalte) + (Q+1)*((Q+1) >= maxspalte);
666     f=0;               // transformiertes Polynom wird nicht mehr gebraucht
667     Abbruch=1;
668  }
669  else {if (M<N) {"Something has gone wrong: M<N";}}
670  if(defined(Protokoll)) {"new M,N:",M,N;}
671
672 //----------------- Abbruchbedingungen fuer die Schleife: --------------------
673  if ((N==1) and (Abbruch!=1) and ((M > maxspalte) or (minimalHNE==1))
674      and (size(a(zeile))>0))
675 //----------------------------------------------------------------------------
676 // Abbruch, wenn die Matrix so voll ist, dass eine neue Spalte angefangen
677 // werden muesste und die letzte Zeile nicht nur Nullen enthaelt
678 // oder wenn die Matrix nicht voll gemacht werden soll (minimale Information)
679 //----------------------------------------------------------------------------
680   { Abbruch=1; hqs[zeile+1]=-1;
681     if (maxspalte < ncols(a(zeile))) { maxspalte=ncols(a(zeile));}
682     if ((minimalHNE==1) and (M <= maxspalte)) {
683 // teile param mit, dass Eintraege der letzten Zeile nur teilw. richtig sind:-
684       hqs[zeile+1]=-M;
685 //------------- markiere den Rest der Zeile als unbekannt: -------------------
686       for (R=M; R <= maxspalte; R++) { a(zeile)[R]=y;}
687     }                  // R wird nicht mehr gebraucht
688   }
689 //========================= Ende der Schleife ================================
690
691 }
692 setring altring;
693 if (Ausgabe == 0) {
694 //-------------------- Ergebnis in den alten Ring transferieren: -------------
695   map zurueck=guenstig,maxideal(1)[1],maxideal(1)[2];
696   matrix amat[zeile+1][maxspalte];
697   ideal uebergabe;
698   for (e=0; e<=zeile; e=e+1) {
699     uebergabe=zurueck(a(e));
700     if (ncols(uebergabe) > 1) {
701      amat[e+1,1..ncols(uebergabe)]=uebergabe;}
702     else {amat[e+1,1]=uebergabe[1];}
703   }
704   if (ringwechsel) { f=zurueck(f); }
705 }
706
707 kill guenstig;
708 if (Ausgabe == 0) { return(list(amat,hqs,getauscht,f));}
709 if (Ausgabe == 1) { return(return_error);}             // error has occurred
710 if (Ausgabe == 2) { return(list(matrix(ideal(0,x)),intvec(1),getauscht,
711                                 poly(0)));}            // HNE is x=0 or y=0
712}
713example
714{ "EXAMPLE:"; echo = 2;
715  ring exring = 7,(x,y),ds;
716  list hne=develop(4x98+2x49y7+x11y14+2y14);
717  print(hne[1]);
718  // therefore the HNE is:
719  // z(-1)= 3*z(0)^7 + z(0)^7*z(1),
720  // z(0) = z(1)*z(2),   (note that there is  1 zero   in the 2nd row before x)
721  // z(1) = z(2)^3*z(3), (note that there are 3 zeroes in the 3rd row before x)
722  // z(2) = z(3)*z(4),
723  // z(3) = -z(4)^2 + 0*z(4)^3 +...+ 0*z(4)^8 + ?*z(4)^9 + ...
724  //   (the missing x in the matrix indicates that this line is not complete.
725  //    It can only occur in the last line of the HNE, and normally does.)
726  hne[2]; pause;
727  param(hne);
728  // returns the parametrisation x(t)= -t14+O(t21), y(t)= -3t98+O(t105)
729  // (the term -t109 in y may have a wrong coefficient)
730  displayHNE(hne);
731}
732
733///////////////////////////////////////////////////////////////////////////////
734
735proc param
736"USAGE:  param(l) takes the output of develop(f)
737        (list l (matrix m, intvec v, int s[,poly g]))
738        and gives a parametrisation for f; the first variable of the active
739        ring is chosen as indefinite. If the ring contains more than
740        two variables, the 3rd variable is chosen (remember that develop takes
741        the first two variables and therefore the other vars should be unused)
742RETURN: ideal of two polynomials: if the HNE was finite, f(param[1],param[2])
743        will be  zero. If not, the real parametrisation will be
744        two power series; then param will return a truncation of these series.
745EXAMPLE: example param;      shows an example
746         example develop;   shows another example
747
748"{
749 //-------------------------- Initialisierungen -------------------------------
750 matrix m=#[1];
751 intvec v=#[2];
752 int switch=#[3];
753 if (switch==-1) {
754   "An error has occurred in develop, so there is no HNE.";
755   return(ideal(0,0));
756 }
757 int fehler,fehlervor,untergrad,untervor,beginn,i,zeile,hilf;
758
759 if (nvars(basering) > 2) { poly z(size(v)+1)=var(3); }
760 else                     { poly z(size(v)+1)=var(1); }
761 poly z(size(v));
762 zeile=size(v);
763 //------------- Parametrisierung der untersten Zeile der HNE -----------------
764 if (v[zeile] > 0) {
765   fehler=0;           // die Parametrisierung wird exakt werden
766   for (i=1; i<=v[zeile]; i++) {
767     z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
768 }}
769 else {
770   untervor=1;         // = Untergrad der vorhergehenden Zeile
771   if (v[zeile]==-1) {
772     fehler=ncols(m)+1;
773     for (i=1; i<=ncols(m); i++) {
774       z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
775       if ((untergrad==0) and (m[zeile,i]!=0)) {untergrad=i;}
776                       // = Untergrad der aktuellen Zeile
777   }}
778   else {
779     fehler= -v[zeile];
780     for (i=1; i<-v[zeile]; i++) {
781       z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
782       if ((untergrad==0) and (m[zeile,i]!=0)) {untergrad=i;}
783   }}
784 }
785 //------------- Parametrisierung der restlichen Zeilen der HNE ---------------
786 for (zeile=size(v)-1; zeile>0; zeile--) {
787   poly z(zeile);
788   beginn=0;             // Beginn der aktuellen Zeile
789   for (i=1; i<=v[zeile]; i++) {
790     z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i;
791     if ((beginn==0) and (m[zeile,i]!=0)) { beginn=i;}
792   }
793   z(zeile)=z(zeile) + z(zeile+1)^v[zeile] * z(zeile+2);
794   if (beginn==0) {
795     if (fehler>0) {     // damit fehler=0 bleibt bei exakter Param.
796     fehlervor=fehler;   // Fehler der letzten Zeile
797     fehler=fehler+untergrad*(v[zeile]-1)+untervor;   // Fehler dieser Zeile
798     hilf=untergrad;
799     untergrad=untergrad*v[zeile]+untervor;
800     untervor=hilf;}}    // untervor = altes untergrad
801   else {
802     fehlervor=fehler;
803     fehler=fehler+untergrad*(beginn-1);
804     untervor=untergrad;
805     untergrad=untergrad*beginn;
806 }}
807 //--------------------- Ausgabe der Fehlerabschaetzung -----------------------
808 if (switch==0) {
809   if (fehler>0) {
810     if (fehlervor>0)
811      { "// ** Warning: result is exact up to order",fehlervor-1,"in",
812        maxideal(1)[1],"and",fehler-1,"in",maxideal(1)[2],"!";}
813     else { "// ** Warning: result is exact up to order",fehler-1,"in",
814        maxideal(1)[2],"!";}}
815   return(ideal(z(2),z(1)));}
816 else {
817   if (fehler>0) {
818     if (fehlervor>0)
819      { "// ** Warning: result is exact up to order",fehler-1,"in",
820        maxideal(1)[1],"and",fehlervor-1,"in",maxideal(1)[2],"!";}
821     else { "// ** Warning: result is exact up to order",fehler-1,"in",
822        maxideal(1)[1],"!";}}
823   return(ideal(z(1),z(2)));}
824}
825example
826{ "EXAMPLE:"; echo = 2;
827 ring exring=0,(x,y,t),ds;
828 poly f=x3+2xy2+y2;
829 list hne=develop(f);
830 list hne_extended1=develop(f,6);
831 list hne_extended2=develop(f,10);
832 pause;
833                                   //   compare the different matrices ...
834 print(hne[1]);
835 print(hne_extended1[1]);
836 print(hne_extended2[1]);
837                                   // ... and the resulting parametrizations:
838 param(hne);
839 param(hne_extended1);
840 param(hne_extended2);
841}
842///////////////////////////////////////////////////////////////////////////////
843
844proc invariants
845"USAGE:  invariants(l) takes the output of develop(f)
846        (list l (matrix m, intvec v, int s[,poly g]))
847        and computes some invariants:
848RETURN:  a list, if l contains a valid HNE:
849    - invariants(l)[1]: intvec         (characteristic exponents)
850    - invariants(1)[2],[3]: 2 x intvec (puiseux pairs, 1st and 2nd components)
851    - invariants(l)[4]: int            (degree of the conductor)
852    - invariants(l)[5]: intvec         (sequence of multiplicities)
853         an empty list, if an error occurred in the procedure develop
854EXAMPLE: example invariants; shows an example
855"{
856 //-------------------------- Initialisierungen -------------------------------
857 matrix m=#[1];
858 intvec v=#[2];
859 int switch=#[3];
860 list ergebnis;
861 if (switch==-1) {
862   "An error has occurred in develop, so there is no HNE.";
863   return(ergebnis);
864 }
865 intvec beta,s,svorl,ordnung,multips,multseq,mpuiseux,npuiseux;
866 int genus,zeile,i,j,k,summe,conductor,ggT;
867 string Ausgabe;
868 int nc=ncols(m); int nr=nrows(m);
869 ordnung[nr]=1;
870         // alle Indizes muessen (gegenueber [Ca]) um 1 erhoeht werden,
871         // weil 0..r nicht als Wertebereich erlaubt ist (aber nrows(m)==r+1)
872
873 //---------------- Bestimme den Untergrad der einzelnen Zeilen ---------------
874 for (zeile=nr; zeile>1; zeile--) {
875   if ((size(ideal(m[zeile,1..nc])) > 1) or (zeile==nr)) { // keine Nullzeile
876      k=1;
877      while (m[zeile,k]==0) {k++;}
878      ordnung[zeile-1]=k*ordnung[zeile]; // vgl. auch Def. von untergrad in
879      genus++;                           // proc param
880      svorl[genus]=zeile;} // werden gerade in umgekehrter Reihenfolge abgelegt
881   else {
882      ordnung[zeile-1]=v[zeile]*ordnung[zeile]+ordnung[zeile+1];
883 }}
884 //----------------- charakteristische Exponenten (beta) ----------------------
885 s[1]=1;
886 for (k=1; k <= genus; k++) { s[k+1]=svorl[genus-k+1];} // s[2]==s(1), u.s.w.
887 beta[1]=ordnung[1]; //charakt. Exponenten: Index wieder verschoben
888 for (k=1; k <= genus; k++) {
889   summe=0;
890   for (i=1; i <= s[k]; i++) {summe=summe+v[i]*ordnung[i];}
891   beta[k+1]=summe+ordnung[s[k]]+ordnung[s[k]+1]-ordnung[1];
892 }
893 //--------------------------- Puiseuxpaare -----------------------------------
894 int produkt=1;
895 for (i=1; i<=genus; i++) {
896   ggT=gcd(beta[1],beta[i+1]*produkt);
897   mpuiseux[i]=beta[i+1]*produkt / ggT;
898   npuiseux[i]=beta[1] / ggT;
899   produkt=produkt*npuiseux[i];
900 }
901 //---------------------- Grad des Konduktors ---------------------------------
902 summe=1-ordnung[1];
903 if (genus > 0) {
904   for (i=2; i <= genus+1; i++) {
905     summe=summe + beta[i] * (ordnung[s[i-1]] - ordnung[s[i]]);
906   }                              // n.b.: Indizierung wieder um 1 verschoben
907 }
908 conductor=summe;
909 //------------------- Multiplizitaetensequenz: -------------------------------
910 k=1;
911 multips=multiplicities(#);
912 for (i=1; i<size(v); i++) {
913   for (j=1; j<=v[i]; j++) {
914     multseq[k]=multips[i];
915     k++;
916 }}
917 multseq[k]=1;
918 //------------------------- Rueckgabe ----------------------------------------
919 ergebnis=beta,mpuiseux,npuiseux,conductor,multseq;
920 return(ergebnis);
921}
922example
923{ "EXAMPLE:"; echo = 2;
924 ring exring=0,(x,y),dp;
925 list hne=develop(y4+2x3y2+x6+x5y);
926 list erg=invariants(hne);
927 erg[1];                   // the characteristic exponents
928 erg[2],erg[3];            // the puiseux pairs in packed form
929 erg[4] / 2;               // the delta-invariant
930 erg[5];                   // the sequence of multiplicities
931                           // To display the invariants more 'nicely':
932 displayInvariants(hne);
933}
934///////////////////////////////////////////////////////////////////////////////
935
936proc displayInvariants
937"USAGE:   displayInvariants(l); takes the output both of
938         develop(f) and reddevelop(f)
939          ( list l (matrix m, intvec v, int s[,poly g])
940            or list of lists in the form l            )
941DISPLAY: invariants of the corresponding branch, resp. of all branches,
942         in a better readable form
943RETURN:  nothing
944EXAMPLE: example displayInvariants;  shows an example
945"{
946 int i,j,k,mul;
947 string Ausgabe;
948 list ergebnis;
949 //-- entferne ueberfluessige Daten zur Erhoehung der Rechengeschwindigkeit: --
950 #=stripHNE(#);
951 //-------------------- Ausgabe eines Zweiges ---------------------------------
952 if (typeof(#[1])=="matrix") {
953   ergebnis=invariants(#);
954   if (size(ergebnis)!=0) {
955    " characteristic exponents  :",ergebnis[1];
956    if (size(ergebnis[1])>1) {
957     for (i=1; i<=size(ergebnis[2]); i++) {
958       Ausgabe=Ausgabe+"("+string(ergebnis[2][i])+","
959       +string(ergebnis[3][i])+")";
960    }}
961    " puiseux pairs             :",Ausgabe;
962    " degree of the conductor   :",ergebnis[4];
963    " delta invariant           :",ergebnis[4]/2;
964    " generators of semigroup   :",puiseux2generators(ergebnis[2],ergebnis[3]);
965    " sequence of multiplicities:",ergebnis[5];
966 }}
967 //-------------------- Ausgabe aller Zweige ----------------------------------
968 else {
969  for (j=1; j<=size(#); j++) {
970    ergebnis=invariants(#[j]);
971    " --- invariants of branch number",j,": ---";
972    " characteristic exponents  :",ergebnis[1];
973    Ausgabe="";
974    if (size(ergebnis[1])>1) {
975     for (i=1; i<=size(ergebnis[2]); i++) {
976       Ausgabe=Ausgabe+"("+string(ergebnis[2][i])+","
977       +string(ergebnis[3][i])+")";
978    }}
979    " puiseux pairs             :",Ausgabe;
980    " degree of the conductor   :",ergebnis[4];
981    " delta invariant           :",ergebnis[4]/2;
982    " generators of semigroup   :",puiseux2generators(ergebnis[2],ergebnis[3]);
983    " sequence of multiplicities:",ergebnis[5];
984    "";
985  }
986  if (size(#)>1) {
987    " -------------- intersection multiplicities : -------------- ";"";
988    Ausgabe="branch |   ";
989    for (j=size(#); j>1; j--) {
990      if (size(string(j))==1) { Ausgabe=Ausgabe+" "+string(j)+"    "; }
991      else                    { Ausgabe=Ausgabe+string(j)+"    "; }
992    }
993    Ausgabe;
994    Ausgabe="-------+";
995    for (j=2; j<size(#); j++) { Ausgabe=Ausgabe+"------"; }
996    Ausgabe=Ausgabe+"-----";
997    Ausgabe;
998  }
999  for (j=1; j<size(#); j++) {
1000    if (size(string(j))==1) { Ausgabe="    "+string(j)+"  |"; }
1001    else                    { Ausgabe="   " +string(j)+"  |"; }
1002    for (k=size(#); k>j; k--) {
1003      mul=intersection(#[j],#[k]);
1004      for (i=1; i<=5-size(string(mul)); i++) { Ausgabe=Ausgabe+" "; }
1005      Ausgabe=Ausgabe+string(mul);
1006      if (k>j+1) { Ausgabe=Ausgabe+","; }
1007    }
1008    Ausgabe;
1009  }
1010 }
1011 return();
1012}
1013example
1014{ "EXAMPLE:"; echo = 2;
1015 ring exring=0,(x,y),dp;
1016 list hne=develop(y4+2x3y2+x6+x5y);
1017 displayInvariants(hne);
1018}
1019///////////////////////////////////////////////////////////////////////////////
1020
1021proc multiplicities
1022"USAGE:   multiplicities(l) takes the output of develop(f)
1023         (list l (matrix m, intvec v, int s[,poly g]))
1024RETURN:  intvec of the different multiplicities that occur during the succes-
1025         sive blowing up of the curve corresponding to f
1026EXAMPLE: example multiplicities;  shows an example
1027"{
1028 matrix m=#[1];
1029 intvec v=#[2];
1030 int switch=#[3];
1031 list ergebnis;
1032 if (switch==-1) {
1033   "An error has occurred in develop, so there is no HNE.";
1034   return(intvec(0));
1035 }
1036 intvec ordnung;
1037 int zeile,k;
1038 int nc=ncols(m); int nr=nrows(m);
1039 ordnung[nr]=1;
1040 //---------------- Bestimme den Untergrad der einzelnen Zeilen ---------------
1041 for (zeile=nr; zeile>1; zeile--) {
1042   if ((size(ideal(m[zeile,1..nc])) > 1) or (zeile==nr)) { // keine Nullzeile
1043      k=1;
1044      while (m[zeile,k]==0) {k++;}
1045      ordnung[zeile-1]=k*ordnung[zeile];
1046   }
1047   else {
1048      ordnung[zeile-1]=v[zeile]*ordnung[zeile]+ordnung[zeile+1];
1049 }}
1050 return(ordnung);
1051}
1052example
1053{ "EXAMPLE:"; echo = 2;
1054  ring r=0,(x,y),dp;
1055  multiplicities(develop(x5+y7));
1056  // The first value is the multiplicity of the curve itself, here it's 5
1057}
1058///////////////////////////////////////////////////////////////////////////////
1059
1060proc puiseux2generators (intvec m, intvec n)
1061"USAGE:   puiseux2generators (m,n);
1062         m,n intvecs with 1st resp. 2nd components of puiseux pairs
1063RETURN:  intvec of the generators of the semigroup of values
1064EXAMPLE: example puiseux2generators;  shows an example
1065"{
1066 intvec beta;
1067 int q=1;
1068 //------------ glatte Kurve (eigentl. waeren m,n leer): ----------------------
1069 if (m==0) { return(intvec(1)); }
1070 //------------------- singulaere Kurve: --------------------------------------
1071 for (int i=1; i<=size(n); i++) { q=q*n[i]; }
1072 beta[1]=q; // == q_0
1073 m=1,m; n=1,n; // m[1] ist damit m_0 usw., genau wie beta[1]==beta_0
1074 for (i=2; i<=size(n); i++) {
1075  beta[i]=m[i]*q/n[i] - m[i-1]*q + n[i-1]*beta[i-1];
1076  q=q/n[i]; // == q_i
1077 }
1078 return(beta);
1079}
1080example
1081{ "EXAMPLE:"; echo = 2;
1082  // take (3,2),(7,2),(15,2),(31,2),(63,2),(127,2) as puiseux pairs:
1083  puiseux2generators(intvec(3,7,15,31,63,127),intvec(2,2,2,2,2,2));
1084}
1085///////////////////////////////////////////////////////////////////////////////
1086
1087proc generators
1088"USAGE:   generators(l); takes the output of develop(f)
1089         (list l (matrix m, intvec v, int s[,poly g]))
1090RETURN:  intvec of the generators of the semigroup of values
1091EXAMPLE: example generators;  shows an example
1092"
1093{
1094 list inv=invariants(#);
1095 return(puiseux2generators(inv[2..3]));
1096}
1097example
1098{ "EXAMPLE:"; echo = 2;
1099  ring r=0,(x,y),dp;
1100  poly f=x15-21x14+8x13y-6x13-16x12y+20x11y2-1x12+8x11y-36x10y2+24x9y3+4x9y2
1101         -16x8y3+26x7y4-6x6y4+8x5y5+4x3y6-1y8;
1102  list hne=develop(f);
1103  generators(hne);
1104}
1105///////////////////////////////////////////////////////////////////////////////
1106
1107proc intersection (list hn1, list hn2)
1108"USAGE:   intersection(hne1,hne2);
1109         hne1,hne2 two lists representing a HNE (normally two entries out of
1110         the output of reddevelop), i.e. list(matrix,intvec,int[,poly])
1111RETURN:  intersection multiplicity of the branches corresponding to hne1 & hne2
1112         (of type int)
1113EXAMPLE: example intersection;  shows an example
1114"
1115{
1116 //------------------ `intersect' ist schon reserviert ... --------------------
1117 int i,j,s,sum,schnitt,unterschied;
1118 matrix a1=hn1[1];
1119 matrix a2=hn2[1];
1120 intvec h1=hn1[2];
1121 intvec h2=hn2[2];
1122 intvec n1=multiplicities(hn1);
1123 intvec n2=multiplicities(hn2);
1124 if (hn1[3]!=hn2[3]) {
1125 //-- die jeweils erste Zeile von hn1,hn2 gehoert zu verschiedenen Parametern -
1126 //---------------- d.h. beide Kurven schneiden sich transversal --------------
1127   schnitt=n1[1]*n2[1];        // = mult(hn1)*mult(hn2)
1128 }
1129 else {
1130 //--------- die jeweils erste Zeile gehoert zum gleichen Parameter -----------
1131   unterschied=0;
1132 //h1[1],h2[1];size(h1),size(h2);
1133   for (s=1; (h1[s]==h2[s]) && (s<size(h1)) && (s<size(h2))
1134              && (unterschied==0); s++) {
1135     for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]); i++) {;}
1136     if (i<=h1[s]) {
1137       unterschied=1;
1138       s--;                    // um s++ am Schleifenende wieder auszugleichen
1139     }
1140   }
1141   if (unterschied==0) {
1142     if ((s<size(h1)) && (s<size(h2))) {
1143       for (i=1; (a1[s,i]==a2[s,i]) && (i<=h1[s]) && (i<=h2[s]); i++) {;}
1144     }
1145     else {
1146 //-------------- Sonderfall: Unterschied in letzter Zeile suchen -------------
1147 // Beachte: Es koennen undefinierte Stellen auftreten, bei abbrechender HNE
1148 // muss die Ende-Markierung weg, h_[r] ist unendlich, die Matrix muss mit
1149 // Nullen fortgesetzt gedacht werden
1150 //----------------------------------------------------------------------------
1151       if (ncols(a1)>ncols(a2)) { j=ncols(a1); }
1152       else                     { j=ncols(a2); }
1153       unterschied=0;
1154       if ((h1[s]>0) && (s==size(h1))) {
1155         a1[s,h1[s]+1]=0;
1156         if (ncols(a1)<=ncols(a2)) { unterschied=1; }
1157       }
1158       if ((h2[s]>0) && (s==size(h2))) {
1159         a2[s,h2[s]+1]=0;
1160         if (ncols(a2)<=ncols(a1)) { unterschied=1; }
1161       }
1162       if (unterschied==1) {                   // mind. eine HNE war endlich
1163         matrix ma1[1][j]=a1[s,1..ncols(a1)];  // und bedarf der Fortsetzung
1164         matrix ma2[1][j]=a2[s,1..ncols(a2)];  // mit Nullen
1165       }
1166       else {
1167         if (ncols(a1)>ncols(a2)) { j=ncols(a2); }
1168         else                     { j=ncols(a1); }
1169         matrix ma1[1][j]=a1[s,1..j];          // Beschr. auf vergleichbaren
1170         matrix ma2[1][j]=a2[s,1..j];          // Teil (der evtl. y's enth.)
1171       }
1172 // print(ma1);print(ma2);
1173       for (i=1; (ma1[1,i]==ma2[1,i]) && (i<j) && (ma1[1,i]!=var(2)); i++) {;}
1174       if (ma1[1,i]==ma2[1,i]) {
1175         "//** The two HNE's are identical!";
1176         "//** You have either tried to intersect a branch with itself,";
1177         "//** or the two branches have been developed separately.";
1178         "//   In the latter case use `extdevelop' to extend the HNE's until",
1179         "they differ.";
1180         return(-1);
1181       }
1182       if ((ma1[1,i]==var(2)) || (ma2[1,i]==var(2))) {
1183         "//** The two HNE's are (so far) identical. This is because they",
1184         "have been";
1185         "//** computed separately. I need more data; use `extdevelop' to",
1186         "extend them,";
1187         if (ma1[1,i]==var(2)) {"//** at least the first one.";}
1188         else                  {"//** at least the second one.";}
1189         return(-1);
1190       }
1191     }
1192   }
1193   sum=0;
1194   h1[size(h1)]=ncols(a1)+42;        // Ersatz fuer h1[r]=infinity
1195   h2[size(h2)]=ncols(a2)+42;
1196   for (j=1; j<s; j++) {sum=sum+h1[j]*n1[j]*n2[j];}
1197   if ((i<=h1[s]) && (i<=h2[s]))    { schnitt=sum+i*n1[s]*n2[s]; }
1198   if (i==h2[s]+1) { schnitt=sum+h2[s]*n1[s]*n2[s]+n2[s+1]*n1[s]; }
1199   if (i==h1[s]+1) { schnitt=sum+h1[s]*n2[s]*n1[s]+n1[s+1]*n2[s]; }
1200 // "s:",s-1,"i:",i,"S:",sum;
1201 }
1202 return(schnitt);
1203}
1204example
1205{
1206  if (nameof(basering)=="HNEring") {
1207   def rettering=HNEring;
1208   kill HNEring;
1209  }
1210  "EXAMPLE:"; echo = 2;
1211  ring r=0,(x,y),dp;
1212  list hne=reddevelop((x2-y3)*(x2+y3));
1213  intersection(hne[1],hne[2]);
1214  kill HNEring,r;
1215  echo = 0;
1216  if (defined(rettering)) { // wenn aktueller Ring vor Aufruf von example
1217   setring rettering;       // HNEring war, muss dieser erst wieder restauriert
1218   def HNEring=rettering;   // werden
1219   export HNEring;
1220  }
1221}
1222///////////////////////////////////////////////////////////////////////////////
1223
1224proc displayHNE(list ldev,list #)
1225"USAGE:   displayHNE(ldev,[,n]); ldev=list
1226         (output list of develop(f) or reddevelop(f)), n=int
1227RETURN:  - if only one argument is given, no return value, but
1228         display an ideal HNE of the following form:
1229            HNE[1]=-y+[]*z(0)^1+[]*z(0)^2+...+z(0)^<>*z(1)
1230            HNE[2]=-x+          []*z(1)^2+...+z(1)^<>*z(2)
1231            HNE[3]=             []*z(2)^2+...+z(2)^<>*z(3)
1232            .......             ..........................
1233            HNE[r+1]=           []*z(r)^2+[]*z(r)^3+......
1234     where x,y are the indeterminates of the basering. The values of [] are
1235     the coefficients of the Hamburger-Noether matrix, the values of <> are
1236     represented in the HN-matrix as 'x'
1237         the 1st line (HNE[1]) means that y==[]*z(0)^1+... ,
1238         the 2nd line (HNE[2]) means that x==[]*z(1)^1+...
1239     so you can see which indeterminate corresponds to which line
1240     (it's also possible that x corresponds to the 1st line and y to the 2nd)
1241         - if a second argument is given, create and export a new ring with
1242         name 'displayring' containing an ideal 'HNE' as described above
1243
1244     If ldev contains the output of reddevelop(f), displayHNE shows the HNE's
1245     of all branches of f in the form described above.
1246     The optional parameter is then ignored.
1247EXAMPLE: example displayHNE; shows an example
1248"
1249{
1250 if ((typeof(ldev[1])=="list") || (typeof(ldev[1])=="none")) {
1251   for (int i=1; i<=size(ldev); i++) {
1252     "// Hamburger-Noether development of branch nr."+string(i)+":";
1253     displayHNE(ldev[i]);"";
1254   }
1255   return();
1256 }
1257 //--------------------- Initialisierungen und Ringwechsel --------------------
1258 matrix m=ldev[1];
1259 intvec v=ldev[2];
1260 int switch=ldev[3];
1261 if (switch==-1) {
1262   "An error has occurred in develop, so there is no HNE.";
1263   return(ideal(0));
1264 }
1265 if (parstr(basering)!="") {
1266   if (charstr(basering)!=string(char(basering))+","+parstr(basering)) {
1267     "Sorry -- not implemented for rings of type (p^k,a),...";
1268     return(ideal(0));
1269   }
1270 }
1271 def altring=basering;
1272 ring dazu=char(altring),z(0..size(v)-1),ls;
1273 def displayring=dazu+altring;
1274 setring displayring;
1275 if (size(#) != 0) {
1276    export displayring;
1277  }
1278 map holematrix=altring,0;        // mappt nur die Monome vom Grad Null
1279 matrix m=holematrix(m);
1280 //--------------------- Erzeuge Matrix n mit n[i,j]=z(j-1)^i -----------------
1281 int i;
1282 matrix n[ncols(m)][nrows(m)];
1283 for (int j=1; j<=nrows(m); j++) {
1284    for (i=1; i<=ncols(m); i++) { n[i,j]=z(j-1)^i; }
1285 }
1286 matrix displaymat=m*n;
1287 ideal HNE;
1288 for (i=1; i<nrows(m); i++) { HNE[i]=displaymat[i,i]+z(i)*z(i-1)^v[i]; }
1289 HNE[nrows(m)]=displaymat[nrows(m),nrows(m)];
1290 if (nrows(m)<2) { HNE[2]=z(0); }
1291 if (switch==0) {
1292    HNE[1] = HNE[1]-var(size(v)+2);
1293    HNE[2] = HNE[2]-var(size(v)+1);
1294 }
1295 else {
1296    HNE[1] = HNE[1]-var(size(v)+1);
1297    HNE[2] = HNE[2]-var(size(v)+2);
1298 }
1299if (size(#) == 0) {
1300   HNE;
1301   return();
1302 }
1303if (size(#) != 0) {
1304   "// basering is now 'displayring' containing ideal 'HNE'";
1305   keepring(displayring);
1306   export(HNE);
1307   return(HNE);
1308 }
1309}
1310example
1311{ "EXAMPLE:"; echo = 2;
1312  ring r=0,(x,y),dp;
1313  poly f=x3+2xy2+y2;
1314  list hn=develop(f);
1315  displayHNE(hn);
1316}
1317///////////////////////////////////////////////////////////////////////////////
1318//                      procedures for reducible curves                      //
1319///////////////////////////////////////////////////////////////////////////////
1320
1321
1322proc newtonhoehne (poly f)
1323"USAGE:   newtonhoehne(f);   f poly
1324ASSUME:  basering = ...,(x,y),ds  or ls
1325RETURN:  list of intvec(x,y) of coordinates of the newtonpolygon of f
1326NOTE:    This proc is only available in versions of Singular that know the
1327         command  system(\"newton\",f);  f poly
1328"
1329{
1330 intvec nm = getnm(f);
1331 if ((nm[1]>0) && (nm[2]>0)) { f=jet(f,nm[1]*nm[2],nm); }
1332 list erg=system("newton",f);
1333 int i; list Ausgabe;
1334 for (i=1; i<=size(erg); i++) { Ausgabe[i]=leadexp(erg[i]); }
1335 return(Ausgabe);
1336}
1337///////////////////////////////////////////////////////////////////////////////
1338
1339proc newtonpoly (poly f)
1340"USAGE:   newtonpoly(f);   f poly
1341RETURN:  list of intvec(x,y) of coordinates of the newtonpolygon of f
1342NOTE: There are many assumptions: (to improve speed)
1343      - The basering must be ...,(x,y),ls
1344      - f(x,0) != 0 != f(0,y), f(0,0) = 0
1345EXAMPLE: example newtonpoly;  shows an example
1346"
1347{
1348 intvec A=(0,ord(subst(f,x,0)));
1349 intvec B=(ord(subst(f,y,0)),0);
1350 intvec C,D,H; list L;
1351 int abbruch,i;
1352 poly hilf;
1353
1354 L[1]=A;
1355 //-------- wirf alle Monome auf oder oberhalb der Geraden AB raus: -----------
1356 f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]));
1357 map xytausch=basering,y,x;
1358 for (i=2; f!=0; i++) {
1359   abbruch=0;
1360   while (abbruch==0) {
1361 // finde den Punkt aus {verbliebene Pkte (a,b) mit a minimal} mit b minimal: -
1362
1363     C=leadexp(f);           // Ordnung ls ist wesentlich!
1364
1365     if (jet(f,A[2]*C[1]-A[1]*C[2]-1,intvec(A[2]-C[2],C[1]-A[1]))==0)
1366       { abbruch=1; }        // keine Monome unterhalb der Geraden AC
1367
1368 // ----- alle Monome auf der Parallelen zur y-Achse durch C wegwerfen: -------
1369 // ------------------ (links von C gibt es sowieso keine mehr) ---------------
1370     else { f=jet(f,-C[1]-1,intvec(-1,0)); }
1371   }
1372 //- finde alle Monome auf der Geraden durch A und C (unterhalb gibt's keine) -
1373   hilf=jet(f,A[2]*C[1]-A[1]*C[2],intvec(A[2]-C[2],C[1]-A[1]));
1374
1375   H=leadexp(xytausch(hilf));
1376   D=H[2],H[1];
1377
1378 // die Alternative waere ein Ringwechsel nach ..,(y,x),ds gewesen
1379 // D ist der naechste Eckpunkt (unterster Punkt auf Geraden durch A,C)
1380
1381   L[i]=D;
1382   A=D;
1383 //----------------- alle Monome auf oder unterhalb DB raus -------------------
1384   f=jet(f,D[2]*B[1]-1,intvec(D[2],B[1]-D[1]));
1385 }
1386 L[i]=B;
1387 return(L);
1388}
1389example
1390{ "EXAMPLE:";
1391  ring @exring_Newt=0,(x,y),ls;
1392  export @exring_Newt;
1393  "  ring exring=0,(x,y),ls;";
1394  echo = 2;
1395  poly f=x5+2x3y-x2y2+3xy5+y6-y7;
1396  newtonpoly(f);
1397  echo = 0;
1398  kill @exring_Newt;
1399}
1400///////////////////////////////////////////////////////////////////////////////
1401
1402proc charPoly(poly f, int M, int N)
1403"USAGE:  charPoly(f,M,N);  f poly in x,y,  M,N int: length and height
1404                          of newtonpolygon of f, which has to be only one line
1405RETURN:  the characteristic polynomial of f
1406EXAMPLE: example charPoly;  shows an example
1407"
1408{
1409 poly charp;
1410 int Np=N/ gcd(M,N);
1411 f=subst(f,x,1);
1412 for(charp=0; f<>0; f=f-lead(f))
1413  { charp=charp+leadcoef(f)*y^(leadexp(f)[2]/ Np);}
1414 return(charp);
1415}
1416example
1417{ "EXAMPLE:"; echo = 2;
1418  ring exring=0,(x,y),dp;
1419  charPoly(y4+2y3x2-yx6+x8,8,4);
1420  charPoly(y6+3y3x2-x4,4,6);
1421}
1422///////////////////////////////////////////////////////////////////////////////
1423
1424proc find_in_list(list L,int p)
1425"USAGE:  find_in_list(L,p); L: list of intvec(x,y)
1426        (sorted in y: L[1][2]>=L[2][2]), int p >= 0
1427RETURN: int i: L[i][2]=p if existent; otherwise i with L[i][2]<p if existent;
1428        otherwise i = size(L)+1;
1429"
1430{
1431 int i;
1432 L[size(L)+1]=intvec(0,-1);          // falls p nicht in L[.][2] vorkommt
1433 for (i=1; L[i][2]>p; i++) {;}
1434 return(i);
1435}
1436///////////////////////////////////////////////////////////////////////////////
1437proc set_list
1438"USAGE: set_list(L,v,l); L list, v intvec, l element of L (arbitrary type)
1439RETURN: L with L[v[1]][v[2]][...][v[n]] set to l
1440WARNING: Make sure there is no conflict of types!
1441         size(v) must be bigger than 1! (i.e. v must NOT be simply an integer)
1442EXAMPLE:                 L=set_list(L,intvec(a,b,c),q);
1443         does exactly what you would expect when typing
1444                         L[a][b][c]=q;
1445         (but the last command will only produce an error message)
1446"
1447{
1448 list L=#[1];
1449 intvec v=#[2];
1450 def ersetze=#[3];
1451 //---------------- Bilde den Verzweigungsbaum nach ---------------------------
1452 def hilf(1)=L[v[1]];
1453 for (int tiefe=2; tiefe<size(v); tiefe++) {
1454   def hilf(tiefe)=hilf(tiefe-1)[v[tiefe]];
1455 }
1456 //---------------- Fuege Aenderung in die Liste ein --------------------------
1457 hilf(size(v)-1)[v[size(v)]]=ersetze;
1458 for (tiefe=size(v)-1; tiefe>1; tiefe--) {
1459   hilf(tiefe-1)[v[tiefe]]=hilf(tiefe);
1460 }
1461 L[v[1]]=hilf(1);
1462 return(L);
1463}
1464///////////////////////////////////////////////////////////////////////////////
1465proc get_last_divisor(int M, int N)
1466"USAGE:   get_last_divisor(M,N); int M,N
1467RETURN:  int Q: M=q1*N+r1, N=q2*r1+r2, ..., ri=Q*r(i+1) (Euclidean alg.)
1468EXAMPLE: example get_last_divisor; shows an example
1469"
1470{
1471 int R=M%N; int Q=M / N;
1472 while (R!=0) {M=N; N=R; R=M%N; Q=M / N;}
1473 return(Q)
1474}
1475example
1476{ "EXAMPLE"; echo = 2;
1477  ring r=0,(x,y),dp;
1478  get_last_divisor(12,10);
1479}
1480///////////////////////////////////////////////////////////////////////////////
1481proc redleit (poly f,intvec S, intvec E)
1482"USAGE:   redleit(f,S,E); f poly,
1483         S,E intvec(x,y): two different points on a line in the newtonpolygon
1484         of f (e.g. the starting and the ending point)
1485RETURN:  poly g: all monomials of f which lay on that line
1486NOTE:    It is not checked whether S and E really belong to the newtonpolygon!
1487EXAMPLE: example redleit;  shows an example
1488"
1489{
1490 if (E[1]<S[1]) { intvec H=E; E=S; S=H; } // S,E verkehrt herum eingegeben
1491 return(jet(f,E[1]*S[2]-E[2]*S[1],intvec(S[2]-E[2],E[1]-S[1])));
1492}
1493example
1494{ "EXAMPLE"; echo = 2;
1495  ring exring=0,(x,y),dp;
1496  redleit(y6+xy4-2x3y2+x4y+x6,intvec(3,2),intvec(4,1));
1497}
1498///////////////////////////////////////////////////////////////////////////////
1499
1500
1501proc extdevelop (list l, int Exaktheit)
1502"USAGE:   extdevelop(l,n);  list l (matrix m, intvec v, int s, poly g), int n
1503         takes the output l of develop(f) ( or extdevelop(L,n),
1504         or one entry of the output of reddevelop(f))  and
1505RETURN:  an extension of the Hamburger-Noether development of f in the form
1506         of develop, i.e. a list L (matrix m,intvec v,...)
1507         The new HN-matrix will have at least n columns
1508         (if the HNE isn't finite).
1509         Thus if f is irreducible, extdevelop(develop(f),n);
1510         (in most cases) will produce the same result as develop(f,n).
1511         Type  help develop;  for more details
1512EXAMPLE: example extdevelop;  shows an example
1513"
1514{
1515 //------------ Initialisierungen und Abfangen unzulaessiger Aufrufe ----------
1516 matrix m=l[1];
1517 intvec v=l[2];
1518 int switch=l[3];
1519 if (nvars(basering) < 2) {
1520   " Sorry. I need two variables in the ring.";
1521   return(list(matrix(maxideal(1)[1]),intvec(0),-1,poly(0)));}
1522 if (switch==-1) {
1523   "An error has occurred in develop, so there is no HNE and no extension.";
1524   return(l);
1525 }
1526 poly f=l[4];
1527 if (f==0) {
1528   " No extension is possible";
1529   return(l);
1530 }
1531 int Q=v[size(v)];
1532 if (Q>0) {
1533   " The HNE was already exact";
1534   return(l);
1535 }
1536 else {
1537   if (Q==-1) { Q=ncols(m); }
1538   else { Q=-Q-1; }
1539 }
1540 int zeile=nrows(m);
1541 int spalten,i,lastside;
1542 ideal lastrow=m[zeile,1..Q];
1543 int ringwechsel=(varstr(basering)!="x,y") or (ordstr(basering)!="ls(2),C");
1544
1545 //------------------------- Ringwechsel, falls noetig ------------------------
1546 if (ringwechsel) {
1547  def altring = basering;
1548  int p = char(basering); // Ringcharakteristik
1549  if (charstr(basering)!=string(p)) {
1550     string tststr=charstr(basering);
1551     tststr=tststr[1..find(tststr,",")-1];     //-> "p^k" bzw. "p"
1552     if (tststr==string(p)) {
1553       if (size(parstr(basering))>1) {         // ring (p,a,..),...
1554        execute "ring extdguenstig=("+charstr(basering)+"),(x,y),ls;";
1555       }
1556       else {                                  // ring (p,a),...
1557        string mipl=string(minpoly);
1558        ring extdguenstig=(p,`parstr(basering)`),(x,y),ls;
1559        if (mipl!="0") { execute "minpoly="+mipl+";"; }
1560       }
1561     }
1562     else {
1563       execute "ring extdguenstig=("+charstr(basering)+"),(x,y),ls;";
1564     }
1565  }
1566  else {                               // charstr(basering)== p : no parameter
1567     ring extdguenstig=p,(x,y),ls;
1568  }
1569  export extdguenstig;
1570  map hole=altring,x,y;
1571  poly f=hole(f);
1572  ideal a=hole(lastrow);
1573  //----- map hat sich als sehr zeitaufwendig bei f gross herausgestellt ------
1574  //-------- (fetch,imap nur 5% schneller), daher nach Moeglichkeit -----------
1575  //-----------------Beibehaltung des alten Ringes: ---------------------------
1576 }
1577 else { ideal a=lastrow; }
1578 list Newton;
1579 number delta;
1580 //--------------------- Fortsetzung der HNE ----------------------------------
1581 while (Q<Exaktheit) {
1582  Newton=newtonpoly(f);
1583  lastside=size(Newton)-1;
1584  if (Newton[lastside][2]!=1) {
1585    " *** The transformed polynomial was not valid!!";}
1586  Q=Newton[lastside+1][1]-Newton[lastside][1];      // Laenge der letzten Seite
1587
1588 //--- quasihomogenes Leitpolynom ist c*x^a*y+d*x^(a+Q) => delta=-d/c: --------
1589  delta=-koeff(f,Newton[lastside+1][1],0)/koeff(f,Newton[lastside][1],1);
1590  a[Q]=delta;
1591  "a("+string(zeile-1)+","+string(Q)+") =",delta;
1592  if (Q<Exaktheit) {
1593   f=T1_Transform(f,delta,Q);
1594   if (defined(Protokoll)) { "transformed polynomial:",f; }
1595   if (subst(f,y,0)==0) {
1596     "The HNE is finite!";
1597     a[Q+1]=x; Exaktheit=Q;
1598     f=0;                          // Speicherersparnis: f nicht mehr gebraucht
1599   }
1600  }
1601 }
1602 //------- Wechsel in alten Ring, Zusammensetzung alte HNE + Erweiterung ------
1603 if (ringwechsel) {
1604  setring altring;
1605  map zurueck=extdguenstig,var(1),var(2);
1606  f=zurueck(f);
1607  lastrow=zurueck(a);
1608 }
1609 else { lastrow=a; }
1610 if (ncols(lastrow)>ncols(m)) { spalten=ncols(lastrow); }
1611 else { spalten=ncols(m); }
1612 matrix mneu[zeile][spalten];
1613 for (i=1; i<nrows(m); i++) {
1614  mneu[i,1..ncols(m)]=m[i,1..ncols(m)];
1615 }
1616 mneu[zeile,1..ncols(lastrow)]=lastrow;
1617 if (lastrow[ncols(lastrow)]!=var(1)) {
1618  if (ncols(lastrow)==spalten) { v[zeile]=-1; }  // keine undefinierten Stellen
1619  else {
1620   v[zeile]=-Q-1;
1621   for (i=ncols(lastrow)+1; i<=spalten; i++) {
1622    mneu[zeile,i]=var(2);           // fuelle nicht def. Stellen der Matrix auf
1623 }}}
1624 else { v[zeile]=Q; }               // HNE war exakt
1625 if (ringwechsel) { kill extdguenstig; }
1626
1627 return(list(mneu,v,switch,f));
1628}
1629example
1630{
1631  if (nameof(basering)=="HNEring") {
1632   def rettering=HNEring;
1633   kill HNEring;
1634  }
1635  "EXAMPLE:"; echo = 2;
1636  ring exring=0,(x,y),dp;
1637  list hne=reddevelop(x14-3y2x11-y3x10-y2x9+3y4x8+y5x7+3y4x6+x5*(-y6+y5)-3y6x3-y7x2+y8);
1638  print(hne[1][1]); // finite HNE
1639  print(extdevelop(hne[1],5)[1]); pause;
1640  print(hne[2][1]); // HNE that can be extended
1641  list ehne=extdevelop(hne[2],5);
1642  print(ehne[1]); // new HN-matrix has 5 columns
1643  param(hne[2]);
1644  param(ehne);
1645  param(extdevelop(ehne,7)); pause;
1646  //////////////////////////
1647  print(develop(x-y2-2y3-3y4)[1]);
1648  print(extdevelop(develop(x-y2-2y3-3y4),4)[1]);
1649  print(extdevelop(develop(x-y2-2y3-3y4),10)[1]);
1650  kill HNEring,exring;
1651  echo = 0;
1652  if (defined(rettering)) { // wenn aktueller Ring vor Aufruf von example
1653   setring rettering;       // HNEring war, muss dieser erst wieder restauriert
1654   def HNEring=rettering;   // werden
1655   export HNEring;
1656  }
1657}
1658///////////////////////////////////////////////////////////////////////////////
1659
1660proc stripHNE (list l)
1661"USAGE:   stripHNE(l);    takes the output both of develop(f) and reddevelop(f)
1662           ( list l (matrix m, intvec v, int s[,poly g])
1663             or list of lists in the form l            )
1664RETURN:  a list in the same format as l, but all polynomials g are set to zero
1665NOTE:    The purpose of this procedure is to remove huge amounts of data
1666         no longer needed. It is useful, if one or more of the polynomials
1667         in l consumes much memory. It is still possible to compute invariants,
1668         parametrizations etc. with the stripped HNE(s), but it is not possible
1669         to use extdevelop with them.
1670EXAMPLE: example stripHNE;  shows an example
1671"
1672{
1673 list h;
1674 if (typeof(l[1])=="matrix") { l[4]=poly(0); }
1675 else {
1676  for (int i=1; i<=size(l); i++) {
1677    h=l[i];
1678    h[4]=poly(0);
1679    l[i]=h;
1680  }
1681 }
1682 return(l);
1683}
1684example
1685{ "EXAMPLE:"; echo = 2;
1686  ring r=0,(x,y),dp;
1687  list hne=develop(x2+y3+y4);
1688  hne;
1689  stripHNE(hne);
1690}
1691///////////////////////////////////////////////////////////////////////////////
1692
1693proc extractHNEs(list HNEs, int transvers)
1694"USAGE:  extractHNEs(HNEs,transvers);  list HNEs (output from HN),
1695        int transvers: 1 if x,y were exchanged, 0 else
1696RETURN: list of Hamburger-Noether-Extensions in the form of reddevelop
1697NOTE:   This procedure is only for internal purpose; examples don't make sense
1698"
1699{
1700 int i,maxspalte,hspalte,hnezaehler;
1701 list HNEaktu,Ergebnis;
1702 for (hnezaehler=1; hnezaehler<=size(HNEs); hnezaehler++) {
1703  maxspalte=0;
1704  HNEaktu=HNEs[hnezaehler];
1705  if (defined(Protokoll)) {"To process:";HNEaktu;}
1706  if (size(HNEs[hnezaehler])!=size(HNEs[hnezaehler][1])+2) {
1707     "The ideals and the hqs in HNEs[",hnezaehler,"] don't match!!";
1708     HNEs[hnezaehler];
1709  }
1710 //------------ ermittle  maximale Anzahl benoetigter Spalten: ----------------
1711  for (i=2; i<size(HNEaktu); i++) {
1712    hspalte=ncols(HNEaktu[i]);
1713    maxspalte=maxspalte*(hspalte < maxspalte)+hspalte*(hspalte >= maxspalte);
1714  }
1715 //------------- schreibe Ausgabe fuer hnezaehler-ten Zweig: ------------------
1716  matrix ma[size(HNEaktu)-2][maxspalte];
1717  for (i=1; i<=(size(HNEaktu)-2); i++) {
1718    if (ncols(HNEaktu[i+1]) > 1) {
1719      ma[i,1..ncols(HNEaktu[i+1])]=HNEaktu[i+1]; }
1720    else { ma[i,1]=HNEaktu[i+1][1];}
1721  }
1722  Ergebnis[hnezaehler]=list(ma,HNEaktu[1],transvers,HNEaktu[size(HNEaktu)]);
1723  kill ma;
1724 }
1725 return(Ergebnis);
1726}
1727///////////////////////////////////////////////////////////////////////////////
1728
1729proc factorfirst(poly f, int M, int N)
1730"USAGE : factorfirst(f,M,N); f poly, M,N int
1731RETURN: number d: f=c*(y^(N/e) - d*x^(M/e))^e with e=gcd(M,N), number c fitting
1732        0 if d does not exist
1733EXAMPLE: example factorfirst;  shows an example
1734"
1735{
1736 number c = koeff(f,0,N);
1737 number delta;
1738 int eps,l;
1739 int p=char(basering);
1740 string ringchar=charstr(basering);
1741
1742 if (c == 0) {"Something has gone wrong! I didn't get N correctly!"; exit;}
1743 int e = gcd(M,N);
1744
1745 if (p==0) { delta = koeff(f,M/ e,N - N/ e) / (-1*e*c); }
1746 else {
1747   if (e%p != 0) { delta = koeff(f,M/ e,N - N/ e) / (-1*e*c); }
1748   else {
1749     eps = e;
1750     for (l = 0; eps%p == 0; l=l+1) { eps=eps/ p;}
1751     if (defined(Protokoll)) {e," -> ",eps,"*",p,"^",l;}
1752     delta = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c);
1753
1754     if ((charstr(basering) != string(p)) and (delta != 0)) {
1755 //------ coefficient field is not Z/pZ => (p^l)th root is not identity -------
1756       delta=0;
1757       if (defined(Protokoll)) {
1758         "trivial factorization not implemented for",
1759         "parameters---I've to use 'factorize'";
1760       }
1761     }
1762   }
1763 }
1764 if (defined(Protokoll)) {"quasihomogeneous leading form:",f," = ",c,
1765        "* (y^"+string(N/ e),"-",delta,"* x^"+string(M/ e)+")^",e," ?";}
1766 if (f != c*(y^(N/ e) - delta*x^(M/ e))^e) {return(0);}
1767 else {return(delta);}
1768}
1769example
1770{ "EXAMPLE:"; echo = 2;
1771  ring exring=7,(x,y),dp;
1772  factorfirst(2*(y3-3x4)^5,20,15);
1773  factorfirst(x14+y7,14,7);
1774  factorfirst(x14+x8y3+y7,14,7);
1775}
1776///////////////////////////////////////////////////////////////////////////////
1777
1778
1779proc reddevelop (poly f)
1780"USAGE:   reddevelop(f); f poly
1781RETURN:  Hamburger-Noether development of f :
1782         A list of lists in the form of develop(f); each entry contains the
1783         data for one of the branches of f.
1784         For more details type 'help develop;'
1785NOTE:    as the Hamburger-Noether development of a reducible curve normally
1786         exists only in a ring extension, reddevelop automatically changes
1787         the basering to  ring HNEring=...,(x,y),ls; !
1788EXAMPLE: example reddevelop;  shows an example
1789"
1790{
1791 def altring = basering;
1792 int p = char(basering);                 // Ringcharakteristik
1793
1794 //-------------------- Tests auf Zulaessigkeit von basering ------------------
1795 if (charstr(basering)=="real") {
1796   " Singular cannot factorize over 'real' as ground field";
1797   return(list());
1798 }
1799 if (size(maxideal(1))<2) {
1800   " A univariate polynomial ring makes no sense !";
1801   return(list());
1802 }
1803 if (size(maxideal(1))>2) {
1804   " Warning: all but the first two variables are ignored!";
1805 }
1806 if (find(charstr(basering),string(char(basering)))!=1) {
1807   " ring of type (p^k,a) not implemented";
1808 //----------------------------------------------------------------------------
1809 // weder primitives Element noch factorize noch map "char p^k" -> "char -p"
1810 // [(p^k,a)->(p,a) ground field] noch fetch
1811 //----------------------------------------------------------------------------
1812   return(list());
1813 }
1814 //----------------- Definition eines neuen Ringes: HNEring -------------------
1815 string namex=varstr(1); string namey=varstr(2);
1816 if (string(char(altring))==charstr(altring)) {       // kein Parameter
1817   ring HNEring = char(altring),(x,y),ls;
1818   map m=altring,x,y;
1819   poly f=m(f);
1820   kill m;
1821 }
1822 else {
1823   string mipl=string(minpoly);
1824   if (mipl=="0") {
1825     " ** WARNING: No algebraic extension of this ground field is possible!";
1826     " ** We try to develop this polynomial, but if the need for an extension";
1827     " ** occurs during the calculation, we cannot proceed with the";
1828     " ** corresponding branches ...";
1829     execute "ring HNEring=("+charstr(basering)+"),(x,y),ls;";
1830 //--- ring ...=(char(.),`parstr()`),... geht nicht, wenn mehr als 1 Param. ---
1831   }
1832   else {
1833    string pa=parstr(altring);
1834    ring HNhelpring=p,`pa`,dp;
1835    execute "poly mipo="+mipl+";";   // Minimalpolynom in Polynom umgewandelt
1836    ring HNEring=(p,a),(x,y),ls;
1837    map getminpol=HNhelpring,a;
1838    mipl=string(getminpol(mipo));    // String umgewandelt mit 'a' als Param.
1839    execute "minpoly="+mipl+";";     // "minpoly=poly is not supported"
1840    kill HNhelpring, getminpol;
1841   }
1842   if (nvars(altring)==2) { poly f=fetch(altring,f); }
1843   else {
1844     map m=altring,x,y;
1845     poly f=m(f);
1846     kill m;
1847   }
1848 }
1849 export HNEring;
1850
1851 if (defined(Protokoll))
1852 {"received polynomial: ",f,", with x =",namex,", y =",namey;}
1853
1854 //----------------------- Variablendefinitionen ------------------------------
1855 int Abbruch,i,NullHNEx,NullHNEy;
1856 string str;
1857 list Newton,Ergebnis,hilflist;
1858
1859 //====================== Tests auf Zulaessigkeit des Polynoms ================
1860
1861 //-------------------------- Test, ob Einheit oder Null ----------------------
1862 if (subst(subst(f,x,0),y,0)!=0) {
1863   "The given polynomial is a unit in the power series ring!";
1864   keepring HNEring;
1865   return(list());                   // there are no HNEs
1866 }
1867 if (f==0) {
1868   "The given polynomial is zero!";
1869   keepring HNEring;
1870   return(list());                   // there are no HNEs
1871 }
1872
1873 //-----------------------  Test auf Quadratfreiheit --------------------------
1874
1875 if ((p==0) and (size(charstr(basering))==1)) {
1876
1877 //-------- Fall basering==0,... : Wechsel in Ring mit char >0 ----------------
1878 // weil squarefree eine Standardbasis berechnen muss (verwendet Syzygien)
1879 // -- wenn f in diesem Ring quadratfrei ist, dann erst recht im Ring HNEring
1880 //----------------------------------------------------------------------------
1881  int testerg=(polytest(f)==0);
1882  ring zweitring = 32003,(x,y),dp;
1883
1884  map polyhinueber=HNEring,x,y;         // fetch geht nicht
1885  poly f=polyhinueber(f);
1886  poly test_sqr=squarefree(f);
1887  if (test_sqr != f) {
1888   "Most probably the given polynomial is not squarefree. But the test was";
1889   "made in characteristic 32003 and not 0 to improve speed. You can";
1890   "(r) redo the test in char 0 (but this may take some time)";
1891   "(c) continue the development, if you're sure that the polynomial IS",
1892   "squarefree";
1893   if (testerg==1) {
1894     "(s) continue the development with a squarefree factor (*)";}
1895   "(q) or just quit the algorithm (default action)";
1896   "";"Please enter the letter of your choice:";
1897   str=read("")[1];                     // : liest 1 Zeichen von der Tastatur
1898   setring HNEring;
1899   map polyhinueber=zweitring,x,y;
1900   if (str=="r") {
1901     poly test_sqr=squarefree(f);
1902     if (test_sqr != f) {
1903      "The given polynomial is in fact not squarefree.";
1904      "I'll continue with the radical.";
1905      f=test_sqr;
1906     }
1907     else {"everything is ok -- the polynomial is squarefree in",
1908           "characteristic 0";}
1909   }
1910   else {
1911     if ((str=="s") and (testerg==1)) {
1912       "(*)attention: it could be that the factor is only one in char 32003!";
1913       f=polyhinueber(test_sqr);
1914     }
1915     else {
1916       if (str<>"c") {
1917         setring altring;kill HNEring;kill zweitring;
1918         return(list());}
1919       else { "if the algorithm doesn't terminate, you were wrong...";}
1920   }}
1921   kill zweitring;
1922   if (defined(Protokoll)) {"I continue with the polynomial",f; }
1923  }
1924  else {
1925    setring HNEring;
1926    kill zweitring;
1927  }
1928 }
1929 //------------------ Fall Char > 0 oder Ring hat Parameter -------------------
1930 else {
1931  poly test_sqr=squarefree(f);
1932  if (test_sqr != f) {
1933   if (test_sqr == 1) {
1934    "The given polynomial is of the form g^"+string(p)+",";
1935    "therefore not squarefree.  You can:";
1936    " (q) quit the algorithm (recommended) or";
1937    " (f) continue with the full radical (using a factorization of the";
1938    "     pure power part; this could take much time)";
1939    "";"Please enter the letter of your choice:";
1940    str=read("")[1];
1941    if (str<>"f") { str="q"; }
1942   }
1943   else {
1944    "The given polynomial is not squarefree.";
1945    if (p != 0)
1946     {
1947      " You can:";
1948      " (c) continue with a squarefree divisor (but factors of the form g^"
1949      +string(p);
1950      "     are lost; this is recommended, takes no more time)";
1951      " (f) continue with the full radical (using a factorization of the";
1952      "     pure power part; this could take much time)";
1953      " (q) quit the algorithm";
1954      "";"Please enter the letter of your choice:";
1955      str=read("")[1];
1956      if ((str<>"f") && (str<>"q")) { str="c"; }
1957     }
1958    else { "I'll continue with the radical."; str="c"; }
1959   }                                // endelse (test_sqr!=1)
1960   if (str=="q") {
1961    setring altring;kill HNEring;
1962    return(list());
1963   }
1964   if (str=="c") { f=test_sqr; }
1965   if (str=="f") { f=allsquarefree(f,test_sqr); }
1966  }
1967  if (defined(Protokoll)) {"I continue with the polynomial",f; }
1968
1969 }
1970 //====================== Ende Test auf Quadratfreiheit =======================
1971 if (subst(subst(f,x,0),y,0)!=0) {
1972   "Sorry. The remaining polynomial is a unit in the power series ring...";
1973   keepring HNEring;
1974   return(list());
1975 }
1976 //---------------------- Test, ob f teilbar durch x oder y -------------------
1977 if (subst(f,y,0)==0) {
1978   f=f/y; NullHNEy=1; }             // y=0 is a solution
1979 if (subst(f,x,0)==0) {
1980   f=f/x; NullHNEx=1; }             // x=0 is a solution
1981
1982 Newton=newtonpoly(f);
1983 i=1; Abbruch=0;
1984 //----------------------------------------------------------------------------
1985 // finde Eckpkt. des Newtonpolys, der den Teil abgrenzt, fuer den x transvers:
1986 // Annahme: Newton ist sortiert, s.d. Newton[1]=Punkt auf der y-Achse,
1987 // Newton[letzt]=Punkt auf der x-Achse
1988 //----------------------------------------------------------------------------
1989 while ((i<size(Newton)) and (Abbruch==0)) {
1990  if ((Newton[i+1][1]-Newton[i][1])>=(Newton[i][2]-Newton[i+1][2]))
1991   {Abbruch=1;}
1992  else {i=i+1;}
1993 }
1994 int grenze1=Newton[i][2];
1995 int grenze2=Newton[i][1];
1996 //----------------------------------------------------------------------------
1997 // Stelle Ring bereit zur Uebertragung der Daten im Fall einer Koerperer-
1998 // weiterung. Definiere Objekte, die spaeter uebertragen werden.
1999 // Binde die Listen (azeilen,...) an den Ring (um sie nicht zu ueberschreiben
2000 // bei Def. in einem anderen Ring).
2001 // Exportiere Objekte, damit sie auch in der proc HN noch da sind
2002 //----------------------------------------------------------------------------
2003 ring HNE_noparam = char(altring),(a,x,y),ls;
2004 export HNE_noparam;
2005 poly f;
2006 list azeilen=ideal(0);
2007 list HNEs=ideal(0);
2008 list aneu=ideal(0);
2009 ideal deltais;
2010 poly delta;                   // nicht number, weil delta von a abhaengen kann
2011 export f,azeilen,HNEs,aneu,deltais,delta;
2012 //----- hier steht die Anzahl bisher benoetigter Ringerweiterungen drin: -----
2013 int EXTHNEnumber=0; export EXTHNEnumber;
2014 setring HNEring;
2015
2016 // -------- Berechne HNE von allen Zweigen, fuer die x transvers ist: --------
2017 if (defined(Protokoll))
2018   {"1st step: Treat newton polygon until height",grenze1;}
2019 if (grenze1>0) {
2020  hilflist=HN(f,grenze1,1);
2021  if (typeof(hilflist[1][1])=="ideal") { hilflist[1]=list(); }
2022 //- fuer den Fall, dass keine Zweige in transz. Erw. berechnet werden konnten-
2023  Ergebnis=extractHNEs(hilflist[1],0);
2024  if (hilflist[2]!=-1) {
2025    if (defined(Protokoll)) {" ring change in HN(",1,") detected";}
2026    poly transfproc=hilflist[2];
2027    map hole=HNE_noparam,transfproc,x,y;
2028    setring HNE_noparam;
2029    f=imap(HNEring,f);
2030    setring EXTHNEring(EXTHNEnumber);
2031    poly f=hole(f);
2032  }
2033 }
2034 if (NullHNEy==1) {
2035  Ergebnis=Ergebnis+list(list(matrix(ideal(0,x)),intvec(1),int(0),poly(0)));
2036 }
2037 // --------------- Berechne HNE von allen verbliebenen Zweigen: --------------
2038 if (defined(Protokoll))
2039    {"2nd step: Treat newton polygon until height",grenze2;}
2040 if (grenze2>0) {
2041  map xytausch=basering,y,x;
2042  kill hilflist;
2043  def letztring=basering;
2044  if (EXTHNEnumber==0) { setring HNEring; }
2045  else                 { setring EXTHNEring(EXTHNEnumber); }
2046  list hilflist=HN(xytausch(f),grenze2,1);
2047  if (typeof(hilflist[1][1])=="ideal") { hilflist[1]=list(); }
2048  if (not(defined(Ergebnis))) {
2049 //-- HN wurde schon mal ausgefuehrt; Ringwechsel beim zweiten Aufruf von HN --
2050    if (defined(Protokoll)) {" ring change in HN(",1,") detected";}
2051    poly transfproc=hilflist[2];
2052    map hole=HNE_noparam,transfproc,x,y;
2053    setring HNE_noparam;
2054    list Ergebnis=imap(letztring,Ergebnis);
2055    setring EXTHNEring(EXTHNEnumber);
2056    list Ergebnis=hole(Ergebnis);
2057  }
2058  Ergebnis=Ergebnis+extractHNEs(hilflist[1],1);
2059 }
2060 if (NullHNEx==1) {
2061  Ergebnis=Ergebnis+list(list(matrix(ideal(0,x)),intvec(1),int(1),poly(0)));
2062 }
2063 //------------------- Loesche globale, nicht mehr benoetigte Objekte: --------
2064 if (EXTHNEnumber>0) {
2065  kill HNEring;
2066  def HNEring=EXTHNEring(EXTHNEnumber);
2067  setring HNEring;
2068  export HNEring;
2069  kill EXTHNEring(1..EXTHNEnumber);
2070 }
2071 kill HNE_noparam;
2072 kill EXTHNEnumber;
2073 keepring basering;
2074 "  ~~~ result:",size(Ergebnis)," branch(es) successfully computed ~~~";
2075 return(Ergebnis);
2076}
2077example
2078{
2079  if (nameof(basering)=="HNEring") {
2080   def rettering=HNEring;
2081   kill HNEring;
2082  }
2083  "EXAMPLE:"; echo = 2;
2084  ring r=0,(x,y),dp;
2085  list hne=reddevelop(x4-y6);
2086  size(hne);           // i.e. x4-y6 has two branches
2087  print(hne[1][1]);    // HN-matrix of 1st branch
2088  param(hne[1]);
2089  param(hne[2]);
2090  displayInvariants(hne);
2091  kill HNEring,r;
2092  pause;
2093  // a more interesting example:
2094  ring r = 32003,(x,y),dp;
2095  poly f = x25+x24-4x23-1x22y+4x22+8x21y-2x21-12x20y-4x19y2+4x20+10x19y+12x18y2
2096          -24x18y-20x17y2-4x16y3+x18+60x16y2+20x15y3-9x16y-80x14y3-10x13y4
2097          +36x14y2+60x12y4+2x11y5-84x12y3-24x10y5+126x10y4+4x8y6-126x8y5+84x6y6
2098          -36x4y7+9x2y8-1y9;
2099  list hne=reddevelop(f);
2100  size(hne);
2101  print(hne[1][1]);
2102  print(hne[4][1]);
2103  // a ring change was necessary, a is a parameter
2104  HNEring;
2105  kill HNEring,r;
2106  echo = 0;
2107  if (defined(rettering)) { // wenn aktueller Ring vor Aufruf von example
2108   setring rettering;       // HNEring war, muss dieser erst wieder restauriert
2109   def HNEring=rettering;   // werden
2110   export HNEring;
2111  }
2112}
2113///////////////////////////////////////////////////////////////////////////////
2114static
2115proc HN (poly f,int grenze, int Aufruf_Ebene)
2116"NOTE: This procedure is only for internal use, it is called via reddevelop
2117"
2118{
2119 //---------- Variablendefinitionen fuer den unverzweigten Teil: --------------
2120 if (defined(Protokoll)) {"procedure HN",Aufruf_Ebene;}
2121 int Abbruch,ende,i,j,e,M,N,Q,R,zeiger,zeile,zeilevorher;
2122 intvec hqs;
2123 poly fvorher;
2124 list erg=ideal(0); list HNEs=ideal(0); // um die Listen an den Ring zu binden
2125
2126 //-------------------- Bedeutung von Abbruch: --------------------------------
2127 //------- 0:keine Verzweigung | 1:Verzweigung,nicht fertig | 2:fertig --------
2128 //
2129 // Struktur von HNEs : Liste von Listen L (fuer jeden Zweig) der Form
2130 // L[1]=intvec (hqs), L[2],L[3],... ideal (die Zeilen (0,1,...) der HNE)
2131 // L[letztes]=poly (transformiertes f)
2132 //----------------------------------------------------------------------------
2133 list Newton;
2134 number delta;
2135 int p = char(basering);                // Ringcharakteristik
2136 list azeilen=ideal(0);
2137 ideal hilfid; list hilflist=ideal(0); intvec hilfvec;
2138
2139 // ======================= der unverzweigte Teil: ============================
2140 while (Abbruch==0) {
2141  Newton=newtonpoly(f);
2142  zeiger=find_in_list(Newton,grenze);
2143  if (Newton[zeiger][2] != grenze)
2144    {"Didn't find an edge in the newton polygon!";}
2145  if (zeiger==size(Newton)-1) {
2146    if (defined(Protokoll)) {"only one relevant side in newton polygon";}
2147    M=Newton[zeiger+1][1]-Newton[zeiger][1];
2148    N=Newton[zeiger][2]-Newton[zeiger+1][2];
2149    R = M%N;
2150    Q = M / N;
2151
2152 //-------- 1. Versuch: ist der quasihomogene Leitterm reine Potenz ? ---------
2153 //              (dann geht alles wie im irreduziblen Fall)
2154 //----------------------------------------------------------------------------
2155    e = gcd(M,N);
2156    delta=factorfirst(redleit(f,Newton[zeiger],Newton[zeiger+1])
2157                      /x^Newton[zeiger][1],M,N);
2158    if (delta==0) {
2159      if (defined(Protokoll)) {" The given polynomial is reducible !";}
2160      Abbruch=1;
2161    }
2162    if (Abbruch==0) {
2163 //-------------- f,zeile retten fuer den Spezialfall (###): ------------------
2164      fvorher=f;zeilevorher=zeile;
2165      if (R==0) {
2166 //------------- transformiere f mit T1, wenn kein Abbruch nachher: -----------
2167        if (N>1) { f = T1_Transform(f,delta,M/ e); }
2168        else     { ende=1; }
2169        if (defined(Protokoll)) {"a("+string(zeile)+","+string(Q)+") =",delta;}
2170        azeilen=set_list(azeilen,intvec(zeile+1,Q),delta);
2171      }
2172      else {
2173 //------------- R > 0 : transformiere f mit T2 -------------------------------
2174        erg=T2_Transform(f,delta,M,N);
2175        f=erg[1];delta=erg[2];
2176 //------- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: ---------
2177        while (R!=0) {
2178         if (defined(Protokoll)) { "h("+string(zeile)+") =",Q; }
2179         hqs[zeile+1]=Q;  //denn zeile beginnt mit dem Wert 0
2180 //------------------ markiere das Zeilenende der HNE: ------------------------
2181         azeilen=set_list(azeilen,intvec(zeile+1,Q+1),x);
2182         zeile=zeile+1;
2183 //----------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ---------
2184         azeilen[zeile+1]=ideal(0);
2185         M=N; N=R; R=M%N; Q=M / N;
2186        }
2187        if (defined(Protokoll)) {"a("+string(zeile)+","+string(Q)+") =",delta;}
2188        azeilen=set_list(azeilen,intvec(zeile+1,Q),delta);
2189      }
2190      if (defined(Protokoll)) {"transformed polynomial: ",f;}
2191      grenze=e;
2192 //----------------------- teste Abbruchbedingungen: --------------------------
2193      if (subst(f,y,0)==0) {              // <==> y|f
2194        "finite HNE of one branch found";
2195        azeilen=set_list(azeilen,intvec(zeile+1,Q+1),x);
2196 //- Q wird nur in hqs eingetragen, wenn der Spezialfall nicht eintritt (s.u.)-
2197        Abbruch=2;
2198        if (grenze>1) {
2199         if (jet(f,1,intvec(0,1))==0) {
2200 //------ jet(...)=alle Monome von f, die nicht durch y2 teilbar sind ---------
2201 "THE TEST FOR SQUAREFREENESS WAS BAD!! The polynomial was NOT squarefree!!!";}
2202         else {
2203 //-------------------------- Spezialfall (###): ------------------------------
2204 // Wir haben das Problem, dass die HNE eines Zweiges hier abbricht, aber ein
2205 // anderer Zweig bis hierher genau die gleiche HNE hat, die noch weiter geht
2206 // Loesung: mache Transform. rueckgaengig und behandle f im Verzweigungsteil
2207 //----------------------------------------------------------------------------
2208          Abbruch=1;
2209          f=fvorher;zeile=zeilevorher;grenze=Newton[zeiger][2];
2210        }}
2211        else {f=0;}     // f nicht mehr gebraucht - spare Speicher
2212        if (Abbruch==2) { hqs[zeile+1]=Q; }
2213      }                 // Spezialfall nicht eingetreten
2214      else {
2215        if (ende==1) {
2216          "HNE of one branch found"; Abbruch=2; hqs[zeile+1]=-Q-1;}
2217      }
2218    }                   // end(if Abbruch==0)
2219  }                     // end(if zeiger...)
2220  else { Abbruch=1;}
2221 }                      // end(while Abbruch==0)
2222
2223 // ===================== der Teil bei Verzweigung: ===========================
2224
2225 if (Abbruch==1) {
2226 //---------- Variablendefinitionen fuer den verzweigten Teil: ----------------
2227  poly leitf,teiler,transformiert;
2228  list faktoren=ideal(0); list aneu=ideal(0);
2229  ideal deltais;
2230  intvec eis;
2231  int zaehler,hnezaehler,zl,zl1,M1,N1,R1,Q1,needext;
2232  int numberofRingchanges,lastRingnumber,ringischanged,flag;
2233  string letztringname;
2234
2235  zeiger=find_in_list(Newton,grenze);
2236  if (defined(Protokoll)) {
2237    "Branching part reached---newton polygon :",Newton;
2238    "relevant part until height",grenze,", from",Newton[zeiger],"on";
2239  }
2240  azeilen=list(hqs)+azeilen; // hat jetzt Struktur von HNEs: hqs in der 1.Zeile
2241
2242 //======= Schleife fuer jede zu betrachtende Seite des Newtonpolygons: =======
2243  for(i=zeiger; i<size(Newton); i++) {
2244   if (defined(Protokoll)) { "we consider side",Newton[i],Newton[i+1]; }
2245   M=Newton[i+1][1]-Newton[i][1];
2246   N=Newton[i][2]-Newton[i+1][2];
2247   R = M%N;
2248   Q = M / N;
2249   e=gcd(M,N);
2250   needext=1;
2251   letztringname=nameof(basering);
2252   lastRingnumber=EXTHNEnumber;
2253
2254 //-- wechsle so lange in Ringerw., bis Leitform vollst. in Linearfakt. zerf.:-
2255   for (numberofRingchanges=1; needext==1; numberofRingchanges++) {
2256    leitf=redleit(f,Newton[i],Newton[i+1])/(x^Newton[i][1]*y^Newton[i+1][2]);
2257    delta=factorfirst(leitf,M,N);
2258    needext=0;
2259    if (delta==0) {
2260
2261 //---------- Sonderbehandlung: faktorisere einige Polynome ueber Q(a): -------
2262     if (charstr(basering)=="0,a") {
2263       faktoren,flag=extrafactor(leitf,M,N);  // damit funktion. Bsp. Baladi 5
2264       if (flag==0)
2265         { " ** I cannot proceed, will produce some error messages...";}
2266     }
2267     else
2268      {
2269 //------------------ faktorisiere das charakt. Polynom: ----------------------
2270      faktoren=factorize(charPoly(leitf,M,N));
2271     }
2272
2273     zaehler=1; eis=0;
2274     for (j=1; j<=size(faktoren[2]); j++) {
2275      teiler=faktoren[1][j];
2276      if (teiler/y != 0) {         // sonst war's eine Einheit --> wegwerfen!
2277        if (defined(Protokoll)) {"factor of leading form found:",teiler;}
2278        if (teiler/y2 == 0) {      // --> Faktor hat die Form cy+d
2279          deltais[zaehler]=-subst(teiler,y,0)/koeff(teiler,0,1); //=-d/c
2280          eis[zaehler]=faktoren[2][j];
2281          zaehler++;
2282        }
2283        else {
2284          " Change of basering necessary!!";
2285          if (defined(Protokoll)) { teiler,"is not properly factored!"; }
2286          if (needext==0) { poly zerlege=teiler; }
2287          needext=1;
2288        }
2289      }
2290     }                             // end(for j)
2291    }
2292    else { deltais=ideal(delta); eis=e;}
2293    if (defined(Protokoll)) {"roots of char. poly:";deltais;
2294                             "with multiplicities:",eis;}
2295    if (needext==1) {
2296 //--------------------- fuehre den Ringwechsel aus: --------------------------
2297      ringischanged=1;
2298      if ((size(parstr(basering))>0) && string(minpoly)=="0") {
2299        " ** We've had bad luck! The HNE cannot completely be calculated!";
2300                                   // HNE in transzendenter Erw. fehlgeschlagen
2301        kill zerlege;
2302        ringischanged=0; break;    // weiter mit gefundenen Faktoren
2303      }
2304      if (parstr(basering)=="") {
2305        EXTHNEnumber++;
2306        splitring(zerlege,"EXTHNEring("+string(EXTHNEnumber)+")");
2307        poly transf=0;
2308        poly transfproc=0;
2309      }
2310      else {
2311        if (defined(translist)) { kill translist; } // Vermeidung einer Warnung
2312        if (numberofRingchanges>1) {  // ein Ringwechsel hat nicht gereicht
2313         list translist=splitring(zerlege,"",list(transf,transfproc));
2314         poly transf=translist[1]; poly transfproc=translist[2];
2315        }
2316        else {
2317         if (defined(transfproc)) { // in dieser proc geschah schon Ringwechsel
2318          EXTHNEnumber++;
2319          list translist=splitring(zerlege,"EXTHNEring("
2320               +string(EXTHNEnumber)+")",list(a,transfproc));
2321          poly transf=translist[1];
2322          poly transfproc=translist[2];
2323         }
2324         else {
2325          EXTHNEnumber++;
2326          list translist=splitring(zerlege,"EXTHNEring("
2327               +string(EXTHNEnumber)+")",a);
2328          poly transf=translist[1];
2329          poly transfproc=transf;
2330        }}
2331      }
2332 //----------------------------------------------------------------------------
2333 // transf enthaelt jetzt den alten Parameter des Ringes, der aktiv war vor
2334 // Beginn der Schleife (evtl. also ueber mehrere Ringwechsel weitergereicht),
2335 // transfproc enthaelt den alten Parm. des R., der aktiv war zu Beginn der
2336 // Prozedur, und der an die aufrufende Prozedur zurueckgegeben werden muss
2337 // transf ist Null, falls der alte Ring keinen Parameter hatte,
2338 // das gleiche gilt fuer transfproc
2339 //----------------------------------------------------------------------------
2340
2341 //------ Neudef. von Variablen, Uebertragung bisher errechneter Daten: -------
2342      poly leitf,teiler,transformiert;
2343      list faktoren; list aneu=ideal(0);
2344      ideal deltais;
2345      number delta;
2346      setring HNE_noparam;
2347      if (defined(letztring)) { kill letztring; }
2348      if (lastRingnumber>0) { def letztring=EXTHNEring(lastRingnumber); }
2349      else                  { def letztring=HNEring; }
2350      f=imap(letztring,f);
2351      setring EXTHNEring(EXTHNEnumber);
2352      map hole=HNE_noparam,transf,x,y;
2353      poly f=hole(f);
2354    }
2355   }    // end (while needext==1) bzw. for (numberofRingchanges)
2356
2357   if (eis==0) { i++; continue; }
2358   if (ringischanged==1) {
2359    list erg,hilflist;              // dienen nur zum Sp. von Zwi.erg.
2360    ideal hilfid;
2361    erg=ideal(0); hilflist=erg;
2362
2363    hole=HNE_noparam,transf,x,y;
2364    setring HNE_noparam;
2365    azeilen=imap(letztring,azeilen);
2366    HNEs=imap(letztring,HNEs);
2367
2368    setring EXTHNEring(EXTHNEnumber);
2369    list azeilen=hole(azeilen);
2370    list HNEs=hole(HNEs);
2371    kill letztring;
2372    ringischanged=0;
2373   }
2374
2375 //============ Schleife fuer jeden gefundenen Faktor der Leitform: ===========
2376   for (j=1; j<=size(eis); j++) {
2377 //-- Mache Transf. T1 oder T2, trage Daten in HNEs ein, falls HNE abbricht: --
2378
2379 //------------------------ Fall R==0: ----------------------------------------
2380    if (R==0) {
2381      transformiert = T1_Transform(f,number(deltais[j]),M/ e);
2382      if (defined(Protokoll)) {
2383        "a("+string(zeile)+","+string(Q)+") =",deltais[j];
2384        "transformed polynomial: ",transformiert;
2385      }
2386      if (subst(transformiert,y,0)==0) {
2387       "finite HNE found";
2388       hnezaehler++;
2389 //------------ trage deltais[j],x ein in letzte Zeile, fertig: ---------------
2390       HNEs[hnezaehler]=azeilen+list(poly(0));
2391       hilfid=HNEs[hnezaehler][zeile+2]; hilfid[Q]=deltais[j]; hilfid[Q+1]=x;
2392       HNEs=set_list(HNEs,intvec(hnezaehler,zeile+2),hilfid);
2393       HNEs=set_list(HNEs,intvec(hnezaehler,1,zeile+1),Q);
2394                                  // aktualisiere Vektor mit den hqs
2395       if (eis[j]>1) {
2396        transformiert=transformiert/y;
2397        if (subst(transformiert,y,0)==0) {
2398 "THE TEST FOR SQUAREFREENESS WAS BAD!! The polynomial was NOT squarefree!!!";}
2399        else {
2400 //------ Spezialfall (###) eingetreten: Noch weitere Zweige vorhanden --------
2401          eis[j]=eis[j]-1;
2402        }
2403       }
2404      }
2405    }
2406    else {
2407 //------------------------ Fall R <> 0: --------------------------------------
2408      erg=T2_Transform(f,number(deltais[j]),M,N);
2409      transformiert=erg[1];delta=erg[2];
2410      if (defined(Protokoll)) {"transformed polynomial: ",transformiert;}
2411      if (subst(transformiert,y,0)==0) {
2412       "finite HNE found";
2413       hnezaehler++;
2414 //---------------- trage endliche HNE in HNEs ein: ---------------------------
2415       HNEs[hnezaehler]=azeilen;  // dupliziere den gemeins. Anfang der HNE's
2416       zl=2;
2417 //----------------------------------------------------------------------------
2418 // Werte von:  zeile: aktuelle Zeilennummer der HNE (gemeinsamer Teil)
2419 //             zl   : die HNE spaltet auf; zeile+zl ist der Index fuer die
2420 // Zeile des aktuellen Zweigs; (zeile+zl-2) ist die tatsaechl. Zeilennr.
2421 // (bei 0 angefangen) der HNE  ([1] <- intvec(hqs), [2] <- 0. Zeile usw.)
2422 //----------------------------------------------------------------------------
2423
2424 //---------- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: ------
2425       M1=M;N1=N;R1=R;Q1=M1/ N1;
2426       while (R1!=0) {
2427        if (defined(Protokoll)) { "h("+string(zeile+zl-2)+") =",Q1; }
2428        HNEs=set_list(HNEs,intvec(hnezaehler,1,zeile+zl-1),Q1);
2429        HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl,Q1+1),x);
2430                                  // markiere das Zeilenende der HNE
2431        zl=zl+1;
2432 //-------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ------------
2433        HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl),ideal(0));
2434
2435        M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1;
2436       }
2437       if (defined(Protokoll)) {
2438         "a("+string(zeile+zl-2)+","+string(Q1)+") =",delta;
2439       }
2440       hilfid=ideal(0);           // = HNEs[hnezaehler][zeile+zl];
2441       hilfid[Q1]=delta;hilfid[Q1+1]=x;
2442       HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl),hilfid);
2443       HNEs=set_list(HNEs,intvec(hnezaehler,1,zeile+zl-1),Q1);
2444                                  // aktualisiere Vektor mit hqs
2445       HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl+1),poly(0));
2446 //-------------------- Ende der Eintragungen in HNEs -------------------------
2447
2448       if (eis[j]>1) {
2449        transformiert=transformiert/y;
2450        if (subst(transformiert,y,0)==0) {
2451 "THE TEST FOR SQUAREFREENESS WAS BAD!! The polynomial was NOT squarefree!!!";}
2452         else {
2453 //--------- Spezialfall (###) eingetreten: Noch weitere Zweige vorhanden -----
2454          eis[j]=eis[j]-1;
2455       }}
2456      }                           // endif (subst()==0)
2457    }                             // endelse (R<>0)
2458
2459 //========== Falls HNE nicht abbricht: Rekursiver Aufruf von HN: =============
2460 //------------------- Berechne HNE von transformiert -------------------------
2461    if (subst(transformiert,y,0)!=0) {
2462     lastRingnumber=EXTHNEnumber;
2463     list HNerg=HN(transformiert,eis[j],Aufruf_Ebene+1);
2464     if (HNerg[2]==-1) {          // kein Ringwechsel in HN aufgetreten
2465       aneu=HNerg[1];  }
2466     else {
2467       if (defined(Protokoll))
2468          {" ring change in HN(",Aufruf_Ebene+1,") detected";}
2469       list aneu=HNerg[1];
2470       poly transfproc=HNerg[2];
2471
2472 //- stelle lokale Var. im neuen Ring wieder her und rette ggf. ihren Inhalt: -
2473       list erg,hilflist,faktoren;
2474       ideal hilfid;
2475       erg=ideal(0); hilflist=erg; faktoren=erg;
2476       poly leitf,teiler,transformiert;
2477
2478       map hole=HNE_noparam,transfproc,x,y;
2479       setring HNE_noparam;
2480       if (lastRingnumber>0) { def letztring=EXTHNEring(lastRingnumber); }
2481       else                  { def letztring=HNEring; }
2482       HNEs=imap(letztring,HNEs);
2483       azeilen=imap(letztring,azeilen);
2484       deltais=imap(letztring,deltais);
2485       delta=imap(letztring,delta);
2486       f=imap(letztring,f);
2487
2488       setring EXTHNEring(EXTHNEnumber);
2489       list HNEs=hole(HNEs);
2490       list azeilen=hole(azeilen);
2491       ideal deltais=hole(deltais);
2492       number delta=number(hole(delta));
2493       poly f=hole(f);
2494     }
2495     kill HNerg;
2496 //----------------------------------------------------------------------------
2497 // HNerg muss jedesmal mit "list" neu definiert werden, weil vorher noch nicht
2498 // ------- klar ist, ob der Ring nach Aufruf von HN noch derselbe ist --------
2499
2500 //============= Verknuepfe bisherige HNE mit von HN gelieferten HNEs: ========
2501     if (R==0) {
2502       HNEs,hnezaehler=constructHNEs(HNEs,hnezaehler,aneu,azeilen,zeile,
2503                       deltais,Q,j);
2504     }
2505     else {
2506      for (zaehler=1; zaehler<=size(aneu); zaehler++) {
2507       hnezaehler++;
2508       HNEs[hnezaehler]=azeilen; // dupliziere den gemeinsamen Anfang der HNE's
2509       zl=2;
2510 //---------------- Trage Beitrag dieser Transformation T2 ein: ---------------
2511 //--------- Zur Bedeutung von zeile, zl: siehe Kommentar weiter oben ---------
2512
2513 //--------- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: -------
2514       M1=M;N1=N;R1=R;Q1=M1/ N1;
2515       while (R1!=0) {
2516        if (defined(Protokoll)) { "h("+string(zeile+zl-2)+") =",Q1; }
2517        HNEs=set_list(HNEs,intvec(hnezaehler,1,zeile+zl-1),Q1);
2518        HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl,Q1+1),x);
2519                                 // Markierung des Zeilenendes der HNE
2520        zl=zl+1;
2521 //-------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ------------
2522        HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl),ideal(0));
2523        M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1;
2524       }
2525       if (defined(Protokoll)) {
2526         "a("+string(zeile+zl-2)+","+string(Q1)+") =",delta;
2527       }
2528       HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl,Q1),delta);
2529
2530 //--- Daten aus T2_Transform sind eingetragen; haenge Daten von HN an: -------
2531       hilfid=HNEs[hnezaehler][zeile+zl];
2532       for (zl1=Q1+1; zl1<=ncols(aneu[zaehler][2]); zl1++) {
2533        hilfid[zl1]=aneu[zaehler][2][zl1];
2534       }
2535       HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl),hilfid);
2536 //--- vorher HNEs[.][zeile+zl]<-aneu[.][2], jetzt [zeile+zl+1] <- [3] usw.: --
2537       for (zl1=3; zl1<=size(aneu[zaehler]); zl1++) {
2538        HNEs=set_list(HNEs,intvec(hnezaehler,zeile+zl+zl1-2),
2539             aneu[zaehler][zl1]);
2540       }
2541 //--- setze die hqs zusammen: HNEs[hnezaehler][1]=HNEs[..][1],aneu[..][1] ----
2542       hilfvec=HNEs[hnezaehler][1],aneu[zaehler][1];
2543       HNEs=set_list(HNEs,intvec(hnezaehler,1),hilfvec);
2544 //----------- weil nicht geht: liste[1]=liste[1],aneu[zaehler][1] ------------
2545      }                     // end(for zaehler)
2546     }                      // endelse (R<>0)
2547    }                       // endif (subst()!=0)  (weiteres Aufblasen mit HN)
2548
2549   }                        // end(for j) (Behandlung der einzelnen delta_i)
2550
2551  }
2552  keepring basering;
2553  if (defined(transfproc)) { return(list(HNEs,transfproc)); }
2554  else                     { return(list(HNEs,poly(-1))); }
2555 // -1 als 2. Rueckgabewert zeigt an, dass kein Ringwechsel stattgefunden hat -
2556 }
2557 else {
2558  HNEs[1]=list(hqs)+azeilen+list(f); // f ist das transform. Poly oder Null
2559  keepring basering;
2560  return(list(HNEs,poly(-1)));
2561 //-- in dieser proc trat keine Verzweigung auf, also auch kein Ringwechsel ---
2562 }
2563}
2564///////////////////////////////////////////////////////////////////////////////
2565
2566static
2567proc constructHNEs (list HNEs,int hnezaehler,list aneu,list azeilen,int zeile,ideal deltais,int Q,int j)
2568"NOTE: This procedure is only for internal use, it is called via HN
2569"
2570{
2571  int zaehler,zl;
2572  ideal hilfid;
2573  list hilflist;
2574  intvec hilfvec;
2575  for (zaehler=1; zaehler<=size(aneu); zaehler++) {
2576     hnezaehler++;
2577     HNEs[hnezaehler]=azeilen;            // dupliziere gemeins. Anfang
2578 //----------------------- trage neu berechnete Daten ein ---------------------
2579     hilfid=HNEs[hnezaehler][zeile+2];
2580     hilfid[Q]=deltais[j];
2581     for (zl=Q+1; zl<=ncols(aneu[zaehler][2]); zl++) {
2582      hilfid[zl]=aneu[zaehler][2][zl];
2583     }
2584     hilflist=HNEs[hnezaehler];hilflist[zeile+2]=hilfid;
2585 //------------------ haenge uebrige Zeilen von aneu[] an: --------------------
2586     for (zl=3; zl<=size(aneu[zaehler]); zl++) {
2587      hilflist[zeile+zl]=aneu[zaehler][zl];
2588     }
2589 //--- setze die hqs zusammen: HNEs[hnezaehler][1]=HNEs[..][1],aneu[..][1] ----
2590     if (hilflist[1]==0) { hilflist[1]=aneu[zaehler][1]; }
2591     else { hilfvec=hilflist[1],aneu[zaehler][1];hilflist[1]=hilfvec; }
2592     HNEs[hnezaehler]=hilflist;
2593  }
2594  return(HNEs,hnezaehler);
2595}
2596///////////////////////////////////////////////////////////////////////////////
2597
2598proc extrafactor (poly leitf, int M, int N)
2599"USAGE:   factors,flag=extrafactor(f,M,N);
2600         list factors, int flag,M,N, poly f in x,y
2601ASSUME:  basering is (0,a),(x,y),ds or ls
2602         Newtonpolygon of f is one side with height N, length M
2603RETURN:  for some special f, factors is a list of the factors of
2604         charPoly(f,M,N), same format as factorize
2605         (see help charPoly; help factorize)
2606         In this case, flag!=0 (TRUE). If extrafactor cannot factorize f,
2607         flag will be 0 (FALSE), factors will be the empty list.
2608NOTE:    This procedure is designed to make reddevelop able to compute some
2609         special cases that need a ring extension in char 0.
2610         It becomes obsolete as soon as factorize works also in such rings.
2611EXAMPLE: example extrafactor;  shows an example
2612"
2613{
2614 list faktoren;
2615
2616      if (a2==-1) {
2617       poly testpol=charPoly(leitf,M,N);
2618       if (testpol==1+2y2+y4) {
2619         faktoren=list(ideal(y+a,y-a),intvec(2,2));
2620         testpol=0;
2621       }
2622 //------------ ist Poly von der Form q*i+r*y^n, n in N, q,r in Q ?: ----------
2623       if ((size(testpol)==2) && (find(string(lead(testpol)),"a")>0)
2624           && (find(string(testpol-lead(testpol)),"a")==0)) {
2625        faktoren=list(ideal(testpol),intvec(1));
2626        testpol=0;
2627       }
2628      }
2629      if (a4==-625) {
2630        poly testpol=charPoly(leitf,M,N);
2631        if (testpol==4a2-4y2)
2632         { faktoren=list(ideal(-4,y+a,y-a),intvec(1,1,1)); testpol=0;}
2633        if (testpol==-4a2-4y2)
2634         { faktoren=list(ideal(-4,y+1/25*a3,y-1/25*a3),intvec(1,1,1));
2635           testpol=0;}
2636      }
2637      if (defined(testpol)==0) { poly testpol=1; }
2638      if (testpol!=0) {
2639        "factorize not implemented in char (0,a)!";
2640        "could not factorize:",charPoly(leitf,M,N); pause;
2641      }
2642 return(faktoren,(testpol==0)); // Test: faktoren==list() geht leider nicht
2643}
2644example
2645{ "EXAMPLE:"; echo=2;
2646  ring r=(0,a),(x,y),ds;
2647  minpoly=a2+1;
2648  poly f=x4+2x2y2+y4;
2649  charPoly(f,4,4);
2650  list factors;
2651  int flag;
2652  factors,flag=extrafactor(f,4,4);
2653  if (flag!=0) {factors;}
2654}
Note: See TracBrowser for help on using the repository browser.