/*
 * $Id$
 */

package edu.jas.ufd;


import java.util.ArrayList;
import java.util.List;

import org.apache.log4j.Logger;

import edu.jas.poly.GenPolynomial;
import edu.jas.poly.GenPolynomialRing;
import edu.jas.poly.PolyUtil;
import edu.jas.structure.GcdRingElem;
import edu.jas.structure.RingFactory;


/**
 * Greatest common divisor algorithms.
 * @author Heinz Kredel
 */

public abstract class GreatestCommonDivisorAbstract> implements
                GreatestCommonDivisor {


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


    private final boolean debug = logger.isDebugEnabled();


    /**
     * Get the String representation.
     * @see java.lang.Object#toString()
     */
    @Override
    public String toString() {
        return getClass().getName();
    }


    /**
     * GenPolynomial base coefficient content.
     * @param P GenPolynomial.
     * @return cont(P).
     */
    public C baseContent(GenPolynomial P) {
        if (P == null) {
            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
        }
        if (P.isZERO()) {
            return P.ring.getZEROCoefficient();
        }
        C d = null;
        for (C c : P.getMap().values()) {
            if (d == null) {
                d = c;
            } else {
                d = d.gcd(c);
            }
            if (d.isONE()) {
                return d;
            }
        }
        if (d.signum() < 0) {
            d = d.negate();
        }
        return d;
    }


    /**
     * GenPolynomial base coefficient primitive part.
     * @param P GenPolynomial.
     * @return pp(P).
     */
    public GenPolynomial basePrimitivePart(GenPolynomial P) {
        if (P == null) {
            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
        }
        if (P.isZERO()) {
            return P;
        }
        C d = baseContent(P);
        if (d.isONE()) {
            return P;
        }
        GenPolynomial pp = P.divide(d);
        if (debug) {
            GenPolynomial p = pp.multiply(d);
            if (!p.equals(P)) {
                throw new ArithmeticException("pp(p)*cont(p) != p: ");
            }
        }
        return pp;
    }


    /**
     * Univariate GenPolynomial greatest common divisor. Uses sparse
     * pseudoRemainder for remainder.
     * @param P univariate GenPolynomial.
     * @param S univariate GenPolynomial.
     * @return gcd(P,S).
     */
    public abstract GenPolynomial baseGcd(GenPolynomial P, GenPolynomial S);


    /**
     * GenPolynomial recursive content.
     * @param P recursive GenPolynomial.
     * @return cont(P).
     */
    public GenPolynomial recursiveContent(GenPolynomial> P) {
        if (P == null) {
            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
        }
        if (P.isZERO()) {
            return P.ring.getZEROCoefficient();
        }
        GenPolynomial d = null;
        for (GenPolynomial c : P.getMap().values()) {
            if (d == null) {
                d = c;
            } else {
                d = gcd(d, c); // go to recursion
            }
            if (d.isONE()) {
                return d;
            }
        }
        return d.abs();
    }


    /**
     * GenPolynomial recursive primitive part.
     * @param P recursive GenPolynomial.
     * @return pp(P).
     */
    public GenPolynomial> recursivePrimitivePart(GenPolynomial> P) {
        if (P == null) {
            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
        }
        if (P.isZERO()) {
            return P;
        }
        GenPolynomial d = recursiveContent(P);
        if (d.isONE()) {
            return P;
        }
        GenPolynomial> pp = PolyUtil. recursiveDivide(P, d);
        return pp;
    }


    /**
     * GenPolynomial base recursive content.
     * @param P recursive GenPolynomial.
     * @return baseCont(P).
     */
    public C baseRecursiveContent(GenPolynomial> P) {
        if (P == null) {
            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
        }
        if (P.isZERO()) {
            GenPolynomialRing cf = (GenPolynomialRing) P.ring.coFac;
            return cf.coFac.getZERO();
        }
        C d = null;
        for (GenPolynomial c : P.getMap().values()) {
            C cc = baseContent(c);
            if (d == null) {
                d = cc;
            } else {
                d = gcd(d, cc);
            }
            if (d.isONE()) {
                return d;
            }
        }
        if (d.signum() < 0) {
            d = d.negate();
        }
        return d;
    }


    /**
     * GenPolynomial base recursive primitive part.
     * @param P recursive GenPolynomial.
     * @return basePP(P).
     */
    public GenPolynomial> baseRecursivePrimitivePart(GenPolynomial> P) {
        if (P == null) {
            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
        }
        if (P.isZERO()) {
            return P;
        }
        C d = baseRecursiveContent(P);
        if (d.isONE()) {
            return P;
        }
        GenPolynomial> pp = PolyUtil. baseRecursiveDivide(P, d);
        return pp;
    }


    /**
     * GenPolynomial recursive greatest common divisor. Uses pseudoRemainder for
     * remainder.
     * @param P recursive GenPolynomial.
     * @param S recursive GenPolynomial.
     * @return gcd(P,S).
     */
    public GenPolynomial> recursiveGcd(GenPolynomial> P,
                    GenPolynomial> S) {
        if (S == null || S.isZERO()) {
            return P;
        }
        if (P == null || P.isZERO()) {
            return S;
        }
        if (P.ring.nvar <= 1) {
            return recursiveUnivariateGcd(P, S);
        }
        // distributed polynomials gcd
        GenPolynomialRing> rfac = P.ring;
        RingFactory> rrfac = rfac.coFac;
        GenPolynomialRing cfac = (GenPolynomialRing) rrfac;
        GenPolynomialRing dfac = cfac.extend(rfac.nvar);
        GenPolynomial Pd = PolyUtil. distribute(dfac, P);
        GenPolynomial Sd = PolyUtil. distribute(dfac, S);
        GenPolynomial Dd = gcd(Pd, Sd);
        // convert to recursive
        GenPolynomial> C = PolyUtil. recursive(rfac, Dd);
        return C;
    }


    /**
     * Univariate GenPolynomial recursive greatest common divisor. Uses
     * pseudoRemainder for remainder.
     * @param P univariate recursive GenPolynomial.
     * @param S univariate recursive GenPolynomial.
     * @return gcd(P,S).
     */
    public abstract GenPolynomial> recursiveUnivariateGcd(GenPolynomial> P,
                    GenPolynomial> S);


    /**
     * GenPolynomial content.
     * @param P GenPolynomial.
     * @return cont(P).
     */
    public GenPolynomial content(GenPolynomial P) {
        if (P == null) {
            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
        }
        GenPolynomialRing pfac = P.ring;
        if (pfac.nvar <= 1) {
            // baseContent not possible by return type
            throw new IllegalArgumentException(this.getClass().getName()
                            + " use baseContent for univariate polynomials");

        }
        GenPolynomialRing cfac = pfac.contract(1);
        GenPolynomialRing> rfac = new GenPolynomialRing>(cfac, 1);

        GenPolynomial> Pr = PolyUtil. recursive(rfac, P);
        GenPolynomial D = recursiveContent(Pr);
        return D;
    }


    /**
     * GenPolynomial primitive part.
     * @param P GenPolynomial.
     * @return pp(P).
     */
    public GenPolynomial primitivePart(GenPolynomial P) {
        if (P == null) {
            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
        }
        if (P.isZERO()) {
            return P;
        }
        GenPolynomialRing pfac = P.ring;
        if (pfac.nvar <= 1) {
            return basePrimitivePart(P);
        }
        GenPolynomialRing cfac = pfac.contract(1);
        GenPolynomialRing> rfac = new GenPolynomialRing>(cfac, 1);

        GenPolynomial> Pr = PolyUtil. recursive(rfac, P);
        GenPolynomial> PP = recursivePrimitivePart(Pr);

        GenPolynomial D = PolyUtil. distribute(pfac, PP);
        return D;
    }


    /**
     * GenPolynomial division. Indirection to GenPolynomial method.
     * @param a GenPolynomial.
     * @param b coefficient.
     * @return a/b.
     */
    public GenPolynomial divide(GenPolynomial a, C b) {
        if (b == null || b.isZERO()) {
            throw new IllegalArgumentException("division by zero");

        }
        if (a == null || a.isZERO()) {
            return a;
        }
        return a.divide(b);
    }


    /**
     * Coefficient greatest common divisor. Indirection to coefficient method.
     * @param a coefficient.
     * @param b coefficient.
     * @return gcd(a,b).
     */
    public C gcd(C a, C b) {
        if (b == null || b.isZERO()) {
            return a;
        }
        if (a == null || a.isZERO()) {
            return b;
        }
        return a.gcd(b);
    }


    /**
     * GenPolynomial greatest common divisor.
     * @param P GenPolynomial.
     * @param S GenPolynomial.
     * @return gcd(P,S).
     */
    public GenPolynomial gcd(GenPolynomial P, GenPolynomial S) {
        if (S == null || S.isZERO()) {
            return P;
        }
        if (P == null || P.isZERO()) {
            return S;
        }
        GenPolynomialRing pfac = P.ring;
        if (pfac.nvar <= 1) {
            GenPolynomial T = baseGcd(P, S);
            return T;
        }
        GenPolynomialRing cfac = pfac.contract(1);
        GenPolynomialRing> rfac;
        if (pfac.getVars() != null && pfac.getVars().length > 0) {
            String[] v = new String[] { pfac.getVars()[pfac.nvar - 1] };
            rfac = new GenPolynomialRing>(cfac, 1, v);
        } else {
            rfac = new GenPolynomialRing>(cfac, 1);
        }
        GenPolynomial> Pr = PolyUtil. recursive(rfac, P);
        GenPolynomial> Sr = PolyUtil. recursive(rfac, S);
        GenPolynomial> Dr = recursiveUnivariateGcd(Pr, Sr);
        GenPolynomial D = PolyUtil. distribute(pfac, Dr);
        return D;
    }


    /**
     * GenPolynomial least common multiple.
     * @param P GenPolynomial.
     * @param S GenPolynomial.
     * @return lcm(P,S).
     */
    public GenPolynomial lcm(GenPolynomial P, GenPolynomial S) {
        if (S == null || S.isZERO()) {
            return S;
        }
        if (P == null || P.isZERO()) {
            return P;
        }
        GenPolynomial C = gcd(P, S);
        GenPolynomial A = P.multiply(S);
        return PolyUtil. basePseudoDivide(A, C);
    }


    /**
     * List of GenPolynomials greatest common divisor.
     * @param A non empty list of GenPolynomials.
     * @return gcd(A_i).
     */
    public GenPolynomial gcd(List> A) {
        if (A == null || A.isEmpty()) {
            throw new IllegalArgumentException("A may not be empty");
        }
        GenPolynomial g = A.get(0);
        for (int i = 1; i < A.size(); i++) {
            GenPolynomial f = A.get(i);
            g = gcd(g, f);
        }
        return g;
    }


    /**
     * Univariate GenPolynomial resultant.
     * @param P univariate GenPolynomial.
     * @param S univariate GenPolynomial.
     * @return res(P,S).
     * @throws UnsupportedOperationException if there is no implementation in
     *             the sub-class.
     */
    @SuppressWarnings("unused")
    public GenPolynomial baseResultant(GenPolynomial P, GenPolynomial S) {
        // can not be abstract
        throw new UnsupportedOperationException("not implmented");
    }


    /**
     * Univariate GenPolynomial recursive resultant.
     * @param P univariate recursive GenPolynomial.
     * @param S univariate recursive GenPolynomial.
     * @return res(P,S).
     * @throws UnsupportedOperationException if there is no implementation in
     *             the sub-class.
     */
    @SuppressWarnings("unused")
    public GenPolynomial> recursiveUnivariateResultant(GenPolynomial> P,
                    GenPolynomial> S) {
        // can not be abstract
        throw new UnsupportedOperationException("not implmented");
    }


    /**
     * GenPolynomial recursive resultant.
     * @param P univariate recursive GenPolynomial.
     * @param S univariate recursive GenPolynomial.
     * @return res(P,S).
     * @throws UnsupportedOperationException if there is no implementation in
     *             the sub-class.
     */
    public GenPolynomial> recursiveResultant(GenPolynomial> P,
                    GenPolynomial> S) {
        if (S == null || S.isZERO()) {
            return S;
        }
        if (P == null || P.isZERO()) {
            return P;
        }
        GenPolynomialRing> rfac = P.ring;
        GenPolynomialRing cfac = (GenPolynomialRing) rfac.coFac;
        GenPolynomialRing dfac = cfac.extend(rfac.getVars());
        GenPolynomial Pp = PolyUtil. distribute(dfac, P);
        GenPolynomial Sp = PolyUtil. distribute(dfac, S);
        GenPolynomial res = resultant(Pp, Sp);
        GenPolynomial> Rr = PolyUtil. recursive(rfac, res);
        return Rr;
    }


    /**
     * GenPolynomial resultant. The input polynomials are considered as
     * univariate polynomials in the main variable.
     * @param P GenPolynomial.
     * @param S GenPolynomial.
     * @return res(P,S).
     * @see edu.jas.ufd.GreatestCommonDivisorSubres#recursiveResultant
     * @throws UnsupportedOperationException if there is no implementation in
     *             the sub-class.
     */
    public GenPolynomial resultant(GenPolynomial P, GenPolynomial S) {
        if (S == null || S.isZERO()) {
            return S;
        }
        if (P == null || P.isZERO()) {
            return P;
        }
        // no more hacked: GreatestCommonDivisorSubres ufd_sr = new GreatestCommonDivisorSubres();
        GenPolynomialRing pfac = P.ring;
        if (pfac.nvar <= 1) {
            return baseResultant(P, S);
        }
        GenPolynomialRing cfac = pfac.contract(1);
        GenPolynomialRing> rfac = new GenPolynomialRing>(cfac, 1);

        GenPolynomial> Pr = PolyUtil. recursive(rfac, P);
        GenPolynomial> Sr = PolyUtil. recursive(rfac, S);

        GenPolynomial> Dr = recursiveUnivariateResultant(Pr, Sr);
        GenPolynomial D = PolyUtil. distribute(pfac, Dr);
        return D;
    }


    /**
     * GenPolynomial co-prime list.
     * @param A list of GenPolynomials.
     * @return B with gcd(b,c) = 1 for all b != c in B and for all non-constant
     *         a in A there exists b in B with b|a. B does not contain zero or
     *         constant polynomials.
     */
    public List> coPrime(List> A) {
        if (A == null || A.isEmpty()) {
            return A;
        }
        List> B = new ArrayList>(A.size());
        // make a coprime to rest of list
        GenPolynomial a = A.get(0);
        //System.out.println("a = " + a);
        if (!a.isZERO() && !a.isConstant()) {
            for (int i = 1; i < A.size(); i++) {
                GenPolynomial b = A.get(i);
                GenPolynomial g = gcd(a, b).abs();
                if (!g.isONE()) {
                    a = PolyUtil. basePseudoDivide(a, g);
                    b = PolyUtil. basePseudoDivide(b, g);
                    GenPolynomial gp = gcd(a, g).abs();
                    while (!gp.isONE()) {
                        a = PolyUtil. basePseudoDivide(a, gp);
                        g = PolyUtil. basePseudoDivide(g, gp);
                        B.add(g); // gcd(a,g) == 1
                        g = gp;
                        gp = gcd(a, gp).abs();
                    }
                    if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
                        B.add(g); // gcd(a,g) == 1
                    }
                }
                if (!b.isZERO() && !b.isConstant()) {
                    B.add(b); // gcd(a,b) == 1
                }
            }
        } else {
            B.addAll(A.subList(1, A.size()));
        }
        // make rest coprime
        B = coPrime(B);
        //System.out.println("B = " + B);
        if (!a.isZERO() && !a.isConstant() /*&& !B.contains(a)*/) {
            a = a.abs();
            B.add(a);
        }
        return B;
    }


    /**
     * GenPolynomial co-prime list.
     * @param A list of GenPolynomials.
     * @return B with gcd(b,c) = 1 for all b != c in B and for all non-constant
     *         a in A there exists b in B with b|a. B does not contain zero or
     *         constant polynomials.
     */
    public List> coPrimeRec(List> A) {
        if (A == null || A.isEmpty()) {
            return A;
        }
        List> B = new ArrayList>();
        // make a co-prime to rest of list
        for (GenPolynomial a : A) {
            //System.out.println("a = " + a);
            B = coPrime(a, B);
            //System.out.println("B = " + B);
        }
        return B;
    }


    /**
     * GenPolynomial co-prime list.
     * @param a GenPolynomial.
     * @param P co-prime list of GenPolynomials.
     * @return B with gcd(b,c) = 1 for all b != c in B and for non-constant a
     *         there exists b in P with b|a. B does not contain zero or constant
     *         polynomials.
     */
    public List> coPrime(GenPolynomial a, List> P) {
        if (a == null || a.isZERO() || a.isConstant()) {
            return P;
        }
        List> B = new ArrayList>(P.size() + 1);
        // make a coprime to elements of the list P
        for (int i = 0; i < P.size(); i++) {
            GenPolynomial b = P.get(i);
            GenPolynomial g = gcd(a, b).abs();
            if (!g.isONE()) {
                a = PolyUtil. basePseudoDivide(a, g);
                b = PolyUtil. basePseudoDivide(b, g);
                // make g co-prime to new a, g is co-prime to c != b in P, B
                GenPolynomial gp = gcd(a, g).abs();
                while (!gp.isONE()) {
                    a = PolyUtil. basePseudoDivide(a, gp);
                    g = PolyUtil. basePseudoDivide(g, gp);
                    if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
                        B.add(g); // gcd(a,g) == 1 and gcd(g,c) == 1 for c != b in P, B
                    }
                    g = gp;
                    gp = gcd(a, gp).abs();
                }
                // make new g co-prime to new b
                gp = gcd(b, g).abs();
                while (!gp.isONE()) {
                    b = PolyUtil. basePseudoDivide(b, gp);
                    g = PolyUtil. basePseudoDivide(g, gp);
                    if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
                        B.add(g); // gcd(a,g) == 1 and gcd(g,c) == 1 for c != b in P, B
                    }
                    g = gp;
                    gp = gcd(b, gp).abs();
                }
                if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
                    B.add(g); // gcd(a,g) == 1 and gcd(g,c) == 1 for c != b in P, B
                }
            }
            if (!b.isZERO() && !b.isConstant() /*&& !B.contains(b)*/) {
                B.add(b); // gcd(a,b) == 1 and gcd(b,c) == 1 for c != b in P, B
            }
        }
        if (!a.isZERO() && !a.isConstant() /*&& !B.contains(a)*/) {
            B.add(a);
        }
        return B;
    }


    /**
     * GenPolynomial test for co-prime list.
     * @param A list of GenPolynomials.
     * @return true if gcd(b,c) = 1 for all b != c in B, else false.
     */
    public boolean isCoPrime(List> A) {
        if (A == null || A.isEmpty()) {
            return true;
        }
        if (A.size() == 1) {
            return true;
        }
        for (int i = 0; i < A.size(); i++) {
            GenPolynomial a = A.get(i);
            for (int j = i + 1; j < A.size(); j++) {
                GenPolynomial b = A.get(j);
                GenPolynomial g = gcd(a, b);
                if (!g.isONE()) {
                    System.out.println("not co-prime, a: " + a);
                    System.out.println("not co-prime, b: " + b);
                    System.out.println("not co-prime, g: " + g);
                    return false;
                }
            }
        }
        return true;
    }


    /**
     * GenPolynomial test for co-prime list of given list.
     * @param A list of GenPolynomials.
     * @param P list of co-prime GenPolynomials.
     * @return true if isCoPrime(P) and for all a in A exists p in P with p | a,
     *         else false.
     */
    public boolean isCoPrime(List> P, List> A) {
        if (!isCoPrime(P)) {
            return false;
        }
        if (A == null || A.isEmpty()) {
            return true;
        }
        for (GenPolynomial q : A) {
            if (q.isZERO() || q.isConstant()) {
                continue;
            }
            boolean divides = false;
            for (GenPolynomial p : P) {
                GenPolynomial a = PolyUtil. baseSparsePseudoRemainder(q, p);
                if (a.isZERO()) { // p divides q
                    divides = true;
                    break;
                }
            }
            if (!divides) {
                System.out.println("no divisor for: " + q);
                return false;
            }
        }
        return true;
    }


    /**
     * Univariate GenPolynomial extended greatest common divisor. Uses sparse
     * pseudoRemainder for remainder.
     * @param P univariate GenPolynomial.
     * @param S univariate GenPolynomial.
     * @return [ gcd(P,S), a, b ] with a*P + b*S = gcd(P,S).
     */
    @SuppressWarnings("unchecked")
    public GenPolynomial[] baseExtendedGcd(GenPolynomial P, GenPolynomial S) {
        //return P.egcd(S);
        GenPolynomial[] hegcd = baseHalfExtendedGcd(P, S);
        GenPolynomial[] ret = (GenPolynomial[]) new GenPolynomial[3];
        ret[0] = hegcd[0];
        ret[1] = hegcd[1];
        GenPolynomial x = hegcd[0].subtract(hegcd[1].multiply(P));
        GenPolynomial[] qr = PolyUtil. basePseudoQuotientRemainder(x, S);
        // assert qr[1].isZERO() 
        ret[2] = qr[0];
        return ret;
    }


    /**
     * Univariate GenPolynomial half extended greatest comon divisor. Uses
     * sparse pseudoRemainder for remainder.
     * @param S GenPolynomial.
     * @return [ gcd(P,S), a ] with a*P + b*S = gcd(P,S).
     */
    @SuppressWarnings("unchecked")
    public GenPolynomial[] baseHalfExtendedGcd(GenPolynomial P, GenPolynomial S) {
        //if ( P == null ) {
        //    throw new IllegalArgumentException("null P not allowed");
        //}
        GenPolynomial[] ret = (GenPolynomial[]) new GenPolynomial[2];
        ret[0] = null;
        ret[1] = null;
        if (S == null || S.isZERO()) {
            ret[0] = P;
            ret[1] = P.ring.getONE();
            return ret;
        }
        if (P == null || P.isZERO()) {
            ret[0] = S;
            ret[1] = S.ring.getZERO();
            return ret;
        }
        if (P.ring.nvar != 1) {
            throw new IllegalArgumentException(this.getClass().getName() + " not univariate polynomials "
                            + P.ring);
        }
        GenPolynomial q = P;
        GenPolynomial r = S;
        GenPolynomial c1 = P.ring.getONE().copy();
        GenPolynomial d1 = P.ring.getZERO().copy();
        while (!r.isZERO()) {
            GenPolynomial[] qr = PolyUtil. basePseudoQuotientRemainder(q, r);
            //q.divideAndRemainder(r);
            q = qr[0];
            GenPolynomial x = c1.subtract(q.multiply(d1));
            c1 = d1;
            d1 = x;
            q = r;
            r = qr[1];
        }
        // normalize ldcf(q) to 1, i.e. make monic
        C g = q.leadingBaseCoefficient();
        if (g.isUnit()) {
            C h = g.inverse();
            q = q.multiply(h);
            c1 = c1.multiply(h);
        }
        //assert ( ((c1.multiply(P)).remainder(S).equals(q) )); 
        ret[0] = q;
        ret[1] = c1;
        return ret;
    }


    /**
     * Univariate GenPolynomial greatest common divisor diophantine version.
     * @param P univariate GenPolynomial.
     * @param S univariate GenPolynomial.
     * @param c univariate GenPolynomial.
     * @return [ a, b ] with a*P + b*S = c and deg(a) < deg(S).
     */
    @SuppressWarnings("unchecked")
    public GenPolynomial[] baseGcdDiophant(GenPolynomial P, GenPolynomial S, GenPolynomial c) {
        GenPolynomial[] egcd = baseExtendedGcd(P, S);
        GenPolynomial g = egcd[0];
        GenPolynomial[] qr = PolyUtil. basePseudoQuotientRemainder(c, g);
        if (!qr[1].isZERO()) {
            throw new ArithmeticException("not solvable, r = " + qr[1] + ", c = " + c + ", g = " + g);
        }
        GenPolynomial q = qr[0];
        GenPolynomial a = egcd[1].multiply(q);
        GenPolynomial b = egcd[2].multiply(q);
        if (!a.isZERO() && a.degree(0) >= S.degree(0)) {
            qr = PolyUtil. basePseudoQuotientRemainder(a, S);
            a = qr[1];
            b = b.sum(P.multiply(qr[0]));
        }
        GenPolynomial[] ret = (GenPolynomial[]) new GenPolynomial[2];
        ret[0] = a;
        ret[1] = b;

        if (debug) {
            GenPolynomial y = ret[0].multiply(P).sum(ret[1].multiply(S));
            if (!y.equals(c)) {
                System.out.println("P  = " + P);
                System.out.println("S  = " + S);
                System.out.println("c  = " + c);
                System.out.println("a  = " + a);
                System.out.println("b  = " + b);
                System.out.println("y  = " + y);
                throw new ArithmeticException("not diophant, x = " + y.subtract(c));
            }
        }
        return ret;
    }


    /**
     * Univariate GenPolynomial partial fraction decomposition.
     * @param A univariate GenPolynomial.
     * @param P univariate GenPolynomial.
     * @param S univariate GenPolynomial.
     * @return [ A0, Ap, As ] with A/(P*S) = A0 + Ap/P + As/S with deg(Ap) <
     *         deg(P) and deg(As) < deg(S).
     */
    @SuppressWarnings("unchecked")
    public GenPolynomial[] basePartialFraction(GenPolynomial A, GenPolynomial P, GenPolynomial S) {
        GenPolynomial[] ret = (GenPolynomial[]) new GenPolynomial[3];
        ret[0] = null;
        ret[1] = null;
        ret[2] = null;
        GenPolynomial ps = P.multiply(S);
        GenPolynomial[] qr = PolyUtil. basePseudoQuotientRemainder(A, ps);
        ret[0] = qr[0];
        GenPolynomial r = qr[1];
        GenPolynomial[] diop = baseGcdDiophant(S, P, r); // switch arguments

        //         GenPolynomial x = diop[0].multiply(S).sum( diop[1].multiply(P) );
        //         if ( !x.equals(r) ) {
        //             System.out.println("r  = " + r);
        //             System.out.println("x  = " + x);
        //             throw new RuntimeException("not partial fraction, x = " + x);
        //         }

        ret[1] = diop[0];
        ret[2] = diop[1];
        if (ret[1].degree(0) >= P.degree(0)) {
            qr = PolyUtil. basePseudoQuotientRemainder(ret[1], P);
            ret[0] = ret[0].sum(qr[0]);
            ret[1] = qr[1];
        }
        if (ret[2].degree(0) >= S.degree(0)) {
            qr = PolyUtil. basePseudoQuotientRemainder(ret[2], S);
            ret[0] = ret[0].sum(qr[0]);
            ret[2] = qr[1];
        }
        return ret;
    }


    /**
     * Univariate GenPolynomial partial fraction decomposition.
     * @param A univariate GenPolynomial.
     * @param P univariate GenPolynomial.
     * @param e exponent for P.
     * @return [ F0, F1, ..., Fe ] with A/(P^e) = sum( Fi / P^i ) with deg(Fi) <
     *         deg(P).
     */
    public List> basePartialFraction(GenPolynomial A, GenPolynomial P, int e) {
        if (A == null || P == null || e == 0) {
            throw new IllegalArgumentException("null A, P or e = 0 not allowed");
        }
        List> pf = new ArrayList>(e);
        if (A.isZERO()) {
            for (int i = 0; i < e; i++) {
                pf.add(A);
            }
            return pf;
        }
        if (e == 1) {
            GenPolynomial[] qr = PolyUtil. basePseudoQuotientRemainder(A, P);
            pf.add(qr[0]);
            pf.add(qr[1]);
            return pf;
        }
        GenPolynomial a = A;
        for (int j = e; j > 0; j--) {
            GenPolynomial[] qr = PolyUtil. basePseudoQuotientRemainder(a, P);
            a = qr[0];
            pf.add(0, qr[1]);
        }
        pf.add(0, a);
        return pf;
    }


    /**
     * Univariate GenPolynomial partial fraction decomposition.
     * @param A univariate GenPolynomial.
     * @param D list of co-prime univariate GenPolynomials.
     * @return [ A0, A1,..., An ] with A/prod(D) = A0 + sum( Ai/Di ) with
     *         deg(Ai) < deg(Di).
     */
    public List> basePartialFraction(GenPolynomial A, List> D) {
        if (D == null || A == null) {
            throw new IllegalArgumentException("null A or D not allowed");
        }
        List> pf = new ArrayList>(D.size() + 1);
        if (A.isZERO() || D.size() == 0) {
            pf.add(A);
            for (int i = 0; i < D.size(); i++) {
                pf.add(A);
            }
            return pf;
        }
        List> Dp = new ArrayList>(D.size() - 1);
        GenPolynomial P = A.ring.getONE();
        GenPolynomial d1 = null;
        for (GenPolynomial d : D) {
            if (d1 == null) {
                d1 = d;
            } else {
                P = P.multiply(d);
                Dp.add(d);
            }
        }
        GenPolynomial[] qr = PolyUtil. basePseudoQuotientRemainder(A, P.multiply(d1));
        GenPolynomial A0 = qr[0];
        GenPolynomial r = qr[1];
        if (D.size() == 1) {
            pf.add(A0);
            pf.add(r);
            return pf;
        }
        GenPolynomial[] diop = baseGcdDiophant(P, d1, r); // switch arguments
        GenPolynomial A1 = diop[0];
        GenPolynomial S = diop[1];
        List> Fr = basePartialFraction(S, Dp);
        A0 = A0.sum(Fr.remove(0));
        pf.add(A0);
        pf.add(A1);
        pf.addAll(Fr);
        return pf;
    }


    /**
     * Test for Univariate GenPolynomial partial fraction decomposition.
     * @param A univariate GenPolynomial.
     * @param D list of (co-prime) univariate GenPolynomials.
     * @param F list of univariate GenPolynomials from a partial fraction
     *            computation.
     * @return true if A/prod(D) = F0 + sum( Fi/Di ) with deg(Fi) < deg(Di), Fi
     *         in F, else false.
     */
    public boolean isBasePartialFraction(GenPolynomial A, List> D,
                    List> F) {
        if (D == null || A == null || F == null) {
            throw new IllegalArgumentException("null A, F or D not allowed");
        }
        if (D.size() != F.size() - 1) {
            return false;
        }
        // A0*prod(D) + sum( Ai * Dip ), Dip = prod(D,j!=i)
        GenPolynomial P = A.ring.getONE();
        for (GenPolynomial d : D) {
            P = P.multiply(d);
        }
        List> Fp = new ArrayList>(F);
        GenPolynomial A0 = Fp.remove(0).multiply(P);
        //System.out.println("A0 = " + A0);
        int j = 0;
        for (GenPolynomial Fi : Fp) {
            P = A.ring.getONE();
            int i = 0;
            for (GenPolynomial d : D) {
                if (i != j) {
                    P = P.multiply(d);
                }
                i++;
            }
            //System.out.println("Fi = " + Fi);
            //System.out.println("P  = " + P);
            A0 = A0.sum(Fi.multiply(P));
            //System.out.println("A0 = " + A0);
            j++;
        }
        boolean t = A.equals(A0);
        if (!t) {
            System.out.println("not isPartFrac = " + A0);
        }
        return t;
    }


    /**
     * Test for Univariate GenPolynomial partial fraction decomposition.
     * @param A univariate GenPolynomial.
     * @param P univariate GenPolynomial.
     * @param e exponent for P.
     * @param F list of univariate GenPolynomials from a partial fraction
     *            computation.
     * @return true if A/(P^e) = F0 + sum( Fi / P^i ) with deg(Fi) < deg(P), Fi
     *         in F, else false.
     */
    public boolean isBasePartialFraction(GenPolynomial A, GenPolynomial P, int e,
                    List> F) {
        if (A == null || P == null || F == null || e == 0) {
            throw new IllegalArgumentException("null A, P, F or e = 0 not allowed");
        }
        GenPolynomial A0 = basePartialFractionValue(P, e, F);
        //         A.ring.getZERO();
        //         for ( GenPolynomial Fi : F ) {
        //             A0 = A0.multiply(P);
        //             A0 = A0.sum(Fi);
        //             //System.out.println("A0 = " + A0);
        //         }
        boolean t = A.equals(A0);
        if (!t) {
            System.out.println("not isPartFrac = " + A0);
        }
        return t;
    }


    /**
     * Test for Univariate GenPolynomial partial fraction decomposition.
     * @param P univariate GenPolynomial.
     * @param e exponent for P.
     * @param F list of univariate GenPolynomials from a partial fraction
     *            computation.
     * @return (F0 + sum( Fi / P^i )) * P^e.
     */
    public GenPolynomial basePartialFractionValue(GenPolynomial P, int e, List> F) {
        if (P == null || F == null || e == 0) {
            throw new IllegalArgumentException("null P, F or e = 0 not allowed");
        }
        GenPolynomial A0 = P.ring.getZERO();
        for (GenPolynomial Fi : F) {
            A0 = A0.multiply(P);
            A0 = A0.sum(Fi);
            //System.out.println("A0 = " + A0);
        }
        return A0;
    }

}