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