/*
 * $Id$
 */

package edu.jas.gbufd;


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

import org.apache.log4j.Logger;

import edu.jas.gb.GroebnerBaseAbstract;
import edu.jas.gb.SolvableGroebnerBaseAbstract;
import edu.jas.gb.SolvableGroebnerBaseSeq;
import edu.jas.gb.SolvableReductionAbstract;
import edu.jas.gb.SolvableReductionSeq;
import edu.jas.poly.ExpVector;
import edu.jas.poly.GenPolynomial;
import edu.jas.poly.GenPolynomialRing;
import edu.jas.poly.GenSolvablePolynomial;
import edu.jas.poly.GenSolvablePolynomialRing;
import edu.jas.poly.PolyUtil;
import edu.jas.structure.GcdRingElem;
import edu.jas.structure.RingElem;


/**
 * Package gbufd utilities.
 * @author Heinz Kredel
 */

public class PolyGBUtil {


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


    private static boolean debug = logger.isDebugEnabled();


    /**
     * Test for resultant.
     * @param A generic polynomial.
     * @param B generic polynomial.
     * @param r generic polynomial.
     * @return true if res(A,B) isContained in ideal(A,B), else false.
     */
    public static > boolean isResultant(GenPolynomial A, GenPolynomial B,
                    GenPolynomial r) {
        if (r == null || r.isZERO()) {
            return true;
        }
        GroebnerBaseAbstract bb = GBFactory. getImplementation(r.ring.coFac);
        List> F = new ArrayList>(2);
        F.add(A);
        F.add(B);
        List> G = bb.GB(F);
        //System.out.println("G = " + G);
        GenPolynomial n = bb.red.normalform(G, r);
        //System.out.println("n = " + n);
        return n.isZERO();
    }


    /**
     * Top pseudo reduction wrt the main variables.
     * @param P generic polynomial.
     * @param A list of generic polynomials sorted according to appearing main
     *            variables.
     * @return top pseudo remainder of P wrt. A for the appearing variables.
     */
    public static > GenPolynomial topPseudoRemainder(List> A,
                    GenPolynomial P) {
        if (A == null || A.isEmpty()) {
            return P.monic();
        }
        if (P.isZERO()) {
            return P;
        }
        //System.out.println("remainder, P = " + P);
        GenPolynomialRing pfac = A.get(0).ring;
        if (pfac.nvar <= 1) { // recursion base 
            GenPolynomial R = PolyUtil. baseSparsePseudoRemainder(P, A.get(0));
            return R.monic();
        }
        // select polynomials according to the main variable
        GenPolynomialRing> rfac = pfac.recursive(1);
        GenPolynomial Q = A.get(0); // wrong, must eventually search polynomial
        GenPolynomial> qr = PolyUtil. recursive(rfac, Q);
        GenPolynomial> pr = PolyUtil. recursive(rfac, P);
        GenPolynomial> rr;
        if (qr.isONE()) {
            return P.ring.getZERO();
        }
        if (qr.degree(0) > 0) {
            rr = PolyUtil. recursiveSparsePseudoRemainder(pr, qr);
            //System.out.println("remainder, pr = " + pr);
            //System.out.println("remainder, qr = " + qr);
            //System.out.println("remainder, rr = " + rr);
        } else {
            rr = pr;
        }
        if (rr.degree(0) > 0) {
            GenPolynomial R = PolyUtil. distribute(pfac, rr);
            return R.monic();
            // not further reduced wrt. other variables = top-reduction only
        }
        List> zeroDeg = zeroDegrees(A);
        GenPolynomial R = topPseudoRemainder(zeroDeg, rr.leadingBaseCoefficient());
        R = R.extend(pfac, 0, 0L);
        return R.monic();
    }


    /**
     * Top coefficient pseudo remainder of the leading coefficient of P wrt A in
     * the main variables.
     * @param P generic polynomial in n+1 variables.
     * @param A list of generic polynomials in n variables sorted according to
     *            appearing main variables.
     * @return pseudo remainder of the leading coefficient of P wrt A.
     */
    public static > GenPolynomial topCoefficientPseudoRemainder(
                    List> A, GenPolynomial P) {
        if (A == null || A.isEmpty()) {
            return P.monic();
        }
        if (P.isZERO()) {
            return P;
        }
        GenPolynomialRing pfac = P.ring;
        GenPolynomialRing pfac1 = A.get(0).ring;
        if (pfac1.nvar <= 1) { // recursion base 
            GenPolynomial a = A.get(0);
            GenPolynomialRing> rfac = pfac.recursive(pfac.nvar - 1);
            GenPolynomial> pr = PolyUtil. recursive(rfac, P);
            // ldcf(P,x_m) = q a + r 
            GenPolynomial> rr = PolyGBUtil. coefficientPseudoRemainderBase(pr, a);
            GenPolynomial R = PolyUtil. distribute(pfac, rr);
            return R.monic();
        }
        // select polynomials according to the main variable
        GenPolynomialRing> rfac1 = pfac1.recursive(1);
        int nv = pfac.nvar - pfac1.nvar;
        GenPolynomialRing> rfac = pfac.recursive(1 + nv);
        GenPolynomialRing>> rfac2 = rfac.recursive(nv);
        if (debug) {
            logger.info("rfac =" + rfac);
        }
        GenPolynomial> pr = PolyUtil. recursive(rfac, P);
        GenPolynomial>> pr2 = PolyUtil.> recursive(rfac2, pr);
        //System.out.println("recursion, pr2 = " + pr2);
        GenPolynomial Q = A.get(0);
        GenPolynomial> qr = PolyUtil. recursive(rfac1, Q);
        GenPolynomial>> rr;
        if (qr.isONE()) {
            return P.ring.getZERO();
        }
        if (qr.degree(0) > 0) {
            // pseudo remainder:  ldcf(P,x_m) = a q + r 
            rr = PolyGBUtil. coefficientPseudoRemainder(pr2, qr);
            //System.out.println("recursion, qr  = " + qr);
            //System.out.println("recursion, pr  = " + pr2);
            //System.out.println("recursion, rr  = " + rr);
        } else {
            rr = pr2;
        }
        // reduction wrt. the other variables
        List> zeroDeg = zeroDegrees(A);
        GenPolynomial> Rr = PolyUtil.> distribute(rfac, rr);
        GenPolynomial R = PolyUtil. distribute(pfac, Rr);
        R = topCoefficientPseudoRemainder(zeroDeg, R);
        return R.monic();
    }


    /**
     * Polynomial leading coefficient pseudo remainder.
     * @param P generic polynomial in n+1 variables.
     * @param A generic polynomial in n variables.
     * @return pseudo remainder of the leading coefficient of P wrt A, with
     *         ldcf(A)m' P = quotient * A + remainder.
     */
    public static > GenPolynomial>> coefficientPseudoRemainder(
                    GenPolynomial>> P, GenPolynomial> A) {
        if (A == null || A.isZERO()) { // findbugs
            throw new ArithmeticException(P + " division by zero " + A);
        }
        if (A.isONE()) {
            return P.ring.getZERO();
        }
        if (P.isZERO() || P.isONE()) {
            return P;
        }
        GenPolynomialRing>> pfac = P.ring;
        GenPolynomialRing> afac = A.ring; // == pfac.coFac
        GenPolynomial>> r = P;
        GenPolynomial> h;
        GenPolynomial>> hr;
        GenPolynomial> ldcf = P.leadingBaseCoefficient();
        long m = ldcf.degree(0);
        long n = A.degree(0);
        GenPolynomial c = A.leadingBaseCoefficient();
        GenPolynomial> cc = afac.getZERO().sum(c);
        //System.out.println("cc = " + cc);
        ExpVector e = A.leadingExpVector();
        for (long i = m; i >= n; i--) {
            if (r.isZERO()) {
                return r;
            }
            GenPolynomial> p = r.leadingBaseCoefficient();
            ExpVector g = r.leadingExpVector();
            long k = p.degree(0);
            if (i == k) {
                GenPolynomial pl = p.leadingBaseCoefficient();
                ExpVector f = p.leadingExpVector();
                f = f.subtract(e);
                r = r.multiply(cc); // coeff cc
                h = A.multiply(pl, f); // coeff ac
                hr = new GenPolynomial>>(pfac, h, g);
                r = r.subtract(hr);
            } else {
                r = r.multiply(cc);
            }
            //System.out.println("r = " + r);
        }
        if (r.degree(0) < P.degree(0)) { // recursion for degree
            r = coefficientPseudoRemainder(r, A);
        }
        return r;
    }


    /**
     * Polynomial leading coefficient pseudo remainder, base case.
     * @param P generic polynomial in 1+1 variables.
     * @param A generic polynomial in 1 variable.
     * @return pseudo remainder of the leading coefficient of P wrt. A, with
     *         ldcf(A)m' P = quotient * A + remainder.
     */
    public static > GenPolynomial> coefficientPseudoRemainderBase(
                    GenPolynomial> P, GenPolynomial A) {
        if (A == null || A.isZERO()) { // findbugs
            throw new ArithmeticException(P + " division by zero " + A);
        }
        if (A.isONE()) {
            return P.ring.getZERO();
        }
        if (P.isZERO() || P.isONE()) {
            return P;
        }
        GenPolynomialRing> pfac = P.ring;
        GenPolynomialRing afac = A.ring; // == pfac.coFac
        GenPolynomial> r = P;
        GenPolynomial h;
        GenPolynomial> hr;
        GenPolynomial ldcf = P.leadingBaseCoefficient();
        long m = ldcf.degree(0);
        long n = A.degree(0);
        C c = A.leadingBaseCoefficient();
        GenPolynomial cc = afac.getZERO().sum(c);
        //System.out.println("cc = " + cc);
        ExpVector e = A.leadingExpVector();
        for (long i = m; i >= n; i--) {
            if (r.isZERO()) {
                return r;
            }
            GenPolynomial p = r.leadingBaseCoefficient();
            ExpVector g = r.leadingExpVector();
            long k = p.degree(0);
            if (i == k) {
                C pl = p.leadingBaseCoefficient();
                ExpVector f = p.leadingExpVector();
                f = f.subtract(e);
                r = r.multiply(cc); // coeff cc
                h = A.multiply(pl, f); // coeff ac
                hr = new GenPolynomial>(pfac, h, g);
                r = r.subtract(hr);
            } else {
                r = r.multiply(cc);
            }
            //System.out.println("r = " + r);
        }
        if (r.degree(0) < P.degree(0)) { // recursion for degree
            r = coefficientPseudoRemainderBase(r, A);
        }
        return r;
    }


    /**
     * Extract polynomials with degree zero in the main variable.
     * @param A list of generic polynomials in n variables.
     * @return Z = [a_i] with deg(a_i,x_n) = 0 and in n-1 variables.
     */
    public static > List> zeroDegrees(List> A) {
        if (A == null || A.isEmpty()) {
            return A;
        }
        GenPolynomialRing pfac = A.get(0).ring;
        GenPolynomialRing> rfac = pfac.recursive(1);
        List> zeroDeg = new ArrayList>(A.size());
        for (int i = 0; i < A.size(); i++) {
            GenPolynomial q = A.get(i);
            GenPolynomial> fr = PolyUtil. recursive(rfac, q);
            if (fr.degree(0) == 0) {
                zeroDeg.add(fr.leadingBaseCoefficient());
            }
        }
        return zeroDeg;
    }


    /**
     * Intersection. Generators for the intersection of ideals.
     * @param pfac polynomial ring
     * @param A list of polynomials
     * @param B list of polynomials
     * @return generators for (A \cap B)
     */
    public static > List> intersect(GenPolynomialRing pfac,
                    List> A, List> B) {
        if (A == null || A.isEmpty()) { // (0)
            return B;
        }
        if (B == null || B.isEmpty()) { // (0)
            return A;
        }
        int s = A.size() + B.size();
        List> c = new ArrayList>(s);
        GenPolynomialRing tfac = pfac.extend(1);
        // term order is also adjusted
        for (GenPolynomial p : A) {
            p = p.extend(tfac, 0, 1L); // t*p
            c.add(p);
        }
        for (GenPolynomial p : B) {
            GenPolynomial q = p.extend(tfac, 0, 1L);
            GenPolynomial r = p.extend(tfac, 0, 0L);
            p = r.subtract(q); // (1-t)*p
            c.add(p);
        }
        GroebnerBaseAbstract bb = GBFactory. getImplementation(tfac.coFac);
        logger.warn("intersect computing GB");
        List> G = bb.GB(c);
        if (debug) {
            logger.debug("intersect GB = " + G);
        }
        List> I = PolyUtil. intersect(pfac, G);
        return I;
    }


    /**
     * Intersection. Generators for the intersection of ideals.
     * @param pfac solvable polynomial ring
     * @param A list of polynomials
     * @param B list of polynomials
     * @return generators for (A \cap B)
     */
    public static > List> intersect(
                    GenSolvablePolynomialRing pfac, List> A,
                    List> B) {
        if (A == null || A.isEmpty()) { // (0)
            return B;
        }
        if (B == null || B.isEmpty()) { // (0)
            return A;
        }
        int s = A.size() + B.size();
        List> c = new ArrayList>(s);
        GenSolvablePolynomialRing tfac = pfac.extend(1);
        // term order is also adjusted
        for (GenSolvablePolynomial p : A) {
            p = (GenSolvablePolynomial) p.extend(tfac, 0, 1L); // t*p
            c.add(p);
        }
        for (GenSolvablePolynomial p : B) {
            GenSolvablePolynomial q = (GenSolvablePolynomial) p.extend(tfac, 0, 1L);
            GenSolvablePolynomial r = (GenSolvablePolynomial) p.extend(tfac, 0, 0L);
            p = (GenSolvablePolynomial) r.subtract(q); // (1-t)*p
            c.add(p);
        }
        SolvableGroebnerBaseAbstract sbb = new SolvableGroebnerBaseSeq();
        //GBFactory. getImplementation(tfac.coFac);
        logger.warn("intersect computing GB");
        List> g = sbb.leftGB(c);
        if (debug) {
            logger.debug("intersect GB = " + g);
        }
        List> I = PolyUtil. intersect(pfac, g);
        return I;
    }


    /**
     * Least common multiple via ideal intersection.
     * @param r solvable polynomial ring.
     * @param n first solvable polynomial.
     * @param d second solvable polynomial.
     * @return lcm(n,d)
     */
    public static > GenSolvablePolynomial syzLcm(GenSolvablePolynomialRing r,
                    GenSolvablePolynomial n, GenSolvablePolynomial d) {
        if (n.isZERO()) {
            return n;
        }
        if (d.isZERO()) {
            return d;
        }
        if (n.isONE()) {
            return d;
        }
        if (d.isONE()) {
            return n;
        }
        List> A = new ArrayList>(1);
        A.add(n);
        List> B = new ArrayList>(1);
        B.add(d);
        List> c = PolyGBUtil. intersect(r, A, B);
        //if (c.size() != 1) {
        // SolvableSyzygyAbstract sz = new SolvableSyzygyAbstract();
        // GenSolvablePolynomial[] oc = sz.leftOreCond(n,d);
        // GenSolvablePolynomial nc = oc[0].multiply(n);
        // System.out.println("nc = " + nc);
        // return nc;
        //}
        GenSolvablePolynomial lcm = null;
        for (GenSolvablePolynomial p : c) {
            if (p == null || p.isZERO()) {
                continue;
            }
            //System.out.println("p = " + p);
            if (lcm == null) {
                lcm = p;
                continue;
            }
            if (lcm.compareTo(p) > 0) {
                lcm = p;
            }
        }
        if (lcm == null) {
            throw new RuntimeException("this cannot happen: lcm == null: " + c);
        }
        return lcm;
    }


    /**
     * Greatest common divisor via least common multiple.
     * @param r solvable polynomial ring.
     * @param n first solvable polynomial.
     * @param d second solvable polynomial.
     * @return gcd(n,d)
     */
    public static > GenSolvablePolynomial syzGcd(GenSolvablePolynomialRing r,
                    GenSolvablePolynomial n, GenSolvablePolynomial d) {
        if (n.isZERO()) {
            return d;
        }
        if (d.isZERO()) {
            return n;
        }
        if (n.isONE()) {
            return n;
        }
        if (d.isONE()) {
            return d;
        }
        if (n.totalDegree() > 3 || d.totalDegree() > 3) { // how avoid too long running GBs ?
        //if (n.totalDegree() + d.totalDegree() > 6) { // how avoid too long running GBs ?
            // && n.length() < 10 && d.length() < 10
            logger.warn("skipping GB computation: degs = " + n.totalDegree() + ", " + d.totalDegree());
            return r.getONE();
        }
        List> A = new ArrayList>(2);
        A.add(n);
        A.add(d);
        SolvableGroebnerBaseAbstract sbb = new SolvableGroebnerBaseSeq();
        logger.warn("syzGcd computing GB: " + A);
        List> G = sbb.rightGB(A); //leftGB, not: sbb.twosidedGB(A);
        if ( logger.isInfoEnabled() ) { 
            logger.info("G = " + G);
        }
        if (G.size() == 1) {
            return G.get(0);
        }
        logger.warn("gcd not determined, set to 1: " + G); // + ", A = " + A);
        return r.getONE();
    }


    /**
     * Solvable quotient and remainder via reduction.
     * @param n first solvable polynomial.
     * @param d second solvable polynomial.
     * @return [ n/d, n - (n/d)*d ]
     */
    public static > GenSolvablePolynomial[] quotientRemainder(
                    GenSolvablePolynomial n, GenSolvablePolynomial d) {
        GenSolvablePolynomial[] res = (GenSolvablePolynomial[]) new GenSolvablePolynomial[2];
        if (d.isZERO()) {
            throw new RuntimeException("division by zero: " + n + "/" + d);
        }
        if (n.isZERO()) {
            res[0] = n;
            res[1] = n;
            return res;
        }
        GenSolvablePolynomialRing r = n.ring;
        if (d.isONE()) {
            res[0] = n;
            res[1] = r.getZERO();
            return res;
        }
        // divide
        List> Q = new ArrayList>(1);
        Q.add(r.getZERO());
        List> D = new ArrayList>(1);
        D.add(d);
        SolvableReductionAbstract sred = new SolvableReductionSeq();
        res[1] = sred.rightNormalform(Q, D, n); // left
        res[0] = Q.get(0);
        return res;
    }


    /**
     * Greatest common divisor and cofactors via least common multiple and
     * reduction.
     * @param r solvable polynomial ring.
     * @param n first solvable polynomial.
     * @param d second solvable polynomial.
     * @return [ g=gcd(n,d), n/g, d/g ]
     */
    public static > GenSolvablePolynomial[] syzGcdCofactors(
                    GenSolvablePolynomialRing r, GenSolvablePolynomial n, GenSolvablePolynomial d) {
        GenSolvablePolynomial[] res = (GenSolvablePolynomial[]) new GenSolvablePolynomial[3];
        res[0] = PolyGBUtil. syzGcd(r, n, d);
        res[1] = n;
        res[2] = d;
        if ( res[0].isONE() ) {
            return res;
        }
        GenSolvablePolynomial[] nqr = PolyGBUtil. quotientRemainder(n, res[0]);
        if ( !nqr[1].isZERO() ) {
            res[0] = r.getONE();
            return res;
	}
        GenSolvablePolynomial[] dqr = PolyGBUtil. quotientRemainder(d, res[0]);
        if ( !dqr[1].isZERO() ) {
            res[0] = r.getONE();
            return res;
	}
        res[1] = nqr[0];
        res[2] = dqr[0];
        return res;
    }


    /**
     * Least common multiple. Just for fun, is not efficient.
     * @param r polynomial ring.
     * @param n first polynomial.
     * @param d second polynomial.
     * @return lcm(n,d)
     */
    public static > GenPolynomial syzLcm(GenPolynomialRing r,
                    GenPolynomial n, GenPolynomial d) {
        List> A = new ArrayList>(1);
        A.add(n);
        List> B = new ArrayList>(1);
        B.add(d);
        List> c = PolyGBUtil. intersect(r, A, B);
        if (c.size() != 1) {
            logger.warn("lcm not uniqe: " + c);
            //throw new RuntimeException("lcm not uniqe: " + c);
        }
        GenPolynomial lcm = c.get(0);
        return lcm;
    }


    /**
     * Greatest common divisor. Just for fun, is not efficient.
     * @param r polynomial ring.
     * @param n first polynomial.
     * @param d second polynomial.
     * @return gcd(n,d)
     */
    public static > GenPolynomial syzGcd(GenPolynomialRing r,
                    GenPolynomial n, GenPolynomial d) {
        if (n.isZERO()) {
            return d;
        }
        if (d.isZERO()) {
            return n;
        }
        if (n.isONE()) {
            return n;
        }
        if (d.isONE()) {
            return d;
        }
        GenPolynomial p = n.multiply(d);
        GenPolynomial lcm = syzLcm(r, n, d);
        GenPolynomial gcd = PolyUtil. basePseudoDivide(p, lcm);
        return gcd;
    }

}