source: git/Singular/LIB/triang.lib @ c27111a

spielwiese
Last change on this file since c27111a was c27111a, checked in by Hans Schönemann <hannes@…>, 16 years ago
*hannes: 7bit clean git-svn-id: file:///usr/local/Singular/svn/trunk@11086 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 35.6 KB
Line 
1//last change: 13.02.2001 (Eric Westenberger)
2//////////////////////////////////////////////////////////////////////////////
3version="$Id: triang.lib,v 1.12 2008-10-01 14:20:22 Singular Exp $";
4category="Symbolic-numerical solving";
5info="
6LIBRARY: triang.lib   Decompose Zero-dimensional Ideals into Triangular Sets
7AUTHOR: D. Hillebrand
8
9PROCEDURES:
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
16LIB "general.lib"; // Muss geladen werden fuer den Befehl sort().
17LIB "elim.lib";    // Muss geladen werden fuer sat().
18
19
20///////////////////////////////////////////////////////////////////////////////
21//
22//  Der Lazard-Algorithmus
23//
24///////////////////////////////////////////////////////////////////////////////
25
26proc triangL (ideal G)
27"USAGE:   triangL(G); G=ideal
28ASSUME:  G is the reduced lexicographical Groebner basis of the
29         zero-dimensional ideal (G), sorted by increasing leading terms.
30RETURN:  a list of finitely many triangular systems, such that
31         the union of their varieties equals the variety of (G).
32NOTE:    Algorithm of Lazard (see: Lazard, D.: Solving zero-dimensional
33         algebraic systems, J. Symb. Comp. 13, 117 - 132, 1992).
34EXAMPLE: 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 ----------------------------------
166example
167{  "EXAMPLE:"; echo = 2;
168    ring rC5 = 0,(e,d,c,b,a),lp;
169    triangL(stdfglm(cyclic(5)));
170}
171
172///////////////////////////////////////////////////////////////////////////////
173
174proc triangLfak (ideal G)
175"USAGE:   triangLfak(G); G=ideal
176ASSUME:  G is the reduced lexicographical Groebner basis of the
177         zero-dimensional ideal (G), sorted by increasing leading terms.
178RETURN:  a list of finitely many triangular systems, such that
179         the union of their varieties equals the variety of (G).
180NOTE:    Algorithm of Lazard with factorization (see: Lazard, D.: Solving
181         zero-dimensional algebraic systems, J. Symb. Comp. 13, 117 - 132, 1992).
182REMARK:  each polynomial of the triangular systems is factorized.
183EXAMPLE: example triangLfak; shows an example
184"
185{
186    return(triangLbas(G,2));
187}
188//-------------------------------- example ----------------------------------
189example
190{  "EXAMPLE:"; echo = 2;
191    ring rC5 = 0,(e,d,c,b,a),lp;
192    triangLfak(stdfglm(cyclic(5)));
193}
194
195///////////////////////////////////////////////////////////////////////////////
196
197static 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).
202RETURN:  Triangulierung von (G).
203         Ist i == 2, dann wird jedes Polynom der Dreiecksbasen
204         der Triangulierung zusaetzlich faktorisiert.
205NOTE:    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
337static proc invertieren (poly p, ideal T)
338"USAGE:   invertieren(p,T); p Polynom, T reduzierte Dreiecksbasis.
339RETURN:  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
407static proc invertieren_oT (poly p, ideal T)
408"USAGE:   invertieren_oT(p,T); T reduzierte Dreiecksbasis,
409                              p <> 0 Polynom, irreduzibel modulo (T).
410RETURN:  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
494static 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).
498RETURN:  Liste aus drei Polynomen d,a,b, wobei
499         a * f + b * g = d mod (T) und
500         (d + (T)) = ((f + (T), g + (T)).
501NOTE:    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
591static 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.
594RETURN:  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
623proc triangM (ideal G, list #)
624"USAGE:   triangM(G[,i]); G=ideal, i=integer,@*
625ASSUME:  G is the reduced lexicographical Groebner basis of the
626         zero-dimensional ideal (G), sorted by increasing leading terms.
627RETURN:  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.
631NOTE:    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).
635EXAMPLE: 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----------------------------------
794example
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
802proc triangMH (ideal G, list #)
803"USAGE:   triangMH(G[,i]); G=ideal, i=integer
804ASSUME:  G is the reduced lexicographical Groebner basis of the
805         zero-dimensional ideal (G), sorted by increasing leading terms.
806RETURN:  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.
809NOTE:    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).
816EXAMPLE: 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 ----------------------------------
993example
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
1006static proc sort_red (ideal G)
1007"USAGE:   sort_red(G);  G lexikographische Groebnerbasis.
1008RETURN:  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
1032static proc Listenrest (list Liste)
1033"USAGE:   Listenrest(Liste); Liste muss eine mindestens einelementige
1034                            Liste sein (leere Liste ergibt Fehler).
1035RETURN:  Liste ohne das erste Element.
1036"
1037{
1038    return(delete(Liste,1));
1039}
1040
1041///////////////////////////////////////////////////////////////////////////////
1042
1043static proc Idealrest (ideal i)
1044"USAGE:   Idealrest(i); i Ideal, i <> (0).
1045RETURN:  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
1059static proc H_anhaengen (list L, int h)
1060"USAGE:   H_anhaengen(L,h); L Liste aus Idealen, h natuerliche Zahl.
1061RETURN:  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
1078static proc faktorisiere_letzten (list L)
1079"USAGE:   faktorisiere_letzten(L); L Liste aus Idealen.
1080RETURN:  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
1106static proc faktorisiere_DB (ideal db)
1107"USAGE:   faktorisiere_DB(db); db reduzierte Dreiecksbasis.
1108RETURN:  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
1131static proc degv (poly f, int v)
1132"USAGE:   degv(f,v);  f Polynom in var(v).
1133RETURN:  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
1142static proc pdiv (poly f, poly g, int v)
1143"USAGE:   pdiv(f,g);  f,g Polynome in var(v), f > g.
1144RETURN:  Liste, wobei [1] = prem(f,g) (Pseudoremainder),
1145                      [2] = pquo(f,g) (Pseudoquotient).
1146NOTE:    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
1171static proc lvar (poly f)
1172"USAGE:   lvar(f);  f nicht-konstantes Polynom.
1173RETURN:  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
1188static proc lcoef (poly f, int v)
1189"USAGE:   lcoef(f,v);  f Polynom, var(v) Ringvariable.
1190RETURN:  Leitkoeffizient von f bzgl. var(v).
1191"
1192{
1193    matrix m = coeffs(f,var(v));
1194    return(m[size(m),1]);
1195}
1196///////////////////////////////////////////////////////////////////////////////
1197
Note: See TracBrowser for help on using the repository browser.