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