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