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