1 | //last change: 13.02.2001 (Eric Westenberger) |
---|
2 | ////////////////////////////////////////////////////////////////////////////// |
---|
3 | version="$Id: triang.lib,v 1.12 2008-10-01 14:20:22 Singular Exp $"; |
---|
4 | category="Symbolic-numerical solving"; |
---|
5 | info=" |
---|
6 | LIBRARY: triang.lib Decompose Zero-dimensional Ideals into Triangular Sets |
---|
7 | AUTHOR: D. Hillebrand |
---|
8 | |
---|
9 | PROCEDURES: |
---|
10 | triangL(G); Decomposition of (G) into triangular systems (Lazard). |
---|
11 | triangLfak(G); Decomp. of (G) into tri. systems plus factorization. |
---|
12 | triangM(G[,.]); Decomposition of (G) into triangular systems (Moeller). |
---|
13 | triangMH(G[,.]); Decomp. of (G) into tri. syst. with disjoint varieties. |
---|
14 | "; |
---|
15 | |
---|
16 | LIB "general.lib"; // Muss geladen werden fuer den Befehl sort(). |
---|
17 | LIB "elim.lib"; // Muss geladen werden fuer sat(). |
---|
18 | |
---|
19 | |
---|
20 | /////////////////////////////////////////////////////////////////////////////// |
---|
21 | // |
---|
22 | // Der Lazard-Algorithmus |
---|
23 | // |
---|
24 | /////////////////////////////////////////////////////////////////////////////// |
---|
25 | |
---|
26 | proc triangL (ideal G) |
---|
27 | "USAGE: triangL(G); G=ideal |
---|
28 | ASSUME: G is the reduced lexicographical Groebner basis of the |
---|
29 | zero-dimensional ideal (G), sorted by increasing leading terms. |
---|
30 | RETURN: a list of finitely many triangular systems, such that |
---|
31 | the union of their varieties equals the variety of (G). |
---|
32 | NOTE: Algorithm of Lazard (see: Lazard, D.: Solving zero-dimensional |
---|
33 | algebraic systems, J. Symb. Comp. 13, 117 - 132, 1992). |
---|
34 | EXAMPLE: example triangL; shows an example |
---|
35 | " |
---|
36 | { |
---|
37 | |
---|
38 | // Test, ob G Groebnerbasis eines nulldimensionalen Ideals ist. |
---|
39 | if (attrib(G,"isSB") <> 1) |
---|
40 | { |
---|
41 | ERROR("The input is no groebner basis."); |
---|
42 | } |
---|
43 | else |
---|
44 | { |
---|
45 | if (dim(G) <> 0) |
---|
46 | { |
---|
47 | ERROR("ideal is not zero-dimensional."); |
---|
48 | } |
---|
49 | } |
---|
50 | |
---|
51 | if (size(G) <= 2) // schon Dreiecksbasis |
---|
52 | { |
---|
53 | return(list(G)); |
---|
54 | } |
---|
55 | |
---|
56 | // Noetige Optionen setzen. |
---|
57 | // obachman: save options so that tey can be reset on exit |
---|
58 | intvec ovec = option(get); |
---|
59 | option(redSB); |
---|
60 | option(returnSB); |
---|
61 | list B; // Bearbeitungsliste |
---|
62 | list F; // Ergebnisliste = Triangulierung |
---|
63 | ideal T, T_neu; |
---|
64 | poly p,p_neu,r,s,lt,lc,lcr; |
---|
65 | list inv; |
---|
66 | int u,v; |
---|
67 | int l,h; |
---|
68 | string st; |
---|
69 | |
---|
70 | T = ideal(G[1]); |
---|
71 | attrib(T,"isSB",1); |
---|
72 | B = list(list(T,2)); |
---|
73 | |
---|
74 | while (size(B) > 0) |
---|
75 | { |
---|
76 | T = B[1][1]; |
---|
77 | h = B[1][2]; |
---|
78 | B = Listenrest(B); |
---|
79 | p = G[h]; |
---|
80 | v = lvar(p); |
---|
81 | lc = lcoef(p,v); |
---|
82 | st=string(nvars(basering)-v+1); |
---|
83 | dbprint(string(timer)+" S"+st+": invertiere Leitkoeffizient von G["+string(h)+"]."); |
---|
84 | |
---|
85 | // Invertiere den Leitkoeffizienten von G[h] bzgl. var(v) modulo (T). |
---|
86 | inv = invertieren(lc,T); |
---|
87 | |
---|
88 | dbprint(string(timer)+" S"+st+": Anzahl Dreiecksbasen: "+string(size(inv))+"."); |
---|
89 | |
---|
90 | while (size(inv) > 0) |
---|
91 | { |
---|
92 | r = inv[1][1]; |
---|
93 | s = inv[1][2]; |
---|
94 | T = inv[1][3]; |
---|
95 | inv = Listenrest(inv); |
---|
96 | if (r == 0) // r == 0?, dann p nicht der ggt der Stufe, |
---|
97 | { |
---|
98 | dbprint(string(timer)+" S"+st+": Leitkoeffizient == 0."); |
---|
99 | B = list(list(T,h+1)) + B; |
---|
100 | } |
---|
101 | else // ansonsten ggt der Stufe gefunden. |
---|
102 | { |
---|
103 | dbprint(string(timer)+" S"+st+": ggt gefunden."); |
---|
104 | lt = var(v)**degv(p,v); |
---|
105 | p_neu = cleardenom(s*lt + reduce(r*(p-lc*lt),T)); |
---|
106 | T_neu = T,p_neu; |
---|
107 | attrib(T_neu,"isSB",1); |
---|
108 | |
---|
109 | // Restlichen Polynome der gleichen Stufe uebergehen. |
---|
110 | for (l = h+1; l <= size(G); l++) |
---|
111 | { |
---|
112 | u = lvar(G[l]); |
---|
113 | if (u <> v) {break;} |
---|
114 | } |
---|
115 | |
---|
116 | // Durchsuche restliche Stufen nach trivialen ggt's. |
---|
117 | while (l <= size(G)) |
---|
118 | { |
---|
119 | lc = lcoef(G[l],u); |
---|
120 | lcr = reduce(lc,T_neu); |
---|
121 | |
---|
122 | while (lcr == 0) // Gehe bis zum ersten Polynom <> 0 |
---|
123 | { // modulo (T_neu). |
---|
124 | l++; |
---|
125 | lc = lcoef(G[l],u); |
---|
126 | lcr = reduce(lc,T_neu); |
---|
127 | } |
---|
128 | |
---|
129 | if (deg(lcr) == 0) // Polynom <> 0 normiert?, dann |
---|
130 | { // an Dreiecksbasis anhaengen, |
---|
131 | dbprint(string(timer)+" S"+string(nvars(basering)-u+1)+": Schnellerkennung: ggt gefunden."); |
---|
132 | lt = var(u)**degv(G[l],u); |
---|
133 | p_neu = cleardenom(lcr*lt + |
---|
134 | reduce(G[l]-lc*lt,T_neu)); |
---|
135 | T_neu = T_neu,p_neu; |
---|
136 | attrib(T_neu,"isSB",1); |
---|
137 | u = lvar(G[l]); |
---|
138 | l++; |
---|
139 | while (l <= size(G)) // Rest der Stufe uebergehen. |
---|
140 | { |
---|
141 | if (lvar(G[l]) <> u) {u = lvar(G[l]); break;} |
---|
142 | l++; |
---|
143 | } |
---|
144 | } |
---|
145 | else // sonst nicht normierte Polynome auf der Stufe. |
---|
146 | { |
---|
147 | break; |
---|
148 | } |
---|
149 | } |
---|
150 | if (l > size(G)) // Ende von G erreicht, dann Dreiecks- |
---|
151 | { // basis der Triangulierung gefunden, |
---|
152 | dbprint(string(timer)+" S"+st+": Dreiecksbasis der Triangulierung gefunden."); |
---|
153 | F = F + list(T_neu); |
---|
154 | } |
---|
155 | else // sonst T_neu,l erneut bearbeiten. |
---|
156 | { |
---|
157 | B = list(list(T_neu,l)) + B; |
---|
158 | } |
---|
159 | } |
---|
160 | } |
---|
161 | } |
---|
162 | option(set, ovec); |
---|
163 | return(F); |
---|
164 | } |
---|
165 | //-------------------------------- example ---------------------------------- |
---|
166 | example |
---|
167 | { "EXAMPLE:"; echo = 2; |
---|
168 | ring rC5 = 0,(e,d,c,b,a),lp; |
---|
169 | triangL(stdfglm(cyclic(5))); |
---|
170 | } |
---|
171 | |
---|
172 | /////////////////////////////////////////////////////////////////////////////// |
---|
173 | |
---|
174 | proc triangLfak (ideal G) |
---|
175 | "USAGE: triangLfak(G); G=ideal |
---|
176 | ASSUME: G is the reduced lexicographical Groebner basis of the |
---|
177 | zero-dimensional ideal (G), sorted by increasing leading terms. |
---|
178 | RETURN: a list of finitely many triangular systems, such that |
---|
179 | the union of their varieties equals the variety of (G). |
---|
180 | NOTE: Algorithm of Lazard with factorization (see: Lazard, D.: Solving |
---|
181 | zero-dimensional algebraic systems, J. Symb. Comp. 13, 117 - 132, 1992). |
---|
182 | REMARK: each polynomial of the triangular systems is factorized. |
---|
183 | EXAMPLE: example triangLfak; shows an example |
---|
184 | " |
---|
185 | { |
---|
186 | return(triangLbas(G,2)); |
---|
187 | } |
---|
188 | //-------------------------------- example ---------------------------------- |
---|
189 | example |
---|
190 | { "EXAMPLE:"; echo = 2; |
---|
191 | ring rC5 = 0,(e,d,c,b,a),lp; |
---|
192 | triangLfak(stdfglm(cyclic(5))); |
---|
193 | } |
---|
194 | |
---|
195 | /////////////////////////////////////////////////////////////////////////////// |
---|
196 | |
---|
197 | static proc triangLbas (ideal G, list #) |
---|
198 | "USAGE: triangLbas(G[,i]); G reduzierte lexikographische Groebnerbasis |
---|
199 | des nulldimensionalen Ideals (G), nach |
---|
200 | Leittermen aufsteigend sortiert, |
---|
201 | i = 1 oder 2 (ohne Angabe i = 1). |
---|
202 | RETURN: Triangulierung von (G). |
---|
203 | Ist i == 2, dann wird jedes Polynom der Dreiecksbasen |
---|
204 | der Triangulierung zusaetzlich faktorisiert. |
---|
205 | NOTE: Algorithmus von Lazard (siehe: Lazard, D.: Solving |
---|
206 | zero-dimensional algebraic systems, |
---|
207 | J. Symb. Comp., Bd. 13, S. 117 - 132, 1992). |
---|
208 | " |
---|
209 | { |
---|
210 | |
---|
211 | // Test, ob G Groebnerbasis eines nulldimensionalen Ideals ist. |
---|
212 | if (attrib(G,"isSB") <> 1) |
---|
213 | { |
---|
214 | ERROR("The input is no groebner basis."); |
---|
215 | } |
---|
216 | else |
---|
217 | { |
---|
218 | if (dim(G) <> 0) |
---|
219 | { |
---|
220 | ERROR("ideal is not zero-dimensional."); |
---|
221 | } |
---|
222 | } |
---|
223 | |
---|
224 | // Noetige Optionen setzen. |
---|
225 | // obachman: save options so that tey can be reset on exit |
---|
226 | intvec ovec = option(get); |
---|
227 | option(redSB); |
---|
228 | option(returnSB); |
---|
229 | // Faktorisierungsschalter setzen. |
---|
230 | int fak; |
---|
231 | if (size(#) > 0 && typeof(#[1]) == "int") |
---|
232 | { |
---|
233 | if (#[1] == 2) {fak = 1;} |
---|
234 | } |
---|
235 | |
---|
236 | if (size(G) <= 2) // schon Dreiecksbasis |
---|
237 | { |
---|
238 | list E; |
---|
239 | if (fak) |
---|
240 | { |
---|
241 | E = faktorisiere_DB(G); |
---|
242 | } |
---|
243 | else |
---|
244 | { |
---|
245 | E = list(G); |
---|
246 | } |
---|
247 | option(set, ovec); |
---|
248 | return(E); |
---|
249 | } |
---|
250 | |
---|
251 | list B; // Bearbeitungsliste |
---|
252 | list F; // Ergebnisliste = Triangulierung |
---|
253 | list T_neu; |
---|
254 | ideal T; |
---|
255 | poly p,p_neu,r,s,lt,lc; |
---|
256 | list inv; |
---|
257 | int v; |
---|
258 | int l,h; |
---|
259 | string st; |
---|
260 | |
---|
261 | // B initialisieren |
---|
262 | if (fak) |
---|
263 | { |
---|
264 | B = H_anhaengen(faktorisiere_letzten(G[1]),2); |
---|
265 | } |
---|
266 | else |
---|
267 | { |
---|
268 | T = ideal(G[1]); |
---|
269 | attrib(T,"isSB",1); |
---|
270 | B = list(list(T,2)); |
---|
271 | } |
---|
272 | |
---|
273 | while (size(B) > 0) |
---|
274 | { |
---|
275 | T = B[1][1]; |
---|
276 | h = B[1][2]; |
---|
277 | B = Listenrest(B); |
---|
278 | p = G[h]; |
---|
279 | v = lvar(p); |
---|
280 | lc = lcoef(p,v); |
---|
281 | st=string(nvars(basering)-v+1); |
---|
282 | dbprint(string(timer)+" S"+st+": invertiere Leitkoeffizient von G["+string(h)+"]."); |
---|
283 | |
---|
284 | // invertiere den Leitkoeffizienten von G[h] bzgl. var(v) modulo (T). |
---|
285 | inv = invertieren(lc,T); |
---|
286 | dbprint(string(timer)+" S"+st+": Anzahl Dreiecksbasen: "+string(size(inv))+"."); |
---|
287 | |
---|
288 | while (size(inv) > 0) |
---|
289 | { |
---|
290 | r = inv[1][1]; |
---|
291 | s = inv[1][2]; |
---|
292 | T = inv[1][3]; |
---|
293 | inv = Listenrest(inv); |
---|
294 | if (r == 0) // r == 0?, dann p nicht der ggt der Stufe, |
---|
295 | { |
---|
296 | dbprint(string(timer)+" S"+st+": Leitkoeffizient == 0."); |
---|
297 | B = list(list(T,h+1)) + B; |
---|
298 | } |
---|
299 | else // ansonsten ggt der Stufe gefunden. |
---|
300 | { |
---|
301 | dbprint(string(timer)+" S"+st+": ggt gefunden."); |
---|
302 | lt = var(v)**degv(p,v); |
---|
303 | p_neu = cleardenom(s*lt + reduce(r*(p-lc*lt),T)); |
---|
304 | if (fak) |
---|
305 | { |
---|
306 | T_neu = faktorisiere_letzten(list(T + p_neu)); |
---|
307 | } |
---|
308 | else |
---|
309 | { |
---|
310 | T_neu = list(T + p_neu); |
---|
311 | } |
---|
312 | |
---|
313 | // Restlichen Polynome der gleichen Stufe uebergehen. |
---|
314 | for (l = h+1; l <= size(G); l++) |
---|
315 | { |
---|
316 | if (lvar(G[l]) <> v) {break;} |
---|
317 | } |
---|
318 | |
---|
319 | if (l > size(G)) // Ende von G erreicht, dann Dreiecks- |
---|
320 | { // basis der Triangulierung gefunden, |
---|
321 | dbprint(string(timer)+" S"+st+": Dreiecksbasis der Triangulierung gefunden."); |
---|
322 | F = F + T_neu; |
---|
323 | } |
---|
324 | else // sonst T_neu,l erneut bearbeiten. |
---|
325 | { |
---|
326 | B = H_anhaengen(T_neu,l) + B; |
---|
327 | } |
---|
328 | } |
---|
329 | } |
---|
330 | } |
---|
331 | option(set, ovec); |
---|
332 | return(F); |
---|
333 | } |
---|
334 | |
---|
335 | /////////////////////////////////////////////////////////////////////////////// |
---|
336 | |
---|
337 | static proc invertieren (poly p, ideal T) |
---|
338 | "USAGE: invertieren(p,T); p Polynom, T reduzierte Dreiecksbasis. |
---|
339 | RETURN: Liste von Tripeln (Listen) bestehend aus einem Polynom ri, |
---|
340 | einem Element aus dem Grundring si, und einer |
---|
341 | reduzierten Dreiecksbasis Ti, i = 1,...,m, so dass |
---|
342 | [T1,...,Tm] eine Triangulierung von T ist und |
---|
343 | p = 0 mod (Ti) falls ri = 0, ansonsten |
---|
344 | ri*p = si mod (Ti) fuer i = 1,...,m. |
---|
345 | " |
---|
346 | { |
---|
347 | // Triviale Faelle vorweg behandeln. |
---|
348 | p = reduce(p,T); |
---|
349 | if (p == 0) {return(list(list(0,0,T)));} |
---|
350 | if (deg(p) == 0) {return(list(list(1,p,T)));} |
---|
351 | |
---|
352 | list inv; |
---|
353 | int zerlegt; // zerlegt ist Schalter dafuer, ob T zerlegt wurde, |
---|
354 | export zerlegt; // einzige globale Variable. |
---|
355 | list toSee = list(list(p,T)); // zu bearbeitende Paare |
---|
356 | list toSave; // Ergebnisliste |
---|
357 | poly pr; |
---|
358 | |
---|
359 | while (size(toSee) > 0) |
---|
360 | { |
---|
361 | zerlegt = 0; // zerlegt = FALSE |
---|
362 | p = toSee[1][1]; |
---|
363 | T = toSee[1][2]; |
---|
364 | toSee = Listenrest(toSee); |
---|
365 | |
---|
366 | // invertieren_oT aufrufen, um p+(T) zu invertieren. |
---|
367 | inv = invertieren_oT(p,T); |
---|
368 | |
---|
369 | if (zerlegt) // T zerlegt?, dann pro neuer Dreiecksbasis |
---|
370 | { // weiteruntersuchen, |
---|
371 | pr = reduce(p,inv[1]); |
---|
372 | if (pr == 0) // es sei denn p ist reduziert gleich 0, |
---|
373 | { |
---|
374 | toSave = list(list(0,0,inv[1])) + toSave; |
---|
375 | attrib(toSave[1][3],"isSB",1); |
---|
376 | } |
---|
377 | else |
---|
378 | { |
---|
379 | toSee = list(list(pr,inv[1])) + toSee; |
---|
380 | attrib(toSee[1][2],"isSB",1); |
---|
381 | } |
---|
382 | pr = reduce(p,inv[2]); |
---|
383 | if (pr == 0) |
---|
384 | { |
---|
385 | toSave = list(list(0,0,inv[2])) + toSave; |
---|
386 | attrib(toSave[1][3],"isSB",1); |
---|
387 | } |
---|
388 | else |
---|
389 | { |
---|
390 | toSee = list(list(pr,inv[2])) + toSee; |
---|
391 | attrib(toSee[1][2],"isSB",1); |
---|
392 | } |
---|
393 | |
---|
394 | } |
---|
395 | else // ansonsten Quasi-Inverses gefunden. |
---|
396 | { |
---|
397 | toSave = list(list(inv[1],inv[2],T)) + toSave; |
---|
398 | attrib(toSave[1][3],"isSB",1); |
---|
399 | } |
---|
400 | } |
---|
401 | kill zerlegt; |
---|
402 | return(toSave); |
---|
403 | } |
---|
404 | |
---|
405 | /////////////////////////////////////////////////////////////////////////////// |
---|
406 | |
---|
407 | static proc invertieren_oT (poly p, ideal T) |
---|
408 | "USAGE: invertieren_oT(p,T); T reduzierte Dreiecksbasis, |
---|
409 | p <> 0 Polynom, irreduzibel modulo (T). |
---|
410 | RETURN: Entweder ein Quasi-Inverses (r,s) von p modulo (T), d.h. |
---|
411 | r*p = s mod (T) und s im Grundring, |
---|
412 | oder eine Triangulierung [T1,T2] von T |
---|
413 | mit (Ti) <> (T) fuer i = 1,2. |
---|
414 | " |
---|
415 | { |
---|
416 | // Quasi-Inverses von Konstante klar. |
---|
417 | if (deg(p) == 0) {return(list(1,p));} |
---|
418 | |
---|
419 | ideal T_bis_k,T1,T2; |
---|
420 | poly g,a,b; |
---|
421 | list rq; |
---|
422 | list gab; |
---|
423 | int i; |
---|
424 | int v = lvar(p); |
---|
425 | int k = nvars(basering)-v+1; |
---|
426 | |
---|
427 | // Zu p passende Dreiecksbasis auswaehlen. |
---|
428 | if (k == 1) |
---|
429 | { |
---|
430 | T_bis_k = 0; |
---|
431 | } |
---|
432 | else |
---|
433 | { |
---|
434 | T_bis_k = T[1..(k-1)]; |
---|
435 | } |
---|
436 | attrib(T_bis_k,"isSB",1); |
---|
437 | |
---|
438 | // Um p+(T) zu invertieren, erw. eukl. Algorithmus anwenden. |
---|
439 | gab = Erw_ggt_oT(T[k],p,v,T_bis_k); |
---|
440 | |
---|
441 | // Entweder Triangulierung von T_bis_k erhalten, |
---|
442 | if (zerlegt) |
---|
443 | { |
---|
444 | // T[k..size(T)] an gab[i] anhaengen und reduzieren. |
---|
445 | T1 = gab[1]; |
---|
446 | T2 = gab[2]; |
---|
447 | for (i = k; i <= size(T); i++) |
---|
448 | { |
---|
449 | T1 = T1 + cleardenom(reduce(T[i],T1)); |
---|
450 | T2 = T2 + cleardenom(reduce(T[i],T2)); |
---|
451 | attrib(T1,"isSB",1); |
---|
452 | attrib(T2,"isSB",1); |
---|
453 | } |
---|
454 | return(list(T1,T2)); |
---|
455 | } |
---|
456 | |
---|
457 | // ansonsten gilt: a*T[k] + b*p = g mod (T_bis_k) |
---|
458 | g = gab[1]; |
---|
459 | a = gab[2]; |
---|
460 | b = gab[3]; |
---|
461 | |
---|
462 | if (degv(g,v) > 0) // Entweder echten Teiler von T[k] gefunden, |
---|
463 | { // dann splitte T, |
---|
464 | rq = pdiv(T[k],g,v); |
---|
465 | |
---|
466 | T1 = T_bis_k; |
---|
467 | T1[k] = cleardenom(reduce(g,T1)); |
---|
468 | attrib(T1,"isSB",1); |
---|
469 | T2 = T_bis_k; |
---|
470 | T2[k] = cleardenom(reduce(rq[2],T2)); |
---|
471 | attrib(T2,"isSB",1); |
---|
472 | |
---|
473 | // T[k..size(T)] an T1, T2 anhaengen und reduzieren. |
---|
474 | for (i = k + 1; i <= size(T); i++) |
---|
475 | { |
---|
476 | T1 = T1 + cleardenom(reduce(T[i],T1)); |
---|
477 | T2 = T2 + cleardenom(reduce(T[i],T2)); |
---|
478 | attrib(T1,"isSB",1); |
---|
479 | attrib(T2,"isSB",1); |
---|
480 | } |
---|
481 | |
---|
482 | // zerlegt = TRUE, hier ist einzige Stelle, wo T zerlegt wird. |
---|
483 | zerlegt = 1; |
---|
484 | return(list(T1,T2)); |
---|
485 | } |
---|
486 | else // ansonsten Quasi-Inverses gefunden. |
---|
487 | { |
---|
488 | return(list(b,g)); |
---|
489 | } |
---|
490 | } |
---|
491 | |
---|
492 | /////////////////////////////////////////////////////////////////////////////// |
---|
493 | |
---|
494 | static proc Erw_ggt_oT (poly f, poly g, int v, ideal T) |
---|
495 | "USAGE: Erw_ggt_oT(f,g,v,T); f,g Polynome in var(v), f > g, |
---|
496 | T Dreiecksbasis mit groesster |
---|
497 | Variable var(v+1). |
---|
498 | RETURN: Liste aus drei Polynomen d,a,b, wobei |
---|
499 | a * f + b * g = d mod (T) und |
---|
500 | (d + (T)) = ((f + (T), g + (T)). |
---|
501 | NOTE: Dies ist der erweiterte euklidische Algorithmus im Polynom- |
---|
502 | ring ueber dem Restklassenring modulo (T) in der Unbestimmten |
---|
503 | var(v). Zusaetzlich wird jeweils der Dividend so normiert, |
---|
504 | dass der Leitkoeffizient bzgl. var(v) im Grundring ist |
---|
505 | und die Leikoeffizienten werden durch Anwendung der |
---|
506 | Subresultant-remainder-sequence Variante des euklidischen |
---|
507 | Algorithmus moeglichst klein gehalten. |
---|
508 | " |
---|
509 | { |
---|
510 | // Initialisierung fuer erw. eukl. Algorithmus. |
---|
511 | poly p1 = f; |
---|
512 | poly p2 = g; |
---|
513 | poly a1 = 1; |
---|
514 | poly a2 = 0; |
---|
515 | poly b1 = 0; |
---|
516 | poly b2 = 1; |
---|
517 | |
---|
518 | // Normierung des Dividenden. |
---|
519 | list pab = normieren_oT(p2,a2,b2,v,T); |
---|
520 | if (zerlegt) {return(pab);} |
---|
521 | p2 = pab[1]; a2 = pab[2]; b2 = pab[3]; |
---|
522 | |
---|
523 | poly p3,a3,b3,q; // Hilfspolynome |
---|
524 | list rq; // Rueckgabelisten der Fkten pdiv und normieren_oT |
---|
525 | vector clden; |
---|
526 | |
---|
527 | // Variablen fuer Subresultanten-Variante initialisieren. |
---|
528 | int delta = degv(p1,v)-degv(p2,v); |
---|
529 | number alpha; |
---|
530 | number beta = (-1)**(delta + 1); |
---|
531 | number psi = -1; |
---|
532 | number m; |
---|
533 | |
---|
534 | // Pseudodivision |
---|
535 | rq = pdiv(p1,p2,v); |
---|
536 | q = rq[2]; |
---|
537 | p3 = reduce(rq[1],T); |
---|
538 | |
---|
539 | while (p3 <> 0) // Hauptschleife des eukl. Algorithmus. |
---|
540 | { |
---|
541 | alpha = number(lcoef(p2,v)); // return von lcoef ist type poly |
---|
542 | m = alpha**(delta + 1); //m*p1 = q*p2 + p3 mod (T) |
---|
543 | p3 = p3 / beta; //m*p1 = q*p2 + beta * p3 mod (T) |
---|
544 | |
---|
545 | // a und b anpassen. |
---|
546 | q = reduce(q,T); |
---|
547 | a3 = reduce(m*a1 - q*a2,T); |
---|
548 | b3 = reduce(m*b1 - q*b2,T); |
---|
549 | a3 = a3 / beta; |
---|
550 | b3 = b3 / beta; |
---|
551 | |
---|
552 | // Variablen tauschen. |
---|
553 | p1 = p2; |
---|
554 | p2 = p3; |
---|
555 | a1 = a2; |
---|
556 | a2 = a3; |
---|
557 | b1 = b2; |
---|
558 | b2 = b3; |
---|
559 | |
---|
560 | // Normierung des Dividenden. |
---|
561 | pab = normieren_oT(p2,a2,b2,v,T); |
---|
562 | if (zerlegt) {return(pab);} |
---|
563 | p2 = pab[1]; a2 = pab[2]; b2 = pab[3]; |
---|
564 | |
---|
565 | if (degv(p2,v) > 0) // Dividend aus Grundring?, dann |
---|
566 | { // Subresultanten-Variabeln anpassen |
---|
567 | if (delta > 0) |
---|
568 | { |
---|
569 | psi = ((-alpha)**delta) / (psi**(delta-1)); |
---|
570 | } |
---|
571 | delta = degv(p1,v)-degv(p2,v); |
---|
572 | beta = (-alpha) * (psi**delta); |
---|
573 | |
---|
574 | rq = pdiv(p1,p2,v); // und Pseudodivision, |
---|
575 | q = rq[2]; |
---|
576 | p3 = reduce(rq[1], T); |
---|
577 | } |
---|
578 | else // ansonsten Ende. |
---|
579 | { |
---|
580 | p3 = 0; |
---|
581 | } |
---|
582 | } |
---|
583 | |
---|
584 | // Gemeinsame Teiler herausdividieren. |
---|
585 | clden = cleardenom([p2,a2,b2]); |
---|
586 | return(list(clden[1],clden[2],clden[3])); |
---|
587 | } |
---|
588 | |
---|
589 | /////////////////////////////////////////////////////////////////////////////// |
---|
590 | |
---|
591 | static proc normieren_oT (poly p, poly a, poly b, int v, ideal T) |
---|
592 | "USAGE: normieren_oT(p,a,b,v,T); p,a,b Polynome in var(v), |
---|
593 | T Dreiecksbasis. |
---|
594 | RETURN: Entweder Triangulierung [T1,T2] von (T) oder |
---|
595 | Liste,wobei [1] = reduce(r*p,T), |
---|
596 | [2] = reduce(r*a,T), |
---|
597 | [3] = reduce(r*b,T), |
---|
598 | und reduce(r*lcoef(p,v),T) ist im Grundring. |
---|
599 | " |
---|
600 | { |
---|
601 | poly lc = lcoef(p,v); |
---|
602 | if (deg(lc) == 0) {return(list(p,a,b));} // lc im Grundring |
---|
603 | lc = cleardenom(lc); |
---|
604 | |
---|
605 | // invertieren_oT aufrufen, um lc+(T) zu invertieren. |
---|
606 | list inv = invertieren_oT(lc,T); |
---|
607 | if (zerlegt) {return(inv);} |
---|
608 | |
---|
609 | // ansonsten Polynom zurueck gekommen. |
---|
610 | p = reduce(inv[1] * p,T); |
---|
611 | a = reduce(inv[1] * a,T); |
---|
612 | b = reduce(inv[1] * b,T); |
---|
613 | return(list(p,a,b)); |
---|
614 | } |
---|
615 | |
---|
616 | |
---|
617 | /////////////////////////////////////////////////////////////////////////////// |
---|
618 | // |
---|
619 | // Der Moeller-Algorithmus |
---|
620 | // |
---|
621 | /////////////////////////////////////////////////////////////////////////////// |
---|
622 | |
---|
623 | proc triangM (ideal G, list #) |
---|
624 | "USAGE: triangM(G[,i]); G=ideal, i=integer,@* |
---|
625 | ASSUME: G is the reduced lexicographical Groebner basis of the |
---|
626 | zero-dimensional ideal (G), sorted by increasing leading terms. |
---|
627 | RETURN: a list of finitely many triangular systems, such that |
---|
628 | the union of their varieties equals the variety of (G). |
---|
629 | If i = 2, then each polynomial of the triangular systems |
---|
630 | is factorized. |
---|
631 | NOTE: Algorithm of Moeller (see: Moeller, H.M.: |
---|
632 | On decomposing systems of polynomial equations with |
---|
633 | finitely many solutions, Appl. Algebra Eng. Commun. Comput. 4, |
---|
634 | 217 - 230, 1993). |
---|
635 | EXAMPLE: example triangM; shows an example |
---|
636 | " |
---|
637 | { |
---|
638 | // Test, ob G Groebnerbasis ist. |
---|
639 | if (attrib(G,"isSB") <> 1) |
---|
640 | { |
---|
641 | ERROR("The input is no groebner basis."); |
---|
642 | } |
---|
643 | |
---|
644 | // Faktorisierungsschalter setzen. |
---|
645 | int fak; |
---|
646 | if (size(#) > 0 && typeof(#[1]) == "int") |
---|
647 | { |
---|
648 | if (#[1] == 2) {fak = 1;} |
---|
649 | } |
---|
650 | |
---|
651 | // Noetige Optionen setzen. |
---|
652 | // obachman: save options so that they can be reset on exit |
---|
653 | intvec ovec = option(get); |
---|
654 | option(redSB); |
---|
655 | option(returnSB); |
---|
656 | |
---|
657 | list F,Fh; // Ergebnislisten |
---|
658 | int nG = size(G); |
---|
659 | if (nG <= 2) // G schon Dreiecksbasis? |
---|
660 | { |
---|
661 | if (fak) |
---|
662 | { |
---|
663 | F = faktorisiere_DB(G); |
---|
664 | } |
---|
665 | else |
---|
666 | { |
---|
667 | F = list(G); |
---|
668 | } |
---|
669 | option(set, ovec); |
---|
670 | return(F); |
---|
671 | } |
---|
672 | |
---|
673 | int u = lvar(G[nG-1]); |
---|
674 | int v = lvar(G[nG]); |
---|
675 | int m = nG; |
---|
676 | int l; |
---|
677 | ideal lc,G1,G2,H; |
---|
678 | poly lcr; |
---|
679 | int i,j,k; |
---|
680 | string s = string(nvars(basering)-v+1); // fuer Protokoll |
---|
681 | |
---|
682 | // Oberes Dreieckssytem abschneiden. |
---|
683 | while(u <> v) |
---|
684 | { |
---|
685 | m = m-1; |
---|
686 | if (m == 2) // G ist bereits Dreiecksbasis. |
---|
687 | { |
---|
688 | if (fak) |
---|
689 | { |
---|
690 | F = faktorisiere_DB(G); |
---|
691 | } |
---|
692 | else |
---|
693 | { |
---|
694 | F = list(G); |
---|
695 | } |
---|
696 | option(set, ovec); |
---|
697 | return (F); |
---|
698 | } |
---|
699 | v = u; |
---|
700 | u = lvar(G[m-1]); |
---|
701 | } |
---|
702 | |
---|
703 | // Leitkoeffizienten der Polynome bzgl. var(u) bestimmen. |
---|
704 | l = m-1; |
---|
705 | while(u == v) |
---|
706 | { |
---|
707 | lc = lc + lcoef(G[l],u); |
---|
708 | l = l - 1; |
---|
709 | u = lvar(G[l]); |
---|
710 | } |
---|
711 | |
---|
712 | // Erster Teil des Splitts. |
---|
713 | G1 = sort_red(ideal(G[1..l]) + lc); // Sortiere und reduziere GB. |
---|
714 | |
---|
715 | s = string(nvars(basering)-v+1); |
---|
716 | dbprint(string(timer)+" S"+s+": Erster Teil des Splitts."); |
---|
717 | |
---|
718 | // Rekursiver Aufruf, um G1 zu triangulieren. |
---|
719 | F = triangM(G1,#); |
---|
720 | |
---|
721 | dbprint(string(timer)+" S"+s+": Oberes Dreieckssystem anhaengen."); |
---|
722 | |
---|
723 | // Oberes Dreieckssystem an jede berechnete Dreiecksbasis anhaengen. |
---|
724 | for (j = m; j <= nG; j++) // meist nur j = m = nG |
---|
725 | { |
---|
726 | for (i = 1; i <= size(F); i++) |
---|
727 | { |
---|
728 | F[i] = F[i] + cleardenom(G[j][1] + reduce(G[j]-G[j][1],F[i])); |
---|
729 | attrib(F[i],"isSB",1); |
---|
730 | } |
---|
731 | if (fak) |
---|
732 | { |
---|
733 | F = faktorisiere_letzten(F); |
---|
734 | } |
---|
735 | } |
---|
736 | |
---|
737 | dbprint(string(timer)+" S"+s+": Zweiter Teil des Splitts, "+string(m-l-1)+" Leitkoeffizient(en)."); |
---|
738 | |
---|
739 | // Zweiten Teil des Splitts berechnen. |
---|
740 | H = ideal(G[1..m]); |
---|
741 | attrib(H,"isSB",1); |
---|
742 | |
---|
743 | for (k = m-l-1; k >= 1; k--) // Leitkoeffizientenliste durchgehen, |
---|
744 | { // angefangen mit lc des kleinsten Polynoms. |
---|
745 | lcr = reduce(lc[k],H); |
---|
746 | if (lcr <> 0) // d.h. (H):lcr <> (1) |
---|
747 | { |
---|
748 | dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Berechnung Idealquotient."); |
---|
749 | G2 = quotient(H,lcr); // Wg. Option returnSB schon Standardbasis. |
---|
750 | attrib(G2,"isSB",1); |
---|
751 | |
---|
752 | dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Rekursiver Aufruf."); |
---|
753 | // Rekursiver Aufruf, um G2 zu triangulieren. |
---|
754 | Fh = triangM(G2,#); |
---|
755 | |
---|
756 | dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Oberes Dreieckssystem, falls vorhanden, anhaengen."); |
---|
757 | // Oberes Dreieckssystem, falls vorhanden, |
---|
758 | // an jede berechnete Dreiecksbasis anhaengen. |
---|
759 | for (j = m + 1; j <= nG; j++) |
---|
760 | { |
---|
761 | for (i = 1; i <= size(Fh); i++) |
---|
762 | { |
---|
763 | Fh[i] = Fh[i] + cleardenom(G[j][1] + |
---|
764 | reduce(G[j]-G[j][1],Fh[i])); |
---|
765 | attrib(Fh[i],"isSB",1); |
---|
766 | } |
---|
767 | if (fak) |
---|
768 | { |
---|
769 | Fh = faktorisiere_letzten(Fh); |
---|
770 | } |
---|
771 | } |
---|
772 | F = F + Fh; |
---|
773 | if (k > 1) |
---|
774 | { |
---|
775 | dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Lex. GB von (H + lc) berechnen, naechsten lc reduzieren."); |
---|
776 | H = std(H + lcr); // Wg. reduce(...,H) oben notwendig. |
---|
777 | } |
---|
778 | } |
---|
779 | else |
---|
780 | { |
---|
781 | dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Leere Varietaet."); |
---|
782 | } |
---|
783 | } |
---|
784 | |
---|
785 | // Singular loescht die Attribute beim Listenverketten, deswegen neu setzen. |
---|
786 | for (i = 1; i <= size(F); i++) |
---|
787 | { |
---|
788 | attrib(F[i],"isSB",1); |
---|
789 | } |
---|
790 | option(set, ovec); |
---|
791 | return(F); |
---|
792 | } |
---|
793 | //-------------------------------- example---------------------------------- |
---|
794 | example |
---|
795 | { "EXAMPLE:"; echo = 2; |
---|
796 | ring rC5 = 0,(e,d,c,b,a),lp; |
---|
797 | triangM(stdfglm(cyclic(5))); //oder: triangM(stdfglm(cyclic(5)),2); |
---|
798 | } |
---|
799 | |
---|
800 | /////////////////////////////////////////////////////////////////////////////// |
---|
801 | |
---|
802 | proc triangMH (ideal G, list #) |
---|
803 | "USAGE: triangMH(G[,i]); G=ideal, i=integer |
---|
804 | ASSUME: G is the reduced lexicographical Groebner basis of the |
---|
805 | zero-dimensional ideal (G), sorted by increasing leading terms. |
---|
806 | RETURN: a list of finitely many triangular systems, such that |
---|
807 | the disjoint union of their varieties equals the variety of (G). |
---|
808 | If i = 2, then each polynomial of the triangular systems is factorized. |
---|
809 | NOTE: Algorithm of Moeller and Hillebrand (see: Moeller, H.M.: |
---|
810 | On decomposing systems of polynomial equations with finitely many |
---|
811 | solutions, Appl. Algebra Eng. Commun. Comput. 4, 217 - 230, 1993 and |
---|
812 | Hillebrand, D.: Triangulierung nulldimensionaler Ideale - |
---|
813 | Implementierung und Vergleich zweier Algorithmen, master thesis, |
---|
814 | Universitaet Dortmund, Fachbereich Mathematik, Prof. Dr. H.M. Moeller, |
---|
815 | 1999). |
---|
816 | EXAMPLE: example triangMH; shows an example |
---|
817 | " |
---|
818 | { |
---|
819 | // Test, ob G Groebnerbasis ist. |
---|
820 | if (attrib(G,"isSB") <> 1) |
---|
821 | { |
---|
822 | ERROR("The input is no groebner basis."); |
---|
823 | } |
---|
824 | if(npars(basering)<>0) |
---|
825 | { |
---|
826 | ERROR("basering has parameters"); |
---|
827 | } |
---|
828 | // Faktorisierungsschalter setzen. |
---|
829 | int fak; |
---|
830 | if (size(#) > 0 && typeof(#[1]) == "int") |
---|
831 | { |
---|
832 | if (#[1] == 2) {fak = 1;} |
---|
833 | } |
---|
834 | |
---|
835 | // Noetige Optionen setzen. |
---|
836 | // obachman: save options so that tey can be reset on exit |
---|
837 | intvec ovec = option(get); |
---|
838 | option(redSB); |
---|
839 | option(redTail); |
---|
840 | option(returnSB); |
---|
841 | |
---|
842 | list F,Fh; // Ergebnislisten |
---|
843 | int nG = size(G); |
---|
844 | if (nG <= 2) // G schon Dreiecksbasis? |
---|
845 | { |
---|
846 | if (fak) |
---|
847 | { |
---|
848 | F = faktorisiere_DB(G); |
---|
849 | } |
---|
850 | else |
---|
851 | { |
---|
852 | F = list(G); |
---|
853 | } |
---|
854 | option(set, ovec); |
---|
855 | return(F); |
---|
856 | } |
---|
857 | |
---|
858 | int u = lvar(G[nG-1]); |
---|
859 | int v = lvar(G[nG]); |
---|
860 | int m = nG; |
---|
861 | int l; |
---|
862 | poly lcr; |
---|
863 | ideal lc,G1,G2,H; |
---|
864 | int i,j,k; |
---|
865 | string s = string(nvars(basering)-v+1); // fuer Protokoll |
---|
866 | |
---|
867 | // Oberes Dreieckssytem abschneiden. |
---|
868 | while(u <> v) |
---|
869 | { |
---|
870 | m = m-1; |
---|
871 | if (m == 2) // G ist Dreiecksbasis. |
---|
872 | { |
---|
873 | if (fak) |
---|
874 | { |
---|
875 | F = faktorisiere_DB(G); |
---|
876 | } |
---|
877 | else |
---|
878 | { |
---|
879 | F = list(G); |
---|
880 | } |
---|
881 | option(set, ovec); |
---|
882 | return(F); |
---|
883 | } |
---|
884 | v = u; |
---|
885 | u = lvar(G[m-1]); |
---|
886 | } |
---|
887 | |
---|
888 | // Leitkoeffizienten der Polynome bzgl. var(u) bestimmen. |
---|
889 | l = m-1; |
---|
890 | while(u == v) |
---|
891 | { |
---|
892 | lc = lc + lcoef(G[l],u); |
---|
893 | l = l - 1; |
---|
894 | u = lvar(G[l]); |
---|
895 | } |
---|
896 | |
---|
897 | // Erster Teil des Splitts. |
---|
898 | G1 = sort_red(ideal(G[1..l]) + lc); // Sortiere und reduziere GB. |
---|
899 | |
---|
900 | s = string(nvars(basering)-v+1); |
---|
901 | dbprint(string(timer)+" S"+s+": Erster Teil des Splitts."); |
---|
902 | |
---|
903 | // Rekursiver Aufruf, um G1 zu triangulieren. |
---|
904 | F = triangMH(G1,#); |
---|
905 | |
---|
906 | dbprint(string(timer)+" S"+s+": Oberes Dreieckssystem anhaengen."); |
---|
907 | |
---|
908 | // Oberes Dreieckssystem an jede berechnete Dreiecksbasis anhaengen. |
---|
909 | for (j = m; j <= nG; j++) // meist nur j = m = nG |
---|
910 | { |
---|
911 | for (i = 1; i <= size(F); i++) |
---|
912 | { |
---|
913 | F[i] = F[i] + cleardenom(G[j][1] + reduce(G[j]-G[j][1],F[i])); |
---|
914 | attrib(F[i],"isSB",1); |
---|
915 | } |
---|
916 | if (fak) |
---|
917 | { |
---|
918 | F = faktorisiere_letzten(F); |
---|
919 | } |
---|
920 | } |
---|
921 | |
---|
922 | dbprint(string(timer)+" S"+s+": Zweiter Teil des Splitts, "+string(m-l-1)+" Leitkoeffizient(en)."); |
---|
923 | |
---|
924 | // Zweiten Teil des Splitts berechnen. |
---|
925 | H = ideal(G[1..l]); // Nur die Polynome in var(u). |
---|
926 | attrib(H,"isSB",1); |
---|
927 | |
---|
928 | for (k = m-l-1; k >= 1; k--) // Leitkoeffizientenliste durchgehen, |
---|
929 | { // angefangen mit lc des kleinsten Polynoms. |
---|
930 | lcr = reduce(lc[k],H); |
---|
931 | if (lcr <> 0) // d.h. (H):lcr <> (1) |
---|
932 | { |
---|
933 | dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Berechnung Saturation."); |
---|
934 | G2 = sat(H,lcr)[1]; // Saturation statt Idealquotient. |
---|
935 | attrib(G2,"isSB",1); |
---|
936 | |
---|
937 | if (G2[1] <> 1) |
---|
938 | { |
---|
939 | dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Rekursiver Aufruf."); |
---|
940 | // Rekursiver Aufruf, um G2 zu triangulieren. |
---|
941 | Fh = triangMH(G2,#); |
---|
942 | |
---|
943 | dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Ggt und oberes Dreieckssystem anhaengen."); |
---|
944 | // Ggt an jede Dreiecksbasis anhaengen. |
---|
945 | for (i = 1; i <= size(Fh); i++) |
---|
946 | { |
---|
947 | Fh[i] = std(Fh[i] + G[m-k]); |
---|
948 | } |
---|
949 | if (fak) |
---|
950 | { |
---|
951 | Fh = faktorisiere_letzten(Fh); |
---|
952 | } |
---|
953 | for (j = m+1; j <= nG; j++) // oberes Dreieckssystem an jede |
---|
954 | { // berechnete Dreiecksbasis anhaengen |
---|
955 | for (i = 1; i <= size(Fh); i++) |
---|
956 | { |
---|
957 | Fh[i] = Fh[i] + cleardenom(G[j][1] |
---|
958 | + reduce(G[j]-G[j][1],Fh[i])); |
---|
959 | attrib(Fh[i],"isSB",1); |
---|
960 | } |
---|
961 | if (fak) |
---|
962 | { |
---|
963 | Fh = faktorisiere_letzten(Fh); |
---|
964 | } |
---|
965 | } |
---|
966 | F = F + Fh; |
---|
967 | } |
---|
968 | else |
---|
969 | { |
---|
970 | dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Leere Varietaet (G2 == (1))."); |
---|
971 | } |
---|
972 | if (k > 1) |
---|
973 | { |
---|
974 | dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Lex. GB von (H + lc) berechnen, naechsten lc reduzieren."); |
---|
975 | H = std(H + lcr); // wegen reduce(...,H) oben |
---|
976 | } |
---|
977 | } |
---|
978 | else |
---|
979 | { |
---|
980 | dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Leere Varietaet (lcr == 0)."); |
---|
981 | } |
---|
982 | } |
---|
983 | |
---|
984 | // Singular loescht die Attribute beim Listenverketten, deswegen neu setzen. |
---|
985 | for (i = 1; i <= size(F); i++) |
---|
986 | { |
---|
987 | attrib(F[i],"isSB",1); |
---|
988 | } |
---|
989 | option(set, ovec); |
---|
990 | return(F); |
---|
991 | } |
---|
992 | //-------------------------------- example ---------------------------------- |
---|
993 | example |
---|
994 | { "EXAMPLE:"; echo = 2; |
---|
995 | ring rC5 = 0,(e,d,c,b,a),lp; |
---|
996 | triangMH(stdfglm(cyclic(5))); |
---|
997 | } |
---|
998 | |
---|
999 | |
---|
1000 | /////////////////////////////////////////////////////////////////////////////// |
---|
1001 | // |
---|
1002 | // Hilfsfunktionen |
---|
1003 | // |
---|
1004 | /////////////////////////////////////////////////////////////////////////////// |
---|
1005 | |
---|
1006 | static proc sort_red (ideal G) |
---|
1007 | "USAGE: sort_red(G); G lexikographische Groebnerbasis. |
---|
1008 | RETURN: reduzierte lex. GB G, deren Elemente nach aufsteigenden |
---|
1009 | Leittermen sortiert sind. |
---|
1010 | " |
---|
1011 | { |
---|
1012 | int i; |
---|
1013 | |
---|
1014 | // sortieren |
---|
1015 | G = sort(G)[1]; |
---|
1016 | |
---|
1017 | // reduzieren |
---|
1018 | ideal Gred = cleardenom(G[1]); |
---|
1019 | attrib(Gred,"isSB",1); |
---|
1020 | |
---|
1021 | for (i = 2; i <= size(G); i++) |
---|
1022 | { |
---|
1023 | Gred = Gred + cleardenom(reduce(G[i],Gred)); |
---|
1024 | attrib(Gred,"isSB",1); |
---|
1025 | } |
---|
1026 | attrib(Gred,"isSB",1); |
---|
1027 | return(Gred); |
---|
1028 | } |
---|
1029 | |
---|
1030 | /////////////////////////////////////////////////////////////////////////////// |
---|
1031 | |
---|
1032 | static proc Listenrest (list Liste) |
---|
1033 | "USAGE: Listenrest(Liste); Liste muss eine mindestens einelementige |
---|
1034 | Liste sein (leere Liste ergibt Fehler). |
---|
1035 | RETURN: Liste ohne das erste Element. |
---|
1036 | " |
---|
1037 | { |
---|
1038 | return(delete(Liste,1)); |
---|
1039 | } |
---|
1040 | |
---|
1041 | /////////////////////////////////////////////////////////////////////////////// |
---|
1042 | |
---|
1043 | static proc Idealrest (ideal i) |
---|
1044 | "USAGE: Idealrest(i); i Ideal, i <> (0). |
---|
1045 | RETURN: Ideal i ohne das erste Element bzw das Nullideal, |
---|
1046 | falls i nur ein Erzeugendes besass. |
---|
1047 | " |
---|
1048 | { |
---|
1049 | int ni = size(i); |
---|
1050 | if (ni == 1) |
---|
1051 | { |
---|
1052 | return(ideal(0)); // Nullideal |
---|
1053 | } |
---|
1054 | return(ideal(i[2..ni])); |
---|
1055 | } |
---|
1056 | |
---|
1057 | /////////////////////////////////////////////////////////////////////////////// |
---|
1058 | |
---|
1059 | static proc H_anhaengen (list L, int h) |
---|
1060 | "USAGE: H_anhaengen(L,h); L Liste aus Idealen, h natuerliche Zahl. |
---|
1061 | RETURN: Macht aus dem Listenelement i ein Paar [i,h] und setzt |
---|
1062 | Attribut "isSB" fuer i. |
---|
1063 | " |
---|
1064 | { |
---|
1065 | int i; |
---|
1066 | list E; // Ergebnisliste |
---|
1067 | |
---|
1068 | for (i = 1; i <= size(L); i++) |
---|
1069 | { |
---|
1070 | E = E + list(list(L[i],h)); |
---|
1071 | attrib(E[size(E)][1],"isSB",1); |
---|
1072 | } |
---|
1073 | return(E); |
---|
1074 | } |
---|
1075 | |
---|
1076 | /////////////////////////////////////////////////////////////////////////////// |
---|
1077 | |
---|
1078 | static proc faktorisiere_letzten (list L) |
---|
1079 | "USAGE: faktorisiere_letzten(L); L Liste aus Idealen. |
---|
1080 | RETURN: Faktorisiert letztes Polynom aus jedem Ideal der Liste und |
---|
1081 | zerlegt das Ideal entsprechend. Attribut "isSB" wird gesetzt. |
---|
1082 | " |
---|
1083 | { |
---|
1084 | int i,j; |
---|
1085 | ideal h; |
---|
1086 | int nh; |
---|
1087 | list factors; |
---|
1088 | list E; // Ergebnisliste |
---|
1089 | for (i = 1; i <= size(L); i++) |
---|
1090 | { |
---|
1091 | h = L[i]; |
---|
1092 | nh = size(h); |
---|
1093 | factors = factorize(h[nh],1); |
---|
1094 | for (j = 1; j <= size(factors[1]); j++) |
---|
1095 | { |
---|
1096 | h[nh] = factors[1][j]; |
---|
1097 | E = insert(E,h,size(E)); |
---|
1098 | attrib(E[size(E)],"isSB",1); |
---|
1099 | } |
---|
1100 | } |
---|
1101 | return(E); |
---|
1102 | } |
---|
1103 | |
---|
1104 | /////////////////////////////////////////////////////////////////////////////// |
---|
1105 | |
---|
1106 | static proc faktorisiere_DB (ideal db) |
---|
1107 | "USAGE: faktorisiere_DB(db); db reduzierte Dreiecksbasis. |
---|
1108 | RETURN: Liste aus reduzierten Dreiecksbasen, |
---|
1109 | die ensteht, indem, mit dem ersten beginnend, |
---|
1110 | jedes Polynom der Dreiecksbasis faktorisiert wird. |
---|
1111 | " |
---|
1112 | { |
---|
1113 | int i,j; |
---|
1114 | poly h; |
---|
1115 | list E = faktorisiere_letzten(db[1]); |
---|
1116 | |
---|
1117 | for (j = 2; j <= size(db); j++) |
---|
1118 | { |
---|
1119 | for (i = 1; i <= size(E); i++) |
---|
1120 | { |
---|
1121 | h = db[j][1] + reduce(db[j]-db[j][1],E[i]); |
---|
1122 | E[i] = E[i] + h; |
---|
1123 | } |
---|
1124 | E = faktorisiere_letzten(E); |
---|
1125 | } |
---|
1126 | return(E); |
---|
1127 | } |
---|
1128 | |
---|
1129 | /////////////////////////////////////////////////////////////////////////////// |
---|
1130 | |
---|
1131 | static proc degv (poly f, int v) |
---|
1132 | "USAGE: degv(f,v); f Polynom in var(v). |
---|
1133 | RETURN: Liefert den Grad von f bzgl. var(v) bzw -1, falls f == 0. |
---|
1134 | " |
---|
1135 | { |
---|
1136 | if (f == 0) {return(-1);} |
---|
1137 | return(size(coeffs(f,var(v))) - 1); |
---|
1138 | } |
---|
1139 | |
---|
1140 | ////////////////////////////////////////////////////////////////////////////// |
---|
1141 | |
---|
1142 | static proc pdiv (poly f, poly g, int v) |
---|
1143 | "USAGE: pdiv(f,g); f,g Polynome in var(v), f > g. |
---|
1144 | RETURN: Liste, wobei [1] = prem(f,g) (Pseudoremainder), |
---|
1145 | [2] = pquo(f,g) (Pseudoquotient). |
---|
1146 | NOTE: Pseudodivision von f durch g. |
---|
1147 | " |
---|
1148 | { |
---|
1149 | // Initialisierung |
---|
1150 | poly q = 0; |
---|
1151 | poly rem = f; |
---|
1152 | int deg_g = degv(g,v); |
---|
1153 | int deg_rem = degv(rem,v); |
---|
1154 | poly lcg = lcoef(g,v); // Im Lazard-Algor. aus dem Grundring. |
---|
1155 | poly h; |
---|
1156 | rem = rem*lcg**(deg_rem - deg_g + 1); // f durch g teilbar machen. |
---|
1157 | |
---|
1158 | // Gewoehnliche Polynomdivision |
---|
1159 | while (rem <> 0 and deg_rem >= deg_g) |
---|
1160 | { |
---|
1161 | h = (lcoef(rem,v)/lcg)*var(v)**(deg_rem - deg_g); |
---|
1162 | q = q + h; |
---|
1163 | rem = rem - h*g; |
---|
1164 | deg_rem = degv(rem,v); |
---|
1165 | } |
---|
1166 | return(list(rem,q)); |
---|
1167 | } |
---|
1168 | |
---|
1169 | ////////////////////////////////////////////////////////////////////////////// |
---|
1170 | |
---|
1171 | static proc lvar (poly f) |
---|
1172 | "USAGE: lvar(f); f nicht-konstantes Polynom. |
---|
1173 | RETURN: groesste Variable von f bzgl. der aktuellen Termordnung. |
---|
1174 | " |
---|
1175 | { |
---|
1176 | int i = 1; |
---|
1177 | intvec exp = leadexp(f); |
---|
1178 | |
---|
1179 | while (1) |
---|
1180 | { |
---|
1181 | if (exp[i] <> 0) {return(i);} |
---|
1182 | i++; |
---|
1183 | } |
---|
1184 | } |
---|
1185 | |
---|
1186 | ////////////////////////////////////////////////////////////////////////////// |
---|
1187 | |
---|
1188 | static proc lcoef (poly f, int v) |
---|
1189 | "USAGE: lcoef(f,v); f Polynom, var(v) Ringvariable. |
---|
1190 | RETURN: Leitkoeffizient von f bzgl. var(v). |
---|
1191 | " |
---|
1192 | { |
---|
1193 | matrix m = coeffs(f,var(v)); |
---|
1194 | return(m[size(m),1]); |
---|
1195 | } |
---|
1196 | /////////////////////////////////////////////////////////////////////////////// |
---|
1197 | |
---|