/*
 * $Id$
 */

package edu.jas.application;


import java.io.Serializable;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashSet;
import java.util.Iterator;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.SortedMap;

import org.apache.log4j.Logger;

import edu.jas.gb.ExtendedGB;
import edu.jas.gb.GroebnerBaseAbstract;
import edu.jas.gb.Reduction;
import edu.jas.gbufd.GBFactory;
import edu.jas.gbufd.PolyGBUtil;
import edu.jas.gbufd.GroebnerBasePartial;
import edu.jas.gbmod.SyzygyAbstract;
import edu.jas.poly.AlgebraicNumber;
import edu.jas.poly.AlgebraicNumberRing;
import edu.jas.poly.ExpVector;
import edu.jas.poly.GenPolynomial;
import edu.jas.poly.GenPolynomialRing;
import edu.jas.poly.OptimizedPolynomialList;
import edu.jas.poly.PolyUtil;
import edu.jas.poly.PolynomialList;
import edu.jas.poly.TermOrder;
import edu.jas.poly.TermOrderOptimization;
import edu.jas.structure.GcdRingElem;
import edu.jas.structure.NotInvertibleException;
import edu.jas.structure.Power;
import edu.jas.structure.RingFactory;
import edu.jas.ufd.FactorAbstract;
import edu.jas.ufd.GCDFactory;
import edu.jas.ufd.GreatestCommonDivisor;
import edu.jas.ufd.PolyUfdUtil;
import edu.jas.ufd.Quotient;
import edu.jas.ufd.QuotientRing;
import edu.jas.ufd.SquarefreeAbstract;
import edu.jas.ufd.SquarefreeFactory;


/**
 * Ideal implements some methods for ideal arithmetic, for example intersection,
 * quotient and zero and positive dimensional ideal decomposition.
 * @author Heinz Kredel
 */
public class Ideal> implements Comparable>, Serializable {


    private static final Logger logger = Logger.getLogger(Ideal.class);


    private final boolean debug = logger.isDebugEnabled();


    /**
     * The data structure is a PolynomialList.
     */
    protected PolynomialList list;


    /**
     * Indicator if list is a Groebner Base.
     */
    protected boolean isGB;


    /**
     * Indicator if test has been performed if this is a Groebner Base.
     */
    protected boolean testGB;


    /**
     * Indicator if list has optimized term order.
     */
    protected boolean isTopt;


    /**
     * Groebner base engine.
     */
    protected final GroebnerBaseAbstract bb;


    /**
     * Reduction engine.
     */
    protected final Reduction red;


    /**
     * Squarefree decomposition engine.
     */
    protected final SquarefreeAbstract engine;


    /**
     * Constructor.
     * @param ring polynomial ring
     */
    public Ideal(GenPolynomialRing ring) {
        this(ring, new ArrayList>());
    }


    /**
     * Constructor.
     * @param ring polynomial ring
     * @param F list of polynomials
     */
    public Ideal(GenPolynomialRing ring, List> F) {
        this(new PolynomialList(ring, F));
    }


    /**
     * Constructor.
     * @param ring polynomial ring
     * @param F list of polynomials
     * @param gb true if F is known to be a Groebner Base, else false
     */
    public Ideal(GenPolynomialRing ring, List> F, boolean gb) {
        this(new PolynomialList(ring, F), gb);
    }


    /**
     * Constructor.
     * @param ring polynomial ring
     * @param F list of polynomials
     * @param gb true if F is known to be a Groebner Base, else false
     * @param topt true if term order is optimized, else false
     */
    public Ideal(GenPolynomialRing ring, List> F, boolean gb, boolean topt) {
        this(new PolynomialList(ring, F), gb, topt);
    }


    /**
     * Constructor.
     * @param list polynomial list
     */
    public Ideal(PolynomialList list) {
        this(list, false);
    }


    /**
     * Constructor.
     * @param list polynomial list
     * @param bb Groebner Base engine
     * @param red Reduction engine
     */
    public Ideal(PolynomialList list, GroebnerBaseAbstract bb, Reduction red) {
        this(list, false, bb, red);
    }


    /**
     * Constructor.
     * @param list polynomial list
     * @param gb true if list is known to be a Groebner Base, else false
     */
    public Ideal(PolynomialList list, boolean gb) {
        //this(list, gb, new GroebnerBaseSeqPairSeq(), new ReductionSeq());
        this(list, gb, GBFactory.getImplementation(list.ring.coFac));
    }


    /**
     * Constructor.
     * @param list polynomial list
     * @param gb true if list is known to be a Groebner Base, else false
     * @param topt true if term order is optimized, else false
     */
    public Ideal(PolynomialList list, boolean gb, boolean topt) {
        //this(list, gb, topt, new GroebnerBaseSeqPairSeq(), new ReductionSeq());
        this(list, gb, topt, GBFactory.getImplementation(list.ring.coFac));
    }


    /**
     * Constructor.
     * @param list polynomial list
     * @param gb true if list is known to be a Groebner Base, else false
     * @param bb Groebner Base engine
     * @param red Reduction engine
     */
    public Ideal(PolynomialList list, boolean gb, GroebnerBaseAbstract bb, Reduction red) {
        this(list, gb, false, bb, red);
    }


    /**
     * Constructor.
     * @param list polynomial list
     * @param gb true if list is known to be a Groebner Base, else false
     * @param bb Groebner Base engine
     */
    public Ideal(PolynomialList list, boolean gb, GroebnerBaseAbstract bb) {
        this(list, gb, false, bb, bb.red);
    }


    /**
     * Constructor.
     * @param list polynomial list
     * @param gb true if list is known to be a Groebner Base, else false
     * @param topt true if term order is optimized, else false
     * @param bb Groebner Base engine
     */
    public Ideal(PolynomialList list, boolean gb, boolean topt, GroebnerBaseAbstract bb) {
        this(list, gb, topt, bb, bb.red);
    }


    /**
     * Constructor.
     * @param list polynomial list
     * @param gb true if list is known to be a Groebner Base, else false
     * @param topt true if term order is optimized, else false
     * @param bb Groebner Base engine
     * @param red Reduction engine
     */
    public Ideal(PolynomialList list, boolean gb, boolean topt, GroebnerBaseAbstract bb,
                    Reduction red) {
        if (list == null || list.list == null) {
            throw new IllegalArgumentException("list and list.list may not be null");
        }
        this.list = list;
        this.isGB = gb;
        this.isTopt = topt;
        this.testGB = (gb ? true : false); // ??
        this.bb = bb;
        this.red = red;
        this.engine = SquarefreeFactory. getImplementation(list.ring.coFac);
    }


    /**
     * Clone this.
     * @return a copy of this.
     */
    public Ideal copy() {
        return new Ideal(list.copy(), isGB, isTopt, bb, red);
    }


    /**
     * Get the List of GenPolynomials.
     * @return list.list
     */
    public List> getList() {
        return list.list;
    }


    /**
     * Get the GenPolynomialRing.
     * @return list.ring
     */
    public GenPolynomialRing getRing() {
        return list.ring;
    }


    /**
     * Get the zero ideal.
     * @return ideal(0)
     */
    public Ideal getZERO() {
        List> z = new ArrayList>(0);
        PolynomialList pl = new PolynomialList(getRing(), z);
        return new Ideal(pl, true, isTopt, bb, red);
    }


    /**
     * Get the one ideal.
     * @return ideal(1)
     */
    public Ideal getONE() {
        List> one = new ArrayList>(1);
        one.add(list.ring.getONE());
        PolynomialList pl = new PolynomialList(getRing(), one);
        return new Ideal(pl, true, isTopt, bb, red);
    }


    /**
     * String representation of the ideal.
     * @see java.lang.Object#toString()
     */
    @Override
    public String toString() {
        return list.toString();
    }


    /**
     * Get a scripting compatible string representation.
     * @return script compatible representation for this Element.
     * @see edu.jas.structure.Element#toScript()
     */
    public String toScript() {
        // Python case
        return list.toScript();
    }


    /**
     * Comparison with any other object. Note: If both ideals are not Groebner
     * Bases, then false may be returned even the ideals are equal.
     * @see java.lang.Object#equals(java.lang.Object)
     */
    @Override
    @SuppressWarnings("unchecked")
    public boolean equals(Object b) {
        if (!(b instanceof Ideal)) {
            logger.warn("equals no Ideal");
            return false;
        }
        Ideal B = null;
        try {
            B = (Ideal) b;
        } catch (ClassCastException ignored) {
            return false;
        }
        //if ( isGB && B.isGB ) {
        //   return list.equals( B.list ); requires also monic polys
        //} else { // compute GBs ?
        return this.contains(B) && B.contains(this);
        //}
    }


    /**
     * Ideal list comparison.
     * @param L other Ideal.
     * @return compareTo() of polynomial lists.
     */
    public int compareTo(Ideal L) {
        return list.compareTo(L.list);
    }


    /**
     * Hash code for this ideal.
     * @see java.lang.Object#hashCode()
     */
    @Override
    public int hashCode() {
        int h;
        h = list.hashCode();
        if (isGB) {
            h = h << 1;
        }
        if (testGB) {
            h += 1;
        }
        return h;
    }


    /**
     * Test if ZERO ideal.
     * @return true, if this is the 0 ideal, else false
     */
    public boolean isZERO() {
        return list.isZERO();
    }


    /**
     * Test if ONE is contained in the ideal. To test for a proper ideal use
     * ! id.isONE().
     * @return true, if this is the 1 ideal, else false
     */
    public boolean isONE() {
        return list.isONE();
    }


    /**
     * Optimize the term order.
     */
    public void doToptimize() {
        if (isTopt) {
            return;
        }
        list = TermOrderOptimization. optimizeTermOrder(list);
        isTopt = true;
        if (isGB) {
            isGB = false;
            doGB();
        }
        return;
    }


    /**
     * Test if this is a Groebner base.
     * @return true, if this is a Groebner base, else false
     */
    public boolean isGB() {
        if (testGB) {
            return isGB;
        }
        logger.warn("isGB computing");
        isGB = bb.isGB(getList());
        testGB = true;
        return isGB;
    }


    /**
     * Do Groebner Base. compute the Groebner Base for this ideal.
     */
    @SuppressWarnings("unchecked")
    public void doGB() {
        if (isGB && testGB) {
            return;
        }
        //logger.warn("GB computing");
        List> G = getList();
        logger.info("GB computing = " + G);
        G = bb.GB(G);
        if (isTopt) {
            List perm = ((OptimizedPolynomialList) list).perm;
            list = new OptimizedPolynomialList(perm, getRing(), G);
        } else {
            list = new PolynomialList(getRing(), G);
        }
        isGB = true;
        testGB = true;
        return;
    }


    /**
     * Groebner Base. Get a Groebner Base for this ideal.
     * @return GB(this)
     */
    public Ideal GB() {
        if (isGB) {
            return this;
        }
        doGB();
        return this;
    }


    /**
     * Ideal containment. Test if B is contained in this ideal. Note: this is
     * eventually modified to become a Groebner Base.
     * @param B ideal
     * @return true, if B is contained in this, else false
     */
    public boolean contains(Ideal B) {
        if (B == null || B.isZERO()) {
            return true;
        }
        return contains(B.getList());
    }


    /**
     * Ideal containment. Test if b is contained in this ideal. Note: this is
     * eventually modified to become a Groebner Base.
     * @param b polynomial
     * @return true, if b is contained in this, else false
     */
    public boolean contains(GenPolynomial b) {
        if (b == null || b.isZERO()) {
            return true;
        }
        if (this.isONE()) {
            return true;
        }
        if (this.isZERO()) {
            return false;
        }
        if (!isGB) {
            doGB();
        }
        GenPolynomial z;
        z = red.normalform(getList(), b);
        if (z == null || z.isZERO()) {
            return true;
        }
        return false;
    }


    /**
     * Ideal containment. Test if each b in B is contained in this ideal. Note:
     * this is eventually modified to become a Groebner Base.
     * @param B list of polynomials
     * @return true, if each b in B is contained in this, else false
     */
    public boolean contains(List> B) {
        if (B == null || B.size() == 0) {
            return true;
        }
        if (this.isONE()) {
            return true;
        }
        if (!isGB) {
            doGB();
        }
        for (GenPolynomial b : B) {
            if (b == null) {
                continue;
            }
            GenPolynomial z = red.normalform(getList(), b);
            if (!z.isZERO()) {
                //System.out.println("contains nf(b) != 0: " + b);
                return false;
            }
        }
        return true;
    }


    /**
     * Summation. Generators for the sum of ideals. Note: if both ideals are
     * Groebner bases, a Groebner base is returned.
     * @param B ideal
     * @return ideal(this+B)
     */
    public Ideal sum(Ideal B) {
        if (B == null || B.isZERO()) {
            return this;
        }
        if (this.isZERO()) {
            return B;
        }
        int s = getList().size() + B.getList().size();
        List> c;
        c = new ArrayList>(s);
        c.addAll(getList());
        c.addAll(B.getList());
        Ideal I = new Ideal(getRing(), c, false);
        if (isGB && B.isGB) {
            I.doGB();
        }
        return I;
    }


    /**
     * Summation. Generators for the sum of ideal and a polynomial. Note: if
     * this ideal is a Groebner base, a Groebner base is returned.
     * @param b polynomial
     * @return ideal(this+{b})
     */
    public Ideal sum(GenPolynomial b) {
        if (b == null || b.isZERO()) {
            return this;
        }
        int s = getList().size() + 1;
        List> c;
        c = new ArrayList>(s);
        c.addAll(getList());
        c.add(b);
        Ideal I = new Ideal(getRing(), c, false);
        if (isGB) {
            I.doGB();
        }
        return I;
    }


    /**
     * Summation. Generators for the sum of this ideal and a list of
     * polynomials. Note: if this ideal is a Groebner base, a Groebner base is
     * returned.
     * @param L list of polynomials
     * @return ideal(this+L)
     */
    public Ideal sum(List> L) {
        if (L == null || L.isEmpty()) {
            return this;
        }
        int s = getList().size() + L.size();
        List> c = new ArrayList>(s);
        c.addAll(getList());
        c.addAll(L);
        Ideal I = new Ideal(getRing(), c, false);
        if (isGB) {
            I.doGB();
        }
        return I;
    }


    /**
     * Product. Generators for the product of ideals. Note: if both ideals are
     * Groebner bases, a Groebner base is returned.
     * @param B ideal
     * @return ideal(this*B)
     */
    public Ideal product(Ideal B) {
        if (B == null || B.isZERO()) {
            return B;
        }
        if (this.isZERO()) {
            return this;
        }
        int s = getList().size() * B.getList().size();
        List> c;
        c = new ArrayList>(s);
        for (GenPolynomial p : getList()) {
            for (GenPolynomial q : B.getList()) {
                q = p.multiply(q);
                c.add(q);
            }
        }
        Ideal I = new Ideal(getRing(), c, false);
        if (isGB && B.isGB) {
            I.doGB();
        }
        return I;
    }


    /**
     * Product. Generators for the product this ideal by a polynomial. 
     * Note: if this ideal is a
     * Groebner base, a Groebner base is returned.
     * @param b polynomial
     * @return ideal(this*b)
     */
    public Ideal product(GenPolynomial b) {
        if (b == null || b.isZERO()) {
            return getZERO();
        }
        if (this.isZERO()) {
            return this;
        }
        List> c;
        c = new ArrayList>(getList().size());
        for (GenPolynomial p : getList()) {
             GenPolynomial q = p.multiply(b);
             c.add(q);
        }
        Ideal I = new Ideal(getRing(), c, false);
        if (isGB) {
            I.doGB();
        }
        return I;
    }


    /**
     * Intersection. Generators for the intersection of ideals. Using an
     * iterative algorithm.
     * @param Bl list of ideals
     * @return ideal(cap_i B_i), a Groebner base
     */
    public Ideal intersect(List> Bl) {
        if (Bl == null || Bl.size() == 0) {
            return getZERO();
        }
        Ideal I = null;
        for (Ideal B : Bl) {
            if (I == null) {
                I = B;
                continue;
            }
            if (I.isONE()) {
                return I;
            }
            I = I.intersect(B);
        }
        return I;
    }


    /**
     * Intersection. Generators for the intersection of ideals.
     * @param B ideal
     * @return ideal(this \cap B), a Groebner base
     */
    public Ideal intersect(Ideal B) {
        if (B == null || B.isZERO()) { // (0)
            return B;
        }
        if (this.isZERO()) {
            return this;
        }
        List> c = PolyGBUtil. intersect(getRing(),getList(),B.getList());
        Ideal I = new Ideal(getRing(), c, true);
        return I;
    }


    /**
     * Intersection. Generators for the intersection of a ideal with a
     * polynomial ring. The polynomial ring of this ideal must be a
     * contraction of R and the TermOrder must be an elimination
     * order.
     * @param R polynomial ring
     * @return ideal(this \cap R)
     */
    public Ideal intersect(GenPolynomialRing R) {
        if (R == null) {
            throw new IllegalArgumentException("R may not be null");
        }
        List> H = PolyUtil. intersect(R,getList());
        return new Ideal(R, H, isGB, isTopt);
    }


    /**
     * Eliminate. Generators for the intersection of a ideal with a polynomial
     * ring. The polynomial rings must have variable names.
     * @param R polynomial ring
     * @return ideal(this \cap R)
     */
    public Ideal eliminate(GenPolynomialRing R) {
        if (R == null) {
            throw new IllegalArgumentException("R may not be null");
        }
        if (list.ring.equals(R)) {
            return this;
        }
        String[] ename = R.getVars();
        Ideal I = eliminate(ename);
        return I.intersect(R);
    }


    /**
     * Eliminate. Preparation of generators for the intersection of a ideal with
     * a polynomial ring.
     * @param ename variables for the elimination ring.
     * @return ideal(this) in K[ename,{vars \ ename}])
     */
    public Ideal eliminate(String... ename) {
        //System.out.println("ename = " + Arrays.toString(ename));
        if (ename == null) {
            throw new IllegalArgumentException("ename may not be null");
        }
        String[] aname = getRing().getVars();
        //System.out.println("aname = " + Arrays.toString(aname));
        if (aname == null) {
            throw new IllegalArgumentException("aname may not be null");
        }

        GroebnerBasePartial bbp = new GroebnerBasePartial(bb, null);
        String[] rname = GroebnerBasePartial.remainingVars(aname, ename);
        //System.out.println("rname = " + Arrays.toString(rname));
        PolynomialList Pl = null;
        if (rname.length == 0) {
            if (Arrays.equals(aname, ename)) {
                return this;
            }
            Pl = bbp.partialGB(getList(), ename); // normal GB
        } else {
            Pl = bbp.elimPartialGB(getList(), rname, ename); // reversed!
        }
        //System.out.println("Pl = " + Pl);
        if (debug) {
            logger.debug("elimination GB = " + Pl);
        }
        Ideal I = new Ideal(Pl, true);
        return I;
    }


    /**
     * Quotient. Generators for the ideal quotient.
     * @param h polynomial
     * @return ideal(this : h), a Groebner base
     */
    public Ideal quotient(GenPolynomial h) {
        if (h == null) { // == (0)
            return this;
        }
        if (h.isZERO()) {
            return this;
        }
        if (this.isZERO()) {
            return this;
        }
        List> H;
        H = new ArrayList>(1);
        H.add(h);
        Ideal Hi = new Ideal(getRing(), H, true);

        Ideal I = this.intersect(Hi);

        List> Q;
        Q = new ArrayList>(I.getList().size());
        for (GenPolynomial q : I.getList()) {
            q = q.divide(h); // remainder == 0
            Q.add(q);
        }
        return new Ideal(getRing(), Q, true /*false?*/);
    }


    /**
     * Quotient. Generators for the ideal quotient.
     * @param H ideal
     * @return ideal(this : H), a Groebner base
     */
    public Ideal quotient(Ideal H) {
        if (H == null) { // == (0)
            return this;
        }
        if (H.isZERO()) {
            return this;
        }
        if (this.isZERO()) {
            return this;
        }
        Ideal Q = null;
        for (GenPolynomial h : H.getList()) {
            Ideal Hi = this.quotient(h);
            if (Q == null) {
                Q = Hi;
            } else {
                Q = Q.intersect(Hi);
            }
        }
        return Q;
    }


    /**
     * Infinite quotient. Generators for the infinite ideal quotient.
     * @param h polynomial
     * @return ideal(this : hs), a Groebner base
     */
    public Ideal infiniteQuotientRab(GenPolynomial h) {
        if (h == null || h.isZERO()) { // == (0)
            return getONE();
        }
        if (h.isONE()) {
            return this;
        }
        if (this.isZERO()) {
            return this;
        }
        Ideal I = this.GB(); // should be already
        List> a = I.getList();
        List> c;
        c = new ArrayList>(a.size() + 1);

        GenPolynomialRing tfac = getRing().extend(1);
        // term order is also adjusted
        for (GenPolynomial p : a) {
            p = p.extend(tfac, 0, 0L); // p
            c.add(p);
        }
        GenPolynomial q = h.extend(tfac, 0, 1L);
        GenPolynomial r = tfac.getONE(); // h.extend( tfac, 0, 0L );
        GenPolynomial hs = q.subtract(r); // 1 - t*h // (1-t)*h
        c.add(hs);
        logger.warn("infiniteQuotientRab computing GB ");
        List> g = bb.GB(c);
        if (debug) {
            logger.info("infiniteQuotientRab    = " + tfac + ", c = " + c);
            logger.info("infiniteQuotientRab GB = " + g);
        }
        Ideal E = new Ideal(tfac, g, true);
        Ideal Is = E.intersect(getRing());
        return Is;
    }


    /**
     * Infinite quotient exponent.
     * @param h polynomial
     * @param Q quotient this : h^\infinity
     * @return s with Q = this : hs
     */
    public int infiniteQuotientExponent(GenPolynomial h, Ideal Q) {
        int s = 0;
        if (h == null) { // == 0
            return s;
        }
        if (h.isZERO() || h.isONE()) {
            return s;
        }
        if (this.isZERO() || this.isONE()) {
            return s;
        }
        //see below: if (this.contains(Q)) {
        //    return s;
        //}
        GenPolynomial p = getRing().getONE();
        for (GenPolynomial q : Q.getList()) {
            if (this.contains(q)) {
                continue;
            }
            //System.out.println("q = " + q + ", p = " + p + ", s = " + s);
            GenPolynomial qp = q.multiply(p);
            while (!this.contains(qp)) {
                p = p.multiply(h);
                s++;
                qp = q.multiply(p);
            }
        }
        return s;
    }


    /**
     * Infinite quotient. Generators for the infinite ideal quotient.
     * @param h polynomial
     * @return ideal(this : hs), a Groebner base
     */
    public Ideal infiniteQuotient(GenPolynomial h) {
        if (h == null) { // == (0)
            return this;
        }
        if (h.isZERO()) {
            return this;
        }
        if (this.isZERO()) {
            return this;
        }
        int s = 0;
        Ideal I = this.GB(); // should be already
        GenPolynomial hs = h;
        Ideal Is = I;

        boolean eq = false;
        while (!eq) {
            Is = I.quotient(hs);
            Is = Is.GB(); // should be already
            logger.info("infiniteQuotient s = " + s);
            eq = Is.contains(I); // I.contains(Is) always
            if (!eq) {
                I = Is;
                s++;
                // hs = hs.multiply( h );
            }
        }
        return Is;
    }


    /**
     * Radical membership test.
     * @param h polynomial
     * @return true if h is contained in the radical of ideal(this), else false.
     */
    public boolean isRadicalMember(GenPolynomial h) {
        if (h == null) { // == (0)
            return true;
        }
        if (h.isZERO()) {
            return true;
        }
        if (this.isZERO()) {
            return true;
        }
        Ideal x = infiniteQuotientRab(h);
        if (debug) {
            logger.debug("infiniteQuotientRab = " + x);
        }
        return x.isONE();
    }


    /**
     * Infinite quotient. Generators for the infinite ideal quotient.
     * @param h polynomial
     * @return ideal(this : hs), a Groebner base
     */
    public Ideal infiniteQuotientOld(GenPolynomial h) {
        if (h == null) { // == (0)
            return this;
        }
        if (h.isZERO()) {
            return this;
        }
        if (this.isZERO()) {
            return this;
        }
        int s = 0;
        Ideal I = this.GB(); // should be already
        GenPolynomial hs = h;

        boolean eq = false;
        while (!eq) {
            Ideal Is = I.quotient(hs);
            Is = Is.GB(); // should be already
            logger.debug("infiniteQuotient s = " + s);
            eq = Is.contains(I); // I.contains(Is) always
            if (!eq) {
                I = Is;
                s++;
                hs = hs.multiply(h);
            }
        }
        return I;
    }


    /**
     * Infinite Quotient. Generators for the ideal infinite quotient.
     * @param H ideal
     * @return ideal(this : Hs), a Groebner base
     */
    public Ideal infiniteQuotient(Ideal H) {
        if (H == null) { // == (0)
            return this;
        }
        if (H.isZERO()) {
            return this;
        }
        if (this.isZERO()) {
            return this;
        }
        Ideal Q = null;
        for (GenPolynomial h : H.getList()) {
            Ideal Hi = this.infiniteQuotient(h);
            if (Q == null) {
                Q = Hi;
            } else {
                Q = Q.intersect(Hi);
            }
        }
        return Q;
    }


    /**
     * Infinite Quotient. Generators for the ideal infinite quotient.
     * @param H ideal
     * @return ideal(this : Hs), a Groebner base
     */
    public Ideal infiniteQuotientRab(Ideal H) {
        if (H == null) { // == (0)
            return this;
        }
        if (H.isZERO()) {
            return this;
        }
        if (this.isZERO()) {
            return this;
        }
        Ideal Q = null;
        for (GenPolynomial h : H.getList()) {
            Ideal Hi = this.infiniteQuotientRab(h);
            if (Q == null) {
                Q = Hi;
            } else {
                Q = Q.intersect(Hi);
            }
        }
        return Q;
    }


    /**
     * Power. Generators for the power of this ideal. Note: if this ideal is a
     * Groebner base, a Groebner base is returned.
     * @param d integer
     * @return ideal(this^d)
     */
    public Ideal power(int d) {
        if (d <= 0) {
            return getONE();
        }
        if (this.isZERO() || this.isONE()) {
            return this;
        }
        Ideal c = this;
        for (int i = 1; i < d; i++) {
            c = c.product(this);
        }
        return c;
    }


    /**
     * Normalform for element.
     * @param h polynomial
     * @return normalform of h with respect to this
     */
    public GenPolynomial normalform(GenPolynomial h) {
        if (h == null) {
            return h;
        }
        if (h.isZERO()) {
            return h;
        }
        if (this.isZERO()) {
            return h;
        }
        GenPolynomial r;
        r = red.normalform(list.list, h);
        return r;
    }


    /**
     * Normalform for list of elements.
     * @param L polynomial list
     * @return list of normalforms of the elements of L with respect to this
     */
    public List> normalform(List> L) {
        if (L == null) {
            return L;
        }
        if (L.size() == 0) {
            return L;
        }
        if (this.isZERO()) {
            return L;
        }
        List> M = new ArrayList>(L.size());
        for (GenPolynomial h : L) {
            GenPolynomial r = normalform(h);
            if (r != null && !r.isZERO()) {
                M.add(r);
            }
        }
        return M;
    }


    /**
     * Annihilator for element modulo this ideal.
     * @param h polynomial
     * @return annihilator of h with respect to this
     */
    public Ideal annihilator(GenPolynomial h) {
        if (h == null || h.isZERO()) {
            return getZERO();
        }
        if (this.isZERO()) {
            return this;
        }
        doGB();
        List> F = new ArrayList>(1 + getList().size());
        F.add(h);
        F.addAll(getList());
        //System.out.println("F = " + F);
        SyzygyAbstract syz = new SyzygyAbstract();
        List>> S = syz.zeroRelationsArbitrary(F);
        //System.out.println("S = " + S);
        List> gen = new ArrayList>(S.size());
        for (List> rel : S) {
            if (rel == null || rel.isEmpty()) {
                continue;
            }
            GenPolynomial p = rel.get(0);
            if (p == null || p.isZERO()) {
                continue;
            }
            gen.add(p);
        }
        Ideal ann = new Ideal(getRing(),gen);
        //System.out.println("ann = " + ann);
        return ann;
    }


    /**
     * Test for annihilator of element modulo this ideal.
     * @param h polynomial
     * @param A ideal
     * @return true, if A is the annihilator of h with respect to this
     */
    public boolean isAnnihilator(GenPolynomial h, Ideal A) {
        Ideal B = A.product(h);
        return contains(B);
    }


    /**
     * Annihilator for ideal modulo this ideal.
     * @param H ideal
     * @return annihilator of H with respect to this
     */
    public Ideal annihilator(Ideal H) {
        if (H == null || H.isZERO()) {
            return getZERO();
        }
        if (this.isZERO()) {
            return this;
        }
        Ideal ann = null;
        for (GenPolynomial h : H.getList()) {
            Ideal Hi = this.annihilator(h);
            if (ann == null) {
                ann = Hi;
            } else {
                ann = ann.intersect(Hi);
            }
        }
        return ann;
    }


    /**
     * Test for annihilator of ideal modulo this ideal.
     * @param H ideal
     * @param A ideal
     * @return true, if A is the annihilator of H with respect to this
     */
    public boolean isAnnihilator(Ideal H, Ideal A) {
        Ideal B = A.product(H);
        return contains(B);
    }


    /**
     * Inverse for element modulo this ideal.
     * @param h polynomial
     * @return inverse of h with respect to this, if defined
     */
    public GenPolynomial inverse(GenPolynomial h) {
        if (h == null || h.isZERO()) {
            throw new NotInvertibleException("zero not invertible");
        }
        if (this.isZERO()) {
            throw new NotInvertibleException("zero ideal");
        }
        if (h.isUnit()) {
            return h.inverse();
        }
        doGB();
        List> F = new ArrayList>(1 + list.list.size());
        F.add(h);
        F.addAll(list.list);
        //System.out.println("F = " + F);
        ExtendedGB x = bb.extGB(F);
        List> G = x.G;
        //System.out.println("G = " + G);
        GenPolynomial one = null;
        int i = -1;
        for (GenPolynomial p : G) {
            i++;
            if (p == null) {
                continue;
            }
            if (p.isUnit()) {
                one = p;
                break;
            }
        }
        if (one == null) {
            throw new NotInvertibleException(" h = " + h);
        }
        List> row = x.G2F.get(i); // != -1
        GenPolynomial g = row.get(0);
        if (g == null || g.isZERO()) {
            throw new NotInvertibleException(" h = " + h);
        }
        // adjust g to get g*h == 1
        GenPolynomial f = g.multiply(h);
        GenPolynomial k = red.normalform(list.list, f);
        if (!k.isONE()) {
            C lbc = k.leadingBaseCoefficient();
            lbc = lbc.inverse();
            g = g.multiply(lbc);
        }
        if (debug) {
            //logger.info("inv G = " + G);
            //logger.info("inv G2F = " + x.G2F);
            //logger.info("inv row "+i+" = " + row);
            //logger.info("inv h = " + h);
            //logger.info("inv g = " + g);
            //logger.info("inv f = " + f);
            f = g.multiply(h);
            k = red.normalform(list.list, f);
            logger.debug("inv k = " + k);
            if (!k.isUnit()) {
                throw new NotInvertibleException(" k = " + k);
            }
        }
        return g;
    }


    /**
     * Test if element is a unit modulo this ideal.
     * @param h polynomial
     * @return true if h is a unit with respect to this, else false
     */
    public boolean isUnit(GenPolynomial h) {
        if (h == null || h.isZERO()) {
            return false;
        }
        if (this.isZERO()) {
            return false;
        }
        List> F = new ArrayList>(1 + list.list.size());
        F.add(h);
        F.addAll(list.list);
        List> G = bb.GB(F);
        for (GenPolynomial p : G) {
            if (p == null) {
                continue;
            }
            if (p.isUnit()) {
                return true;
            }
        }
        return false;
    }


    /**
     * Radical approximation. Squarefree generators for the ideal.
     * @return squarefree(this), a Groebner base
     */
    public Ideal squarefree() {
        if (this.isZERO()) {
            return this;
        }
        Ideal R = this;
        Ideal Rp = null;
        List> li, ri;
        while (true) {
            li = R.getList();
            ri = new ArrayList>(li); //.size() );
            for (GenPolynomial h : li) {
                GenPolynomial r = engine.squarefreePart(h);
                ri.add(r);
            }
            Rp = new Ideal(R.getRing(), ri, false);
            Rp.doGB();
            if (R.equals(Rp)) {
                break;
            }
            R = Rp;
        }
        return R;
    }


    /**
     * Ideal common zero test.
     * @return -1, 0 or 1 if dimension(this) &eq; -1, 0 or ≥ 1.
     */
    public int commonZeroTest() {
        if (this.isZERO()) {
            return 1;
        }
        if (!isGB) {
            doGB();
        }
        if (this.isONE()) {
            return -1;
        }
        return bb.commonZeroTest(getList());
    }


    /**
     * Test if this ideal is maximal.
     * @return true, if this is maximal and not one, else false.
     */
    public boolean isMaximal() {
        if (commonZeroTest() != 0) {
            return false;
        }
        for (Long d : univariateDegrees()) {
            if (d > 1L) {
                // todo: test if irreducible
                return false;
            }
        }
        return true;
    }


    /**
     * Univariate head term degrees.
     * @return a list of the degrees of univariate head terms.
     */
    public List univariateDegrees() {
        List ud = new ArrayList();
        if (this.isZERO()) {
            return ud;
        }
        if (!isGB) {
            doGB();
        }
        if (this.isONE()) {
            return ud;
        }
        return bb.univariateDegrees(getList());
    }


    /**
     * Ideal dimension.
     * @return a dimension container (dim,maxIndep,list(maxIndep),vars).
     */
    public Dimension dimension() {
        int t = commonZeroTest();
        Set S = new HashSet();
        Set> M = new HashSet>();
        if (t <= 0) {
            return new Dimension(t, S, M, this.list.ring.getVars());
        }
        int d = 0;
        Set U = new HashSet();
        for (int i = 0; i < this.list.ring.nvar; i++) {
            U.add(i);
        }
        M = dimension(S, U, M);
        for (Set m : M) {
            int dp = m.size();
            if (dp > d) {
                d = dp;
                S = m;
            }
        }
        return new Dimension(d, S, M, this.list.ring.getVars());
    }


    /**
     * Ideal dimension.
     * @param S is a set of independent variables.
     * @param U is a set of variables of unknown status.
     * @param M is a list of maximal sets of independent variables.
     * @return a list of maximal sets of independent variables, eventually
     *         containing S.
     */
    protected Set> dimension(Set S, Set U, Set> M) {
        Set> MP = M;
        Set UP = new HashSet(U);
        for (Integer j : U) {
            UP.remove(j);
            Set SP = new HashSet(S);
            SP.add(j);
            if (!containsHT(SP, getList())) {
                MP = dimension(SP, UP, MP);
            }
        }
        boolean contained = false;
        for (Set m : MP) {
            if (m.containsAll(S)) {
                contained = true;
                break;
            }
        }
        if (!contained) {
            MP.add(S);
        }
        return MP;
    }


    /**
     * Ideal head term containment test.
     * @param G list of polynomials.
     * @param H index set.
     * @return true, if the vaiables of the head terms of each polynomial in G
     *         are contained in H, else false.
     */
    protected boolean containsHT(Set H, List> G) {
        Set S = null;
        for (GenPolynomial p : G) {
            if (p == null) {
                continue;
            }
            ExpVector e = p.leadingExpVector();
            if (e == null) {
                continue;
            }
            int[] v = e.dependencyOnVariables();
            if (v == null) {
                continue;
            }
            //System.out.println("v = " + Arrays.toString(v));
            if (S == null) { // revert indices
                S = new HashSet(H.size());
                int r = e.length() - 1;
                for (Integer i : H) {
                    S.add(r - i);
                }
            }
            if (contains(v, S)) { // v \subset S
                return true;
            }
        }
        return false;
    }


    /**
     * Set containment. is v \subset H.
     * @param v index array.
     * @param H index set.
     * @return true, if each element of v is contained in H, else false .
     */
    protected boolean contains(int[] v, Set H) {
        for (int i = 0; i < v.length; i++) {
            if (!H.contains(v[i])) {
                return false;
            }
        }
        return true;
    }


    /**
     * Construct univariate polynomials of minimal degree in all variables in
     * zero dimensional ideal(G).
     * @return list of univariate polynomial of minimal degree in each variable
     *         in ideal(G)
     */
    public List> constructUnivariate() {
        List> univs = new ArrayList>();
        for (int i = list.ring.nvar - 1; i >= 0; i--) {
            GenPolynomial u = constructUnivariate(i);
            univs.add(u);
        }
        return univs;
    }


    /**
     * Construct univariate polynomial of minimal degree in variable i in zero
     * dimensional ideal(G).
     * @param i variable index.
     * @return univariate polynomial of minimal degree in variable i in ideal(G)
     */
    public GenPolynomial constructUnivariate(int i) {
        doGB();
        return bb.constructUnivariate(i, getList());
    }


    /**
     * Zero dimensional radical decompostition. See Seidenbergs lemma 92, and
     * BWK lemma 8.13.
     * @return intersection of radical ideals G_i with ideal(this) subseteq
     *         cap_i( ideal(G_i) )
     */
    public List> zeroDimRadicalDecomposition() {
        List> dec = new ArrayList>();
        if (this.isZERO()) {
            return dec;
        }
        IdealWithUniv iwu = new IdealWithUniv(this, new ArrayList>());
        dec.add(iwu);
        if (this.isONE()) {
            return dec;
        }
        if (list.ring.coFac.characteristic().signum() > 0 && !list.ring.coFac.isFinite()) {
            logger.warn("must use prime decomposition for char p and infinite coefficient rings, found "
                            + list.ring.coFac.toScript());
            return zeroDimPrimeDecomposition();
        }
        for (int i = list.ring.nvar - 1; i >= 0; i--) {
            List> part = new ArrayList>();
            for (IdealWithUniv id : dec) {
                //System.out.println("id = " + id + ", i = " + i);
                GenPolynomial u = id.ideal.constructUnivariate(i);
                SortedMap, Long> facs = engine.baseSquarefreeFactors(u);
                if (facs == null || facs.size() == 0 || (facs.size() == 1 && facs.get(facs.firstKey()) == 1L)) {
                    List> iup = new ArrayList>();
                    iup.addAll(id.upolys);
                    iup.add(u);
                    IdealWithUniv Ipu = new IdealWithUniv(id.ideal, iup);
                    part.add(Ipu);
                    continue; // irreducible
                }
                if (logger.isInfoEnabled()) {
                    logger.info("squarefree facs = " + facs);
                }
                GenPolynomialRing mfac = id.ideal.list.ring;
                int j = mfac.nvar - 1 - i;
                for (GenPolynomial p : facs.keySet()) {
                    // make p multivariate
                    GenPolynomial pm = p.extendUnivariate(mfac, j);
                    // mfac.parse( p.toString() );
                    //stem.out.println("pm = " + pm);
                    Ideal Ip = id.ideal.sum(pm);
                    List> iup = new ArrayList>();
                    iup.addAll(id.upolys);
                    iup.add(p);
                    IdealWithUniv Ipu = new IdealWithUniv(Ip, iup);
                    if (debug) {
                        logger.info("ideal with squarefree facs = " + Ipu);
                    }
                    part.add(Ipu);
                }
            }
            dec = part;
            //part = new ArrayList>();
        }
        return dec;
    }


    /**
     * Test for Zero dimensional radical. See Seidenbergs lemma 92, and BWK
     * lemma 8.13.
     * @return true if this is an zero dimensional radical ideal, else false
     */
    public boolean isZeroDimRadical() {
        if (this.isZERO()) {
            return false;
        }
        if (this.isONE()) {
            return false; // not 0-dim
        }
        if (list.ring.coFac.characteristic().signum() > 0 && !list.ring.coFac.isFinite()) {
            logger.warn("radical only for char 0 or finite coefficient rings, but found "
                            + list.ring.coFac.toScript());
        }
        for (int i = list.ring.nvar - 1; i >= 0; i--) {
            GenPolynomial u = constructUnivariate(i);
            boolean t = engine.isSquarefree(u);
            if (!t) {
                System.out.println("not squarefree " + engine.squarefreePart(u) + ", " + u);
                return false;
            }
        }
        return true;
    }


    /**
     * Zero dimensional ideal irreducible decompostition. See algorithm DIRGZD
     * of BGK 1986 and also PREDEC of the Gröbner bases book 1993.
     * @return intersection H, of ideals G_i with ideal(this) subseteq cap_i(
     *         ideal(G_i) ) and each ideal G_i has only irreducible minimal
     *         univariate polynomials and the G_i are pairwise co-prime.
     */
    public List> zeroDimDecomposition() {
        List> dec = new ArrayList>();
        if (this.isZERO()) {
            return dec;
        }
        IdealWithUniv iwu = new IdealWithUniv(this, new ArrayList>());
        dec.add(iwu);
        if (this.isONE()) {
            return dec;
        }
        FactorAbstract ufd = FactorFactory. getImplementation(list.ring.coFac);
        for (int i = list.ring.nvar - 1; i >= 0; i--) {
            List> part = new ArrayList>();
            for (IdealWithUniv id : dec) {
                //System.out.println("id.ideal = " + id.ideal);
                GenPolynomial u = id.ideal.constructUnivariate(i);
                SortedMap, Long> facs = ufd.baseFactors(u);
                if (facs.size() == 0 || (facs.size() == 1 && facs.get(facs.firstKey()) == 1L)) {
                    List> iup = new ArrayList>();
                    iup.addAll(id.upolys);
                    iup.add(u);
                    IdealWithUniv Ipu = new IdealWithUniv(id.ideal, iup);
                    part.add(Ipu);
                    continue; // irreducible
                }
                if (debug) {
                    logger.info("irreducible facs = " + facs);
                }
                GenPolynomialRing mfac = id.ideal.list.ring;
                int j = mfac.nvar - 1 - i;
                for (GenPolynomial p : facs.keySet()) {
                    // make p multivariate
                    GenPolynomial pm = p.extendUnivariate(mfac, j);
                    // mfac.parse( p.toString() );
                    //System.out.println("pm = " + pm);
                    Ideal Ip = id.ideal.sum(pm);
                    List> iup = new ArrayList>();
                    iup.addAll(id.upolys);
                    iup.add(p);
                    IdealWithUniv Ipu = new IdealWithUniv(Ip, iup);
                    part.add(Ipu);
                }
            }
            dec = part;
            //part = new ArrayList>();
        }
        return dec;
    }


    /**
     * Zero dimensional ideal irreducible decompostition extension. One step
     * decomposition via a minimal univariate polynomial in the lowest variable,
     * used after each normalPosition step.
     * @param upol list of univariate polynomials
     * @param og list of other generators for the ideal
     * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
     *         ideal(G_i) ) and all minimal univariate polynomials of all G_i
     *         are irreducible
     */
    public List> zeroDimDecompositionExtension(List> upol,
                    List> og) {
        if (upol == null || upol.size() + 1 != list.ring.nvar) {
            throw new IllegalArgumentException("univariate polynomial list not correct " + upol);
        }
        List> dec = new ArrayList>();
        if (this.isZERO()) {
            return dec;
        }
        IdealWithUniv iwu = new IdealWithUniv(this, upol);
        if (this.isONE()) {
            dec.add(iwu);
            return dec;
        }
        FactorAbstract ufd = FactorFactory. getImplementation(list.ring.coFac);
        int i = list.ring.nvar - 1;
        //IdealWithUniv id = new IdealWithUniv(this,upol);
        GenPolynomial u = this.constructUnivariate(i);
        SortedMap, Long> facs = ufd.baseFactors(u);
        if (facs.size() == 1 && facs.get(facs.firstKey()) == 1L) {
            List> iup = new ArrayList>();
            iup.add(u); // new polynomial first
            iup.addAll(upol);
            IdealWithUniv Ipu = new IdealWithUniv(this, iup, og);
            dec.add(Ipu);
            return dec;
        }
        if (true) {
            logger.info("irreducible facs = " + facs);
        }
        GenPolynomialRing mfac = list.ring;
        int j = mfac.nvar - 1 - i;
        for (GenPolynomial p : facs.keySet()) {
            // make p multivariate
            GenPolynomial pm = p.extendUnivariate(mfac, j);
            //System.out.println("pm = " + pm);
            Ideal Ip = this.sum(pm);
            List> iup = new ArrayList>();
            iup.add(p); // new polynomial first
            iup.addAll(upol);
            IdealWithUniv Ipu = new IdealWithUniv(Ip, iup, og);
            dec.add(Ipu);
        }
        return dec;
    }


    /**
     * Test for zero dimensional ideal decompostition.
     * @param L intersection of ideals G_i with ideal(G) subseteq cap_i(
     *            ideal(G_i) ) and all minimal univariate polynomials of all G_i
     *            are irreducible
     * @return true if L is a zero dimensional irreducible decomposition of
     *         this, else false
     */
    public boolean isZeroDimDecomposition(List> L) {
        if (L == null || L.size() == 0) {
            if (this.isZERO()) {
                return true;
            }
            return false;
        }
        // add lower variables if L contains ideals from field extensions
        GenPolynomialRing ofac = list.ring;
        int r = ofac.nvar;
        int rp = L.get(0).ideal.list.ring.nvar;
        int d = rp - r;
        //System.out.println("d = " + d);
        Ideal Id = this;
        if (d > 0) {
            GenPolynomialRing nfac = ofac.extendLower(d);
            //System.out.println("nfac = " + nfac);
            List> elist = new ArrayList>(list.list.size());
            for (GenPolynomial p : getList()) {
                //System.out.println("p = " + p);
                GenPolynomial q = p.extendLower(nfac, 0, 0L);
                //System.out.println("q = "  + q);
                elist.add(q);
            }
            Id = new Ideal(nfac, elist, isGB, isTopt);
        }
        // test if this is contained in the intersection
        for (IdealWithUniv I : L) {
            boolean t = I.ideal.contains(Id);
            if (!t) {
                System.out.println("not contained " + this + " in " + I.ideal);
                return false;
            }
        }
        // test if all univariate polynomials are contained in the respective ideal
        //List> upprod = new ArrayList>(rp);
        for (IdealWithUniv I : L) {
            GenPolynomialRing mfac = I.ideal.list.ring;
            int i = 0;
            for (GenPolynomial p : I.upolys) {
                GenPolynomial pm = p.extendUnivariate(mfac, i++);
                //System.out.println("pm = " + pm + ", p = " + p);
                boolean t = I.ideal.contains(pm);
                if (!t) {
                    System.out.println("not contained " + pm + " in " + I.ideal);
                    return false;
                }
            }
        }
        return true;
    }


    /**
     * Compute normal position for variables i and j.
     * @param i first variable index
     * @param j second variable index
     * @param og other generators for the ideal
     * @return this + (z - x_j - t x_i) in the ring C[z, x_1, ..., x_r]
     */
    public IdealWithUniv normalPositionFor(int i, int j, List> og) {
        // extend variables by one
        GenPolynomialRing ofac = list.ring;
        if (ofac.tord.getEvord() != TermOrder.INVLEX) {
            throw new IllegalArgumentException("invalid term order for normalPosition " + ofac.tord);
        }
        if (ofac.characteristic().signum() == 0) {
            return normalPositionForChar0(i, j, og);
        }
        return normalPositionForCharP(i, j, og);
    }


    /**
     * Compute normal position for variables i and j, characteristic zero.
     * @param i first variable index
     * @param j second variable index
     * @param og other generators for the ideal
     * @return this + (z - x_j - t x_i) in the ring C[z, x_1, ..., x_r]
     */
    IdealWithUniv normalPositionForChar0(int i, int j, List> og) {
        // extend variables by one
        GenPolynomialRing ofac = list.ring;
        GenPolynomialRing nfac = ofac.extendLower(1);
        List> elist = new ArrayList>(list.list.size() + 1);
        for (GenPolynomial p : getList()) {
            GenPolynomial q = p.extendLower(nfac, 0, 0L);
            //System.out.println("q = "  + q);
            elist.add(q);
        }
        List> ogen = new ArrayList>();
        if (og != null && og.size() > 0) {
            for (GenPolynomial p : og) {
                GenPolynomial q = p.extendLower(nfac, 0, 0L);
                //System.out.println("q = "  + q);
                ogen.add(q);
            }
        }
        Ideal I = new Ideal(nfac, elist, true);
        //System.out.println("I = "  + I);
        int ip = list.ring.nvar - 1 - i;
        int jp = list.ring.nvar - 1 - j;
        GenPolynomial xi = nfac.univariate(ip);
        GenPolynomial xj = nfac.univariate(jp);
        GenPolynomial z = nfac.univariate(nfac.nvar - 1);
        // compute GBs until value of t is OK
        Ideal Ip;
        GenPolynomial zp;
        int t = 0;
        do {
            t--;
            // zp = z - ( xj - xi * t )
            zp = z.subtract(xj.subtract(xi.multiply(nfac.fromInteger(t))));
            zp = zp.monic();
            Ip = I.sum(zp);
            //System.out.println("Ip = " + Ip);
            if (-t % 5 == 0) {
                logger.info("normal position, t = " + t);
            }
        } while (!Ip.isNormalPositionFor(i + 1, j + 1));
        if (debug) {
            logger.info("normal position = " + Ip);
        }
        ogen.add(zp);
        IdealWithUniv Ips = new IdealWithUniv(Ip, null, ogen);
        return Ips;
    }


    /**
     * Compute normal position for variables i and j, positive characteristic.
     * @param i first variable index
     * @param j second variable index
     * @param og other generators for the ideal
     * @return this + (z - x_j - t x_i) in the ring C[z, x_1, ..., x_r]
     */
    @SuppressWarnings("unchecked")
    IdealWithUniv normalPositionForCharP(int i, int j, List> og) {
        // extend variables by one
        GenPolynomialRing ofac = list.ring;
        GenPolynomialRing nfac = ofac.extendLower(1);
        List> elist = new ArrayList>(list.list.size() + 1);
        for (GenPolynomial p : getList()) {
            GenPolynomial q = p.extendLower(nfac, 0, 0L);
            //System.out.println("q = "  + q);
            elist.add(q);
        }
        List> ogen = new ArrayList>();
        if (og != null && og.size() > 0) {
            for (GenPolynomial p : og) {
                GenPolynomial q = p.extendLower(nfac, 0, 0L);
                //System.out.println("q = "  + q);
                ogen.add(q);
            }
        }
        Ideal I = new Ideal(nfac, elist, true);
        //System.out.println("I = "  + I);
        int ip = list.ring.nvar - 1 - i;
        int jp = list.ring.nvar - 1 - j;
        GenPolynomial xi = nfac.univariate(ip);
        GenPolynomial xj = nfac.univariate(jp);
        GenPolynomial z = nfac.univariate(nfac.nvar - 1);
        // compute GBs until value of t is OK
        Ideal Ip;
        GenPolynomial zp;
        AlgebraicNumberRing afac = null;
        Iterator> aiter = null;
        //String obr = "";
        //String cbr = "";
        int t = 0;
        do {
            t--;
            // zp = z - ( xj - xi * t )
            GenPolynomial tn;
            if (afac == null) {
                tn = nfac.fromInteger(t);
                if (tn.isZERO()) {
                    RingFactory fac = nfac.coFac;
                    //int braces = 2;
                    while (!(fac instanceof AlgebraicNumberRing)) {
                        if (fac instanceof GenPolynomialRing) {
                            GenPolynomialRing pfac = (GenPolynomialRing) (Object) fac;
                            fac = pfac.coFac;
                        } else if (fac instanceof QuotientRing) {
                            QuotientRing pfac = (QuotientRing) (Object) fac;
                            fac = pfac.ring.coFac;
                        } else {
                            throw new ArithmeticException(
                                            "field elements exhausted, need algebraic extension of base ring");
                        }
                        //braces++;
                    }
                    //for (int ii = 0; ii < braces; ii++) {
                    //    obr += "{ ";
                    //    cbr += " }";
                    //}
                    afac = (AlgebraicNumberRing) (Object) fac;
                    logger.info("afac = " + afac.toScript());
                    aiter = afac.iterator();
                    AlgebraicNumber an = aiter.next();
                    for (int kk = 0; kk < afac.characteristic().intValue(); kk++) {
                        an = aiter.next();
                    }
                    //System.out.println("an,iter = " + an);
                    //tn = nfac.parse(obr + an.toString() + cbr);
                    tn = nfac.parse(an.toString());
                    //System.out.println("tn = " + tn);
                    //if (false) {
                    //    throw new RuntimeException("probe");
                    //}
                }
            } else {
                if (!aiter.hasNext()) {
                    throw new ArithmeticException(
                                    "field elements exhausted, normal position not reachable: !aiter.hasNext(): "
                                                    + t);
                }
                AlgebraicNumber an = aiter.next();
                //System.out.println("an,iter = " + an);
                //tn = nfac.parse(obr + an.toString() + cbr);
                tn = nfac.parse(an.toString());
                //System.out.println("tn = " + tn);
            }
            if (tn.isZERO()) {
                throw new ArithmeticException(
                                "field elements exhausted, normal position not reachable: tn == 0: " + t);
            }
            zp = z.subtract(xj.subtract(xi.multiply(tn)));
            zp = zp.monic();
            Ip = I.sum(zp);
            //System.out.println("Ip = " + Ip);
            if (-t % 4 == 0) {
                logger.info("normal position, t = " + t);
                logger.info("normal position, GB = " + Ip);
                if (t < -550) {
                    throw new ArithmeticException("normal position not reached in " + t + " steps");
                }
            }
        } while (!Ip.isNormalPositionFor(i + 1, j + 1));
        if (debug) {
            logger.info("normal position = " + Ip);
        }
        ogen.add(zp);
        IdealWithUniv Ips = new IdealWithUniv(Ip, null, ogen);
        return Ips;
    }


    /**
     * Test if this ideal is in normal position for variables i and j.
     * @param i first variable index
     * @param j second variable index
     * @return true if this is in normal position with respect to i and j
     */
    public boolean isNormalPositionFor(int i, int j) {
        // called in extended ring!
        int ip = list.ring.nvar - 1 - i;
        int jp = list.ring.nvar - 1 - j;
        boolean iOK = false;
        boolean jOK = false;
        for (GenPolynomial p : getList()) {
            ExpVector e = p.leadingExpVector();
            int[] dov = e.dependencyOnVariables();
            //System.out.println("dov = " + Arrays.toString(dov));
            if (dov.length == 0) {
                throw new IllegalArgumentException("ideal dimension is not zero");
            }
            if (dov[0] == ip) {
                if (e.totalDeg() != 1) {
                    return false;
                }
                iOK = true;
            } else if (dov[0] == jp) {
                if (e.totalDeg() != 1) {
                    return false;
                }
                jOK = true;
            }
            if (iOK && jOK) {
                return true;
            }
        }
        return iOK && jOK;
    }


    /**
     * Normal position index, separate for polynomials with more than 2
     * variables. See also mas.masring.DIPDEC0#DIGISR
     * @return (i,j) for non-normal variables
     */
    public int[] normalPositionIndex2Vars() {
        int[] np = null;
        int i = -1;
        int j = -1;
        for (GenPolynomial p : getList()) {
            ExpVector e = p.leadingExpVector();
            int[] dov = e.dependencyOnVariables();
            //System.out.println("dov_head = " + Arrays.toString(dov));
            if (dov.length == 0) {
                throw new IllegalArgumentException("ideal dimension is not zero " + p);
            }
            // search bi-variate head terms
            if (dov.length >= 2) {
                i = dov[0];
                j = dov[1];
                break;
            }
            int n = dov[0];
            GenPolynomial q = p.reductum();
            e = q.degreeVector();
            dov = e.dependencyOnVariables();
            //System.out.println("dov_red  = " + Arrays.toString(dov));
            int k = Arrays.binarySearch(dov, n);
            int len = 2;
            if (k >= 0) {
                len = 3;
            }
            // search bi-variate reductas
            if (dov.length >= len) {
                switch (k) {
                case 0:
                    i = dov[1];
                    j = dov[2];
                    break;
                case 1:
                    i = dov[0];
                    j = dov[2];
                    break;
                case 2:
                    i = dov[0];
                    j = dov[1];
                    break;
                default:
                    i = dov[0];
                    j = dov[1];
                    break;
                }
                break;
            }
        }
        if (i < 0 || j < 0) {
            return np;
        }
        // adjust index
        i = list.ring.nvar - 1 - i;
        j = list.ring.nvar - 1 - j;
        np = new int[] { j, i }; // reverse
        logger.info("normalPositionIndex2Vars, np = " + Arrays.toString(np));
        return np;
    }


    /**
     * Normal position index, separate multiple univariate polynomials. See also
     * 
     * mas.masring.DIPDEC0#DIGISM
     * @return (i,j) for non-normal variables
     */
    public int[] normalPositionIndexUnivars() {
        int[] np = null; //new int[] { -1, -1 };
        int i = -1;
        int j = -1;
        // search multiple univariate polynomials with degree >= 2
        for (GenPolynomial p : getList()) {
            ExpVector e = p.degreeVector();
            int[] dov = e.dependencyOnVariables();
            long t = e.totalDeg(); // lt(p) would be enough
            //System.out.println("dov_univ = " + Arrays.toString(dov) + ", e = " + e);
            if (dov.length == 0) {
                throw new IllegalArgumentException("ideal dimension is not zero");
            }
            if (dov.length == 1 && t >= 2L) {
                if (i == -1) {
                    i = dov[0];
                } else if (j == -1) {
                    j = dov[0];
                    if (i > j) {
                        int x = i;
                        i = j;
                        j = x;
                    }
                }
            }
            if (i >= 0 && j >= 0) {
                break;
            }
        }
        if (i < 0 || j < 0) {
            // search polynomials with univariate head term and degree >= 2
            for (GenPolynomial p : getList()) {
                ExpVector e = p.leadingExpVector();
                long t = e.totalDeg();
                if (t >= 2) {
                    e = p.degreeVector();
                    int[] dov = e.dependencyOnVariables();
                    //System.out.println("dov_univ2 = " + Arrays.toString(dov) + " e = " + e);
                    if (dov.length == 0) {
                        throw new IllegalArgumentException("ideal dimension is not zero");
                    }
                    if (dov.length >= 2) {
                        i = dov[0];
                        j = dov[1];
                    }
                }
                if (i >= 0 && j >= 0) {
                    break;
                }
            }
        }
        if (i < 0 || j < 0) {
            return np;
        }
        // adjust index
        i = list.ring.nvar - 1 - i;
        j = list.ring.nvar - 1 - j;
        np = new int[] { j, i }; // reverse
        logger.info("normalPositionIndexUnivars, np = " + Arrays.toString(np));
        return np;
    }


    /**
     * Zero dimensional ideal decompostition for real roots. See algorithm
     * mas.masring.DIPDEC0#DINTSR.
     * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
     *         ideal(G_i) ) and each G_i contains at most bi-variate polynomials
     *         and all univariate minimal polynomials are irreducible
     */
    public List> zeroDimRootDecomposition() {
        List> dec = zeroDimDecomposition();
        if (this.isZERO()) {
            return dec;
        }
        if (this.isONE()) {
            return dec;
        }
        List> rdec = new ArrayList>();
        while (dec.size() > 0) {
            IdealWithUniv id = dec.remove(0);
            int[] ri = id.ideal.normalPositionIndex2Vars();
            if (ri == null || ri.length != 2) {
                rdec.add(id);
            } else {
                IdealWithUniv I = id.ideal.normalPositionFor(ri[0], ri[1], id.others);
                List> rd = I.ideal.zeroDimDecompositionExtension(id.upolys, I.others);
                //System.out.println("r_rd = " + rd);
                dec.addAll(rd);
            }
        }
        return rdec;
    }


    /**
     * Zero dimensional ideal prime decompostition. See algorithm
     * mas.masring.DIPDEC0#DINTSS.
     * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
     *         ideal(G_i) ) and each G_i is a prime ideal
     */
    public List> zeroDimPrimeDecomposition() {
        List> pdec = zeroDimPrimeDecompositionFE();
        List> dec = new ArrayList>();
        if (pdec.size() == 1) { // already prime
            IdealWithUniv Ip = pdec.get(0);
            int s = Ip.upolys.size() - getRing().nvar; // skip field ext univariate polys
            List> upol = Ip.upolys.subList(s, Ip.upolys.size());
            Ip = new IdealWithUniv(this, upol);
            dec.add(Ip);
            return dec;
        }
        for (IdealWithUniv Ip : pdec) {
            if (Ip.ideal.getRing().nvar == getRing().nvar) { // no field extension
                dec.add(Ip);
            } else { // remove field extension
                // add other generators for performance
                Ideal Id = Ip.ideal;
                if (Ip.others != null) {
                    //System.out.println("adding Ip.others = " + Ip.others);
                    List> pp = new ArrayList>();
                    pp.addAll(Id.getList());
                    pp.addAll(Ip.others);
                    Id = new Ideal(Id.getRing(), pp);
                }
                Ideal Is = Id.eliminate(getRing());
                //System.out.println("Is = " + Is);
                int s = Ip.upolys.size() - getRing().nvar; // skip field ext univariate polys
                List> upol = Ip.upolys.subList(s, Ip.upolys.size());
                IdealWithUniv Iu = new IdealWithUniv(Is, upol);
                //,Ip.others); used above and must be ignored here 
                dec.add(Iu);
            }
        }
        return dec;
    }


    /**
     * Zero dimensional ideal prime decompostition, with field extension. See
     * algorithm mas.masring.DIPDEC0#DINTSS.
     * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
     *         ideal(G_i) ) and each G_i is a prime ideal with eventually
     *         containing field extension variables
     */
    public List> zeroDimPrimeDecompositionFE() {
        List> dec = zeroDimRootDecomposition();
        if (this.isZERO()) {
            return dec;
        }
        if (this.isONE()) {
            return dec;
        }
        List> rdec = new ArrayList>();
        while (dec.size() > 0) {
            IdealWithUniv id = dec.remove(0);
            int[] ri = id.ideal.normalPositionIndexUnivars();
            if (ri == null || ri.length != 2) {
                rdec.add(id);
            } else {
                IdealWithUniv I = id.ideal.normalPositionFor(ri[0], ri[1], id.others);
                List> rd = I.ideal.zeroDimDecompositionExtension(id.upolys, I.others);
                //System.out.println("rd = " + rd);
                dec.addAll(rd);
            }
        }
        return rdec;
    }


    /**
     * Zero dimensional ideal associated primary ideal. See algorithm
     * mas.masring.DIPIDEAL#DIRLPI.
     * @param P prime ideal associated to this
     * @return primary ideal of this with respect to the associated pime ideal P
     */
    public Ideal primaryIdeal(Ideal P) {
        Ideal Qs = P;
        Ideal Q;
        int e = 0;
        do {
            Q = Qs;
            e++;
            Qs = Q.product(P);
        } while (Qs.contains(this));
        boolean t;
        Ideal As;
        do {
            As = this.sum(Qs);
            t = As.contains(Q);
            if (!t) {
                Q = Qs;
                e++;
                Qs = Q.product(P);
            }
        } while (!t);
        logger.info("exponent = " + e);
        return As;
    }


    /**
     * Zero dimensional ideal primary decompostition.
     * @return list of primary components of primary ideals G_i (pairwise
     *         co-prime) with ideal(this) = cap_i( ideal(G_i) ) together with
     *         the associated primes
     */
    public List> zeroDimPrimaryDecomposition() {
        List> pdec = zeroDimPrimeDecomposition();
        if (logger.isInfoEnabled()) {
            logger.info("prim decomp = " + pdec);
        }
        return zeroDimPrimaryDecomposition(pdec);
    }


    /**
     * Zero dimensional ideal elimination to original ring.
     * @param pdec list of prime ideals G_i
     * @return intersection of pairwise co-prime prime ideals G_i in the ring of
     *         this with ideal(this) = cap_i( ideal(G_i) )
     */
    public List> zeroDimElimination(List> pdec) {
        List> dec = new ArrayList>();
        if (this.isZERO()) {
            return dec;
        }
        if (this.isONE()) {
            dec.add(pdec.get(0));
            return dec;
        }
        List> qdec = new ArrayList>();
        for (IdealWithUniv Ip : pdec) {
            //System.out.println("Ip = " + Ip);
            List> epol = new ArrayList>();
            epol.addAll(Ip.ideal.getList());
            GenPolynomialRing mfac = Ip.ideal.list.ring;
            int j = 0;
            // add univariate polynomials for performance
            for (GenPolynomial p : Ip.upolys) {
                GenPolynomial pm = p.extendUnivariate(mfac, j++);
                if (j != 1) { // skip double
                    epol.add(pm);
                }
            }
            // add other generators for performance
            if (Ip.others != null) {
                epol.addAll(Ip.others);
            }
            Ideal Ipp = new Ideal(mfac, epol);
            // logger.info("eliminate_1 = " + Ipp);
            TermOrder to = null;
            if (mfac.tord.getEvord() != TermOrder.IGRLEX) {
                List> epols = new ArrayList>();
                to = new TermOrder(TermOrder.IGRLEX);
                GenPolynomialRing smfac = new GenPolynomialRing(mfac.coFac, mfac.nvar, to,
                                mfac.getVars());
                for (GenPolynomial p : epol) {
                    GenPolynomial pm = smfac.copy(p);
                    epols.add(pm.monic());
                }
                //epol = epols; 
                Ipp = new Ideal(smfac, epols);
            }
            epol = red.irreducibleSet(Ipp.getList());
            Ipp = new Ideal(Ipp.getRing(), epol);
            if (logger.isInfoEnabled()) {
                logger.info("eliminate = " + Ipp);
            }
            Ideal Is = Ipp.eliminate(list.ring);
            //System.out.println("Is = " + Is);
            if (to != null && !Is.list.ring.equals(list.ring)) {
                List> epols = new ArrayList>();
                for (GenPolynomial p : Is.getList()) {
                    GenPolynomial pm = list.ring.copy(p);
                    epols.add(pm);
                }
                Is = new Ideal(list.ring, epols);
                //System.out.println("Is = " + Is);
            }
            int k = Ip.upolys.size() - list.ring.nvar;
            List> up = new ArrayList>();
            for (int i = 0; i < list.ring.nvar; i++) {
                up.add(Ip.upolys.get(i + k));
            }
            IdealWithUniv Ie = new IdealWithUniv(Is, up);
            qdec.add(Ie);
        }
        return qdec;
    }


    /**
     * Zero dimensional ideal primary decompostition.
     * @param pdec list of prime ideals G_i with no field extensions
     * @return list of primary components of primary ideals G_i (pairwise
     *         co-prime) with ideal(this) = cap_i( ideal(G_i) ) together with
     *         the associated primes
     */
    public List> zeroDimPrimaryDecomposition(List> pdec) {
        List> dec = new ArrayList>();
        if (this.isZERO()) {
            return dec;
        }
        if (this.isONE()) {
            PrimaryComponent pc = new PrimaryComponent(pdec.get(0).ideal, pdec.get(0));
            dec.add(pc);
            return dec;
        }
        for (IdealWithUniv Ip : pdec) {
            Ideal Qs = this.primaryIdeal(Ip.ideal);
            PrimaryComponent pc = new PrimaryComponent(Qs, Ip);
            dec.add(pc);
        }
        return dec;
    }


    /**
     * Test for primary ideal decompostition.
     * @param L list of primary components G_i
     * @return true if ideal(this) == cap_i( ideal(G_i) )
     */
    public boolean isPrimaryDecomposition(List> L) {
        // test if this is contained in the intersection
        for (PrimaryComponent I : L) {
            boolean t = I.primary.contains(this);
            if (!t) {
                System.out.println("not contained " + this + " in " + I);
                return false;
            }
        }
        Ideal isec = null;
        for (PrimaryComponent I : L) {
            if (isec == null) {
                isec = I.primary;
            } else {
                isec = isec.intersect(I.primary);
            }
        }
        return this.contains(isec);
    }


    /**
     * Ideal extension.
     * @param vars list of variables for a polynomial ring for extension
     * @return ideal G, with coefficients in
     *         QuotientRing(GenPolynomialRing(vars))
     */
    public IdealWithUniv> extension(String... vars) {
        GenPolynomialRing fac = getRing();
        GenPolynomialRing efac = new GenPolynomialRing(fac.coFac, vars.length, fac.tord, vars);
        IdealWithUniv> ext = extension(efac);
        return ext;
    }


    /**
     * Ideal extension.
     * @param efac polynomial ring for extension
     * @return ideal G, with coefficients in QuotientRing(efac)
     */
    public IdealWithUniv> extension(GenPolynomialRing efac) {
        QuotientRing qfac = new QuotientRing(efac);
        IdealWithUniv> ext = extension(qfac);
        return ext;
    }


    /**
     * Ideal extension.
     * @param qfac quotient polynomial ring for extension
     * @return ideal G, with coefficients in qfac
     */
    public IdealWithUniv> extension(QuotientRing qfac) {
        GenPolynomialRing fac = getRing();
        GenPolynomialRing efac = qfac.ring;
        String[] rvars = GroebnerBasePartial.remainingVars(fac.getVars(), efac.getVars());
        //System.out.println("rvars = " + Arrays.toString(rvars));

        GroebnerBasePartial bbp = new GroebnerBasePartial();
        //wrong: OptimizedPolynomialList pgb = bbp.partialGB(getList(),rvars);
        OptimizedPolynomialList pgb = bbp.elimPartialGB(getList(), rvars, efac.getVars());
        if (logger.isInfoEnabled()) {
            logger.info("rvars = " + Arrays.toString(rvars));
            logger.info("partialGB = " + pgb);
        }

        GenPolynomialRing> rfac = new GenPolynomialRing>(efac,
                        rvars.length, fac.tord, rvars);
        List> plist = pgb.list;
        List>> rpgb = PolyUtil. recursive(rfac, plist);
        //System.out.println("rfac = " + rfac);
        GenPolynomialRing> qpfac = new GenPolynomialRing>(qfac, rfac);
        List>> qpgb = PolyUfdUtil. quotientFromIntegralCoefficients(qpfac, rpgb);
        //System.out.println("qpfac = " + qpfac);

        // compute f 
        GreatestCommonDivisor ufd = GCDFactory.getImplementation(fac.coFac);
        GenPolynomial f = null; // qfac.ring.getONE();
        for (GenPolynomial> p : rpgb) {
            if (f == null) {
                f = p.leadingBaseCoefficient();
            } else {
                f = ufd.lcm(f, p.leadingBaseCoefficient());
            }
        }
        //SquarefreeAbstract sqf = SquarefreeFactory.getImplementation(fac.coFac);
        //not required: f = sqf.squarefreePart(f);
        GenPolynomial> fp = rfac.getONE().multiply(f);
        GenPolynomial> fq = PolyUfdUtil. quotientFromIntegralCoefficients(qpfac, fp);
        if (logger.isInfoEnabled()) {
            logger.info("extension f = " + f);
            logger.info("ext = " + qpgb);
        }
        List>> upols = new ArrayList>>(0);
        List>> opols = new ArrayList>>(1);
        opols.add(fq);

        qpgb = PolyUtil.> monic(qpgb);
        Ideal> ext = new Ideal>(qpfac, qpgb);
        IdealWithUniv> extu = new IdealWithUniv>(ext, upols, opols);
        return extu;
    }


    /**
     * Ideal contraction and permutation.
     * @param eideal extension ideal of this.
     * @return contraction ideal of eideal in this polynomial ring
     */
    public IdealWithUniv permContraction(IdealWithUniv> eideal) {
        return Ideal. permutation(getRing(), Ideal. contraction(eideal));
    }


    /**
     * Ideal contraction.
     * @param eid extension ideal of this.
     * @return contraction ideal of eid in distributed polynomial ring
     */
    public static > IdealWithUniv contraction(IdealWithUniv> eid) {
        Ideal> eideal = eid.ideal;
        List>> qgb = eideal.getList();
        QuotientRing qfac = (QuotientRing) eideal.getRing().coFac;
        GenPolynomialRing> rfac = new GenPolynomialRing>(qfac.ring,
                        eideal.getRing());
        GenPolynomialRing dfac = qfac.ring.extend(eideal.getRing().getVars());
        TermOrder to = new TermOrder(qfac.ring.tord.getEvord());
        dfac = new GenPolynomialRing(dfac.coFac, dfac.nvar, to, dfac.getVars());
        //System.out.println("qfac = " + qfac);
        //System.out.println("rfac = " + rfac);
        //System.out.println("dfac = " + dfac);
        // convert polynomials
        List>> cgb = PolyUfdUtil. integralFromQuotientCoefficients(rfac,
                        qgb);
        List> dgb = PolyUtil. distribute(dfac, cgb);
        Ideal cont = new Ideal(dfac, dgb);
        // convert other polynomials
        List> opols = new ArrayList>();
        if (eid.others != null && eid.others.size() > 0) {
            List>> orpol = PolyUfdUtil. integralFromQuotientCoefficients(
                            rfac, eid.others);
            List> opol = PolyUtil. distribute(dfac, orpol);
            opols.addAll(opol);
        }
        // convert univariate polynomials
        List> upols = new ArrayList>(0);
        int i = 0;
        for (GenPolynomial> p : eid.upolys) {
            GenPolynomial> pm = p.extendUnivariate(eideal.getRing(), i++);
            //System.out.println("pm = " + pm + ", p = " + p);
            GenPolynomial> urpol = PolyUfdUtil
                            . integralFromQuotientCoefficients(rfac, pm);
            GenPolynomial upol = PolyUtil. distribute(dfac, urpol);
            upols.add(upol);
            //System.out.println("upol = " + upol);
        }
        // compute f 
        GreatestCommonDivisor ufd = GCDFactory.getImplementation(qfac.ring.coFac);
        GenPolynomial f = null; // qfac.ring.getONE();
        for (GenPolynomial> p : cgb) {
            if (f == null) {
                f = p.leadingBaseCoefficient();
            } else {
                f = ufd.lcm(f, p.leadingBaseCoefficient());
            }
        }
        GenPolynomial> fp = rfac.getONE().multiply(f);
        f = PolyUtil. distribute(dfac, fp);
        if (logger.isInfoEnabled()) {
            logger.info("contraction f = " + f);
            logger.info("cont = " + cont);
        }
        opols.add(f);
        if (f.isONE()) {
            IdealWithUniv cf = new IdealWithUniv(cont, upols, opols);
            return cf;
        }
        // compute ideal quotient by f
        Ideal Q = cont.infiniteQuotientRab(f);
        IdealWithUniv Qu = new IdealWithUniv(Q, upols, opols);
        return Qu;
    }


    /**
     * Ideal permutation.
     * @param oring polynomial ring to which variables are back permuted.
     * @param Cont ideal to be permuted
     * @return permutation of cont in polynomial ring oring
     */
    public static > IdealWithUniv permutation(GenPolynomialRing oring,
                    IdealWithUniv Cont) {
        Ideal cont = Cont.ideal;
        GenPolynomialRing dfac = cont.getRing();
        // (back) permutation of variables
        String[] ovars = oring.getVars();
        String[] dvars = dfac.getVars();
        //System.out.println("ovars = " + Arrays.toString(ovars));
        //System.out.println("dvars = " + Arrays.toString(dvars));
        if (Arrays.equals(ovars, dvars)) { // nothing to do
            return Cont;
        }
        List perm = GroebnerBasePartial.getPermutation(dvars, ovars);
        //System.out.println("perm  = " + perm);
        GenPolynomialRing pfac = TermOrderOptimization. permutation(perm, cont.getRing());
        if (logger.isInfoEnabled()) {
            logger.info("pfac = " + pfac);
        }
        List> ppolys = TermOrderOptimization. permutation(perm, pfac, cont.getList());
        //System.out.println("ppolys = " + ppolys);
        cont = new Ideal(pfac, ppolys);
        if (logger.isDebugEnabled()) {
            logger.info("perm cont = " + cont);
        }
        List> opolys = TermOrderOptimization. permutation(perm, pfac, Cont.others);
        //System.out.println("opolys = " + opolys);
        List> upolys = TermOrderOptimization. permutation(perm, pfac, Cont.upolys);
        //System.out.println("opolys = " + opolys);
        IdealWithUniv Cu = new IdealWithUniv(cont, upolys, opolys);
        return Cu;
    }


    /**
     * Ideal radical.
     * @return the radical ideal of this
     */
    public Ideal radical() {
        List> rdec = radicalDecomposition();
        List> dec = new ArrayList>(rdec.size());
        for (IdealWithUniv ru : rdec) {
            dec.add(ru.ideal);
        }
        Ideal R = intersect(dec);
        return R;
    }


    /**
     * Ideal radical decompostition.
     * @return intersection of ideals G_i with radical(this) eq cap_i(
     *         ideal(G_i) ) and each G_i is a radical ideal and the G_i are
     *         pairwise co-prime
     */
    public List> radicalDecomposition() {
        // check dimension
        int z = commonZeroTest();
        List> dec = new ArrayList>();
        List> ups = new ArrayList>();
        // dimension -1
        if (z < 0) {
            IdealWithUniv id = new IdealWithUniv(this, ups);
            dec.add(id); // see GB book
            return dec;
        }
        // dimension 0
        if (z == 0) {
            dec = zeroDimRadicalDecomposition();
            return dec;
        }
        // dimension > 0
        if (this.isZERO()) {
            return dec;
        }
        if (list.ring.coFac.characteristic().signum() > 0 && !list.ring.coFac.isFinite()) {
            // must not be the case at this point
            logger.warn("must use prime decomposition for char p and infinite coefficient rings, found "
                            + list.ring.coFac.toScript());
            return primeDecomposition();
        }
        Dimension dim = dimension();
        if (logger.isInfoEnabled()) {
            logger.info("dimension = " + dim);
        }

        // a short maximal independent set with small variables
        Set> M = dim.M;
        Set min = null;
        for (Set m : M) {
            if (min == null) {
                min = m;
                continue;
            }
            if (m.size() < min.size()) {
                min = m;
            }
        }
        int ms = min.size();
        Integer[] ia = new Integer[0];
        int mx = min.toArray(ia)[ms - 1];
        for (Set m : M) {
            if (m.size() == ms) {
                int mxx = m.toArray(ia)[ms - 1];
                if (mxx < mx) {
                    min = m;
                    mx = mxx;
                }
            }
        }
        //System.out.println("min = " + min);
        String[] mvars = new String[min.size()];
        int j = 0;
        for (Integer i : min) {
            mvars[j++] = dim.v[i];
        }
        if (logger.isInfoEnabled()) {
            logger.info("extension for variables = " + Arrays.toString(mvars) + ", indexes = " + min);
        }
        // reduce to dimension zero
        IdealWithUniv> Ext = extension(mvars);
        if (logger.isInfoEnabled()) {
            logger.info("extension = " + Ext);
        }

        List>> edec = Ext.ideal.zeroDimRadicalDecomposition();
        if (logger.isInfoEnabled()) {
            logger.info("0-dim radical decomp = " + edec);
        }
        // remove field extensions are not required
        // reconstruct dimension
        for (IdealWithUniv> ep : edec) {
            IdealWithUniv cont = permContraction(ep);
            //System.out.println("cont = " + cont);
            dec.add(cont);
        }
        IdealWithUniv extcont = permContraction(Ext);
        //System.out.println("extcont = " + extcont);

        // get f
        List> ql = extcont.others;
        if (ql.size() == 0) { // should not happen
            return dec;
        }
        GenPolynomial fx = ql.get(0);
        //System.out.println("cont(Ext) fx = " + fx + ", " + fx.ring);
        if (fx.isONE()) {
            return dec;
        }
        Ideal T = sum(fx);
        //System.out.println("T.rec = " + T.getList());
        if (T.isONE()) {
            logger.info("1 in ideal for " + fx);
            return dec;
        }
        if (logger.isInfoEnabled()) {
            logger.info("radical decomp ext-cont fx = " + fx);
            logger.info("recursion radical decomp T = " + T);
        }
        // recursion:
        List> Tdec = T.radicalDecomposition();
        if (logger.isInfoEnabled()) {
            logger.info("recursion radical decomp = " + Tdec);
        }
        dec.addAll(Tdec);
        return dec;
    }


    /**
     * Ideal irreducible decompostition.
     * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
     *         ideal(G_i) ) and each G_i is an ideal with irreducible univariate
     *         polynomials (after extension to a zero dimensional ideal) and the
     *         G_i are pairwise co-prime
     */
    public List> decomposition() {
        // check dimension
        int z = commonZeroTest();
        List> dec = new ArrayList>();
        // dimension -1
        if (z < 0) {
            //List> ups = new ArrayList>();
            //IdealWithUniv id = new IdealWithUniv(this, ups);
            //dec.add(id); see GB book
            return dec;
        }
        // dimension 0
        if (z == 0) {
            dec = zeroDimDecomposition();
            return dec;
        }
        // dimension > 0
        if (this.isZERO()) {
            return dec;
        }
        Dimension dim = dimension();
        if (logger.isInfoEnabled()) {
            logger.info("dimension = " + dim);
        }

        // shortest maximal independent set
        Set> M = dim.M;
        Set min = null;
        for (Set m : M) {
            if (min == null) {
                min = m;
                continue;
            }
            if (m.size() < min.size()) {
                min = m;
            }
        }
        //System.out.println("min = " + min);
        String[] mvars = new String[min.size()];
        int j = 0;
        for (Integer i : min) {
            mvars[j++] = dim.v[i];
        }
        if (logger.isInfoEnabled()) {
            logger.info("extension for variables = " + Arrays.toString(mvars));
        }
        // reduce to dimension zero
        IdealWithUniv> Ext = extension(mvars);
        if (logger.isInfoEnabled()) {
            logger.info("extension = " + Ext);
        }

        List>> edec = Ext.ideal.zeroDimDecomposition();
        if (logger.isInfoEnabled()) {
            logger.info("0-dim irred decomp = " + edec);
        }
        // remove field extensions are not required
        // reconstruct dimension
        for (IdealWithUniv> ep : edec) {
            IdealWithUniv cont = permContraction(ep);
            //System.out.println("cont = " + cont);
            dec.add(cont);
        }
        IdealWithUniv extcont = permContraction(Ext);
        //System.out.println("extcont = " + extcont);

        // get f
        List> ql = extcont.others;
        if (ql.size() == 0) { // should not happen
            return dec;
        }
        GenPolynomial fx = ql.get(0);
        //System.out.println("cont(Ext) fx = " + fx + ", " + fx.ring);
        if (fx.isONE()) {
            return dec;
        }
        Ideal T = sum(fx);
        //System.out.println("T.rec = " + T.getList());
        if (T.isONE()) {
            logger.info("1 in ideal for " + fx);
            return dec;
        }
        if (logger.isInfoEnabled()) {
            logger.info("irred decomp ext-cont fx = " + fx);
            logger.info("recursion irred decomp T = " + T);
        }
        // recursion:
        List> Tdec = T.decomposition();
        if (logger.isInfoEnabled()) {
            logger.info("recursion irred decomposition = " + Tdec);
        }
        dec.addAll(Tdec);
        return dec;
    }


    /**
     * Ideal prime decompostition.
     * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
     *         ideal(G_i) ) and each G_i is a prime ideal and the G_i are
     *         pairwise co-prime
     */
    public List> primeDecomposition() {
        // check dimension
        int z = commonZeroTest();
        List> dec = new ArrayList>();
        // dimension -1
        if (z < 0) {
            //List> ups = new ArrayList>();
            //IdealWithUniv id = new IdealWithUniv(this, ups);
            //dec.add(id); see GB book
            return dec;
        }
        // dimension 0
        if (z == 0) {
            dec = zeroDimPrimeDecomposition();
            return dec;
        }
        // dimension > 0
        if (this.isZERO()) {
            return dec;
        }
        Dimension dim = dimension();
        if (logger.isInfoEnabled()) {
            logger.info("dimension = " + dim);
        }

        // shortest maximal independent set
        Set> M = dim.M;
        Set min = null;
        for (Set m : M) {
            if (min == null) {
                min = m;
                continue;
            }
            if (m.size() < min.size()) {
                min = m;
            }
        }
        //System.out.println("min = " + min);
        String[] mvars = new String[min.size()];
        int j = 0;
        for (Integer i : min) {
            mvars[j++] = dim.v[i];
        }
        if (logger.isInfoEnabled()) {
            logger.info("extension for variables = " + Arrays.toString(mvars));
        }
        // reduce to dimension zero
        IdealWithUniv> Ext = extension(mvars);
        if (logger.isInfoEnabled()) {
            logger.info("extension = " + Ext);
        }
        List>> edec = Ext.ideal.zeroDimPrimeDecomposition();
        if (logger.isInfoEnabled()) {
            logger.info("0-dim prime decomp = " + edec);
        }
        // remove field extensions, already done
        // reconstruct dimension
        for (IdealWithUniv> ep : edec) {
            IdealWithUniv cont = permContraction(ep);
            //System.out.println("cont = " + cont);
            dec.add(cont);
        }
        // get f
        IdealWithUniv extcont = permContraction(Ext);
        //System.out.println("extcont = " + extcont);
        List> ql = extcont.others;
        if (ql.size() == 0) { // should not happen
            return dec;
        }
        GenPolynomial fx = ql.get(0);
        //System.out.println("cont(Ext) fx = " + fx + ", " + fx.ring);
        if (fx.isONE()) {
            return dec;
        }
        // compute exponent not required
        Ideal T = sum(fx);
        //System.out.println("T.rec = " + T.getList());
        if (T.isONE()) {
            logger.info("1 in ideal for " + fx);
            return dec;
        }
        if (logger.isInfoEnabled()) {
            logger.info("prime decomp ext-cont fx = " + fx);
            logger.info("recursion prime decomp T = " + T);
        }
        // recursion:
        List> Tdec = T.primeDecomposition();
        if (logger.isInfoEnabled()) {
            logger.info("recursion prime decomp = " + Tdec);
        }
        dec.addAll(Tdec);
        return dec;
    }


    /**
     * Test for ideal decompostition.
     * @param L intersection of ideals G_i with ideal(G) eq cap_i(ideal(G_i) )
     * @return true if L is a decomposition of this, else false
     */
    public boolean isDecomposition(List> L) {
        if (L == null || L.size() == 0) {
            if (this.isZERO()) {
                return true;
            }
            return false;
        }
        GenPolynomialRing ofac = list.ring;
        int r = ofac.nvar;
        int rp = L.get(0).ideal.list.ring.nvar;
        int d = rp - r;
        //System.out.println("d = " + d);
        Ideal Id = this;
        if (d > 0) { // add lower variables
            GenPolynomialRing nfac = ofac.extendLower(d);
            //System.out.println("nfac = " + nfac);
            List> elist = new ArrayList>(list.list.size());
            for (GenPolynomial p : getList()) {
                //System.out.println("p = " + p);
                GenPolynomial q = p.extendLower(nfac, 0, 0L);
                //System.out.println("q = "  + q);
                elist.add(q);
            }
            Id = new Ideal(nfac, elist, isGB, isTopt);
        }

        // test if this is contained in the intersection
        for (IdealWithUniv I : L) {
            boolean t = I.ideal.contains(Id);
            if (!t) {
                System.out.println("not contained " + this + " in " + I.ideal);
                return false;
            }
        }
        //         // test if all univariate polynomials are contained in the respective ideal
        //         for (IdealWithUniv I : L) {
        //             GenPolynomialRing mfac = I.ideal.list.ring;
        //             int i = 0;
        //             for (GenPolynomial p : I.upolys) {
        //                 GenPolynomial pm = p.extendUnivariate(mfac, i++);
        //                 //System.out.println("pm = " + pm + ", p = " + p);
        //                 boolean t = I.ideal.contains(pm);
        //                 if (!t) {
        //                     System.out.println("not contained " + pm + " in " + I.ideal);
        //                     return false;
        //                 }
        //             }
        //         }
        return true;
    }


    /**
     * Ideal primary decompostition.
     * @return list of primary components of primary ideals G_i (pairwise
     *         co-prime) with ideal(this) = cap_i( ideal(G_i) ) together with
     *         the associated primes
     */
    public List> primaryDecomposition() {
        // check dimension
        int z = commonZeroTest();
        List> dec = new ArrayList>();
        // dimension -1
        if (z < 0) {
            //List> ups = new ArrayList>();
            //IdealWithUniv id = new IdealWithUniv(this, ups);
            //PrimaryComponent pc = new PrimaryComponent(this, id);
            //dec.add(pc); see GB book
            return dec;
        }
        // dimension 0
        if (z == 0) {
            dec = zeroDimPrimaryDecomposition();
            return dec;
        }
        // dimension > 0
        if (this.isZERO()) {
            return dec;
        }
        Dimension dim = dimension();
        if (logger.isInfoEnabled()) {
            logger.info("dimension = " + dim);
        }

        // shortest maximal independent set
        Set> M = dim.M;
        Set min = null;
        for (Set m : M) {
            if (min == null) {
                min = m;
                continue;
            }
            if (m.size() < min.size()) {
                min = m;
            }
        }
        //System.out.println("min = " + min);
        String[] mvars = new String[min.size()];
        int j = 0;
        for (Integer i : min) {
            mvars[j++] = dim.v[i];
        }
        if (logger.isInfoEnabled()) {
            logger.info("extension for variables = " + Arrays.toString(mvars));
        }
        // reduce to dimension zero
        IdealWithUniv> Ext = extension(mvars);
        if (logger.isInfoEnabled()) {
            logger.info("extension = " + Ext);
        }

        List>> edec = Ext.ideal.zeroDimPrimaryDecomposition();
        if (logger.isInfoEnabled()) {
            logger.info("0-dim primary decomp = " + edec);
        }
        // remove field extensions, already done
        // reconstruct dimension
        List>> upq = new ArrayList>>();
        for (PrimaryComponent> ep : edec) {
            IdealWithUniv> epu = new IdealWithUniv>(ep.primary, upq);
            IdealWithUniv contq = permContraction(epu);
            IdealWithUniv contp = permContraction(ep.prime);
            PrimaryComponent pc = new PrimaryComponent(contq.ideal, contp);
            //System.out.println("pc = " + pc);
            dec.add(pc);
        }

        // get f
        IdealWithUniv extcont = permContraction(Ext);
        if (debug) {
            logger.info("cont(Ext) = " + extcont);
        }
        List> ql = extcont.others;
        if (ql.size() == 0) { // should not happen
            return dec;
        }
        GenPolynomial fx = ql.get(0);
        //System.out.println("cont(Ext) fx = " + fx + ", " + fx.ring);
        if (fx.isONE()) {
            return dec;
        }
        // compute exponent
        int s = this.infiniteQuotientExponent(fx, extcont.ideal);
        if (s == 0) {
            logger.info("exponent is 0 ");
            return dec;
        }
        if (s > 1) {
            fx = Power.> positivePower(fx, s);
        }
        if (debug) {
            logger.info("exponent fx = " + s + ", fx^s = " + fx);
        }

        Ideal T = sum(fx);
        //System.out.println("T.rec = " + T.getList());
        if (T.isONE()) {
            logger.info("1 in ideal for " + fx);
            return dec;
        }
        if (logger.isInfoEnabled()) {
            logger.info("primmary decomp ext-cont fx = " + fx);
            logger.info("recursion primary decomp T = " + T);
        }
        // recursion:
        List> Tdec = T.primaryDecomposition();
        if (logger.isInfoEnabled()) {
            logger.info("recursion primary decomp = " + Tdec);
        }
        dec.addAll(Tdec);
        return dec;
    }

}