//////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Symbolic-numerical solving";
info="
LIBRARY: triang.lib   Decompose Zero-dimensional Ideals into Triangular Sets
AUTHOR: D. Hillebrand

PROCEDURES:
 triangL(G);          Decomposition of (G) into triangular systems (Lazard).
 triangLfak(G);       Decomp. of (G) into tri. systems plus factorization.
 triangM(G[,.]);      Decomposition of (G) into triangular systems (Moeller).
 triangMH(G[,.]);     Decomp. of (G) into tri. syst. with disjoint varieties.
";

LIB "general.lib"; // Muss geladen werden fuer den Befehl sort().
LIB "elim.lib";    // Muss geladen werden fuer sat().


///////////////////////////////////////////////////////////////////////////////
//
//  Der Lazard-Algorithmus
//
///////////////////////////////////////////////////////////////////////////////

proc triangL (ideal G)
"USAGE:   triangL(G); G=ideal
ASSUME:  G is the reduced lexicographical Groebner basis of the
         zero-dimensional ideal (G), sorted by increasing leading terms.
RETURN:  a list of finitely many triangular systems, such that
         the union of their varieties equals the variety of (G).
NOTE:    Algorithm of Lazard (see: Lazard, D.: Solving zero-dimensional
         algebraic systems, J. Symb. Comp. 13, 117 - 132, 1992).
EXAMPLE: example triangL; shows an example
"
{

    // Test, ob G Groebnerbasis eines nulldimensionalen Ideals ist.
    if (attrib(G,"isSB") <> 1)
    {
        ERROR("The input is no groebner basis.");
    }
    else
    {
        if (dim(G) <> 0)
        {
            ERROR("ideal is not zero-dimensional.");
        }
    }

    if (size(G) <= 2)  // schon Dreiecksbasis
    {
        return(list(G));
    }

    // Noetige Optionen setzen.
    // obachman: save options so that tey can be reset on exit
    intvec ovec = option(get);
    option(redSB);
    option(returnSB);
    list B;    // Bearbeitungsliste
    list F;    // Ergebnisliste = Triangulierung
    ideal T, T_neu;
    poly p,p_neu,r,s,lt,lc,lcr;
    list inv;
    int u,v;
    int l,h;
    string st;

    T = ideal(G[1]);
    attrib(T,"isSB",1);
    B = list(list(T,2));

    while (size(B) > 0)
    {
        T   = B[1][1];
        h   = B[1][2];
        B   = Listenrest(B);
        p   = G[h];
        v   = lvar(p);
        lc  = lcoef(p,v);
          st=string(nvars(basering)-v+1);
          dbprint(string(timer)+" S"+st+": invertiere Leitkoeffizient von G["+string(h)+"].");

        // Invertiere den Leitkoeffizienten von G[h] bzgl. var(v) modulo (T).
        inv = invertieren(lc,T);

          dbprint(string(timer)+" S"+st+": Anzahl Dreiecksbasen: "+string(size(inv))+".");

        while (size(inv) > 0)
        {
            r   = inv[1][1];
            s   = inv[1][2];
            T   = inv[1][3];
            inv = Listenrest(inv);
            if (r == 0)  // r == 0?, dann p nicht der ggt der Stufe,
            {
                   dbprint(string(timer)+" S"+st+": Leitkoeffizient == 0.");
                B = list(list(T,h+1)) + B;
            }
            else         // ansonsten ggt der Stufe gefunden.
            {
                   dbprint(string(timer)+" S"+st+": ggt gefunden.");
                lt    = var(v)**degv(p,v);
                p_neu = cleardenom(s*lt + reduce(r*(p-lc*lt),T));
                T_neu = T,p_neu;
                attrib(T_neu,"isSB",1);

                // Restlichen Polynome der gleichen Stufe uebergehen.
                for (l = h+1; l <= size(G); l++)
                {
                    u = lvar(G[l]);
                    if (u <> v) {break;}
                }

                // Durchsuche restliche Stufen nach trivialen ggt's.
                while (l <= size(G))
                {
                    lc = lcoef(G[l],u);
                    lcr = reduce(lc,T_neu);

                    while (lcr == 0) // Gehe bis zum ersten Polynom <> 0
                    {                // modulo (T_neu).
                        l++;
                        lc = lcoef(G[l],u);
                        lcr = reduce(lc,T_neu);
                    }

                    if (deg(lcr) == 0) // Polynom <> 0 normiert?, dann
                    {                  // an Dreiecksbasis anhaengen,
                          dbprint(string(timer)+" S"+string(nvars(basering)-u+1)+": Schnellerkennung: ggt gefunden.");
                        lt = var(u)**degv(G[l],u);
                        p_neu = cleardenom(lcr*lt +
                                  reduce(G[l]-lc*lt,T_neu));
                        T_neu = T_neu,p_neu;
                        attrib(T_neu,"isSB",1);
                        u = lvar(G[l]);
                        l++;
                        while (l <= size(G)) // Rest der Stufe uebergehen.
                        {
                            if (lvar(G[l]) <> u) {u = lvar(G[l]); break;}
                            l++;
                        }
                    }
                    else // sonst nicht normierte Polynome auf der Stufe.
                    {
                        break;
                    }
                }
                if (l > size(G))  // Ende von G erreicht, dann Dreiecks-
                {                 // basis der Triangulierung gefunden,
                       dbprint(string(timer)+" S"+st+": Dreiecksbasis der Triangulierung gefunden.");
                    F = F + list(T_neu);
                }
                else              // sonst T_neu,l erneut bearbeiten.
                {
                    B = list(list(T_neu,l)) + B;
                }
            }
        }
    }
    option(set, ovec);
    return(F);
}
//-------------------------------- example ----------------------------------
example
{  "EXAMPLE:"; echo = 2;
    ring rC5 = 0,(e,d,c,b,a),lp;
    triangL(stdfglm(cyclic(5)));
}

///////////////////////////////////////////////////////////////////////////////

proc triangLfak (ideal G)
"USAGE:   triangLfak(G); G=ideal
ASSUME:  G is the reduced lexicographical Groebner basis of the
         zero-dimensional ideal (G), sorted by increasing leading terms.
RETURN:  a list of finitely many triangular systems, such that
         the union of their varieties equals the variety of (G).
NOTE:    Algorithm of Lazard with factorization (see: Lazard, D.: Solving
         zero-dimensional algebraic systems, J. Symb. Comp. 13, 117 - 132, 1992).
REMARK:  each polynomial of the triangular systems is factorized.
EXAMPLE: example triangLfak; shows an example
"
{
    return(triangLbas(G,2));
}
//-------------------------------- example ----------------------------------
example
{  "EXAMPLE:"; echo = 2;
    ring rC5 = 0,(e,d,c,b,a),lp;
    triangLfak(stdfglm(cyclic(5)));
}

///////////////////////////////////////////////////////////////////////////////

static proc triangLbas (ideal G, list #)
"USAGE:   triangLbas(G[,i]); G reduzierte lexikographische Groebnerbasis
                            des nulldimensionalen Ideals (G), nach
                            Leittermen aufsteigend sortiert,
                            i = 1 oder 2 (ohne Angabe i = 1).
RETURN:  Triangulierung von (G).
         Ist i == 2, dann wird jedes Polynom der Dreiecksbasen
         der Triangulierung zusaetzlich faktorisiert.
NOTE:    Algorithmus von Lazard (siehe: Lazard, D.: Solving
         zero-dimensional algebraic systems,
         J. Symb. Comp., Bd. 13, S. 117 - 132, 1992).
"
{

    // Test, ob G Groebnerbasis eines nulldimensionalen Ideals ist.
    if (attrib(G,"isSB") <> 1)
    {
        ERROR("The input is no groebner basis.");
    }
    else
    {
        if (dim(G) <> 0)
        {
            ERROR("ideal is not zero-dimensional.");
        }
    }

    // Noetige Optionen setzen.
    // obachman: save options so that tey can be reset on exit
    intvec ovec = option(get);
    option(redSB);
    option(returnSB);
    // Faktorisierungsschalter setzen.
    int fak;
    if (size(#) > 0 && typeof(#[1]) == "int")
    {
        if (#[1] == 2) {fak = 1;}
    }

    if (size(G) <= 2)  // schon Dreiecksbasis
    {
        list E;
        if (fak)
        {
            E = faktorisiere_DB(G);
        }
        else
        {
            E = list(G);
        }
        option(set, ovec);
        return(E);
    }

    list B;    // Bearbeitungsliste
    list F;    // Ergebnisliste = Triangulierung
    list T_neu;
    ideal T;
    poly p,p_neu,r,s,lt,lc;
    list inv;
    int v;
    int l,h;
    string st;

    // B initialisieren
    if (fak)
    {
        B = H_anhaengen(faktorisiere_letzten(G[1]),2);
    }
    else
    {
        T = ideal(G[1]);
        attrib(T,"isSB",1);
        B = list(list(T,2));
    }

    while (size(B) > 0)
    {
        T   = B[1][1];
        h   = B[1][2];
        B   = Listenrest(B);
        p   = G[h];
        v   = lvar(p);
        lc  = lcoef(p,v);
          st=string(nvars(basering)-v+1);
          dbprint(string(timer)+" S"+st+": invertiere Leitkoeffizient von G["+string(h)+"].");

        // invertiere den Leitkoeffizienten von G[h] bzgl. var(v) modulo (T).
        inv = invertieren(lc,T);
          dbprint(string(timer)+" S"+st+": Anzahl Dreiecksbasen: "+string(size(inv))+".");

        while (size(inv) > 0)
        {
            r   = inv[1][1];
            s   = inv[1][2];
            T   = inv[1][3];
            inv = Listenrest(inv);
            if (r == 0)  // r == 0?, dann p nicht der ggt der Stufe,
            {
                   dbprint(string(timer)+" S"+st+": Leitkoeffizient == 0.");
                B = list(list(T,h+1)) + B;
            }
            else         // ansonsten ggt der Stufe gefunden.
            {
                   dbprint(string(timer)+" S"+st+": ggt gefunden.");
                lt    = var(v)**degv(p,v);
                p_neu = cleardenom(s*lt + reduce(r*(p-lc*lt),T));
                if (fak)
                {
                    T_neu = faktorisiere_letzten(list(T + p_neu));
                }
                else
                {
                    T_neu = list(T + p_neu);
                }

                // Restlichen Polynome der gleichen Stufe uebergehen.
                for (l = h+1; l <= size(G); l++)
                {
                     if (lvar(G[l]) <> v) {break;}
                }

                if (l > size(G))  // Ende von G erreicht, dann Dreiecks-
                {                 // basis der Triangulierung gefunden,
                       dbprint(string(timer)+" S"+st+": Dreiecksbasis der Triangulierung gefunden.");
                    F = F + T_neu;
                }
                else              // sonst T_neu,l erneut bearbeiten.
                {
                    B  = H_anhaengen(T_neu,l) + B;
                }
            }
        }
    }
    option(set, ovec);
    return(F);
}

///////////////////////////////////////////////////////////////////////////////

static proc invertieren (poly p, ideal T)
"USAGE:   invertieren(p,T); p Polynom, T reduzierte Dreiecksbasis.
RETURN:  Liste von Tripeln (Listen) bestehend aus einem Polynom ri,
         einem Element aus dem Grundring si, und einer
         reduzierten Dreiecksbasis Ti, i = 1,...,m, so dass
         [T1,...,Tm] eine Triangulierung von T ist und
         p = 0 mod (Ti) falls ri = 0, ansonsten
         ri*p = si mod (Ti) fuer i = 1,...,m.
"
{
    // Triviale Faelle vorweg behandeln.
    p = reduce(p,T);
    if (p == 0) {return(list(list(0,0,T)));}
    if (deg(p) == 0) {return(list(list(1,p,T)));}

    list inv;
    int zerlegt;      // zerlegt ist Schalter dafuer, ob T zerlegt wurde,
    export zerlegt;   // einzige globale Variable.
    list toSee = list(list(p,T));  // zu bearbeitende Paare
    list toSave;                   // Ergebnisliste
    poly pr;

    while (size(toSee) > 0)
    {
        zerlegt = 0;             // zerlegt = FALSE
        p = toSee[1][1];
        T = toSee[1][2];
        toSee = Listenrest(toSee);

        // invertieren_oT aufrufen, um p+(T) zu invertieren.
        inv = invertieren_oT(p,T);

        if (zerlegt)   // T zerlegt?, dann pro neuer Dreiecksbasis
        {              // weiteruntersuchen,
            pr = reduce(p,inv[1]);
            if (pr == 0) // es sei denn p ist reduziert gleich 0,
            {
                toSave = list(list(0,0,inv[1])) + toSave;
                attrib(toSave[1][3],"isSB",1);
            }
            else
            {
                toSee  = list(list(pr,inv[1])) + toSee;
                attrib(toSee[1][2],"isSB",1);
            }
            pr = reduce(p,inv[2]);
            if (pr == 0)
            {
                toSave = list(list(0,0,inv[2])) + toSave;
                attrib(toSave[1][3],"isSB",1);
            }
            else
            {
                toSee  = list(list(pr,inv[2])) + toSee;
                attrib(toSee[1][2],"isSB",1);
            }

        }
        else           // ansonsten Quasi-Inverses gefunden.
        {
            toSave = list(list(inv[1],inv[2],T)) + toSave;
            attrib(toSave[1][3],"isSB",1);
        }
    }
    kill zerlegt;
    return(toSave);
}

///////////////////////////////////////////////////////////////////////////////

static proc invertieren_oT (poly p, ideal T)
"USAGE:   invertieren_oT(p,T); T reduzierte Dreiecksbasis,
                              p <> 0 Polynom, irreduzibel modulo (T).
RETURN:  Entweder ein Quasi-Inverses (r,s) von p modulo (T), d.h.
         r*p = s mod (T) und s im Grundring,
         oder eine Triangulierung [T1,T2] von T
         mit (Ti) <> (T) fuer i = 1,2.
"
{
    // Quasi-Inverses von Konstante klar.
    if (deg(p) == 0) {return(list(1,p));}

    ideal T_bis_k,T1,T2;
    poly g,a,b;
    list rq;
    list gab;
    int i;
    int v = lvar(p);
    int k = nvars(basering)-v+1;

    // Zu p passende Dreiecksbasis auswaehlen.
    if (k == 1)
    {
        T_bis_k = 0;
    }
    else
    {
        T_bis_k = T[1..(k-1)];
    }
    attrib(T_bis_k,"isSB",1);

    // Um p+(T) zu invertieren, erw. eukl. Algorithmus anwenden.
    gab = Erw_ggt_oT(T[k],p,v,T_bis_k);

    // Entweder Triangulierung von T_bis_k erhalten,
    if (zerlegt)
    {
        // T[k..size(T)] an gab[i] anhaengen und reduzieren.
        T1 = gab[1];
        T2 = gab[2];
        for (i = k; i <= size(T); i++)
        {
            T1 = T1 + cleardenom(reduce(T[i],T1));
            T2 = T2 + cleardenom(reduce(T[i],T2));
            attrib(T1,"isSB",1);
            attrib(T2,"isSB",1);
        }
        return(list(T1,T2));
    }

    // ansonsten gilt: a*T[k] + b*p = g mod (T_bis_k)
    g = gab[1];
    a = gab[2];
    b = gab[3];

    if (degv(g,v) > 0)        // Entweder echten Teiler von T[k] gefunden,
    {                         // dann splitte T,
        rq = pdiv(T[k],g,v);

        T1 = T_bis_k;
        T1[k] = cleardenom(reduce(g,T1));
        attrib(T1,"isSB",1);
        T2 = T_bis_k;
        T2[k] = cleardenom(reduce(rq[2],T2));
        attrib(T2,"isSB",1);

        // T[k..size(T)] an T1, T2 anhaengen und reduzieren.
        for (i = k + 1; i <= size(T); i++)
        {
            T1 = T1 + cleardenom(reduce(T[i],T1));
            T2 = T2 + cleardenom(reduce(T[i],T2));
            attrib(T1,"isSB",1);
            attrib(T2,"isSB",1);
        }

        // zerlegt = TRUE, hier ist einzige Stelle, wo T zerlegt wird.
        zerlegt = 1;
        return(list(T1,T2));
    }
    else                      // ansonsten Quasi-Inverses gefunden.
    {
        return(list(b,g));
    }
}

///////////////////////////////////////////////////////////////////////////////

static proc Erw_ggt_oT (poly f, poly g, int v, ideal T)
"USAGE:   Erw_ggt_oT(f,g,v,T);  f,g Polynome in var(v), f > g,
                               T Dreiecksbasis mit groesster
                               Variable var(v+1).
RETURN:  Liste aus drei Polynomen d,a,b, wobei
         a * f + b * g = d mod (T) und
         (d + (T)) = ((f + (T), g + (T)).
NOTE:    Dies ist der erweiterte euklidische Algorithmus im Polynom-
         ring ueber dem Restklassenring modulo (T) in der Unbestimmten
         var(v). Zusaetzlich wird jeweils der Dividend so normiert,
         dass der Leitkoeffizient bzgl. var(v) im Grundring ist
         und die Leikoeffizienten werden durch Anwendung der
         Subresultant-remainder-sequence Variante des euklidischen
         Algorithmus moeglichst klein gehalten.
"
{
    // Initialisierung fuer erw. eukl. Algorithmus.
    poly p1 = f;
    poly p2 = g;
    poly a1 = 1;
    poly a2 = 0;
    poly b1 = 0;
    poly b2 = 1;

    // Normierung des Dividenden.
    list pab = normieren_oT(p2,a2,b2,v,T);
    if (zerlegt) {return(pab);}
    p2 = pab[1]; a2 = pab[2]; b2 = pab[3];

    poly p3,a3,b3,q;  // Hilfspolynome
    list rq;          // Rueckgabelisten der Fkten pdiv und normieren_oT
    vector clden;

    // Variablen fuer Subresultanten-Variante initialisieren.
    int delta = degv(p1,v)-degv(p2,v);
    number alpha;
    number beta = (-1)**(delta + 1);
    number psi = -1;
    number m;

    // Pseudodivision
    rq = pdiv(p1,p2,v);
    q  = rq[2];
    p3 = reduce(rq[1],T);

    while (p3 <> 0)   // Hauptschleife des eukl. Algorithmus.
    {
        alpha = number(lcoef(p2,v));  // return von lcoef ist type poly
        m  = alpha**(delta + 1);   //m*p1 = q*p2 + p3  mod (T)
        p3 = p3 / beta;            //m*p1 = q*p2 + beta * p3  mod (T)

        // a und b anpassen.
        q  = reduce(q,T);
        a3 = reduce(m*a1 - q*a2,T);
        b3 = reduce(m*b1 - q*b2,T);
        a3 = a3 / beta;
        b3 = b3 / beta;

        // Variablen tauschen.
        p1 = p2;
        p2 = p3;
        a1 = a2;
        a2 = a3;
        b1 = b2;
        b2 = b3;

        // Normierung des Dividenden.
        pab = normieren_oT(p2,a2,b2,v,T);
        if (zerlegt) {return(pab);}
        p2 = pab[1]; a2 = pab[2]; b2 = pab[3];

        if (degv(p2,v) > 0)     // Dividend aus Grundring?, dann
        {                       // Subresultanten-Variabeln anpassen
            if (delta > 0)
            {
                psi = ((-alpha)**delta) / (psi**(delta-1));
            }
            delta = degv(p1,v)-degv(p2,v);
            beta  = (-alpha) * (psi**delta);

            rq = pdiv(p1,p2,v); // und Pseudodivision,
            q  = rq[2];
            p3 = reduce(rq[1], T);
        }
        else                    // ansonsten Ende.
        {
            p3 = 0;
        }
    }

    // Gemeinsame Teiler herausdividieren.
    clden = cleardenom([p2,a2,b2]);
    return(list(clden[1],clden[2],clden[3]));
}

///////////////////////////////////////////////////////////////////////////////

static proc normieren_oT (poly p, poly a, poly b, int v, ideal T)
"USAGE:   normieren_oT(p,a,b,v,T); p,a,b Polynome in var(v),
                                  T Dreiecksbasis.
RETURN:  Entweder Triangulierung [T1,T2] von (T) oder
         Liste,wobei [1] = reduce(r*p,T),
                     [2] = reduce(r*a,T),
                     [3] = reduce(r*b,T),
         und reduce(r*lcoef(p,v),T) ist im Grundring.
"
{
    poly lc = lcoef(p,v);
    if (deg(lc) == 0) {return(list(p,a,b));} // lc im Grundring
    lc = cleardenom(lc);

    // invertieren_oT aufrufen, um lc+(T) zu invertieren.
    list inv = invertieren_oT(lc,T);
    if (zerlegt) {return(inv);}

    // ansonsten Polynom zurueck gekommen.
    p = reduce(inv[1] * p,T);
    a = reduce(inv[1] * a,T);
    b = reduce(inv[1] * b,T);
    return(list(p,a,b));
}


///////////////////////////////////////////////////////////////////////////////
//
//  Der Moeller-Algorithmus
//
///////////////////////////////////////////////////////////////////////////////

proc triangM (ideal G, list #)
"USAGE:   triangM(G[,i]); G=ideal, i=integer,@*
ASSUME:  G is the reduced lexicographical Groebner basis of the
         zero-dimensional ideal (G), sorted by increasing leading terms.
RETURN:  a list of finitely many triangular systems, such that
         the union of their varieties equals the variety of (G).
         If i = 2, then each polynomial of the triangular systems
         is factorized.
NOTE:    Algorithm of Moeller (see: Moeller, H.M.: On decomposing systems of
         polynomial equations with finitely many solutions, Appl. Algebra Eng.
         Commun. Comput. 4, 217 - 230, 1993).
EXAMPLE: example triangM; shows an example
"
{
    // Test, ob G Groebnerbasis ist.
    if (attrib(G,"isSB") <> 1)
    {
        ERROR("The input is no groebner basis.");
    }

    // Faktorisierungsschalter setzen.
    int fak;
    if (size(#) > 0 && typeof(#[1]) == "int")
    {
        if (#[1] == 2) {fak = 1;}
    }

    // Noetige Optionen setzen.
    // obachman: save options so that they can be reset on exit
    intvec ovec = option(get);
    option(redSB);
    option(returnSB);

    list F,Fh;       // Ergebnislisten
    int nG = size(G);
    if (nG <= 2)     // G schon Dreiecksbasis?
    {
        if (fak)
        {
            F = faktorisiere_DB(G);
        }
        else
        {
            F = list(G);
        }
        option(set, ovec);
        return(F);
    }

    int u = lvar(G[nG-1]);
    int v = lvar(G[nG]);
    int m = nG;
    int l;
    ideal lc,G1,G2,H;
    poly lcr;
    int i,j,k;
      string s = string(nvars(basering)-v+1); // fuer Protokoll

    // Oberes Dreieckssytem abschneiden.
    while(u <> v)
    {
        m = m-1;
        if (m == 2)  // G ist bereits Dreiecksbasis.
        {
            if (fak)
            {
                F = faktorisiere_DB(G);
            }
            else
            {
                F = list(G);
            }
            option(set, ovec);
            return (F);
        }
        v = u;
        u = lvar(G[m-1]);
    }

    // Leitkoeffizienten der Polynome bzgl. var(u) bestimmen.
    l = m-1;
    while(u == v)
    {
        lc = lc + lcoef(G[l],u);
        l  = l - 1;
        u  = lvar(G[l]);
    }

    // Erster Teil des Splitts.
    G1 = sort_red(ideal(G[1..l]) + lc); // Sortiere und reduziere GB.

      s = string(nvars(basering)-v+1);
      dbprint(string(timer)+" S"+s+": Erster Teil des Splitts.");

    // Rekursiver Aufruf, um G1 zu triangulieren.
    F = triangM(G1,#);

      dbprint(string(timer)+" S"+s+": Oberes Dreieckssystem anhaengen.");

    // Oberes Dreieckssystem an jede berechnete Dreiecksbasis anhaengen.
    for (j = m; j <= nG; j++)  // meist nur j = m = nG
    {
        for (i = 1; i <= size(F); i++)
        {
            F[i] = F[i] + cleardenom(G[j][1] + reduce(G[j]-G[j][1],F[i]));
            attrib(F[i],"isSB",1);
        }
        if (fak)
        {
            F = faktorisiere_letzten(F);
        }
    }

        dbprint(string(timer)+" S"+s+": Zweiter Teil des Splitts, "+string(m-l-1)+" Leitkoeffizient(en).");

    // Zweiten Teil des Splitts berechnen.
    H = ideal(G[1..m]);
    attrib(H,"isSB",1);

    for (k = m-l-1; k >= 1; k--)  // Leitkoeffizientenliste durchgehen,
    {                             // angefangen mit lc des kleinsten Polynoms.
        lcr = reduce(lc[k],H);
        if (lcr <> 0)             // d.h. (H):lcr <> (1)
        {
              dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Berechnung Idealquotient.");
            G2 = quotient(H,lcr); // Wg. Option returnSB schon Standardbasis.
            attrib(G2,"isSB",1);

              dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Rekursiver Aufruf.");
            // Rekursiver Aufruf, um G2 zu triangulieren.
            Fh = triangM(G2,#);

              dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Oberes Dreieckssystem, falls vorhanden, anhaengen.");
            // Oberes Dreieckssystem, falls vorhanden,
            // an jede berechnete Dreiecksbasis anhaengen.
            for (j = m + 1; j <= nG; j++)
            {
                for (i = 1; i <= size(Fh); i++)
                {
                    Fh[i] = Fh[i] + cleardenom(G[j][1] +
                                    reduce(G[j]-G[j][1],Fh[i]));
                    attrib(Fh[i],"isSB",1);
                }
                if (fak)
                {
                    Fh = faktorisiere_letzten(Fh);
                }
            }
            F = F + Fh;
            if (k > 1)
            {
                  dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Lex. GB von (H + lc) berechnen, naechsten lc reduzieren.");
                H = std(H + lcr); // Wg. reduce(...,H) oben notwendig.
            }
        }
        else
        {
              dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Leere Varietaet.");
        }
    }

    // Singular loescht die Attribute beim Listenverketten, deswegen neu setzen.
    for (i = 1; i <= size(F); i++)
    {
        attrib(F[i],"isSB",1);
    }
    option(set, ovec);
    return(F);
}
//-------------------------------- example----------------------------------
example
{ "EXAMPLE:"; echo = 2;
   ring rC5 = 0,(e,d,c,b,a),lp;
   triangM(stdfglm(cyclic(5))); //oder: triangM(stdfglm(cyclic(5)),2);
}

///////////////////////////////////////////////////////////////////////////////

proc triangMH (ideal G, list #)
"USAGE:   triangMH(G[,i]); G=ideal, i=integer
ASSUME:  G is the reduced lexicographical Groebner basis of the
         zero-dimensional ideal (G), sorted by increasing leading terms.
RETURN:  a list of finitely many triangular systems, such that
         the disjoint union of their varieties equals the variety of (G).
         If i = 2, then each polynomial of the triangular systems is factorized.
NOTE:    Algorithm of Moeller and Hillebrand (see: Moeller, H.M.:
         On decomposing systems of polynomial equations with finitely many
         solutions, Appl. Algebra Eng. Commun. Comput. 4, 217 - 230, 1993 and
         Hillebrand, D.: Triangulierung nulldimensionaler Ideale -
         Implementierung und Vergleich zweier Algorithmen, master thesis,
         Universitaet Dortmund, Fachbereich Mathematik, Prof. Dr. H.M. Moeller,
         1999).
EXAMPLE: example triangMH; shows an example
"
{
    // Test, ob G Groebnerbasis ist.
    if (attrib(G,"isSB") <> 1)
    {
        ERROR("The input is no groebner basis.");
    }
    if(npars(basering)<>0)
    {
        ERROR("basering has parameters");
    }
    // Faktorisierungsschalter setzen.
    int fak;
    if (size(#) > 0 && typeof(#[1]) == "int")
    {
        if (#[1] == 2) {fak = 1;}
    }

    // Noetige Optionen setzen.
    // obachman: save options so that tey can be reset on exit
    intvec ovec = option(get);
    option(redSB);
    option(redTail);
    option(returnSB);

    list F,Fh;       // Ergebnislisten
    int nG = size(G);
    if (nG <= 2)     // G schon Dreiecksbasis?
    {
        if (fak)
        {
            F = faktorisiere_DB(G);
        }
        else
        {
            F = list(G);
        }
        option(set, ovec);
        return(F);
    }

    int u = lvar(G[nG-1]);
    int v = lvar(G[nG]);
    int m = nG;
    int l;
    poly lcr;
    ideal lc,G1,G2,H;
    int i,j,k;
      string s = string(nvars(basering)-v+1); // fuer Protokoll

    // Oberes Dreieckssytem abschneiden.
    while(u <> v)
    {
        m = m-1;
        if (m == 2)  // G ist Dreiecksbasis.
        {
            if (fak)
            {
                F = faktorisiere_DB(G);
            }
            else
            {
                F = list(G);
            }
            option(set, ovec);
            return(F);
         }
        v = u;
        u = lvar(G[m-1]);
    }

    // Leitkoeffizienten der Polynome bzgl. var(u) bestimmen.
    l = m-1;
    while(u == v)
    {
        lc = lc + lcoef(G[l],u);
        l  = l - 1;
        u  = lvar(G[l]);
    }

    // Erster Teil des Splitts.
    G1 = sort_red(ideal(G[1..l]) + lc); // Sortiere und reduziere GB.

      s = string(nvars(basering)-v+1);
      dbprint(string(timer)+" S"+s+": Erster Teil des Splitts.");

    // Rekursiver Aufruf, um G1 zu triangulieren.
    F = triangMH(G1,#);

      dbprint(string(timer)+" S"+s+": Oberes Dreieckssystem anhaengen.");

    // Oberes Dreieckssystem an jede berechnete Dreiecksbasis anhaengen.
    for (j = m; j <= nG; j++)  // meist nur j = m = nG
    {
        for (i = 1; i <= size(F); i++)
        {
            F[i] = F[i] + cleardenom(G[j][1] + reduce(G[j]-G[j][1],F[i]));
            attrib(F[i],"isSB",1);
        }
        if (fak)
        {
            F = faktorisiere_letzten(F);
        }
    }

        dbprint(string(timer)+" S"+s+": Zweiter Teil des Splitts, "+string(m-l-1)+" Leitkoeffizient(en).");

    // Zweiten Teil des Splitts berechnen.
    H = ideal(G[1..l]);   // Nur die Polynome in var(u).
    attrib(H,"isSB",1);

    for (k = m-l-1; k >= 1; k--)  // Leitkoeffizientenliste durchgehen,
    {                             // angefangen mit lc des kleinsten Polynoms.
        lcr = reduce(lc[k],H);
        if (lcr <> 0)             // d.h. (H):lcr <> (1)
        {
              dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Berechnung Saturation.");
            G2 = sat(H,lcr)[1];     // Saturation statt Idealquotient.
            attrib(G2,"isSB",1);

            if (G2[1] <> 1)
            {
                  dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Rekursiver Aufruf.");
                // Rekursiver Aufruf, um G2 zu triangulieren.
                Fh = triangMH(G2,#);

                  dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Ggt und oberes Dreieckssystem anhaengen.");
                // Ggt an jede Dreiecksbasis anhaengen.
                for (i = 1; i <= size(Fh); i++)
                {
                    Fh[i] = std(Fh[i] + G[m-k]);
                }
                if (fak)
                {
                    Fh = faktorisiere_letzten(Fh);
                }
                for (j = m+1; j <= nG; j++)  // oberes Dreieckssystem an jede
                {                            // berechnete Dreiecksbasis anhaengen
                    for (i = 1; i <= size(Fh); i++)
                    {
                        Fh[i] = Fh[i] + cleardenom(G[j][1]
                                      + reduce(G[j]-G[j][1],Fh[i]));
                        attrib(Fh[i],"isSB",1);
                    }
                    if (fak)
                    {
                        Fh = faktorisiere_letzten(Fh);
                    }
                }
                F = F + Fh;
            }
            else
            {
                  dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Leere Varietaet (G2 == (1)).");
            }
            if (k > 1)
            {
                  dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Lex. GB von (H + lc) berechnen, naechsten lc reduzieren.");
                H = std(H + lcr); // wegen reduce(...,H) oben
            }
        }
        else
        {
              dbprint(string(timer)+" S"+s+"L"+string(m-l-k)+": Leere Varietaet (lcr == 0).");
        }
    }

    // Singular loescht die Attribute beim Listenverketten, deswegen neu setzen.
    for (i = 1; i <= size(F); i++)
    {
        attrib(F[i],"isSB",1);
    }
    option(set, ovec);
    return(F);
}
//-------------------------------- example ----------------------------------
example
{ "EXAMPLE:"; echo = 2;
   ring rC5 = 0,(e,d,c,b,a),lp;
   triangMH(stdfglm(cyclic(5)));
}


///////////////////////////////////////////////////////////////////////////////
//
//  Hilfsfunktionen
//
///////////////////////////////////////////////////////////////////////////////

static proc sort_red (ideal G)
"USAGE:   sort_red(G);  G lexikographische Groebnerbasis.
RETURN:  reduzierte lex. GB G, deren Elemente nach aufsteigenden
         Leittermen sortiert sind.
"
{
    int i;

    // sortieren
    G = sort(G)[1];

    // reduzieren
    ideal Gred = cleardenom(G[1]);
    attrib(Gred,"isSB",1);

    for (i = 2; i <= size(G); i++)
    {
        Gred = Gred + cleardenom(reduce(G[i],Gred));
        attrib(Gred,"isSB",1);
    }
    attrib(Gred,"isSB",1);
    return(Gred);
}

///////////////////////////////////////////////////////////////////////////////

static proc Listenrest (list Liste)
"USAGE:   Listenrest(Liste); Liste muss eine mindestens einelementige
                            Liste sein (leere Liste ergibt Fehler).
RETURN:  Liste ohne das erste Element.
"
{
    return(delete(Liste,1));
}

///////////////////////////////////////////////////////////////////////////////

static proc Idealrest (ideal i)
"USAGE:   Idealrest(i); i Ideal, i <> (0).
RETURN:  Ideal i ohne das erste Element bzw das Nullideal,
         falls i nur ein Erzeugendes besass.
"
{
    int ni = size(i);
    if (ni == 1)
    {
        return(ideal(0));    // Nullideal
    }
    return(ideal(i[2..ni]));
}

///////////////////////////////////////////////////////////////////////////////

static proc H_anhaengen (list L, int h)
"USAGE:   H_anhaengen(L,h); L Liste aus Idealen, h natuerliche Zahl.
RETURN:  Macht aus dem Listenelement i ein Paar [i,h] und setzt
         Attribut "isSB" fuer i.
"
{
    int i;
    list E; // Ergebnisliste

    for (i = 1; i <= size(L); i++)
    {
        E = E + list(list(L[i],h));
        attrib(E[size(E)][1],"isSB",1);
    }
    return(E);
}

///////////////////////////////////////////////////////////////////////////////

static proc faktorisiere_letzten (list L)
"USAGE:   faktorisiere_letzten(L); L Liste aus Idealen.
RETURN:  Faktorisiert letztes Polynom aus jedem Ideal der Liste und
         zerlegt das Ideal entsprechend. Attribut "isSB" wird gesetzt.
"
{
    int i,j;
    ideal h;
    int nh;
    list factors;
    list E;   // Ergebnisliste
    for (i = 1; i <= size(L); i++)
    {
        h = L[i];
        nh = size(h);
        factors = factorize(h[nh],1);
        for (j = 1; j <= size(factors[1]); j++)
        {
            h[nh] = factors[1][j];
            E = insert(E,h,size(E));
            attrib(E[size(E)],"isSB",1);
        }
    }
    return(E);
}

///////////////////////////////////////////////////////////////////////////////

static proc faktorisiere_DB (ideal db)
"USAGE:   faktorisiere_DB(db); db reduzierte Dreiecksbasis.
RETURN:  Liste aus reduzierten Dreiecksbasen,
         die ensteht, indem, mit dem ersten beginnend,
         jedes Polynom der Dreiecksbasis faktorisiert wird.
"
{
    int i,j;
    poly h;
    list E = faktorisiere_letzten(db[1]);

    for (j = 2; j <= size(db); j++)
    {
        for (i = 1; i <= size(E); i++)
        {
            h    = db[j][1] + reduce(db[j]-db[j][1],E[i]);
            E[i] = E[i] + h;
        }
        E = faktorisiere_letzten(E);
    }
    return(E);
}

///////////////////////////////////////////////////////////////////////////////

static proc degv (poly f, int v)
"USAGE:   degv(f,v);  f Polynom in var(v).
RETURN:  Liefert den Grad von f bzgl. var(v) bzw -1, falls f == 0.
"
{
    if (f == 0) {return(-1);}
    return(size(coeffs(f,var(v))) - 1);
}

//////////////////////////////////////////////////////////////////////////////

static proc pdiv (poly f, poly g, int v)
"USAGE:   pdiv(f,g);  f,g Polynome in var(v), f > g.
RETURN:  Liste, wobei [1] = prem(f,g) (Pseudoremainder),
                      [2] = pquo(f,g) (Pseudoquotient).
NOTE:    Pseudodivision von f durch g.
"
{
    // Initialisierung
    poly q      = 0;
    poly rem    = f;
    int deg_g   = degv(g,v);
    int deg_rem = degv(rem,v);
    poly lcg    = lcoef(g,v);     // Im Lazard-Algor. aus dem Grundring.
    poly h;
    rem = rem*lcg**(deg_rem - deg_g + 1); // f durch g teilbar machen.

    // Gewoehnliche Polynomdivision
    while (rem <> 0 and deg_rem >= deg_g)
    {
        h       = (lcoef(rem,v)/lcg)*var(v)**(deg_rem - deg_g);
        q       = q + h;
        rem     = rem - h*g;
        deg_rem = degv(rem,v);
    }
    return(list(rem,q));
}

//////////////////////////////////////////////////////////////////////////////

static proc lvar (poly f)
"USAGE:   lvar(f);  f nicht-konstantes Polynom.
RETURN:  groesste Variable von f bzgl. der aktuellen Termordnung.
"
{
    int i = 1;
    intvec exp = leadexp(f);

    while (1)
    {
        if (exp[i] <> 0) {return(i);}
        i++;
    }
}

//////////////////////////////////////////////////////////////////////////////

static proc lcoef (poly f, int v)
"USAGE:   lcoef(f,v);  f Polynom, var(v) Ringvariable.
RETURN:  Leitkoeffizient von f bzgl. var(v).
"
{
    matrix m = coeffs(f,var(v));
    return(m[size(m),1]);
}
///////////////////////////////////////////////////////////////////////////////