/////////////////////////////////////////////////////////////////////////////// version="$Id$"; category="Singularities"; info=" LIBRARY: hnoether.lib Hamburger-Noether (Puiseux) Expansion AUTHORS: Martin Lamm, lamm@mathematik.uni-kl.de Christoph Lossen, lossen@mathematik.uni-kl.de OVERVIEW: A library for computing the Hamburger-Noether expansion (analogue of Puiseux expansion over fields of arbitrary characteristic) of a reduced plane curve singularity following [Campillo, A.: Algebroid curves in positive characteristic, Springer LNM 813 (1980)]. @* The library contains also procedures for computing the (topological) numerical invariants of plane curve singularities. PROCEDURES: hnexpansion(f [,\"ess\"]); Hamburger-Noether (HN) expansion of f develop(f [,n]); HN expansion of irreducible plane curve germs extdevelop(hne,n); extension of the H-N expansion hne of f param(hne [,s]); parametrization of branches described by HN data displayHNE(hne); display HN expansion as an ideal invariants(hne); invariants of f, e.g. the characteristic exponents displayInvariants(hne); display invariants of f multsequence(hne); sequence of multiplicities displayMultsequence(hne); display sequence of multiplicities intersection(hne1,hne2); intersection multiplicity of two local branches is_irred(f); test whether f is irreducible as power series delta(f); delta invariant of f newtonpoly(f); (local) Newton polygon of f is_NND(f); test whether f is Newton non-degenerate stripHNE(hne); reduce amount of memory consumed by hne puiseux2generators(m,n); convert Puiseux pairs to generators of semigroup separateHNE(hne1,hne2); number of quadratic transf. needed for separation squarefree(f); a squarefree divisor of the polynomial f allsquarefree(f,l); the maximal squarefree divisor of the polynomial f further_hn_proc(); show further procedures useful for interactive use KEYWORDS: Hamburger-Noether expansion; Puiseux expansion; curve singularities "; // essdevelop(f); HN expansion of essential branches // multiplicities(hne); multiplicities of blowed up curves /////////////////////////////////////////////////////////////////////////////// LIB "primitiv.lib"; LIB "inout.lib"; LIB "sing.lib"; /////////////////////////////////////////////////////////////////////////////// proc further_hn_proc() "USAGE: further_hn_proc(); NOTE: The library @code{hnoether.lib} contains some more procedures which are not shown when typing @code{help hnoether.lib;}. They may be useful for interactive use (e.g. if you want to do the calculation of an HN development \"by hand\" to see the intermediate results), and they can be enumerated by calling @code{further_hn_proc()}. @* Use @code{help ;} for detailed information about each of them. " { " The following procedures are also part of `hnoether.lib': getnm(f); intersection pts. of Newton polygon with axes T_Transform(f,Q,N); returns f(y,xy^Q)/y^NQ (f: poly, Q,N: int) T1_Transform(f,d,M); returns f(x,y+d*x^M) (f: poly,d:number,M:int) T2_Transform(f,d,M,N,ref); a composition of T1 & T koeff(f,I,J); gets coefficient of indicated monomial of polynomial f redleit(f,S,E); restriction of monomials of f to line (S-E) leit(f,n,m); special case of redleit (for irred. polynomials) testreducible(f,n,m); tests whether f is reducible charPoly(f,M,N); characteristic polynomial of f find_in_list(L,p); find int p in list L get_last_divisor(M,N); last divisor in Euclid's algorithm factorfirst(f,M,N); try to factor f without `factorize' factorlist(L); factorize a list L of polynomials referencepoly(D); a polynomial f s.t. D is the Newton diagram of f"; // static procedures not useful for interactive use: // polytest(f); tests coefficients and exponents of polynomial f // extractHNEs(H,t); extracts output H of HN to output of hnexpansion // HN(f,grenze); recursive subroutine for hnexpansion // constructHNEs(...); subroutine for HN } example { echo=2; further_hn_proc(); } /////////////////////////////////////////////////////////////////////////////// proc getnm (poly f) "USAGE: getnm(f); f bivariate polynomial RETURN: intvec(n,m) : (0,n) is the intersection point of the Newton polygon of f with the y-axis, n=-1 if it doesn't exist (m,0) is its intersection point with the x-axis, m=-1 if this point doesn't exist ASSUME: ring has ordering `ls' or `ds' EXAMPLE: example getnm; shows an example " { // assume being called by develop ==> ring ordering is ls (ds would also work) return(ord(subst(f,var(1),0)),ord(subst(f,var(2),0))); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),ds; poly f = x5+x4y3-y2+y4; getnm(f); } /////////////////////////////////////////////////////////////////////////////// proc leit (poly f, int n, int m) "USAGE: leit(f,n,m); poly f, int n,m RETURN: all monomials on the line from (0,n) to (m,0) in the Newton diagram EXAMPLE: example leit; shows an example " { return(jet(f,m*n,intvec(n,m))-jet(f,m*n-1,intvec(n,m))) } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),ds; poly f = x5+x4y3-y2+y4; leit(f,2,5); } /////////////////////////////////////////////////////////////////////////////// proc testreducible (poly f, int n, int m) "USAGE: testreducible(f,n,m); f poly, n,m int RETURN: 1 if there are points in the Newton diagram below the line (0,n)-(m,0) 0 else EXAMPLE: example testreducible; shows an example " { return(size(jet(f,m*n-1,intvec(n,m))) != 0) } example { "EXAMPLE:"; echo = 2; ring rg=0,(x,y),ls; testreducible(x2+y3-xy4,3,2); } /////////////////////////////////////////////////////////////////////////////// proc T_Transform (poly f, int Q, int N) "USAGE: T_Transform(f,Q,N); f poly, Q,N int RETURN: f(y,xy^Q)/y^NQ if x,y are the ring variables NOTE: this is intended for irreducible power series f EXAMPLE: example T_Transform; shows an example " { map T = basering,var(2),var(1)*var(2)^Q; return(T(f)/var(2)^(N*Q)); } example { "EXAMPLE:"; echo = 2; ring exrg=0,(x,y),ls; export exrg; T_Transform(x3+y2-xy3,1,2); kill exrg; } /////////////////////////////////////////////////////////////////////////////// proc T1_Transform (poly f, number d, int Q) "USAGE: T1_Transform(f,d,Q); f poly, d number, Q int RETURN: f(x,y+d*x^Q) if x,y are the ring variables EXAMPLE: example T1_Transform; shows an example " { map T1 = basering,var(1),var(2)+d*var(1)^Q; return(T1(f)); } example { "EXAMPLE:"; echo = 2; ring exrg=0,(x,y),ls; export exrg; T1_Transform(y2-2xy+x2+x2y,1,1); kill exrg; } /////////////////////////////////////////////////////////////////////////////// proc T2_Transform (poly f_neu, number d, int M, int N, poly refpoly) "USAGE: T2_Transform(f,d,M,N,ref); f poly, d number; M,N int; ref poly RETURN: list: poly T2(f,d',M,N), number d' in \{ d, 1/d \} ASSUME: ref has the same Newton polygon as f (but can be simpler) for this you can e.g. use the proc `referencepoly' or simply f again COMMENT: T2 is a composition of T_Transform and T1_Transform; the exact definition can be found in Rybowicz: `Sur le calcul des places ...' or in Lamm: `Hamburger-Noether-Entwicklung von Kurvensingularitaeten' SEE ALSO: T_Transform, T1_Transform, referencepoly EXAMPLE: example T2_Transform; shows an example " { //---------------------- compute gcd and extgcd of N,M ----------------------- int ggt=gcd(M,N); M=M/ggt; N=N/ggt; list ts=extgcd(M,N); int tau,sigma=ts[2],-ts[3]; int s,t; poly xp=var(1); poly yp=var(2); poly hilf; if (sigma<0) { tau=-tau; sigma=-sigma;} // es gilt: 0<=tau<=N, 0<=sigma<=M, |N*sigma-M*tau| = 1 = ggT(M,N) if (N*sigma < M*tau) { d = 1/d; } //--------------------------- euklid. Algorithmus ---------------------------- int R; int M1,N1=M,N; for ( R=M1%N1; R!=0; ) { M1=N1; N1=R; R=M1%N1;} int Q=M1 / N1; map T1 = basering,xp,yp+d*xp^Q; map Tstar=basering,xp^(N-Q*tau)*yp^tau,xp^(M-sigma*Q)*yp^sigma; if (defined(HNDebugOn)) { "Trafo. T2: x->x^"+string(N-Q*tau)+"*y^"+string(tau)+", y->x^" +string(M-sigma*Q)+"*y^"+string(sigma); "delt =",d,"Q =",Q,"tau,sigma =",tau,sigma; } //------------------- Durchfuehrung der Transformation T2 -------------------- f_neu=Tstar(f_neu); refpoly=Tstar(refpoly); //--- dividiere f_neu so lange durch x & y, wie die Division aufgeht, // benutze ein Referenzpolynom mit gleichem Newtonpolynom wie f_neu zur // Beschleunigung: --- for (hilf=refpoly/xp; hilf*xp==refpoly; hilf=refpoly/xp) {refpoly=hilf; s++;} for (hilf=refpoly/yp; hilf*yp==refpoly; hilf=refpoly/yp) {refpoly=hilf; t++;} f_neu=f_neu/(xp^s*yp^t); return(list(T1(f_neu),d)); } example { "EXAMPLE:"; echo = 2; ring exrg=0,(x,y),ds; export exrg; poly f=y2-2x2y+x6-x5y+x4y2; T2_Transform(f,1/2,4,1,f); T2_Transform(f,1/2,4,1,referencepoly(newtonpoly(f,1))); // if size(referencepoly) << size(f) the 2nd example would be faster referencepoly(newtonpoly(f,1)); kill exrg; } /////////////////////////////////////////////////////////////////////////////// proc koeff (poly f, int I, int J) "USAGE: koeff(f,I,J); f bivariate polynomial, I,J integers RETURN: if f = sum(a(i,j)*x^i*y^j), then koeff(f,I,J)= a(I,J) (of type number) NOTE: J must be in the range of the exponents of the 2nd ring variable EXAMPLE: example koeff; shows an example " { matrix mat = coeffs(coeffs(f,var(2))[J+1,1],var(1)); if (size(mat) <= I) { return(0);} else { return(leadcoef(mat[I+1,1]));} } example { "EXAMPLE:"; echo = 2; ring r=0,(x,y),dp; koeff(x2+2xy+3xy2-x2y-2y3,1,2); } /////////////////////////////////////////////////////////////////////////////// proc squarefree (poly f) "USAGE: squarefree(f); f poly ASSUME: f is a bivariate polynomial (in the first 2 ring variables). RETURN: poly, a squarefree divisor of f. NOTE: Usually, the return value is the greatest squarefree divisor, but there is one exception: factors with a p-th root, p the characteristic of the basering, are lost. SEE ALSO: allsquarefree EXAMPLE: example squarefree; shows some examples. " { //----------------- Wechsel in geeigneten Ring & Variablendefinition --------- if (nvars(basering)!=2) { ERROR("basering must have exactly 2 variables for Hnoether::squarefree"); } def altring = basering; int e; int gcd_ok=1; string mipl="0"; if (size(parstr(altring))==1) { mipl=string(minpoly); } //---- test: char = (p^k,a) (-> gcd not implemented) or (p,a) (gcd works) ---- if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) { string tststr=charstr(basering); tststr=tststr[1..find(tststr,",")-1]; //-> "p^k" bzw. "p" gcd_ok=(tststr==string(char(basering))); } execute("ring rsqrf = ("+charstr(altring)+"),(x,y),dp;"); if ((gcd_ok!=0) && (mipl!="0")) { execute("minpoly="+mipl+";"); } poly f=fetch(altring,f); poly dif,g,l; if ((char(basering)==0) and (charstr(basering)!=string(char(basering))) and (mipl!="0")) { gcd_ok=0; // since Singular 1.2 gcd no longer implemented } if (gcd_ok!=0) { //--------------------- Berechne f/ggT(f,df/dx,df/dy) ------------------------ dif=diff(f,x); if (dif==0) { g=f; } // zur Beschleunigung else { g=gcd(f,dif); } if (g!=1) { // sonst schon sicher, dass f quadratfrei dif=diff(f,y); if (dif!=0) { g=gcd(g,dif); } } if (g!=1) { e=0; if (g==f) { l=1; } // zur Beschleunigung else { module m=syz(ideal(g,f)); if (deg(m[2,1])>0) { "!! The Singular command 'syz' has returned a wrong result !!"; l=1; // Division f/g muss aufgehen } else { l=m[1,1]; } } } else { e=1; } } else { //------------------- Berechne syz(f,df/dx) oder syz(f,df/dy) ---------------- //-- Achtung: Ist f reduzibel, koennen Faktoren mit Ableitung Null verloren -- //-- gehen! Ist aber nicht weiter schlimm, weil char (p^k,a) nur im irred. -- //-- Fall vorkommen kann. Wenn f nicht g^p ist, wird auf jeden Fall -- //------------------------ ein Faktor gefunden. ------------------------------ dif=diff(f,x); if (dif == 0) { dif=diff(f,y); if (dif==0) { e=2; l=1; } // f is of power divisible by char of basefield else { l=syz(ideal(dif,f))[1,1]; // x^p+y^(p-1) abgedeckt if (subst(f,x,0)==0) { l=l*x; } if (deg(l)==deg(f)) { e=1;} else {e=0;} } } else { l=syz(ideal(dif,f))[1,1]; if (subst(f,y,0)==0) { l=l*y; } if (deg(l)==deg(f)) { e=1;} else {e=0;} } } //--------------- Wechsel in alten Ring und Rueckgabe des Ergebnisses -------- setring altring; if (e==1) { return(f); } // zur Beschleunigung else { poly l=fetch(rsqrf,l); return(l); } } example { "EXAMPLE:"; echo = 2; ring exring=3,(x,y),dp; squarefree((x3+y)^2); squarefree((x+y)^3*(x-y)^2); // Warning: (x+y)^3 is lost squarefree((x+y)^4*(x-y)^2); // result is (x+y)*(x-y) } /////////////////////////////////////////////////////////////////////////////// proc allsquarefree (poly f, poly l) "USAGE : allsquarefree(f,g); f,g poly ASSUME: g is the output of @code{squarefree(f)}. RETURN: the greatest squarefree divisor of f. NOTE : This proc uses factorize to get the missing factors of f not in g and, therefore, may be slow. SEE ALSO: squarefree EXAMPLE: example allsquarefree; shows an example " { //------------------------ Wechsel in geeigneten Ring ------------------------ def altring = basering; string mipl="0"; if (size(parstr(altring))==1) { mipl=string(minpoly); } if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) { string tststr=charstr(basering); tststr=tststr[1..find(tststr,",")-1]; //-> "p^k" bzw. "p" if (tststr!=string(char(basering))) { " Sorry -- not implemented for this ring (gcd doesn't work)"; return(l); } } execute("ring rsqrf = ("+charstr(altring)+"),(x,y),dp;"); if (mipl!="0") { execute("minpoly="+mipl+";"); } poly f=fetch(altring,f); poly l=fetch(altring,l); //---------- eliminiere bereits mit squarefree gefundene Faktoren ------------ poly g=l; while (deg(g)!=0) { f=syz(ideal(g,f))[1,1]; // f=f/g; g=gcd(f,l); } // jetzt f=h^p //--------------- Berechne uebrige Faktoren mit factorize -------------------- if (deg(f)>0) { g=1; //*CL old: ideal factf=factorize(f,1); //* for (int i=1; i<=size(factf); i++) { g=g*factf[i]; } ideal factf=factorize(f)[1]; for (int i=2; i<=size(factf); i++) { g=g*factf[i]; } poly testp=squarefree(g); if (deg(testp) 16001 or exponent divisible by 32003, if there is one 0 else (in this case computing a squarefree divisor in characteristic 32003 could make sense) NOTE: this procedure is only useful in characteristic zero, because otherwise there is no appropriate ordering of the leading coefficients " { poly verbrecher=0; intvec leitexp; for (; (f<>0) and (verbrecher==0); f=f-lead(f)) { if ((leadcoef(f)<-16001) or (leadcoef(f)>16001)) {verbrecher=lead(f);} leitexp=leadexp(f); if (( ((leitexp[1] % 32003) == 0) and (leitexp[1]<>0)) or ( ((leitexp[2] % 32003) == 0) and (leitexp[2]<>0)) ) {verbrecher=lead(f);} } return(verbrecher); } ////////////////////////////////////////////////////////////////////////////// proc develop "USAGE: develop(f [,n]); f poly, n int ASSUME: f is a bivariate polynomial (in the first 2 ring variables) and irreducible as power series (for reducible f use @code{hnexpansion}). RETURN: list @code{L} with: @texinfo @table @asis @item @code{L[1]}; matrix: Each row contains the coefficients of the corresponding line of the Hamburger-Noether expansion (HNE). The end of the line is marked in the matrix by the first ring variable (usually x). @item @code{L[2]}; intvec: indicating the length of lines of the HNE @item @code{L[3]}; int: 0 if the 1st ring variable was transversal (with respect to f), @* 1 if the variables were changed at the beginning of the computation, @* -1 if an error has occurred. @item @code{L[4]}; poly: the transformed polynomial of f to make it possible to extend the Hamburger-Noether development a posteriori without having to do all the previous calculation once again (0 if not needed) @item @code{L[5]}; int: 1 if the curve has exactly one branch (i.e., is irreducible), @* 0 else (i.e., the curve has more than one HNE, or f is not valid). @end table @end texinfo DISPLAY: The (non zero) elements of the HNE (if not called by another proc). NOTE: The optional parameter @code{n} affects only the computation of the LAST line of the HNE. If it is given, the HN-matrix @code{L[1]} will have at least @code{n} columns. @* Otherwise, the number of columns will be chosen minimal such that the matrix contains all necessary information (i.e., all lines of the HNE but the last (which is in general infinite) have place). @* If @code{n} is negative, the algorithm is stopped as soon as the computed information is sufficient for @code{invariants(L)}, but the HN-matrix @code{L[1]} may still contain undetermined elements, which are marked with the 2nd variable (of the basering). @* For time critical computations it is recommended to use @code{ring ...,(x,y),ls} as basering - it increases the algorithm's speed. @* If @code{printlevel>=0} comments are displayed (default is @code{printlevel=0}). SEE ALSO: hnexpansion, extdevelop, displayHNE EXAMPLES: example develop; shows an example example parametrize; shows an example for using the 2nd parameter " { //--------- Abfangen unzulaessiger Ringe: 1) nur eine Unbestimmte ------------ poly f=#[1]; if (size(#) > 1) {int maxspalte=#[2];} else {int maxspalte= 1 ; } if (nvars(basering) < 2) { " Sorry. I need two variables in the ring."; return(list(matrix(maxideal(1)[1]),intvec(0),-1,poly(0),0));} if (nvars(basering) > 2) { dbprint(printlevel-voice+2, " Warning! You have defined too many variables! All variables except the first two will be ignored!" ); } string namex=varstr(1); string namey=varstr(2); list return_error=matrix(maxideal(1)[2]),intvec(0),int(-1),poly(0),int(0); //------------- 2) mehrere Unbestimmte, weitere unzulaessige Ringe ----------- // Wir koennen einheitlichen Rueckgabewert benutzen, aus dem ersichtlich ist, // dass ein Fehler aufgetreten ist: return_error. //---------------------------------------------------------------------------- if (charstr(basering)=="real") { " The algorithm doesn't work with 'real' as coefficient field."; // denn : map from characteristic -1 to -1 not implemented return(return_error); } if ((char(basering)!=0) and (charstr(basering)!=string(char(basering)))) { //-- teste, ob char = (p^k,a) (-> a primitiv; erlaubt) oder (p,a[,b,...]) ---- string tststr=charstr(basering); tststr=tststr[1..find(tststr,",")-1]; //-> "p^k" bzw. "p" int primit=(tststr==string(char(basering))); if (primit!=0) { " Such extensions of Z/p are not implemented."; " Please try (p^k,a) as ground field or use `hnexpansion'."; return(return_error); } } //---- Ende der unzulaessigen Ringe; Ringwechsel in einen guenstigen Ring: --- int ringwechsel=(varstr(basering)!="x,y") or (ordstr(basering)!="ls(2),C"); def altring = basering; if (ringwechsel) { string mipl=string(minpoly); execute("ring guenstig = ("+charstr(altring)+"),(x,y),ls;"); if ((char(basering)==0) && (mipl!="0")) { execute("minpoly="+mipl+";"); }} else { def guenstig=basering; } export guenstig; //-------------------------- Initialisierungen ------------------------------- map m=altring,x,y; if (ringwechsel) { poly f=m(f); } if (defined(HNDebugOn)) {"received polynomial: ",f,", where x =",namex,", y =",namey;} kill m; int M,N,Q,R,l,e,hilf,eps,getauscht,Abbruch,zeile,exponent,Ausgabe; // Werte von Ausgabe: 0 : normale HNE-Matrix, // 1 : Fehler aufgetreten - Matrix (namey) zurueck // 2 : Die HNE ist eine Nullzeile - Matrix (0) zurueck // int maxspalte=1; geaendert: wird jetzt am Anfang gesetzt int minimalHNE=0; // Flag fuer minimale HNE-Berechnung int einzweig=1; // Flag fuer Irreduzibilit"at intvec hqs; // erhaelt die Werte von h(zeile)=Q; if (maxspalte<0) { minimalHNE=1; maxspalte=1; } number c,delt; int p = char(basering); string ringchar=charstr(basering); map xytausch = basering,y,x; if ((p!=0) and (ringchar != string(p))) { // coefficient field is extension of Z/pZ execute("int n_elements="+ ringchar[1,size(ringchar)-size(parstr(basering))-1]+";"); // number of elements of actual ring number generat=par(1); // generator of the coefficient field of the ring } //========= Abfangen von unzulaessigen oder trivialen Eingaben =============== //------------ Nullpolynom oder Einheit im Potenzreihenring: ----------------- if (f == 0) { dbprint(printlevel+1,"The given polynomial is the zero-polynomial !"); Abbruch=1; Ausgabe=1; } else { intvec nm = getnm(f); N = nm[1]; M = nm[2]; // Berechne Schnittpunkte Newtonpolygon mit Achsen if (N == 0) { dbprint(printlevel+1,"The given polynomial is a unit as power series !"); Abbruch=1; Ausgabe=1; } else { if (N == -1) { if ((voice==2) && (printlevel > -1)) { "The HNE is x = 0"; } Abbruch=1; Ausgabe=2; getauscht=1; if (M <> 1) { einzweig=0; } } else { if (M == -1) { if ((voice==2) && (printlevel > -1)) { "The HNE is y = 0"; } Abbruch=1; Ausgabe=2; if (N <> 1) { einzweig=0; } }}} } //--------------------- Test auf Quadratfreiheit ----------------------------- if (Abbruch==0) { //-------- Fall basering==0,... : Wechsel in Ring mit char >0 ---------------- // weil squarefree eine Standardbasis berechnen muss (verwendet Syzygien) // -- wenn f in diesem Ring quadratfrei ist, dann erst recht im Ring guenstig //---------------------------------------------------------------------------- if ((p==0) and (size(charstr(basering))==1)) { int testerg=(polytest(f)==0); ring zweitring = 32003,(x,y),dp; map polyhinueber=guenstig,x,y; // fetch geht nicht poly f=polyhinueber(f); poly test_sqr=squarefree(f); if (test_sqr != f) { if (printlevel>0) { "Most probably the given polynomial is not squarefree. But the test was"; "made in characteristic 32003 and not 0 to improve speed. You can"; "(r) redo the test in char 0 (but this may take some time)"; "(c) continue the development, if you're sure that the polynomial", "IS squarefree"; if (testerg==1) { "(s) continue the development with a squarefree factor (*)";} "(q) or just quit the algorithm (default action)"; "";"Please enter the letter of your choice:"; string str=read("")[1]; } else { string str="r"; } // printlevel <= 0: non-interactive behaviour setring guenstig; map polyhinueber=zweitring,x,y; if (str=="r") { poly test_sqr=squarefree(f); if (test_sqr != f) { if (printlevel>0) { "The given polynomial is in fact not squarefree."; } else { "The given polynomial is not squarefree!"; } "I'll continue with the radical."; if (printlevel>0) { pause("Hit RETURN to continue:"); } f=test_sqr; } else { dbprint(printlevel, "everything is ok -- the polynomial is squarefree in char(k)=0"); } } else { if ((str=="s") and (testerg==1)) { "(*) attention: it could be that the factor is only one in char 32003!"; f=polyhinueber(test_sqr); } else { if (str<>"c") { setring altring;kill guenstig;kill zweitring; return(return_error);} else { "if the algorithm doesn't terminate, you were wrong...";} }} kill zweitring; nm = getnm(f); // N,M haben sich evtl. veraendert N = nm[1]; M = nm[2]; // Berechne Schnittpunkte Newtonpolynom mit Achsen if (defined(HNDebugOn)) {"I continue with the polynomial",f; } } else { setring guenstig; kill zweitring; } } // ------------------- Fall Charakteristik > 0 ------------------------------- else { poly test_sqr=squarefree(f); if (test_sqr == 1) { "The given polynomial is of the form g^"+string(p)+", therefore", "reducible.";"Please try again."; setring altring; kill guenstig; return(return_error);} if (test_sqr != f) { "The given polynomial is not squarefree. I'll continue with the radical."; if (p != 0) {"But if the polynomial contains a factor of the form g^"+string(p)+","; "this factor will be lost.";} if (printlevel>0) { pause("Hit RETURN to continue:"); } f=test_sqr; nm = getnm(f); // N,M haben sich veraendert N = nm[1]; M = nm[2]; // Berechne Schnittpunkte Newtonpolynom mit Achsen if (defined(HNDebugOn)) {"I continue with the polynomial",f; } } } // endelse(p==0) if (N==0) { " Sorry. The remaining polynomial is a unit in the power series ring..."; setring altring;kill guenstig;return(return_error); } //---------------------- gewaehrleiste, dass x transvers ist ----------------- if (M < N) { f = xytausch(f); // Variablentausch : x jetzt transvers getauscht = 1; // den Tausch merken M = M+N; N = M-N; M = M-N; // M, N auch vertauschen } if (defined(HNDebugOn)) { if (getauscht) {"x<->y were exchanged; polynomial is now ",f;} else {"x , y were not exchanged";} "M resp. N are now",M,N; } } // end(if Abbruch==0) ideal a(0); while (Abbruch==0) { //================= Beginn der Schleife (eigentliche Entwicklung) ============ //------------------- ist das Newtonpolygon eine gerade Linie? --------------- if (testreducible(f,N,M)) { dbprint(printlevel+1," The given polynomial is not irreducible"); kill guenstig; setring altring; return(return_error); // Abbruch der Prozedur! } R = M%N; Q = M / N; //-------------------- Fall Rest der Division R = 0 : ------------------------ if (R == 0) { c = koeff(f,0,N); if (c == 0) {"Something has gone wrong! I didn't get N correctly!"; exit;} e = gcd(M,N); //----------------- Test, ob leitf = c*(y^N - delta*x^(m/e))^e ist ----------- if (p==0) { delt = koeff(f,M/ e,N - N/ e) / (-1*e*c); if (defined(HNDebugOn)) {"quasihomogeneous leading form:", leit(f,N,M)," = ",c,"* (y -",delt,"* x^"+string(M/ e)+")^",e," ?";} if (leit(f,N,M) != c*(y^(N/ e) - delt*x^(M/ e))^e) { dbprint(printlevel+1," The given polynomial is reducible !"); Abbruch=1; Ausgabe=1; } } else { // p!=0 if (e%p != 0) { delt = koeff(f,M/ e,N - N/ e) / (-1*e*c); if (defined(HNDebugOn)) {"quasihomogeneous leading form:", leit(f,N,M)," = ",c,"* (y -",delt,"* x^"+string(M/ e)+")^",e," ?";} if (leit(f,N,M) != c*(y^(N/ e) - delt*x^(M/ e))^e) { dbprint(printlevel+1," The given polynomial is reducible !"); Abbruch=1; Ausgabe=1; } } else { // e%p == 0 eps = e; for (l = 0; eps%p == 0; l=l+1) { eps=eps/ p;} if (defined(HNDebugOn)) {e," -> ",eps,"*",p,"^",l;} delt = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c); if ((ringchar != string(p)) and (delt != 0)) { //- coeff. field is not Z/pZ => we`ve to correct delta by taking (p^l)th root- if (delt == generat) {exponent=1;} else { if (delt == 1) {exponent=0;} else { exponent=pardeg(delt); //-- an dieser Stelle kann ein Fehler auftreten, wenn wir eine transzendente - //-- Erweiterung von Z/pZ haben: dann ist das hinzuadjungierte Element kein - //-- Erzeuger der mult. Gruppe, d.h. in Z/pZ (a) gibt es i.allg. keinen - //-- Exponenten mit z.B. a2+a = a^exp - //---------------------------------------------------------------------------- }} delt = generat^(extgcd(n_elements-1,p^l)[3]*exponent); } if (defined(HNDebugOn)) {"quasihomogeneous leading form:", leit(f,N,M)," = ",c,"* (y^"+string(N/ e),"-",delt,"* x^" +string(M/ e)+")^",e," ?";} if (leit(f,N,M) != c*(y^(N/ e) - delt*x^(M/ e))^e) { dbprint(printlevel+1," The given polynomial is reducible !"); Abbruch=1; Ausgabe=1; } } } if (Abbruch == 0) { f = T1_Transform(f,delt,M/ e); dbprint(printlevel-voice+2,"a("+string(zeile)+","+string(Q)+") = " +string(delt)); a(zeile)[Q]=delt; if (defined(HNDebugOn)) {"transformed polynomial: ",f;}} nm=getnm(f); N=nm[1]; M=nm[2]; // Neuberechnung des Newtonpolygons } //--------------------------- Fall R > 0 : ----------------------------------- else { dbprint(printlevel-voice+2, "h("+string(zeile)+ ") ="+string(Q)); hqs[zeile+1]=Q; // denn zeile beginnt mit dem Wert 0 a(zeile)[Q+1]=x; // Markierung des Zeilenendes der HNE maxspalte=maxspalte*((Q+1) < maxspalte) + (Q+1)*((Q+1) >= maxspalte); // Anpassung der Sp.zahl der HNE-Matrix f = T_Transform(f,Q,N); if (defined(HNDebugOn)) {"transformed polynomial: ",f;} zeile=zeile+1; //------------ Bereitstellung von Speicherplatz fuer eine neue Zeile: -------- ideal a(zeile); M=N;N=R; } //--------------- schneidet das Newtonpolygon beide Achsen? ------------------ if (M==-1) { dbprint(printlevel-voice+2,"The HNE is finite!"); a(zeile)[Q+1]=x; // Markiere das Ende der Zeile hqs[zeile+1]=Q; maxspalte=maxspalte*((Q+1) < maxspalte) + (Q+1)*((Q+1) >= maxspalte); if (N <> 1) { einzweig=0; } f=0; // transformiertes Polynom wird nicht mehr gebraucht Abbruch=1; } else {if (M maxspalte) or (minimalHNE==1)) and (size(a(zeile))>0)) //---------------------------------------------------------------------------- // Abbruch, wenn die Matrix so voll ist, dass eine neue Spalte angefangen // werden muesste und die letzte Zeile nicht nur Nullen enthaelt // oder wenn die Matrix nicht voll gemacht werden soll (minimale Information) //---------------------------------------------------------------------------- { Abbruch=1; hqs[zeile+1]=-1; if (maxspalte < ncols(a(zeile))) { maxspalte=ncols(a(zeile));} if ((minimalHNE==1) and (M <= maxspalte)) { // teile param mit, dass Eintraege der letzten Zeile nur teilw. richtig sind:- hqs[zeile+1]=-M; //------------- markiere den Rest der Zeile als unbekannt: ------------------- for (R=M; R <= maxspalte; R++) { a(zeile)[R]=y;} } // R wird nicht mehr gebraucht } //========================= Ende der Schleife ================================ } setring altring; if (Ausgabe == 0) { //-------------------- Ergebnis in den alten Ring transferieren: ------------- map zurueck=guenstig,maxideal(1)[1],maxideal(1)[2]; matrix amat[zeile+1][maxspalte]; ideal uebergabe; for (e=0; e<=zeile; e=e+1) { uebergabe=zurueck(a(e)); if (ncols(uebergabe) > 1) { amat[e+1,1..ncols(uebergabe)]=uebergabe;} else {amat[e+1,1]=uebergabe[1];} } if (ringwechsel) { if (nvars(altring)==2) { f=fetch(guenstig,f); } else { f=zurueck(f); } } } kill guenstig; if ((einzweig==0) && (voice==2) && (printlevel > -1)) { "// Note: The curve is reducible, but we were able to compute a HNE."; "// This means the result is only one of several existing HNE's."; } if (Ausgabe == 0) { return(list(amat,hqs,getauscht,f,einzweig));} if (Ausgabe == 1) { return(return_error);} // error has occurred if (Ausgabe == 2) { return(list(matrix(ideal(0,x)),intvec(1),getauscht, poly(0),einzweig));} // HNE is x=0 or y=0 } example { "EXAMPLE:"; echo = 2; ring exring = 7,(x,y),ds; list Hne=develop(4x98+2x49y7+x11y14+2y14); print(Hne[1]); // therefore the HNE is: // z(-1)= 3*z(0)^7 + z(0)^7*z(1), // z(0) = z(1)*z(2), (there is 1 zero in the 2nd row before x) // z(1) = z(2)^3*z(3), (there are 3 zeroes in the 3rd row) // z(2) = z(3)*z(4), // z(3) = -z(4)^2 + 0*z(4)^3 +...+ 0*z(4)^8 + ?*z(4)^9 + ... // (the missing x in the last line indicates that it is not complete.) Hne[2]; param(Hne); // parametrization: x(t)= -t^14+O(t^21), y(t)= -3t^98+O(t^105) // (the term -t^109 in y may have a wrong coefficient) displayHNE(Hne); } /////////////////////////////////////////////////////////////////////////////// // procedures to extract information out of HNE // /////////////////////////////////////////////////////////////////////////////// proc param (list L, list #) "USAGE: param(L [,s]); L list, s any type (optional) ASSUME: L is the output of @code{develop(f)}, or of @code{extdevelop(develop(f),n)}, or (one entry in) the list of HN data created by @code{hnexpansion(f[,\"ess\"])}. RETURN: If L are the HN data of an irreducible plane curve singularity f: a parametrization for f in the following format: @* - if only the list L is given, the result is an ideal of two polynomials p[1],p[2]: if the HNE was finite then f(p[1],p[2])=0}; if not, the true parametrization will be given by two power series, and p[1],p[2] are truncations of these series.@* - if the optional parameter s is given, the result is a list l: l[1]=param(L) (ideal) and l[2]=intvec with two entries indicating the highest degree up to which the coefficients of the monomials in l[1] are exact (entry -1 means that the corresponding parametrization is exact). If L collects the HN data of a reducible plane curve singularity f, the return value is a list of parametrizations in the respective format. NOTE: If the basering has only 2 variables, the first variable is chosen as indefinite. Otherwise, the 3rd variable is chosen. SEE ALSO: develop, extdevelop KEYWORDS: parametrization EXAMPLE: example param; shows an example example develop; shows another example " { //-------------------------- Initialisierungen ------------------------------- int return_list; if (size(#)>0) { return_list=1; } if (typeof(L[1])=="list") { // output of hnexpansion (> 1 branch) list Ergebnis; for (int i=1; i<=size(L); i++) { dbprint(printlevel-voice+4,"// Parametrization of branch number " +string(i)+" computed."); printlevel=printlevel+1; if (return_list==1) { Ergebnis[i]=param(L[i],1); } else { Ergebnis[i]=param(L[i]); } printlevel=printlevel-1; } return(Ergebnis); } else { matrix m=L[1]; intvec v=L[2]; int switch=L[3]; } if (switch==-1) { "An error has occurred in develop, so there is no HNE."; return(ideal(0,0)); } int fehler,fehlervor,untergrad,untervor,beginn,i,zeile,hilf; if (nvars(basering) > 2) { poly z(size(v)+1)=var(3); } else { poly z(size(v)+1)=var(1); } poly z(size(v)); zeile=size(v); //------------- Parametrisierung der untersten Zeile der HNE ----------------- if (v[zeile] > 0) { fehler=0; // die Parametrisierung wird exakt werden for (i=1; i<=v[zeile]; i++) { z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i; } } else { untervor=1; // = Untergrad der vorhergehenden Zeile if (v[zeile]==-1) { fehler=ncols(m)+1; for (i=1; i<=ncols(m); i++) { z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i; if ((untergrad==0) and (m[zeile,i]!=0)) {untergrad=i;} // = Untergrad der aktuellen Zeile } } else { fehler= -v[zeile]; for (i=1; i<-v[zeile]; i++) { z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i; if ((untergrad==0) and (m[zeile,i]!=0)) {untergrad=i;} } } } //------------- Parametrisierung der restlichen Zeilen der HNE --------------- for (zeile=size(v)-1; zeile>0; zeile--) { poly z(zeile); beginn=0; // Beginn der aktuellen Zeile for (i=1; i<=v[zeile]; i++) { z(zeile)=z(zeile)+m[zeile,i]*z(zeile+1)^i; if ((beginn==0) and (m[zeile,i]!=0)) { beginn=i;} } z(zeile)=z(zeile) + z(zeile+1)^v[zeile] * z(zeile+2); if (beginn==0) { if (fehler>0) { // damit fehler=0 bleibt bei exakter Param. fehlervor=fehler; // Fehler der letzten Zeile fehler=fehler+untergrad*(v[zeile]-1)+untervor; // Fehler dieser Zeile hilf=untergrad; untergrad=untergrad*v[zeile]+untervor; untervor=hilf;} // untervor = altes untergrad } else { fehlervor=fehler; fehler=fehler+untergrad*(beginn-1); untervor=untergrad; untergrad=untergrad*beginn; } } //--------------------- Ausgabe der Fehlerabschaetzung ----------------------- if (switch==0) { if (fehler>0) { if (fehlervor>0) { dbprint(printlevel-voice+4,""+ "// ** Warning: result is exact up to order "+string(fehlervor-1)+ " in "+ string(var(1))+" and "+string(fehler-1)+" in " + string(var(2))+" !"); } else { dbprint(printlevel-voice+4,""+ "// ** Warning: result is exact up to order "+ string(fehler-1)+ " in "+string(var(2))+" !"); } } if (return_list==0) { return(ideal(z(2),z(1))); } else { return(list(ideal(z(2),z(1)),intvec(fehlervor-1,fehler-1))); } } else { if (fehler>0) { if (fehlervor>0) { dbprint(printlevel-voice+4,""+ "// ** Warning: result is exact up to order "+string(fehler-1)+ " in "+ string(var(1))+" and "+string(fehlervor-1)+" in " + string(var(2))+" !"); } else { dbprint(printlevel-voice+4,""+ "// ** Warning: result is exact up to order "+ string(fehler-1)+ " in "+string(var(1))+" !"); } } if (return_list==0) { return(ideal(z(1),z(2))); } else { return(list(ideal(z(1),z(2)),intvec(fehler-1,fehlervor-1))); } } } example { "EXAMPLE:"; echo = 2; ring exring=0,(x,y,t),ds; poly f=x3+2xy2+y2; list Hne=develop(f); list hne_extended=extdevelop(Hne,10); // compare the HNE matrices ... print(Hne[1]); print(hne_extended[1]); // ... and the resulting parametrizations: param(Hne); param(hne_extended); param(hne_extended,0); // An example with more than one branch: list L=hnexpansion(f*(x2+y4)); def HNring = L[1]; setring HNring; param(hne); } /////////////////////////////////////////////////////////////////////////////// proc invariants "USAGE: invariants(INPUT); INPUT list or poly ASSUME: @code{INPUT} is the output of @code{develop(f)}, or of @code{extdevelop(develop(f),n)}, or one entry of the list of HN data computed by @code{hnexpansion(f[,\"ess\"])}. RETURN: list @code{INV} of the following format: @format INV[1]: intvec (characteristic exponents) INV[2]: intvec (generators of the semigroup) INV[3]: intvec (Puiseux pairs, 1st components) INV[4]: intvec (Puiseux pairs, 2nd components) INV[5]: int (degree of the conductor) INV[6]: intvec (sequence of multiplicities) @end format If @code{INPUT} contains no valid HN expansion, the empty list is returned. ASSUME: @code{INPUT} is a bivariate polynomial f, or the output of @code{hnexpansion(f)}, or the list of HN data computed by @code{hnexpansion(f [,\"ess\"])}. RETURN: list @code{INV}, such that @code{INV[i]} coincides with the output of @code{invariants(develop(f[i]))}, where f[i] is the i-th branch of f, and the last entry of @code{INV} contains further invariants of f in the format: @format INV[last][1] : intmat (contact matrix of the branches) INV[last][2] : intmat (intersection multiplicities of the branches) INV[last][3] : int (delta invariant of f) @end format NOTE: In case the Hamburger-Noether expansion of the curve f is needed for other purposes as well it is better to calculate this first with the aid of @code{hnexpansion} and use it as input instead of the polynomial itself. SEE ALSO: hnexpansion, develop, displayInvariants, multsequence, intersection KEYWORDS: characteristic exponents; semigroup of values; Puiseux pairs; conductor, degree; multiplicities, sequence of EXAMPLE: example invariants; shows an example " { //---- INPUT = poly, or HNEring, or hne of reducible curve ----------------- if (typeof(#[1])!="matrix") { if (typeof(#[1])=="poly") { list L=hnexpansion(#[1]); if (typeof(L[1])=="ring") { def altring = basering; def HNring = L[1]; setring HNring; list Ergebnis = invariants(hne); setring altring; kill HNring; return(Ergebnis); } else { return(invariants(L)); } } if (typeof(#[1])=="ring") { def altring = basering; def HNring = #[1]; setring HNring; list Ergebnis = invariants(hne); setring altring; kill HNring; return(Ergebnis); } if (typeof(#[1])=="list") { list hne=#; list Ergebnis; for (int lauf=1;lauf<=size(hne);lauf++) { Ergebnis[lauf]=invariants(hne[lauf]); } // Calculate the intersection matrix and the intersection multiplicities. intmat contact[size(hne)][size(hne)]; intmat intersectionmatrix[size(hne)][size(hne)]; int Lauf; for (lauf=1;lauf<=size(hne);lauf++) { for (Lauf=lauf+1;Lauf<=size(hne);Lauf++) { contact[lauf,Lauf]=separateHNE(hne[lauf],hne[Lauf]); contact[Lauf,lauf]=contact[lauf,Lauf]; intersectionmatrix[lauf,Lauf]=intersection(hne[lauf],hne[Lauf]); intersectionmatrix[Lauf,lauf]=intersectionmatrix[lauf,Lauf]; } } // Calculate the delta invariant. int inters; int del=Ergebnis[size(hne)][5]/2; for(lauf=1;lauf<=size(hne)-1;lauf++) { del=del+Ergebnis[lauf][5]/2; for(Lauf=lauf+1;Lauf<=size(hne);Lauf++) { inters=inters+intersectionmatrix[lauf,Lauf]; } } del=del+inters; list LAST=contact,intersectionmatrix,del; Ergebnis[size(hne)+1]=LAST; return(Ergebnis); } } //-------------------------- Initialisierungen ------------------------------- matrix m=#[1]; intvec v=#[2]; int switch=#[3]; list ergebnis; if (switch==-1) { "An error has occurred in develop, so there is no HNE."; return(ergebnis); } intvec beta,s,svorl,ordnung,multseq,mpuiseux,npuiseux,halbgr; int genus,zeile,i,j,k,summe,conductor,ggT; string Ausgabe; int nc=ncols(m); int nr=nrows(m); ordnung[nr]=1; // alle Indizes muessen (gegenueber [Ca]) um 1 erhoeht werden, // weil 0..r nicht als Wertebereich erlaubt ist (aber nrows(m)==r+1) //---------------- Bestimme den Untergrad der einzelnen Zeilen --------------- for (zeile=nr; zeile>1; zeile--) { if ((size(ideal(m[zeile,1..nc])) > 1) or (zeile==nr)) { // keine Nullzeile k=1; while (m[zeile,k]==0) {k++;} ordnung[zeile-1]=k*ordnung[zeile]; // vgl. auch Def. von untergrad in genus++; // proc param svorl[genus]=zeile;} // werden gerade in umgekehrter Reihenfolge abgelegt else { ordnung[zeile-1]=v[zeile]*ordnung[zeile]+ordnung[zeile+1]; }} //----------------- charakteristische Exponenten (beta) ---------------------- s[1]=1; for (k=1; k <= genus; k++) { s[k+1]=svorl[genus-k+1];} // s[2]==s(1), u.s.w. beta[1]=ordnung[1]; //charakt. Exponenten: Index wieder verschoben for (k=1; k <= genus; k++) { summe=0; for (i=1; i <= s[k]; i++) {summe=summe+v[i]*ordnung[i];} beta[k+1]=summe+ordnung[s[k]]+ordnung[s[k]+1]-ordnung[1]; } //--------------------------- Puiseuxpaare ----------------------------------- int produkt=1; for (i=1; i<=genus; i++) { ggT=gcd(beta[1],beta[i+1]*produkt); mpuiseux[i]=beta[i+1]*produkt / ggT; npuiseux[i]=beta[1] / ggT; produkt=produkt*npuiseux[i]; } //---------------------- Grad des Konduktors --------------------------------- summe=1-ordnung[1]; if (genus > 0) { for (i=2; i <= genus+1; i++) { summe=summe + beta[i] * (ordnung[s[i-1]] - ordnung[s[i]]); } // n.b.: Indizierung wieder um 1 verschoben } conductor=summe; //------------------- Erzeuger der Halbgruppe: ------------------------------- halbgr=puiseux2generators(mpuiseux,npuiseux); //------------------- Multiplizitaetensequenz: ------------------------------- k=1; for (i=1; i1)) { tester=tester-1; } if ((multseq[tester]!=1) and (multseq[tester]!=k-tester)) { for (i=k+1; i<=tester+multseq[tester]; i++) { multseq[i]=1; } } //--- Ende T.Keilen --- 06.05.02 //------------------------- Rueckgabe ---------------------------------------- ergebnis=beta,halbgr,mpuiseux,npuiseux,conductor,multseq; return(ergebnis); } example { "EXAMPLE:"; echo = 2; ring exring=0,(x,y),dp; list Hne=develop(y4+2x3y2+x6+x5y); list INV=invariants(Hne); INV[1]; // the characteristic exponents INV[2]; // the generators of the semigroup of values INV[3],INV[4]; // the Puiseux pairs in packed form INV[5] / 2; // the delta-invariant INV[6]; // the sequence of multiplicities // To display the invariants more 'nicely': displayInvariants(Hne); ///////////////////////////// INV=invariants((x2-y3)*(x3-y5)); INV[1][1]; // the characteristic exponents of the first branch INV[2][6]; // the sequence of multiplicities of the second branch print(INV[size(INV)][1]); // the contact matrix of the branches print(INV[size(INV)][2]); // the intersection numbers of the branches INV[size(INV)][3]; // the delta invariant of the curve } /////////////////////////////////////////////////////////////////////////////// proc displayInvariants "USAGE: displayInvariants(INPUT); INPUT list or poly ASSUME: @code{INPUT} is a bivariate polynomial, or the output of @code{develop(f)}, resp. of @code{extdevelop(develop(f),n)}, or (one entry of) the list of HN data computed by @code{hnexpansion(f[,\"ess\"])}. RETURN: none DISPLAY: invariants of the corresponding branch, resp. of all branches, in a better readable form. NOTE: If the Hamburger-Noether expansion of the curve f is needed for other purposes as well it is better to calculate this first with the aid of @code{hnexpansion} and use it as input instead of the polynomial itself. SEE ALSO: invariants, intersection, develop, hnexpansion EXAMPLE: example displayInvariants; shows an example " { // INPUT = polynomial or ring if (typeof(#[1])=="poly") { list L=hnexpansion(#[1]); if (typeof(L[1])=="ring") { def HNring = L[1]; setring HNring; displayInvariants(hne); return(); } else { displayInvariants(L); return(); } } if (typeof(#[1])=="ring") { def HNring = #[1]; setring HNring; displayInvariants(hne); return(); } // INPUT = hne of a plane curve int i,j,k,mul; string Ausgabe; list ergebnis; //-- entferne ueberfluessige Daten zur Erhoehung der Rechengeschwindigkeit: -- #=stripHNE(#); //-------------------- Ausgabe eines Zweiges --------------------------------- if (typeof(#[1])=="matrix") { ergebnis=invariants(#); if (size(ergebnis)!=0) { " characteristic exponents :",ergebnis[1]; " generators of semigroup :",ergebnis[2]; if (size(ergebnis[1])>1) { for (i=1; i<=size(ergebnis[3]); i++) { Ausgabe=Ausgabe+"("+string(ergebnis[3][i])+"," +string(ergebnis[4][i])+")"; }} " Puiseux pairs :",Ausgabe; " degree of the conductor :",ergebnis[5]; " delta invariant :",ergebnis[5]/2; " sequence of multiplicities:",ergebnis[6]; }} //-------------------- Ausgabe aller Zweige ---------------------------------- else { ergebnis=invariants(#); intmat contact=ergebnis[size(#)+1][1]; intmat intersectionmatrix=ergebnis[size(#)+1][2]; for (j=1; j<=size(#); j++) { " --- invariants of branch number",j,": ---"; " characteristic exponents :",ergebnis[j][1]; " generators of semigroup :",ergebnis[j][2]; Ausgabe=""; if (size(ergebnis[j][1])>1) { for (i=1; i<=size(ergebnis[j][3]); i++) { Ausgabe=Ausgabe+"("+string(ergebnis[j][3][i])+"," +string(ergebnis[j][4][i])+")"; }} " Puiseux pairs :",Ausgabe; " degree of the conductor :",ergebnis[j][5]; " delta invariant :",ergebnis[j][5]/2; " sequence of multiplicities:",ergebnis[j][6]; ""; } if (size(#)>1) { " -------------- contact numbers : -------------- ";""; Ausgabe="branch | "; for (j=size(#); j>1; j--) { if (size(string(j))==1) { Ausgabe=Ausgabe+" "+string(j)+" "; } else { Ausgabe=Ausgabe+string(j)+" "; } } Ausgabe; Ausgabe="-------+"; for (j=2; jj; k--) { mul=contact[j,k];//separateHNE(#[j],#[k]); for (i=1; i<=5-size(string(mul)); i++) { Ausgabe=Ausgabe+" "; } Ausgabe=Ausgabe+string(mul); if (k>j+1) { Ausgabe=Ausgabe+","; } } Ausgabe; } ""; if (size(#)>1) { " -------------- intersection multiplicities : -------------- ";""; Ausgabe="branch | "; for (j=size(#); j>1; j--) { if (size(string(j))==1) { Ausgabe=Ausgabe+" "+string(j)+" "; } else { Ausgabe=Ausgabe+string(j)+" "; } } Ausgabe; Ausgabe="-------+"; for (j=2; jj; k--) { mul=intersectionmatrix[j,k];//intersection(#[j],#[k]); for (i=1; i<=5-size(string(mul)); i++) { Ausgabe=Ausgabe+" "; } Ausgabe=Ausgabe+string(mul); if (k>j+1) { Ausgabe=Ausgabe+","; } } Ausgabe; } ""; " -------------- delta invariant of the curve : ",ergebnis[size(#)+1][3]; } return(); } example { "EXAMPLE:"; echo = 2; ring exring=0,(x,y),dp; list Hne=develop(y4+2x3y2+x6+x5y); displayInvariants(Hne); } /////////////////////////////////////////////////////////////////////////////// proc multiplicities "USAGE: multiplicities(L); L list ASSUME: L is the output of @code{develop(f)}, or of @code{extdevelop(develop(f),n)}, or one entry in the list @code{hne} in the ring created by @code{hnexpansion(f[,\"ess\"])}. RETURN: intvec of the different multiplicities that occur when successively blowing-up the curve singularity corresponding to f. SEE ALSO: multsequence, develop EXAMPLE: example multiplicities; shows an example " { matrix m=#[1]; intvec v=#[2]; int switch=#[3]; list ergebnis; if (switch==-1) { "An error has occurred in develop, so there is no HNE."; return(intvec(0)); } intvec ordnung; int zeile,k; int nc=ncols(m); int nr=nrows(m); ordnung[nr]=1; //---------------- Bestimme den Untergrad der einzelnen Zeilen --------------- for (zeile=nr; zeile>1; zeile--) { if ((size(ideal(m[zeile,1..nc])) > 1) or (zeile==nr)) { // keine Nullzeile k=1; while (m[zeile,k]==0) {k++;} ordnung[zeile-1]=k*ordnung[zeile]; } else { ordnung[zeile-1]=v[zeile]*ordnung[zeile]+ordnung[zeile+1]; }} return(ordnung); } example { "EXAMPLE:"; echo = 2; int p=printlevel; printlevel=-1; ring r=0,(x,y),dp; multiplicities(develop(x5+y7)); // The first value is the multiplicity of the curve itself, here it's 5 printlevel=p; } /////////////////////////////////////////////////////////////////////////////// proc puiseux2generators (intvec m, intvec n) "USAGE: puiseux2generators(m,n); m,n intvec ASSUME: m, resp. n, represent the 1st, resp. 2nd, components of Puiseux pairs (e.g., @code{m=invariants(L)[3]}, @code{n=invariants(L)[4]}). RETURN: intvec of the generators of the semigroup of values. SEE ALSO: invariants EXAMPLE: example puiseux2generators; shows an example " { intvec beta; int q=1; //------------ glatte Kurve (eigentl. waeren m,n leer): ---------------------- if (m==0) { return(intvec(1)); } //------------------- singulaere Kurve: -------------------------------------- for (int i=1; i<=size(n); i++) { q=q*n[i]; } beta[1]=q; // == q_0 m=1,m; n=1,n; // m[1] ist damit m_0 usw., genau wie beta[1]==beta_0 for (i=2; i<=size(n); i++) { beta[i]=m[i]*q/n[i] - m[i-1]*q + n[i-1]*beta[i-1]; q=q/n[i]; // == q_i } return(beta); } example { "EXAMPLE:"; echo = 2; // take (3,2),(7,2),(15,2),(31,2),(63,2),(127,2) as Puiseux pairs: puiseux2generators(intvec(3,7,15,31,63,127),intvec(2,2,2,2,2,2)); } /////////////////////////////////////////////////////////////////////////////// proc intersection (list hn1, list hn2) "USAGE: intersection(hne1,hne2); hne1, hne2 lists ASSUME: @code{hne1, hne2} represent an HN expansion of an irreducible plane curve singularity (that is, are the output of @code{develop(f)}, or of @code{extdevelop(develop(f),n)}, or one entry of the list of HN data computed by @code{hnexpansion(f[,\"ess\"])}). RETURN: int, the intersection multiplicity of the irreducible plane curve singularities corresponding to @code{hne1} and @code{hne2}. SEE ALSO: hnexpansion, displayInvariants KEYWORDS: intersection multiplicity EXAMPLE: example intersection; shows an example " { //------------------ `intersect' ist schon reserviert ... -------------------- int i,j,s,sum,schnitt,unterschied; matrix a1=hn1[1]; matrix a2=hn2[1]; intvec h1=hn1[2]; intvec h2=hn2[2]; intvec n1=multiplicities(hn1); intvec n2=multiplicities(hn2); if (hn1[3]!=hn2[3]) { //-- die jeweils erste Zeile von hn1,hn2 gehoert zu verschiedenen Parametern - //---------------- d.h. beide Kurven schneiden sich transversal -------------- schnitt=n1[1]*n2[1]; // = mult(hn1)*mult(hn2) } else { //--------- die jeweils erste Zeile gehoert zum gleichen Parameter ----------- unterschied=0; for (s=1; (h1[s]==h2[s]) && (sncols(a2)) { j=ncols(a1); } else { j=ncols(a2); } unterschied=0; if ((h1[s]>0) && (s==size(h1))) { a1[s,h1[s]+1]=0; if (ncols(a1)<=ncols(a2)) { unterschied=1; } } if ((h2[s]>0) && (s==size(h2))) { a2[s,h2[s]+1]=0; if (ncols(a2)<=ncols(a1)) { unterschied=1; } } if (unterschied==1) { // mind. eine HNE war endlich matrix ma1[1][j]=a1[s,1..ncols(a1)]; // und bedarf der Fortsetzung matrix ma2[1][j]=a2[s,1..ncols(a2)]; // mit Nullen } else { if (ncols(a1)>ncols(a2)) { j=ncols(a2); } else { j=ncols(a1); } matrix ma1[1][j]=a1[s,1..j]; // Beschr. auf vergleichbaren matrix ma2[1][j]=a2[s,1..j]; // Teil (der evtl. y's enth.) } for (i=1; (ma1[1,i]==ma2[1,i]) && (i1)) { tester=tester-1; } if((multseq[tester]!=1) and (multseq[tester]!=k-tester)) { for (i=k+1; i<=tester+multseq[tester]; i++) { multseq[i]=1; } } //--- Ende T.Keilen --- 06.05.02 return(multseq); } //---------------------------- mehrere Zweige -------------------------------- else { list HNEs=#; int anzahl=size(HNEs); int maxlength=0; int bisher; intvec schnitt,ones; ones[anzahl]=0; ones=ones+1; // = 1,1,...,1 for (i=1; i= maxlength); maxlength=maxlength*(schnitt[i]+1 < maxlength) + (schnitt[i]+1)*(schnitt[i]+1 >= maxlength); } j=size(multsequence(HNEs[anzahl])); maxlength=maxlength*(j < maxlength) + j*(j >= maxlength); //-------------- Konstruktion der ersten zu berechnenden Matrix --------------- intmat allmults[maxlength][anzahl]; for (i=1; i<=maxlength; i++) { allmults[i,1..anzahl]=ones[1..anzahl]; } for (i=1; i<=anzahl; i++) { ones=multsequence(HNEs[i]); allmults[1..size(ones),i]=ones[1..size(ones)]; } //---------------------- Konstruktion der zweiten Matrix ---------------------- intmat separate[maxlength][anzahl]; for (i=1; i<=maxlength; i++) { k=1; bisher=0; if (anzahl==1) { separate[i,1]=1; } for (j=1; jncols(a2)) { j=ncols(a1); } else { j=ncols(a2); } unterschied=0; if ((h1[s]>0) && (s==size(h1))) { a1[s,h1[s]+1]=0; if (ncols(a1)<=ncols(a2)) { unterschied=1; } } if ((h2[s]>0) && (s==size(h2))) { a2[s,h2[s]+1]=0; if (ncols(a2)<=ncols(a1)) { unterschied=1; } } if (unterschied==1) { // mind. eine HNE war endlich matrix ma1[1][j]=a1[s,1..ncols(a1)]; // und bedarf der Fortsetzung matrix ma2[1][j]=a2[s,1..ncols(a2)]; // mit Nullen } else { if (ncols(a1)>ncols(a2)) { j=ncols(a2); } else { j=ncols(a1); } matrix ma1[1][j]=a1[s,1..j]; // Beschr. auf vergleichbaren matrix ma2[1][j]=a2[s,1..j]; // Teil (der evtl. y's enth.) } for (i=1; (ma1[1,i]==ma2[1,i]) && (i*z(1) x = []*z(1)^2+...+z(1)^<>*z(2) z(1) = []*z(2)^2+...+z(2)^<>*z(3) ....... .......................... z(r-1) = []*z(r)^2+[]*z(r)^3+...... @end example where @code{x},@code{y} are the first 2 variables of the basering. The values of @code{[]} are the coefficients of the Hamburger-Noether matrix, the values of @code{<>} are represented by @code{x} in the HN matrix.@* - if a second argument is given and if the input are the HN data of an irreducible plane curve singularity, return a ring containing an ideal @code{HNE} as described above.@* - if L corresponds to the output of @code{hnexpansion(f)} or to the list of HN data computed by @code{hnexpansion(f[,\"ess\"])}, @code{displayHNE(L[,n])} shows the HNE's of all branches of f in the format described above. The optional parameter is then ignored. NOTE: The 1st line of the above ideal (i.e., @code{HNE[1]}) means that @code{y=[]*z(0)^1+...}, the 2nd line (@code{HNE[2]}) means that @code{x=[]*z(1)^2+...}, so you can see which indeterminate corresponds to which line (it's also possible that @code{x} corresponds to the 1st line and @code{y} to the 2nd). SEE ALSO: develop, hnexpansion EXAMPLE: example displayHNE; shows an example " { if ((typeof(ldev[1])=="list") || (typeof(ldev[1])=="none")) { for (int i=1; i<=size(ldev); i++) { "// Hamburger-Noether development of branch nr."+string(i)+":"; displayHNE(ldev[i]);""; } return(); } //--------------------- Initialisierungen und Ringwechsel -------------------- matrix m=ldev[1]; intvec v=ldev[2]; int switch=ldev[3]; if (switch==-1) { "An error has occurred throughout the expansion, so there is no HNE."; return(ideal(0)); } def altring=basering; ///////////////////////////////////////////////////////// // Change by T. Keilen 08.06.2002 // ring + ring does not work if one ring is an algebraic extension /* if (parstr(basering)!="") { if (charstr(basering)!=string(char(basering))+","+parstr(basering)) { execute ("ring dazu=("+charstr(basering)+"),z(0.."+string(size(v)-1)+"),ls;"); } else { ring dazu=char(altring),z(0..size(v)-1),ls; } } else { ring dazu=char(altring),z(0..size(v)-1),ls; } def displayring=dazu+altring; */ execute("ring displayring=("+charstr(basering)+"),(z(0.."+string(size(v)-1)+"),"+varstr(basering)+"),(ls("+string(size(v))+"),"+ordstr(basering)+");"); // End change by T. Keilen ////////////////////////////////////////////////////////////// setring displayring; map holematrix=altring,0; // mappt nur die Monome vom Grad Null matrix m=holematrix(m); int i,j; // lossen: check the last row for finiteness (06/2004) int rowM=nrows(m); int colM=ncols(m); int undef_bd=v[size(v)]; if ( undef_bd<-1 ){ for (j=-undef_bd; j<=colM; j++) { m[rowM,j]=0; } } //--------------------- Erzeuge Matrix n mit n[i,j]=z(j-1)^i ----------------- matrix n[colM][rowM]; for (j=1; j<=rowM; j++) { for (i=1; i<=colM; i++) { n[i,j]=z(j-1)^i; } } matrix displaymat=m*n; ideal HNE; for (i=1; i=" +string(-undef_bd)+")"; } if (undef_bd==-1) { stHNE(size(v))=stHNE(size(v))+" + ..... (terms of degree >=" +string(colM+1)+")"; } if (switch==0) { stHNE(1) = " "+string(var(size(v)+2))+" = "+stHNE(1); } else { stHNE(1) = " "+string(var(size(v)+1))+" = "+stHNE(1); } stHNE(1); if (ncols(HNE)==1) {return();} if (switch==0) { stHNE(2) = " "+string(var(size(v)+1))+" = "+stHNE(2); } else { stHNE(2) = " "+string(var(size(v)+2))+" = "+stHNE(2); } stHNE(2); for (j=3; j<=ncols(HNE); j++){ stHNE(j)= " "+"z(" +string(j-2)+ ") = "+stHNE(j); stHNE(j); } return(); } if (rowM<2) { HNE[2]=z(0); } if (switch==0) { HNE[1] = HNE[1]-var(size(v)+2); HNE[2] = HNE[2]-var(size(v)+1); } else { HNE[1] = HNE[1]-var(size(v)+1); HNE[2] = HNE[2]-var(size(v)+2); } if (size(#) == 0) { HNE; return(); } if (size(#) != 0) { HNE; export(HNE); return(displayring); } } example { "EXAMPLE:"; echo = 2; ring r=0,(x,y),dp; poly f=x3+2xy2+y2; list hn=develop(f); displayHNE(hn); } /////////////////////////////////////////////////////////////////////////////// // procedures for reducible curves // /////////////////////////////////////////////////////////////////////////////// // proc newtonhoehne (poly f) // USAGE: newtonhoehne(f); f poly // ASSUME: basering = ...,(x,y),ds or ls // RETURN: list of intvec(x,y) of coordinates of the newtonpolygon of f // NOTE: This proc is only available in versions of Singular that know the // command system("newton",f); f poly // { // intvec nm = getnm(f); // if ((nm[1]>0) && (nm[2]>0)) { f=jet(f,nm[1]*nm[2],nm); } // list erg=system("newton",f); // int i; list Ausgabe; // for (i=1; i<=size(erg); i++) { Ausgabe[i]=leadexp(erg[i]); } // return(Ausgabe); // } /////////////////////////////////////////////////////////////////////////////// proc newtonpoly (poly f, int #) "USAGE: newtonpoly(f); f poly ASSUME: basering has exactly two variables; @* f is convenient, that is, f(x,0) != 0 != f(0,y). RETURN: list of intvecs (= coordinates x,y of the Newton polygon of f). NOTE: Procedure uses @code{execute}; this can be avoided by calling @code{newtonpoly(f,1)} if the ordering of the basering is @code{ls}. KEYWORDS: Newton polygon EXAMPLE: example newtonpoly; shows an example " { if (size(#)>=1) { if (typeof(#[1])=="int") { // this is done to avoid the "execute" command for procedures in // hnoether.lib def is_ls=#[1]; } } if (defined(is_ls)<=0) { def @Rold=basering; execute("ring @RR=("+charstr(basering)+"),("+varstr(basering)+"),ls;"); poly f=imap(@Rold,f); } intvec A=(0,ord(subst(f,var(1),0))); intvec B=(ord(subst(f,var(2),0)),0); intvec C,H; list L; int abbruch,i; poly hilf; L[1]=A; f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1])); if (defined(is_ls)) { map xytausch=basering,var(2),var(1); } else { map xytausch=@RR,var(2),var(1); } for (i=2; f!=0; i++) { abbruch=0; while (abbruch==0) { C=leadexp(f); if(jet(f,A[2]*C[1]-A[1]*C[2]-1,intvec(A[2]-C[2],C[1]-A[1]))==0) { abbruch=1; } else { f=jet(f,-C[1]-1,intvec(-1,0)); } } hilf=jet(f,A[2]*C[1]-A[1]*C[2],intvec(A[2]-C[2],C[1]-A[1])); H=leadexp(xytausch(hilf)); A=H[2],H[1]; L[i]=A; f=jet(f,A[2]*B[1]-1,intvec(A[2],B[1]-A[1])); } L[i]=B; if (defined(is_ls)) { return(L); } else { setring @Rold; return(L); } } example { "EXAMPLE:"; echo = 2; ring r=0,(x,y),ls; poly f=x5+2x3y-x2y2+3xy5+y6-y7; newtonpoly(f); } /////////////////////////////////////////////////////////////////////////////// proc is_NND (poly f, list #) "USAGE: is_NND(f[,mu,NP]); f poly, mu int, NP list of intvecs ASSUME: f is convenient, that is, f(x,0) != 0 != f(0,y);@* mu (optional) is Milnor number of f.@* NP (optional) is output of @code{newtonpoly(f)}. RETURN: int: 1 if f is Newton non-degenerate, 0 otherwise. SEE ALSO: newtonpoly KEYWORDS: Newton non-degenerate; Newton polygon EXAMPLE: example is_NND; shows examples " { int i; int i_print=printlevel-voice+2; if (size(#)==0) { int mu=milnor(f); list NP=newtonpoly(f); } else { if (typeof(#[1])=="int") { def mu=#[1]; def NP=#[2]; for (i=1;i<=size(NP);i++) { if (typeof(NP[i])!="intvec") { print("third input cannot be Newton polygon ==> ignored ") NP=newtonpoly(f); i=size(NP)+1; } } } else { print("second input cannot be Milnor number ==> ignored ") int mu=milnor(f); NP=newtonpoly(f); } } // computation of the Newton number: int s=size(NP); int nN=-NP[1][2]-NP[s][1]+1; intmat m[2][2]; for(i=1;i<=s-1;i++) { m=NP[i+1],NP[i]; nN=nN+det(m); } if(mu==nN) { // the Newton-polygon is non-degenerate // REFERENCE? (tfuer mehr als 2 Variable gilt nicht, dass mu=nu impliziert, // dass NP nicht ausgeartet ist!, Siehe KOMMENTAR in equising.lib in esIdeal) return(1); } else { return(0); } } example { "EXAMPLE:"; echo = 2; ring r=0,(x,y),ls; poly f=x5+y3; is_NND(f); poly g=(x-y)^5+3xy5+y6-y7; is_NND(g); // if already computed, one should give the Minor number and Newton polygon // as second and third input: int mu=milnor(g); list NP=newtonpoly(g); is_NND(g,mu,NP); } /////////////////////////////////////////////////////////////////////////////// proc charPoly(poly f, int M, int N) "USAGE: charPoly(f,M,N); f bivariate poly, M,N int: length and height of Newton polygon of f, which has to be only one line RETURN: the characteristic polynomial of f EXAMPLE: example charPoly; shows an example " { poly charp; int Np=N/ gcd(M,N); f=subst(f,var(1),1); for(charp=0; f<>0; f=f-lead(f)) { charp=charp+leadcoef(f)*var(2)^(leadexp(f)[2]/ Np);} return(charp); } example { "EXAMPLE:"; echo = 2; ring exring=0,(x,y),dp; charPoly(y4+2y3x2-yx6+x8,8,4); charPoly(y6+3y3x2-x4,4,6); } /////////////////////////////////////////////////////////////////////////////// proc find_in_list(list L,int p) "USAGE: find_in_list(L,p); L: list of intvec(x,y) (sorted in y: L[1][2]>=L[2][2]), int p >= 0 RETURN: int i: L[i][2]=p if existent; otherwise i with L[i][2]

p; i++) {;} return(i); } example { "EXAMPLE:"; echo = 2; list L = intvec(0,4), intvec(1,2), intvec(2,1), intvec(4,0); find_in_list(L,1); L[find_in_list(L,2)]; } /////////////////////////////////////////////////////////////////////////////// proc get_last_divisor(int M, int N) "USAGE: get_last_divisor(M,N); int M,N RETURN: int Q: M=q1*N+r1, N=q2*r1+r2, ..., ri=Q*r(i+1) (Euclidean alg.) EXAMPLE: example get_last_divisor; shows an example " { int R=M%N; int Q=M / N; while (R!=0) {M=N; N=R; R=M%N; Q=M / N;} return(Q) } example { "EXAMPLE"; echo = 2; ring r=0,(x,y),dp; get_last_divisor(12,10); } /////////////////////////////////////////////////////////////////////////////// proc redleit (poly f,intvec S, intvec E) "USAGE: redleit(f,S,E); f poly, S,E intvec(x,y) S,E are two different points on a line in the Newton diagram of f RETURN: poly g: all monomials of f which lie on or below that line NOTE: The main purpose is that if the line defined by S and E is part of the Newton polygon, the result is the quasihomogeneous leading form of f w.r.t. that line. SEE ALSO: newtonpoly EXAMPLE: example redleit; shows an example " { if (E[1]0) { " The HNE was already exact"; return(l); } else { if (Q==-1) { Q=ncols(m); } else { Q=-Q-1; } } int zeile=nrows(m); int spalten,i,M; ideal lastrow=m[zeile,1..Q]; int ringwechsel=(varstr(basering)!="x,y") or (ordstr(basering)!="ls(2),C"); //------------------------- Ringwechsel, falls noetig ------------------------ if (ringwechsel) { def altring = basering; int p = char(basering); if (charstr(basering)!=string(p)) { string tststr=charstr(basering); tststr=tststr[1..find(tststr,",")-1]; //-> "p^k" bzw. "p" if (tststr==string(p)) { if (size(parstr(basering))>1) { // ring (p,a,..),... execute("ring extdguenstig=("+charstr(basering)+"),(x,y),ls;"); } else { // ring (p,a),... string mipl=string(minpoly); ring extdguenstig=(p,`parstr(basering)`),(x,y),ls; if (mipl!="0") { execute("minpoly="+mipl+";"); } } } else { execute("ring extdguenstig=("+charstr(basering)+"),(x,y),ls;"); } } else { // charstr(basering)== p : no parameter ring extdguenstig=p,(x,y),ls; } export extdguenstig; map hole=altring,x,y; //----- map kann sehr zeitaufwendig sein, daher Vermeidung, wo moeglich: ----- if (nvars(altring)==2) { poly f=fetch(altring,f); } else { poly f=hole(f); } ideal a=hole(lastrow); } else { ideal a=lastrow; } list Newton=newtonpoly(f,1); int M1=Newton[size(Newton)-1][1]; // konstant number delt; if (Newton[size(Newton)-1][2]!=1) { " *** The transformed polynomial was not valid!!";} else { //--------------------- Fortsetzung der HNE ---------------------------------- while (Q delta=-d/c: ------- delt=-koeff(f,M,0)/koeff(f,M1,1); a[Q]=delt; dbprint(printlevel-voice+2,"a("+string(zeile-1)+","+string(Q)+") = "+string(delt)); if (Qncols(m)) { spalten=ncols(lastrow); } else { spalten=ncols(m); } matrix mneu[zeile][spalten]; for (i=1; i= maxspalte); } //------------- schreibe Ausgabe fuer hnezaehler-ten Zweig: ------------------ matrix ma[size(HNEaktu)-2][maxspalte]; for (i=1; i<=(size(HNEaktu)-2); i++) { if (ncols(HNEaktu[i+1]) > 1) { ma[i,1..ncols(HNEaktu[i+1])]=HNEaktu[i+1]; } else { ma[i,1]=HNEaktu[i+1][1];} } Ergebnis[hnezaehler]=list(ma,HNEaktu[1],transvers,HNEaktu[size(HNEaktu)]); kill ma; } return(Ergebnis); } /////////////////////////////////////////////////////////////////////////////// proc factorfirst(poly f, int M, int N) "USAGE : factorfirst(f,M,N); f poly, M,N int RETURN: number d such that f=const*(y^(N/e) - d*x^(M/e))^e, where e=gcd(M,N), 0 if such a d does not exist EXAMPLE: example factorfirst; shows an example " { number c = koeff(f,0,N); number delt; int eps,l; int p=char(basering); string ringchar=charstr(basering); if (c == 0) {"Something has gone wrong! I didn't get N correctly!"; exit;} int e = gcd(M,N); if (p==0) { delt = koeff(f,M/ e,N - N/ e) / (-1*e*c); } else { if (e%p != 0) { delt = koeff(f,M/ e,N - N/ e) / (-1*e*c); } else { eps = e; for (l = 0; eps%p == 0; l=l+1) { eps=eps/ p;} if (defined(HNDebugOn)) {e," -> ",eps,"*",p,"^",l;} delt = koeff(f,(M/ e)*p^l,(N/ e)*p^l*(eps-1)) / (-1*eps*c); if ((charstr(basering) != string(p)) and (delt != 0)) { //------ coefficient field is not Z/pZ => (p^l)th root is not identity ------- delt=0; if (defined(HNDebugOn)) { "trivial factorization not implemented for", "parameters---I've to use 'factorize'"; } } } } if (defined(HNDebugOn)) {"quasihomogeneous leading form:",f," = ",c, "* (y^"+string(N/ e),"-",delt,"* x^"+string(M/ e)+")^",e," ?";} if (f != c*(var(2)^(N/ e) - delt*var(1)^(M/ e))^e) {return(0);} else {return(delt);} } example { "EXAMPLE:"; echo = 2; ring exring=7,(x,y),dp; factorfirst(2*(y3-3x4)^5,20,15); factorfirst(x14+y7,14,7); factorfirst(x14+x8y3+y7,14,7); } /////////////////////////////////////////////////////////////////////////// proc hnexpansion(poly f,list #) "USAGE: hnexpansion(f[,\"ess\"]); f poly ASSUME: f is a bivariate polynomial (in the first 2 ring variables) RETURN: list @code{L}, containing Hamburger-Noether data of @code{f}: If the computation of the HNE required no field extension, @code{L} is a list of lists @code{L[i]} (corresponding to the output of @code{develop}, applied to a branch of @code{f}, but the last entry being omitted): @texinfo @table @asis @item @code{L[i][1]}; matrix: Each row contains the coefficients of the corresponding line of the Hamburger-Noether expansion (HNE) for the i-th branch. The end of the line is marked in the matrix by the first ring variable (usually x). @item @code{L[i][2]}; intvec: indicating the length of lines of the HNE @item @code{L[i][3]}; int: 0 if the 1st ring variable was transversal (with respect to the i-th branch), @* 1 if the variables were changed at the beginning of the computation, @* -1 if an error has occurred. @item @code{L[i][4]}; poly: the transformed equation of the i-th branch to make it possible to extend the Hamburger-Noether data a posteriori without having to do all the previous calculation once again (0 if not needed). @end table @end texinfo If the computation of the HNE required a field extension, the first entry @code{L[1]} of the list is a ring, in which a list @code{hne} of lists (the HN data, as above) and a polynomial @code{f} (image of @code{f} over the new field) are stored. @* If called with an additional input parameter, @code{hnexpansion} computes only one representative for each class of conjugate branches (over the ground field active when calling the procedure). In this case, the returned list @code{L} always has only two entries: @code{L[1]} is either a list of lists (the HN data) or a ring (as above), and @code{L[2]} is an integer vector (the number of branches in the respective conjugacy classes). NOTE: If f is known to be irreducible as a power series, @code{develop(f)} could be chosen instead to avoid a change of basering during the computations. @* Increasing @code{printlevel} leads to more and more comments. @* Having defined a variable @code{HNDebugOn} leads to a maximum number of comments. SEE ALSO: develop, extdevelop, param, displayHNE EXAMPLE: example hnexpansion; shows an example " { int essential; if (size(#)==1) { essential=1; } int field_ext; def altring=basering; //--------- Falls Ring (p^k,a),...: Wechsel in (p,a),... + minpoly ----------- if ( (find(charstr(basering),string(char(basering)))!=1) && (charstr(basering)<>"real")&& (charstr(basering)<>"complex") ) { string strmip=string(minpoly); string strf=string(f); execute("ring tempr=("+string(char(basering))+","+parstr(basering)+"),(" +varstr(basering)+"),dp;"); execute("minpoly="+strmip+";"); execute("poly f="+strf+";"); field_ext=1; def L=pre_HN(f,essential); if (size(L)==0) { return(list()); } def HNEring=L[1]; setring HNEring; if ((typeof(hne[1])=="ideal")) { return(list()); } if ((voice==2) && (printlevel > -1)) { "// Attention: The parameter",par(1),"may have changed its meaning!"; "// It needs no longer be a generator of the cyclic group of unities!"; } dbprint(printlevel-voice+2, "// result: "+string(size(hne))+" branch(es) successfully computed."); } else { def L=pre_HN(f,essential); if (size(L)==0) { return(list()); } if (L[2]==1) { field_ext=1; } intvec hne_conj=L[3]; def HNEring=L[1]; setring HNEring; if ((typeof(hne[1])=="ideal")) { return(list()); } dbprint(printlevel-voice+2, "// result: "+string(size(hne))+" branch(es) successfully computed."); } // ----- Lossen 10/02 : the branches have to be resorted to be able to // ----- display the multsequence in a nice way if (size(hne)>2) { int i,j,k,m; list dummy; int nbsave; int no_br = size(hne); intmat nbhd[no_br][no_br]; for (i=1;i= nbhd[i,j]) and (k map data to original ring setring altring; if ((npars(altring)==1) and (minpoly!=0)) { ring HNhelpring=char(altring),(a,x,y),ls; list hne=imap(HNEring,hne); setring altring; map mmm=HNhelpring,par(1),var(1),var(2); list hne=mmm(hne); kill mmm,HNhelpring; } else { list hne=fetch(HNEring,hne); } kill HNEring; if (essential==1) { dbprint(printlevel-voice+3,""+ "// No change of ring necessary, return value is a list: // first entry = list : HN expansion of essential branches. // second entry = intvec: numbers of conjugated branches "); return(list(hne,hne_conj)); } else { dbprint(printlevel-voice+3,""+ "// No change of ring necessary, return value is HN expansion."); return(hne); } } } example { "EXAMPLE:"; echo = 2; ring r=0,(x,y),dp; // First, an example which requires no field extension: list Hne=hnexpansion(x4-y6); size(Hne); // number of branches displayHNE(Hne); // HN expansion of branches param(Hne[1]); // parametrization of 1st branch param(Hne[2]); // parametrization of 2nd branch // An example which requires a field extension: list L=hnexpansion((x4-y6)*(y2+x4)); def R=L[1]; setring R; displayHNE(hne); basering; setring r; kill R; // Computing only one representative per conjugacy class: L=hnexpansion((x4-y6)*(y2+x4),"ess"); def R=L[1]; setring R; displayHNE(hne); L[2]; // number of branches in respective conjugacy classes } /////////////////////////////////////////////////////////////////////////////// static proc pre_HN (poly f, int essential) "NOTE: This procedure is only for internal use, it is called via hnexpansion RETURN: list: first entry = HNEring (containing list hne, poly f) second entry = 0 if no change of base ring necessary 1 if change of base ring necessary third entry = numbers of conjugates ( if essential = 1 ) if some error has occured, the empty list is returned " { def altring = basering; int p = char(basering); int field_ext; intvec hne_conj; //-------------------- Tests auf Zulaessigkeit von basering ------------------ if (charstr(basering)=="real") { " Singular cannot factorize over 'real' as ground field"; return(list()); } if (size(maxideal(1))<2) { " A univariate polynomial ring makes no sense !"; return(list()); } if (size(maxideal(1))>2) { dbprint(printlevel-voice+2, " Warning: all variables except the first two will be ignored!"); } if (find(charstr(basering),string(char(basering)))!=1) { " ring of type (p^k,a) not implemented"; //---------------------------------------------------------------------------- // weder primitives Element noch factorize noch map "char p^k" -> "char -p" // [(p^k,a)->(p,a) ground field] noch fetch //---------------------------------------------------------------------------- return(list()); } //----------------- Definition eines neuen Ringes: HNEring ------------------- string namex=varstr(1); string namey=varstr(2); if (string(char(altring))==charstr(altring)) { // kein Parameter, nicht 'real' ring HNEring = char(altring),(x,y),ls; map m=altring,x,y; poly f=m(f); export f; kill m; } else { string mipl=string(minpoly); if (mipl=="0") { "// ** WARNING: Algebraic extension of given ground field not possible!"; "// ** We try to develop this polynomial, but if the need for a field"; "// ** extension occurs during the calculation, we cannot proceed with"; "// ** the corresponding branches."; execute("ring HNEring=("+charstr(basering)+"),(x,y),ls;"); } else { string pa=parstr(altring); ring HNhelpring=p,`pa`,dp; execute("poly mipo="+mipl+";"); // Minimalpolynom in Polynom umgewandelt ring HNEring=(p,a),(x,y),ls; map getminpol=HNhelpring,a; mipl=string(getminpol(mipo)); // String umgewandelt mit 'a' als Param. execute("minpoly="+mipl+";"); // "minpoly=poly is not supported" kill HNhelpring, getminpol; } if (nvars(altring)==2) { poly f=fetch(altring,f); export f; } else { if (defined(pa)) { // Parameter hatte vorher anderen Namen als 'a' ring HNhelpring=p,(`pa`,x,y),ls; poly f=imap(altring,f); setring HNEring; map m=HNhelpring,a,x,y; poly f=m(f); kill HNhelpring; } else { map m=altring,x,y; poly f=m(f); } export f; kill m; } } if (defined(HNDebugOn)) {"received polynomial: ",f,", with x =",namex,", y =",namey;} //----------------------- Variablendefinitionen ------------------------------ int Abbruch,i,NullHNEx,NullHNEy; string str; list Newton,hne; // --- changed for SINGULAR 3: --- hne=ideal(0); export hne; //====================== Tests auf Zulaessigkeit des Polynoms ================ //-------------------------- Test, ob Einheit oder Null ---------------------- if (subst(subst(f,x,0),y,0)!=0) { dbprint(printlevel+1, "The given polynomial is a unit in the power series ring!"); setring altring; kill HNEring; return(list()); // there are no HNEs } if (f==0) { dbprint(printlevel+1,"The given polynomial is zero!"); setring altring; kill HNEring; return(list()); // there are no HNEs } //----------------------- Test auf Quadratfreiheit -------------------------- if ((p==0) and (size(charstr(basering))==1)) { //-------- Fall basering==0,... : Wechsel in Ring mit char >0 ---------------- // weil squarefree eine Standardbasis berechnen muss (verwendet Syzygien) // -- wenn f in diesem Ring quadratfrei ist, dann erst recht im Ring HNEring //---------------------------------------------------------------------------- int testerg=(polytest(f)==0); ring zweitring = 32003,(x,y),dp; map polyhinueber=HNEring,x,y; // fetch geht nicht poly f=polyhinueber(f); poly test_sqr=squarefree(f); if (test_sqr != f) { if (printlevel>0) { "Most probably the given polynomial is not squarefree. But the test was"; "made in characteristic 32003 and not 0 to improve speed. You can"; "(r) redo the test in char 0 (but this may take some time)"; "(c) continue the development, if you're sure that the polynomial IS", "squarefree"; if (testerg==1) { "(s) continue the development with a squarefree factor (*)";} "(q) or just quit the algorithm (default action)"; "";"Please enter the letter of your choice:"; str=read("")[1]; // reads one character } else { str="r"; } // printlevel <= 0: non-interactive behaviour setring HNEring; map polyhinueber=zweitring,x,y; if (str=="r") { poly test_sqr=squarefree(f); if (test_sqr != f) { if (printlevel>0) { "The given polynomial is in fact not squarefree."; } else { "The given polynomial is not squarefree!"; } "I'll continue with the radical."; f=test_sqr; } else { dbprint(printlevel, "everything is ok -- the polynomial is squarefree in characteristic 0"); } } else { if ((str=="s") and (testerg==1)) { "(*)attention: it could be that the factor is only one in char 32003!"; f=polyhinueber(test_sqr); } else { if (str<>"c") { setring altring; kill HNEring;kill zweitring; return(list());} else { "if the algorithm doesn't terminate, you were wrong...";} }} kill zweitring; if (defined(HNDebugOn)) {"I continue with the polynomial",f; } } else { setring HNEring; kill zweitring; } } //------------------ Fall Char > 0 oder Ring hat Parameter ------------------- else { poly test_sqr=squarefree(f); if (test_sqr != f) { if (printlevel>0) { if (test_sqr == 1) { "The given polynomial is of the form g^"+string(p)+","; "therefore not squarefree. You can:"; " (q) quit the algorithm (recommended) or"; " (f) continue with the full radical (using a factorization of the"; " pure power part; this could take much time)"; "";"Please enter the letter of your choice:"; str=read("")[1]; if (str<>"f") { str="q"; } } else { "The given polynomial is not squarefree."; if (p != 0) { " You can:"; " (c) continue with a squarefree divisor (but factors of the form g^" +string(p); " are lost; this is recommended, takes no extra time)"; " (f) continue with the full radical (using a factorization of the"; " pure power part; this could take some time)"; " (q) quit the algorithm"; "";"Please enter the letter of your choice:"; str=read("")[1]; if ((str<>"f") && (str<>"q")) { str="c"; } } else { "I'll continue with the radical."; str="c"; } } // endelse (test_sqr!=1) } else { "//** Error: The given polynomial is not squarefree!"; "//** Since the global variable `printlevel' has the value",printlevel, "we stop here."; "// Either call me again with a squarefree polynomial f or assign"; " printlevel=1;"; "// before calling me with a non-squarefree f."; "// If printlevel > 0, I present some possibilities how to proceed."; str="q"; } if (str=="q") { setring altring;kill HNEring; return(list()); } if (str=="c") { f=test_sqr; } if (str=="f") { f=allsquarefree(f,test_sqr); } } if (defined(HNDebugOn)) {"I continue with the polynomial",f; } } //====================== Ende Test auf Quadratfreiheit ======================= if (subst(subst(f,x,0),y,0)!=0) { "The polynomial is a unit in the power series ring. No HNE computed."; setring altring;kill HNEring; return(list()); } //---------------------- Test, ob f teilbar durch x oder y ------------------- if (subst(f,y,0)==0) { f=f/y; NullHNEy=1; } // y=0 is a solution if (subst(f,x,0)==0) { f=f/x; NullHNEx=1; } // x=0 is a solution Newton=newtonpoly(f,1); i=1; Abbruch=0; //---------------------------------------------------------------------------- // finde Eckpkt. des Newtonpolys, der den Teil abgrenzt, fuer den x transvers: // Annahme: Newton ist sortiert, s.d. Newton[1]=Punkt auf der y-Achse, // Newton[letzt]=Punkt auf der x-Achse //---------------------------------------------------------------------------- while ((i=(Newton[i][2]-Newton[i+1][2])) {Abbruch=1;} else {i=i+1;} } int grenze1=Newton[i][2]; int grenze2=Newton[i][1]; //---------------------------------------------------------------------------- // Stelle Ring bereit zur Uebertragung der Daten im Fall einer Koerperer- // weiterung. Definiere Objekte, die spaeter uebertragen werden. // Binde die Listen (azeilen,...) an den Ring (um sie nicht zu ueberschreiben // bei Def. in einem anderen Ring). //---------------------------------------------------------------------------- ring HNE_noparam = char(altring),(a,x,y),ls; poly f; list azeilen=ideal(0); list HNEs=ideal(0); list aneu=ideal(0); list faktoren=ideal(0); ideal deltais; poly delt; //----- hier steht die Anzahl bisher benoetigter Ringerweiterungen drin: ----- int EXTHNEnumber=0; list EXTHNEring; list HNE_RingDATA; int number_of_letztring; setring HNEring; number_of_letztring=0; // ================= Die eigentliche Berechnung der HNE: ===================== // ------- Berechne HNE von allen Zweigen, fuer die x transversal ist: ------- if (defined(HNDebugOn)) {"1st step: Treat Newton polygon until height",grenze1;} if (grenze1>0) { if (EXTHNEnumber>0){ EXTHNEring = EXTHNEring(1..EXTHNEnumber); } HNE_RingDATA = list(HNEring, HNE_noparam, EXTHNEnumber, EXTHNEring, number_of_letztring); list hilflist=HN(HNE_RingDATA,f,grenze1,1,essential,0,hne_conj,1); kill HNEring, HNE_noparam; if (EXTHNEnumber>0) { kill EXTHNEring(1..EXTHNEnumber);} def HNEring = hilflist[1][1]; def HNE_noparam = hilflist[1][2]; EXTHNEnumber = hilflist[1][3]; for (i=1; i<=EXTHNEnumber; i++) { def EXTHNEring(i)=hilflist[1][4][i]; } if (hilflist[2]==0) { setring HNEring; number_of_letztring=0; } else { setring EXTHNEring(hilflist[2]);} if (hilflist[3]==1){field_ext=1;} hne_conj=hilflist[5]; if (number_of_letztring != hilflist[2]) { // Ringwechsel in Prozedur HN map hole=HNE_noparam,transfproc,x,y; setring HNE_noparam; if (not(defined(f))) {poly f;} f=imap(HNEring,f); setring EXTHNEring(EXTHNEnumber); if (not(defined(f))) {poly f; f=hole(f); export f;} else {f=hole(f);} } number_of_letztring = hilflist[2]; kill hilflist; } if (NullHNEy==1) { if ((typeof(hne[1])=="ideal")) { hne=list(); } hne=hne+list(list(matrix(ideal(0,x)),intvec(1),int(0),poly(0))); if (hne_conj==0) { hne_conj=1; } else { hne_conj = hne_conj, 1; } } // --------------- Berechne HNE von allen verbliebenen Zweigen: -------------- if (defined(HNDebugOn)) {"2nd step: Treat Newton polygon until height",grenze2;} if (grenze2>0) { if (EXTHNEnumber>0){ EXTHNEring = EXTHNEring(1..EXTHNEnumber); } if (essential==1) { number_of_letztring=0; } if (number_of_letztring==0) { setring HNEring; } else { setring EXTHNEring(number_of_letztring); } map xytausch=basering,y,x; HNE_RingDATA = list(HNEring, HNE_noparam, EXTHNEnumber, EXTHNEring, number_of_letztring); list hilflist=HN(HNE_RingDATA,xytausch(f),grenze2,1,essential,1,hne_conj,1); kill HNEring, HNE_noparam; if (EXTHNEnumber>0){ kill EXTHNEring(1..EXTHNEnumber); } def HNEring = hilflist[1][1]; def HNE_noparam = hilflist[1][2]; EXTHNEnumber = hilflist[1][3]; for (i=1; i<=EXTHNEnumber; i++) { def EXTHNEring(i)=hilflist[1][4][i]; } if (hilflist[2]==0) { setring HNEring; number_of_letztring=0; } else { setring EXTHNEring(hilflist[2]); number_of_letztring=hilflist[2]; } if (hilflist[3]==1){field_ext=1;} hne_conj=hilflist[5]; kill hilflist; } if (NullHNEx==1) { if ((typeof(hne[1])=="ideal")) { hne=list(); } hne=hne+list(list(matrix(ideal(0,x)),intvec(1),int(1),poly(0))); if (hne_conj==0) { hne_conj=1; } else { hne_conj = hne_conj, 1; } } // --- aufraeumen --- if (defined(HNEakut)){ kill HNEakut,faktoren,deltais,transformiert,teiler,leitf; } if (defined(hilflist)) {kill hilflist;} if (defined(erg)) {kill erg;} if (defined(delt)) {kill delt;} if (defined(azeilen)) { kill azeilen;} if (defined(aneu)) { kill aneu;} if (defined(transfproc)) { kill transfproc;} if (defined(transf)) { kill transf;} if (not(defined(f))) { poly f = imap(HNEring,f); export f; } return(list(basering,field_ext,hne_conj)); } ////////////////////////////////////////////////////////////////////////////// proc essdevelop (poly f) "USAGE: essdevelop(f); f poly NOTE: command is obsolete, use hnexpansion(f,\"ess\") instead. SEE ALSO: hnexpansion, develop, extdevelop, param " { printlevel=printlevel+1; list Ergebnis=hnexpansion(f,1); printlevel=printlevel-1; return(Ergebnis); } /////////////////////////////////////////////////////////////////////////////// static proc HN (list HNE_RingDATA,poly fneu,int grenze,Aufruf_Ebene, essential,getauscht,intvec hne_conj,int conj_factor) "NOTE: This procedure is only for internal use, it is called via pre_HN RETURN: list: first entry = list of HNErings, second entry = number of new base ring (0 for HNEring, -1 for HNE_noparam, i for EXTHNEring(i)) third entry = 0 if no field extension necessary 1 if field extension necessary forth entry = HNEs (only if no change of basering) " { //---------- Variablendefinitionen fuer den unverzweigten Teil: -------------- if (defined(HNDebugOn)) {"procedure HN",Aufruf_Ebene;} int Abbruch,ende,i,j,k,e,M,N,Q,R,zeiger,zeile,zeilevorher,dd,ii; intvec hqs; int field_ext; int ring_changed, hneshift; intvec conjugates,conj2,conj1; list EXTHNEring; def HNEring = HNE_RingDATA[1]; def HNE_noparam = HNE_RingDATA[2]; int EXTHNEnumber = HNE_RingDATA[3]; for (i=1; i<=EXTHNEnumber; i++) { def EXTHNEring(i)=HNE_RingDATA[4][i]; } int number_of_letztring = HNE_RingDATA[5]; if (defined(basering)) { if (number_of_letztring==0) { kill HNEring; def HNEring=basering; } else { kill EXTHNEring(number_of_letztring); def EXTHNEring(number_of_letztring)=basering; } } else { if ( number_of_letztring==0) { setring HNEring; } else { setring EXTHNEring(number_of_letztring); } } if (not(defined(hne))) {list hne;} poly fvorher; list erg=ideal(0); list HNEs=ideal(0); // um die Listen an den Ring zu binden //-------------------- Bedeutung von Abbruch: -------------------------------- //------- 0:keine Verzweigung | 1:Verzweigung,nicht fertig | 2:fertig -------- // // Struktur von HNEs : Liste von Listen L (fuer jeden Zweig) der Form // L[1]=intvec (hqs), L[2],L[3],... ideal (die Zeilen (0,1,...) der HNE) // L[letztes]=poly (transformiertes f) //---------------------------------------------------------------------------- list Newton; number delt; int p = char(basering); // Ringcharakteristik list azeilen=ideal(0); ideal hilfid; intvec hilfvec; // ======================= der unverzweigte Teil: ============================ while (Abbruch==0) { Newton=newtonpoly(fneu,1); zeiger=find_in_list(Newton,grenze); if (Newton[zeiger][2] != grenze) {"Didn't find an edge in the Newton polygon!";} if (zeiger==size(Newton)-1) { if (defined(HNDebugOn)) {"only one relevant side in Newton polygon";} M=Newton[zeiger+1][1]-Newton[zeiger][1]; N=Newton[zeiger][2]-Newton[zeiger+1][2]; R = M%N; Q = M / N; //-------- 1. Versuch: ist der quasihomogene Leitterm reine Potenz ? ------ // (dann geht alles wie im irreduziblen Fall) //------------------------------------------------------------------------- e = gcd(M,N); delt=factorfirst(redleit(fneu,Newton[zeiger],Newton[zeiger+1]) /x^Newton[zeiger][1],M,N); if (delt==0) { if (defined(HNDebugOn)) {" The given polynomial is reducible !";} Abbruch=1; } if (Abbruch==0) { //----------- fneu,zeile retten fuer den Spezialfall (###): ------------- fvorher=fneu;zeilevorher=zeile; if (R==0) { //-------- transformiere fneu mit T1, wenn kein Abbruch nachher: ------ if (N>1) { fneu = T1_Transform(fneu,delt,M/ e); } else { ende=1; } if (defined(HNDebugOn)) {"a("+string(zeile)+","+string(Q)+") =",delt;} azeilen[zeile+1][Q]=delt; } else { //------------- R > 0 : transformiere fneu mit T2 --------------------- erg=T2_Transform(fneu,delt,M,N,referencepoly(Newton)); fneu=erg[1];delt=erg[2]; //----- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: ---- while (R!=0) { if (defined(HNDebugOn)) { "h("+string(zeile)+") =",Q; } hqs[zeile+1]=Q; // denn zeile beginnt mit dem Wert 0 //--------------- markiere das Zeilenende der HNE: ------------------- azeilen[zeile+1][Q+1]=x; zeile=zeile+1; //-------- Bereitstellung von Speicherplatz fuer eine neue Zeile: ---- azeilen[zeile+1]=ideal(0); M=N; N=R; R=M%N; Q=M / N; } if (defined(HNDebugOn)) {"a("+string(zeile)+","+string(Q)+") =",delt;} azeilen[zeile+1][Q]=delt; } if (defined(HNDebugOn)) {"transformed polynomial: ",fneu;} grenze=e; //----------------------- teste Abbruchbedingungen: --------------------- if (subst(fneu,y,0)==0) { // <==> y|fneu dbprint(printlevel-voice+3,"finite HNE of one branch found"); // voice abzufragen macht bei rekursiven procs keinen Sinn azeilen[zeile+1][Q+1]=x; //----- Q wird nur in hqs eingetragen, wenn der Spezialfall nicht // eintritt (siehe unten) ----- Abbruch=2; if (grenze>1) { if (jet(fneu,1,intvec(0,1))==0) { //- jet(...)=alle Monome von fneu, die nicht durch y2 teilbar sind - "THE TEST FOR SQUAREFREENESS WAS BAD!!"; " The polynomial was NOT squarefree!!!";} else { //----------------------- Spezialfall (###): ----------------------- // Wir haben das Problem, dass die HNE eines Zweiges hier abbricht, // aber ein anderer Zweig bis hierher genau die gleiche HNE hat, die // noch weiter geht // Loesung: mache Transform. rueckgaengig und behandle fneu im // Verzweigungsteil //------------------------------------------------------------------ Abbruch=1; fneu=fvorher;zeile=zeilevorher;grenze=Newton[zeiger][2]; }} else {fneu=0;} // fneu nicht mehr gebraucht - spare Speicher if (Abbruch==2) { hqs[zeile+1]=Q; } } // Spezialfall nicht eingetreten else { if (ende==1) { dbprint(printlevel-voice+2,"HNE of one branch found"); Abbruch=2; hqs[zeile+1]=-Q-1;} } } // end(if Abbruch==0) } // end(if zeiger...) else { Abbruch=1;} } // end(while Abbruch==0) // ===================== der Teil bei Verzweigung: =========================== if (Abbruch==1) { //---------- Variablendefinitionen fuer den verzweigten Teil: --------------- poly leitf,teiler,transformiert; list aneu=ideal(0); list faktoren; ideal deltais; list HNEakut=ideal(0); intvec eis; int zaehler,hnezaehler,zl,zl1,M1,N1,R1,Q1,needext; int numberofRingchanges,lastRingnumber,ringischanged,flag; string letztringname; zeiger=find_in_list(Newton,grenze); if (defined(HNDebugOn)) { "Branching part reached---Newton polygon :",Newton; "relevant part until height",grenze,", from",Newton[zeiger],"on"; } azeilen=list(hqs)+azeilen; // hat jetzt Struktur von HNEs: hqs in der 1.Zeile //======= Schleife fuer jede zu betrachtende Seite des Newtonpolygons: ====== for(i=zeiger; inumber_of_letztring)) { // ----- setze ring zurueck fuer neue Kante ----- // ---- (damit konjugierte Zweige erkennbar) ----- hneshift=hneshift+hnezaehler; hnezaehler=0; ring_changed=0; def SaveRing = EXTHNEring(EXTHNEnumber); setring SaveRing; if (not(defined(HNEs))) { // HN wurde zum 2.Mal von pre_HN aufgerufen list HNEs=ideal(0); } for (k=number_of_letztring+1; k<=EXTHNEnumber; k++) { kill EXTHNEring(k);} EXTHNEnumber=number_of_letztring; if (EXTHNEnumber==0) { setring HNEring; } else { setring EXTHNEring(EXTHNEnumber); } if (not(defined(HNEs))) { list HNEs; } HNEs=ideal(0); deltais=0; delt=0; if (defined(zerlege)) { kill zerlege; } } if (defined(HNDebugOn)) { "we consider side",Newton[i],Newton[i+1]; } M=Newton[i+1][1]-Newton[i][1]; N=Newton[i][2]-Newton[i+1][2]; R = M%N; Q = M / N; e=gcd(M,N); needext=1; letztringname=nameof(basering); lastRingnumber=EXTHNEnumber; faktoren=list(ideal(charPoly(redleit(fneu,Newton[i],Newton[i+1]) /(x^Newton[i][1]*y^Newton[i+1][2]),M,N) ), intvec(1)); // = (zu faktoriserendes Poly, 1) conjugates=conj_factor; //-- wechsle so lange in Ringerweiterungen, bis Leitform vollstaendig // in Linearfaktoren zerfaellt ----- for (numberofRingchanges=1; needext==1; numberofRingchanges++) { leitf=redleit(fneu,Newton[i],Newton[i+1])/ (x^Newton[i][1]*y^Newton[i+1][2]); delt=factorfirst(leitf,M,N); needext=0; if (delt==0) { //---------- Sonderbehandlung: faktorisiere einige Polynome ueber Q(a): -- if ((charstr(basering)=="0,a") and (essential==0)) { // ==================================================== // neu CL: 06.10.05 poly CHPOLY=charPoly(leitf,M,N); poly tstpoly; if (defined(faktoren)!=0) { // Test, damit kein Fehler eingebaut (vermutlich nicht notwendig) tstpoly = faktoren[1][1]^faktoren[2][1]; for (k=2; k<=size(faktoren[1]); k++) { tstpoly = tstpoly * faktoren[1][k]^faktoren[2][k]; } tstpoly = CHPOLY-tstpoly*(CHPOLY/tstpoly); kill CHPOLY; } if ((numberofRingchanges>1) and (defined(faktoren)!=0) and (tstpoly==0)) { def L_help=factorlist(faktoren,conjugates); faktoren=L_help[1]; conjugates=L_help[2]; kill L_help; } else { faktoren=factorize(charPoly(leitf,M,N)); conjugates=conj_factor; for (k=2;k<=size(faktoren[2]);k++) {conjugates=conjugates,conj_factor;} } kill tstpoly; // Ende neu (vorher nur else Fall) // ==================================================== } else { //------------------ faktorisiere das charakt. Polynom: ---------------- if ((numberofRingchanges==1) or (essential==0)) { def L_help=factorlist(faktoren,conjugates); faktoren=L_help[1]; conjugates=L_help[2]; kill L_help; } else { // eliminiere alle konjugierten Nullstellen bis auf eine: ideal hilf_id; for (zaehler=1; zaehler<=size(faktoren[1]); zaehler++) { hilf_id=factorize(faktoren[1][zaehler],1); if (size(hilf_id)>1) { poly fac=hilf_id[1]; dd=deg(fac); // Zur Sicherheit: if (deg(fac)==0) { fac=hilf_id[2]; } for (k=2;k<=size(hilf_id);k++) { dd=dd+deg(hilf_id[k]); if (deg(hilf_id[k]) wegwerfen! if (defined(HNDebugOn)) {"factor of leading form found:",teiler;} if (teiler/y2 == 0) { // --> Faktor hat die Form cy+d deltais[zaehler]=-subst(teiler,y,0)/koeff(teiler,0,1); //=-d/c eis[zaehler]=faktoren[2][j]; conj2[zaehler]=conjugates[j]; zaehler++; } else { dbprint(printlevel-voice+2, " Change of basering (field extension) necessary!"); if (defined(HNDebugOn)) { teiler,"is not yet properly factorized!"; } if (needext==0) { poly zerlege=teiler; } needext=1; field_ext=1; } } } // end(for j) } else { deltais=ideal(delt); eis=e; conj2=conj_factor; } if (defined(HNDebugOn)) {"roots of char. poly:";deltais; "with multiplicities:",eis;} if (needext==1) { //--------------------- fuehre den Ringwechsel aus: --------------------- ringischanged=1; if ((size(parstr(basering))>0) && string(minpoly)=="0") { " ** We've had bad luck! The HNE cannot be calculated completely!"; // HNE in transzendenter Erweiterung fehlgeschlagen kill zerlege; ringischanged=0; break; // weiter mit gefundenen Faktoren } if (parstr(basering)=="") { EXTHNEnumber++; def EXTHNEring(EXTHNEnumber) = splitring(zerlege); setring EXTHNEring(EXTHNEnumber); poly transf=0; poly transfproc=0; ring_changed=1; export transfproc; } else { if (numberofRingchanges>1) { // ein Ringwechsel hat nicht gereicht def helpring = splitring(zerlege,list(transf,transfproc,faktoren)); kill EXTHNEring(EXTHNEnumber); def EXTHNEring(EXTHNEnumber)=helpring; setring EXTHNEring(EXTHNEnumber); kill helpring; poly transf=erg[1]; poly transfproc=erg[2]; ring_changed=1; list faktoren=erg[3]; export transfproc; kill erg; } else { if (ring_changed==1) { // in dieser proc geschah schon Ringwechsel EXTHNEnumber++; def EXTHNEring(EXTHNEnumber) = splitring(zerlege,list(a,transfproc)); setring EXTHNEring(EXTHNEnumber); poly transf=erg[1]; poly transfproc=erg[2]; export transfproc; kill erg; } else { // parameter war vorher da EXTHNEnumber++; def EXTHNEring(EXTHNEnumber) = splitring(zerlege,a); setring EXTHNEring(EXTHNEnumber); poly transf=erg[1]; poly transfproc=transf; ring_changed=1; export transfproc; kill erg; }} } //----------------------------------------------------------------------- // transf enthaelt jetzt den alten Parameter des Ringes, der aktiv war // vor Beginn der Schleife (evtl. also ueber mehrere Ringwechsel // weitergereicht), // transfproc enthaelt den alten Parameter des Ringes, der aktiv war zu // Beginn der Prozedur, und der an die aufrufende Prozedur zurueckgegeben // werden muss // transf ist Null, falls der alte Ring keinen Parameter hatte, // das gleiche gilt fuer transfproc //----------------------------------------------------------------------- //---- Neudef. von Variablen, Uebertragung bisher errechneter Daten: ---- poly leitf,teiler,transformiert; list aneu=ideal(0); ideal deltais; number delt; setring HNE_noparam; if (defined(letztring)) { kill letztring; } if (EXTHNEnumber>1) { def letztring=EXTHNEring(EXTHNEnumber-1); } else { def letztring=HNEring; } if (not defined(fneu)) {poly fneu;} fneu=imap(letztring,fneu); if (not defined(f)) {poly f;} f=imap(letztring,f); if (not defined(hne)) {list hne;} hne=imap(letztring,hne); if (not defined(faktoren)) {list faktoren; } faktoren=imap(letztring,faktoren); if (not(defined(azeilen))){list azeilen,HNEs;} azeilen=imap(letztring,azeilen); HNEs=imap(letztring,HNEs); setring EXTHNEring(EXTHNEnumber); if (not(defined(hole))) { map hole; } hole=HNE_noparam,transf,x,y; poly fneu=hole(fneu); if (not defined(faktoren)) { list faktoren; faktoren=hole(faktoren); } if (not(defined(f))) { poly f=hole(f); list hne=hole(hne); export f,hne; } } } // end (while needext==1) bzw. for (numberofRingchanges) if (eis==0) { i++; continue; } if (ringischanged==1) { list erg,HNEakut; // dienen nur zum Sp. von Zwi.erg. ideal hilfid; erg=ideal(0); HNEakut=erg; hole=HNE_noparam,transf,x,y; setring HNE_noparam; if (not(defined(azeilen))){list azeilen,HNEs;} azeilen=imap(letztring,azeilen); HNEs=imap(letztring,HNEs); setring EXTHNEring(EXTHNEnumber); list azeilen=hole(azeilen); list HNEs=hole(HNEs); kill letztring; ringischanged=0; } //============ Schleife fuer jeden gefundenen Faktor der Leitform: ========= for (j=1; j<=size(eis); j++) { //---- Mache Transformation T1 oder T2, trage Daten in HNEs ein, // falls HNE abbricht: ---- //------------------------ Fall R==0: ------------------------------------- if (R==0) { transformiert = T1_Transform(fneu,number(deltais[j]),M/ e); if (defined(HNDebugOn)) { "a("+string(zeile)+","+string(Q)+") =",deltais[j]; "transformed polynomial: ",transformiert; } if (subst(transformiert,y,0)==0) { dbprint(printlevel-voice+3,"finite HNE found"); hnezaehler++; //-------- trage deltais[j],x ein in letzte Zeile, fertig: ------------- HNEakut=azeilen+list(poly(0)); // =HNEs[hnezaehler]; hilfid=HNEakut[zeile+2]; hilfid[Q]=deltais[j]; hilfid[Q+1]=x; HNEakut[zeile+2]=hilfid; HNEakut[1][zeile+1]=Q; // aktualisiere Vektor mit den hqs HNEs[hnezaehler]=HNEakut; conj1[hneshift+hnezaehler]=conj2[j]; if (eis[j]>1) { transformiert=transformiert/y; if (subst(transformiert,y,0)==0){"THE TEST FOR SQUAREFREENESS WAS BAD!" +"! The polynomial was NOT squarefree!!!";} else { //--- Spezialfall (###) eingetreten: Noch weitere Zweige vorhanden -- eis[j]=eis[j]-1; } } } } else { //------------------------ Fall R <> 0: --------------------------------- erg=T2_Transform(fneu,number(deltais[j]),M,N,referencepoly(Newton)); transformiert=erg[1];delt=erg[2]; if (defined(HNDebugOn)) {"transformed polynomial: ",transformiert;} if (subst(transformiert,y,0)==0) { dbprint(printlevel-voice+3,"finite HNE found"); hnezaehler++; //---------------- trage endliche HNE in HNEs ein: --------------------- HNEakut=azeilen; // dupliziere den gemeins. Anfang der HNE's zl=2; // (kommt schliesslich nach HNEs[hnezaehler]) //---------------------------------------------------------------------- // Werte von: zeile: aktuelle Zeilennummer der HNE (gemeinsamer Teil) // zl : die HNE spaltet auf; zeile+zl ist der Index fuer die // Zeile des aktuellen Zweigs; (zeile+zl-2) ist die // tatsaechl. Zeilennr. (bei 0 angefangen) der HNE // ([1] <- intvec(hqs), [2] <- 0. Zeile usw.) //---------------------------------------------------------------------- //----- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: ----- M1=M;N1=N;R1=R;Q1=M1/ N1; while (R1!=0) { if (defined(HNDebugOn)) { "h("+string(zeile+zl-2)+") =",Q1; } HNEakut[1][zeile+zl-1]=Q1; HNEakut[zeile+zl][Q1+1]=x; // markiere das Zeilenende der HNE zl=zl+1; //----- Bereitstellung von Speicherplatz fuer eine neue Zeile: -------- HNEakut[zeile+zl]=ideal(0); M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1; } if (defined(HNDebugOn)) { "a("+string(zeile+zl-2)+","+string(Q1)+") =",delt; } HNEakut[zeile+zl][Q1] =delt; HNEakut[zeile+zl][Q1+1]=x; HNEakut[1][zeile+zl-1] =Q1; // aktualisiere Vektor mit hqs HNEakut[zeile+zl+1]=poly(0); HNEs[hnezaehler]=HNEakut; conj1[hneshift+hnezaehler]=conj2[j]; //-------------------- Ende der Eintragungen in HNEs ------------------- if (eis[j]>1) { transformiert=transformiert/y; if (subst(transformiert,y,0)==0){"THE TEST FOR SQUAREFREENESS WAS BAD!" +" The polynomial was NOT squarefree!!!";} else { //--- Spezialfall (###) eingetreten: Noch weitere Zweige vorhanden -- eis[j]=eis[j]-1; }} } // endif (subst()==0) } // endelse (R<>0) //========== Falls HNE nicht abbricht: Rekursiver Aufruf von HN: ========== //------------------- Berechne HNE von transformiert ---------------------- if (subst(transformiert,y,0)!=0) { lastRingnumber=EXTHNEnumber; if (EXTHNEnumber>0){ EXTHNEring = EXTHNEring(1..EXTHNEnumber); } HNE_RingDATA = list( HNEring, HNE_noparam, EXTHNEnumber, EXTHNEring, lastRingnumber); if (defined(HNerg)) {kill HNerg;} list HNerg=HN(HNE_RingDATA,transformiert,eis[j],Aufruf_Ebene+1, essential,getauscht,hne_conj,conj2[j]); HNE_RingDATA = HNerg[1]; if (conj1==0) { conj1=HNerg[5]; } else { conj1 = conj1,HNerg[5]; } if (HNerg[3]==1) { field_ext=1; } if (HNerg[2]==lastRingnumber) { // kein Ringwechsel in HN aufgetreten if (not(defined(aneu))) { list aneu; } aneu = HNerg[4]; } else { // Ringwechsel aufgetreten if (defined(HNDebugOn)) {" ring change in HN(",Aufruf_Ebene+1,") detected";} EXTHNEnumber = HNerg[1][3]; for (ii=lastRingnumber+1; ii<=EXTHNEnumber; ii++) { def EXTHNEring(ii)=HNerg[1][4][ii]; } if (HNerg[2]==0) { setring HNEring; } else { setring EXTHNEring(HNerg[2]); } def tempRing=HNerg[4]; def aneu=imap(tempRing,HNEs); kill tempRing; //--- stelle lokale Variablen im neuen Ring wieder her, und rette // gegebenenfalls ihren Inhalt: ---- list erg,faktoren,HNEakut; ideal hilfid; erg=ideal(0); faktoren=erg; HNEakut=erg; poly leitf,teiler,transformiert; map hole=HNE_noparam,transfproc,x,y; setring HNE_noparam; if (lastRingnumber>0) { def letztring=EXTHNEring(lastRingnumber); } else { def letztring=HNEring; } if (not defined(HNEs)) {list HNEs;} HNEs=imap(letztring,HNEs); if (not defined(azeilen)) {list azeilen;} azeilen=imap(letztring,azeilen); if (not defined(deltais)) {ideal deltais;} deltais=imap(letztring,deltais); if (not defined(delt)) {poly delt;} delt=imap(letztring,delt); if (not defined(fneu)) {poly fneu;} fneu=imap(letztring,fneu); if (not defined(f)) {poly f;} f=imap(letztring,f); if (not defined(hne)) {list hne;} hne=imap(letztring,hne); setring EXTHNEring(EXTHNEnumber); list HNEs=hole(HNEs); list azeilen=hole(azeilen); ideal deltais=hole(deltais); number delt=number(hole(delt)); poly fneu=hole(fneu); if (not(defined(f))) { poly f=hole(f); list hne=hole(hne); export f,hne; } } //========== Verknuepfe bisherige HNE mit von HN gelieferten HNEs: ====== if (R==0) { HNEs,hnezaehler=constructHNEs(HNEs,hnezaehler,aneu,azeilen,zeile, deltais,Q,j); kill aneu; } else { for (zaehler=1; zaehler<=size(aneu); zaehler++) { hnezaehler++; HNEakut=azeilen; // dupliziere den gemeinsamen Anfang der HNE's zl=2; // (kommt schliesslich nach HNEs[hnezaehler]) //------------ Trage Beitrag dieser Transformation T2 ein: ------------- //------ Zur Bedeutung von zeile, zl: siehe Kommentar weiter oben ------ //----- vollziehe Euklid.Alg. nach, um die HN-Matrix zu berechnen: ----- M1=M;N1=N;R1=R;Q1=M1/ N1; while (R1!=0) { if (defined(HNDebugOn)) { "h("+string(zeile+zl-2)+") =",Q1; } HNEakut[1][zeile+zl-1]=Q1; HNEakut[zeile+zl][Q1+1]=x; // Markierung des Zeilenendes der HNE zl=zl+1; //----- Bereitstellung von Speicherplatz fuer eine neue Zeile: -------- HNEakut[zeile+zl]=ideal(0); M1=N1; N1=R1; R1=M1%N1; Q1=M1 / N1; } if (defined(HNDebugOn)) { "a("+string(zeile+zl-2)+","+string(Q1)+") =",delt; } HNEakut[zeile+zl][Q1]=delt; //-- Daten aus T2_Transform sind eingetragen; haenge Daten von HN an: -- hilfid=HNEakut[zeile+zl]; for (zl1=Q1+1; zl1<=ncols(aneu[zaehler][2]); zl1++) { hilfid[zl1]=aneu[zaehler][2][zl1]; } HNEakut[zeile+zl]=hilfid; // ------ vorher HNEs[.][zeile+zl]<-aneu[.][2], // jetzt [zeile+zl+1] <- [3] usw.: -------- for (zl1=3; zl1<=size(aneu[zaehler]); zl1++) { HNEakut[zeile+zl+zl1-2]=aneu[zaehler][zl1]; } //--- setze hqs zusammen: HNEs[hnezaehler][1]=HNEs[..][1],aneu[..][1] -- hilfvec=HNEakut[1],aneu[zaehler][1]; HNEakut[1]=hilfvec; //-------- weil nicht geht: liste[1]=liste[1],aneu[zaehler][1] --------- HNEs[hnezaehler]=HNEakut; } // end(for zaehler) kill aneu; } // endelse (R<>0) } // endif (subst()!=0) (weiteres Aufblasen mit HN) } // end(for j) (Behandlung der einzelnen delta_i) // ------------------------- new for essdevelop ---------------------------- if ((essential==1) and (defined(SaveRing))) { // ----- uebertrage Daten in gemeinsame Koerpererweiterung --------------- if (EXTHNEnumber>number_of_letztring) { // ----- fuer aktuelle Kante war Koerpererweiterung erforderlich ------- EXTHNEnumber++; if (not(defined(minPol))) { poly miniPol; } miniPol=minpoly; setring SaveRing; if (not(defined(miniPol))) { poly miniPol; } miniPol=minpoly; setring HNE_noparam; if (not(defined(a_x))){ map a_x,a_y; poly mp_save, mp_new; } mp_save=imap(SaveRing,miniPol); mp_new=imap(EXTHNEring(EXTHNEnumber-1),miniPol); if (mp_save==mp_new) { // Sonderfall: wieder gleicher Ring def EXTHNEring(EXTHNEnumber)=SaveRing; setring EXTHNEring(EXTHNEnumber); if (not(defined(f))) {poly f; f=hole(f); export f;} list dummyL=imap(EXTHNEring(EXTHNEnumber-1),HNEs); if (not(defined(HNEs))) { def HNEs=list(); } if ((size(HNEs)==1) and (typeof(HNEs[1])=="ideal")) {HNEs=list();} HNEs[size(HNEs)+1..size(HNEs)+size(dummyL)]=dummyL[1..size(dummyL)]; kill dummyL,SaveRing; } else { // verschiedene Ringe a_x=HNE_noparam,x,0,0; a_y=HNE_noparam,y,0,0; mp_save=a_x(mp_save); // minpoly aus SaveRing mit a --> x mp_new=a_y(mp_new); // minpoly aus SaveRing mit a --> y setring SaveRing; poly mp_new=imap(HNE_noparam,mp_new); list Lfac=factorize(mp_new,1); poly fac=Lfac[1][1]; for (k=2;k<=size(Lfac[1]);k++) { if (deg(Lfac[1][k])0) { if ((size(HNEs)>1) or (typeof(HNEs[1])!="ideal") or (size(HNEs[1])>0)) { if ((typeof(hne[1])=="ideal")) { hne=list(); } hne=hne+extractHNEs(HNEs,getauscht); if (hne_conj==0) { hne_conj=conj1; } else { hne_conj = hne_conj, conj1; } } } } else { // HN wurde rekursiv aufgerufen if (number_of_letztring!=EXTHNEnumber) { // Ringwechsel hatte stattgefunden string mipl_alt = string(minpoly); execute("ring tempRing = ("+charstr(basering)+"),("+varstr(basering)+ "),("+ordstr(basering)+");"); execute("minpoly="+ mipl_alt +";"); list HNEs=imap(EXTHNEring(EXTHNEnumber),HNEs); export HNEs; if (defined(HNDebugOn)) {" ! tempRing defined ! ";} } if (conj1!=0) { hne_conj=conj1; } else { hne_conj=conj_factor; } } if (EXTHNEnumber>0){ EXTHNEring = EXTHNEring(1..EXTHNEnumber); } HNE_RingDATA = list(HNEring, HNE_noparam, EXTHNEnumber, EXTHNEring); if (number_of_letztring==EXTHNEnumber) { return(list(HNE_RingDATA,EXTHNEnumber,field_ext,HNEs,hne_conj)); } else { if (defined(tempRing)) { return(list(HNE_RingDATA,EXTHNEnumber,field_ext,tempRing,hne_conj)); } return(list(HNE_RingDATA,EXTHNEnumber,field_ext,0,hne_conj)); } } /////////////////////////////////////////////////////////////////////////////// static proc constructHNEs (list HNEs,int hnezaehler,list aneu,list azeilen, int zeile,ideal deltais,int Q,int j) "NOTE: This procedure is only for internal use, it is called via HN" { int zaehler,zl; ideal hilfid; list hilflist; intvec hilfvec; for (zaehler=1; zaehler<=size(aneu); zaehler++) { hnezaehler++; // HNEs[hnezaehler]=azeilen; // dupliziere gemeins. Anfang //----------------------- trage neu berechnete Daten ein --------------------- hilfid=azeilen[zeile+2]; hilfid[Q]=deltais[j]; for (zl=Q+1; zl<=ncols(aneu[zaehler][2]); zl++) { hilfid[zl]=aneu[zaehler][2][zl]; } hilflist=azeilen; hilflist[zeile+2]=hilfid; //------------------ haenge uebrige Zeilen von aneu[] an: -------------------- for (zl=3; zl<=size(aneu[zaehler]); zl++) { hilflist[zeile+zl]=aneu[zaehler][zl]; } //--- setze die hqs zusammen: HNEs[hnezaehler][1]=HNEs[..][1],aneu[..][1] ---- if (hilflist[1]==0) { hilflist[1]=aneu[zaehler][1]; } else { hilfvec=hilflist[1],aneu[zaehler][1]; hilflist[1]=hilfvec; } HNEs[hnezaehler]=hilflist; } return(HNEs,hnezaehler); } /////////////////////////////////////////////////////////////////////////////// proc referencepoly (list newton) "USAGE: referencepoly(newton); newton is list of intvec(x,y) which represents points in the Newton diagram (e.g. output of the proc newtonpoly) RETURN: a polynomial which has newton as Newton diagram SEE ALSO: newtonpoly EXAMPLE: example referencepoly; shows an example " { poly f; for (int i=1; i<=size(newton); i++) { f=f+var(1)^newton[i][1]*var(2)^newton[i][2]; } return(f); } example { "EXAMPLE:"; echo = 2; ring exring=0,(x,y),ds; referencepoly(list(intvec(0,4),intvec(2,3),intvec(5,1),intvec(7,0))); } /////////////////////////////////////////////////////////////////////////////// proc factorlist (list L, list #) "USAGE: factorlist(L); L a list in the format of `factorize' RETURN: the nonconstant irreducible factors of L[1][1]^L[2][1] * L[1][2]^L[2][2] *...* L[1][size(L[1])]^... with multiplicities (same format as factorize) SEE ALSO: factorize EXAMPLE: example factorlist; shows an example " { int k; if ((size(#)>=1) and ((typeof(#[1])=="intvec") or (typeof(#[1])=="int"))) { int with_conj = 1; intvec C = #[1]; } else { int with_conj = 0; intvec C = L[2]; } // eine Sortierung der Faktoren eruebrigt sich, weil keine zwei versch. // red.Fakt. einen gleichen irred. Fakt. haben koennen (I.3.27 Diplarb.) int i,gross; list faktoren,hilf; intvec conjugates; ideal hil1,hil2; intvec v,w,hilf_conj; for (i=1; (L[1][i] == jet(L[1][i],0)) && (i1) { v=hilf[2]; faktoren=list(ideal(hil1[2..gross]),intvec(v[2..gross])); conjugates=intvec(hilf_conj[2..gross]); } else { faktoren=hilf; conjugates=hilf_conj; } } else { faktoren=L; conjugates=C; } for (i++; i<=size(L[2]); i++) { //------------------------- linearer Term -- irreduzibel --------------------- if (L[1][i] == jet(L[1][i],1)) { if (L[1][i] != jet(L[1][i],0)) { // konst. Faktoren eliminieren hil1=faktoren[1]; hil1[size(hil1)+1]=L[1][i]; faktoren[1]=hil1; v=faktoren[2],L[2][i]; conjugates=conjugates,C[i]; faktoren[2]=v; } } //----------------------- nichtlinearer Term -- faktorisiere ----------------- else { hilf=factorize(L[1][i]); hilf_conj=C[i]; for (k=2;k<=size(hilf[2]);k++){ hilf_conj=hilf_conj,C[i]; } hilf[2]=hilf[2]*L[2][i]; hil1=faktoren[1]; hil2=hilf[1]; gross=size(hil2); // hil2[1] ist konstant, wird weggelassen: hil1[(size(hil1)+1)..(size(hil1)+gross-1)]=hil2[2..gross]; // ideal+ideal does not work due to simplification; // ideal,ideal not allowed faktoren[1]=hil1; w=hilf[2]; v=faktoren[2],w[2..gross]; conjugates=conjugates,hilf_conj[2..gross]; faktoren[2]=v; } } if (with_conj==0) { return(faktoren); } else { return(list(faktoren,conjugates)); } // for essential development } example { "EXAMPLE:"; echo = 2; ring exring=0,(x,y),ds; list L=list(ideal(x,(x-y)^2*(x+y+1),x+y),intvec(2,2,1)); L; factorlist(L); } /////////////////////////////////////////////////////////////////////////////// proc delta "USAGE: delta(INPUT); INPUT a polynomial defining an isolated plane curve singularity at 0, or the Hamburger-Noether expansion thereof, i.e. the output of @code{develop(f)}, or the output of @code{hnexpansion(f)}, or the list of HN data computed by @code{hnexpansion(f)}. RETURN: int, the delta invariant of the singularity at 0, that is, the vector space dimension of R~/R, (R~ the normalization of the local ring of the singularity). NOTE: In case the Hamburger-Noether expansion of the curve f is needed for other purposes as well it is better to calculate this first with the aid of @code{hnexpansion} and use it as input instead of the polynomial itself. SEE ALSO: deltaLoc, invariants KEYWORDS: delta invariant EXAMPLE: example delta; shows an example " { if (typeof(#[1])=="poly") { // INPUT = polynomial defining the singularity list L=hnexpansion(#[1]); if (typeof(L[1])=="ring") { def altring = basering; def HNring = L[1]; setring HNring; return(delta(hne)); } else { return(delta(L)); } } if (typeof(#[1])=="ring") { // INPUT = HNEring of curve def HNring = #[1]; setring HNring; return(delta(hne)); } if (typeof(#[1])=="matrix") { // INPUT = hne of an irreducible curve return(invariants(#)[5]/2); } else { // INPUT = hne of a reducible curve list INV=invariants(#); return(INV[size(INV)][3]); } } example { "EXAMPLE:"; echo = 2; ring r = 32003,(x,y),ds; poly f = x25+x24-4x23-1x22y+4x22+8x21y-2x21-12x20y-4x19y2+4x20+10x19y +12x18y2-24x18y-20x17y2-4x16y3+x18+60x16y2+20x15y3-9x16y -80x14y3-10x13y4+36x14y2+60x12y4+2x11y5-84x12y3-24x10y5 +126x10y4+4x8y6-126x8y5+84x6y6-36x4y7+9x2y8-1y9; delta(f); } ///////////////////////////////////////////////////////////////////////////////