/*
 * $Id$
 */

package edu.jas.ufd;


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

import org.apache.log4j.Logger;

import edu.jas.arith.BigInteger;
import edu.jas.arith.ModInteger;
import edu.jas.arith.ModIntegerRing;
import edu.jas.arith.ModLongRing;
import edu.jas.arith.Modular;
import edu.jas.arith.ModularRingFactory;
import edu.jas.poly.ExpVector;
import edu.jas.poly.GenPolynomial;
import edu.jas.poly.GenPolynomialRing;
import edu.jas.poly.Monomial;
import edu.jas.poly.PolyUtil;
import edu.jas.structure.GcdRingElem;
import edu.jas.structure.RingFactory;
import edu.jas.structure.Power;


/**
 * Hensel utilities for ufd.
 * @author Heinz Kredel
 */

public class HenselUtil {


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


    private static final boolean debug = logger.isDebugEnabled();


    /**
     * Modular Hensel lifting algorithm on coefficients. Let p =
     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
     * with gcd(A,B) == 1 mod p and S A + T B == 1 mod p. See Algorithm 6.1. in
     * Geddes et.al.. Linear version, as it does not lift S A + T B == 1 mod
     * p^{e+1}.
     * @param C GenPolynomial
     * @param A GenPolynomial
     * @param B other GenPolynomial
     * @param S GenPolynomial
     * @param T GenPolynomial
     * @param M bound on the coefficients of A1 and B1 as factors of C.
     * @return [A1,B1,Am,Bm] = lift(C,A,B), with C = A1 * B1 mod p^e, Am = A1
     *         mod p^e, Bm = B1 mod p^e .
     */
    @SuppressWarnings("unchecked")
    public static  & Modular> HenselApprox liftHensel(
            GenPolynomial C, BigInteger M, GenPolynomial A, GenPolynomial B,
            GenPolynomial S, GenPolynomial T) throws NoLiftingException {
        if (C == null || C.isZERO()) {
            return new HenselApprox(C, C, A, B);
        }
        if (A == null || A.isZERO() || B == null || B.isZERO()) {
            throw new IllegalArgumentException("A and B must be nonzero");
        }
        GenPolynomialRing fac = C.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        // setup factories
        GenPolynomialRing pfac = A.ring;
        RingFactory p = pfac.coFac;
        RingFactory q = p;
        ModularRingFactory P = (ModularRingFactory) p;
        ModularRingFactory Q = (ModularRingFactory) q;
        BigInteger Qi = Q.getIntegerModul();
        BigInteger M2 = M.multiply(M.fromInteger(2));
        BigInteger Mq = Qi;

        // normalize c and a, b factors, assert p is prime
        GenPolynomial Ai;
        GenPolynomial Bi;
        BigInteger c = C.leadingBaseCoefficient();
        C = C.multiply(c); // sic
        MOD a = A.leadingBaseCoefficient();
        if (!a.isONE()) { // A = A.monic();
            A = A.divide(a);
            S = S.multiply(a);
        }
        MOD b = B.leadingBaseCoefficient();
        if (!b.isONE()) { // B = B.monic();
            B = B.divide(b);
            T = T.multiply(b);
        }
        MOD cm = P.fromInteger(c.getVal());
        A = A.multiply(cm);
        B = B.multiply(cm);
        T = T.divide(cm);
        S = S.divide(cm);
        Ai = PolyUtil.integerFromModularCoefficients(fac, A);
        Bi = PolyUtil.integerFromModularCoefficients(fac, B);
        // replace leading base coefficients
        ExpVector ea = Ai.leadingExpVector();
        ExpVector eb = Bi.leadingExpVector();
        Ai.doPutToMap(ea, c);
        Bi.doPutToMap(eb, c);

        // polynomials mod p
        GenPolynomial Ap;
        GenPolynomial Bp;
        GenPolynomial A1p = A;
        GenPolynomial B1p = B;
        GenPolynomial Ep;

        // polynomials over the integers
        GenPolynomial E;
        GenPolynomial Ea;
        GenPolynomial Eb;
        GenPolynomial Ea1;
        GenPolynomial Eb1;

        while (Mq.compareTo(M2) < 0) {
            // compute E=(C-AB)/q over the integers
            E = C.subtract(Ai.multiply(Bi));
            if (E.isZERO()) {
                logger.info("leaving on zero E");
                break;
            }
            try {
                E = E.divide(Qi);
            } catch (RuntimeException e) {
                // useful in debuging
                //System.out.println("C  = " + C );
                //System.out.println("Ai = " + Ai );
                //System.out.println("Bi = " + Bi );
                //System.out.println("E  = " + E );
                //System.out.println("Qi = " + Qi );
                throw e;
            }
            // E mod p
            Ep = PolyUtil. fromIntegerCoefficients(pfac, E);
            //logger.info("Ep = " + Ep);

            // construct approximation mod p
            Ap = S.multiply(Ep); // S,T ++ T,S
            Bp = T.multiply(Ep);
            GenPolynomial[] QR;
            QR = Ap.quotientRemainder(B);
            GenPolynomial Qp;
            GenPolynomial Rp;
            Qp = QR[0];
            Rp = QR[1];
            A1p = Rp;
            B1p = Bp.sum(A.multiply(Qp));

            // construct q-adic approximation, convert to integer
            Ea = PolyUtil.integerFromModularCoefficients(fac, A1p);
            Eb = PolyUtil.integerFromModularCoefficients(fac, B1p);
            Ea1 = Ea.multiply(Qi);
            Eb1 = Eb.multiply(Qi);

            Ea = Ai.sum(Eb1); // Eb1 and Ea1 are required
            Eb = Bi.sum(Ea1); //--------------------------
            assert (Ea.degree(0) + Eb.degree(0) <= C.degree(0));
            //if ( Ea.degree(0)+Eb.degree(0) > C.degree(0) ) { // debug
            //   throw new RuntimeException("deg(A)+deg(B) > deg(C)");
            //}

            // prepare for next iteration
            Mq = Qi;
            Qi = Q.getIntegerModul().multiply(P.getIntegerModul());
            // Q = new ModIntegerRing(Qi.getVal());
            if (ModLongRing.MAX_LONG.compareTo(Qi.getVal()) > 0) {
                Q = (ModularRingFactory) new ModLongRing(Qi.getVal());
            } else {
                Q = (ModularRingFactory) new ModIntegerRing(Qi.getVal());
            }
            Ai = Ea;
            Bi = Eb;
        }
        GreatestCommonDivisorAbstract ufd = new GreatestCommonDivisorPrimitive();

        // remove normalization
        BigInteger ai = ufd.baseContent(Ai);
        Ai = Ai.divide(ai);
        BigInteger bi = null;
        try {
            bi = c.divide(ai);
            Bi = Bi.divide(bi); // divide( c/a )
        } catch (RuntimeException e) {
            //System.out.println("C  = " + C );
            //System.out.println("Ai = " + Ai );
            //System.out.println("Bi = " + Bi );
            //System.out.println("c  = " + c );
            //System.out.println("ai = " + ai );
            //System.out.println("bi = " + bi );
            //System.out.println("no exact lifting possible " + e);
            throw new NoLiftingException("no exact lifting possible " + e);
        }
        return new HenselApprox(Ai, Bi, A1p, B1p);
    }


    /**
     * Modular Hensel lifting algorithm on coefficients. Let p =
     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
     * with gcd(A,B) == 1 mod p. See algorithm 6.1. in Geddes et.al. and
     * algorithms 3.5.{5,6} in Cohen.
     * @param C GenPolynomial
     * @param A GenPolynomial
     * @param B other GenPolynomial
     * @param M bound on the coefficients of A1 and B1 as factors of C.
     * @return [A1,B1] = lift(C,A,B), with C = A1 * B1.
     */
    @SuppressWarnings("unchecked")
    public static  & Modular> HenselApprox liftHensel(
            GenPolynomial C, BigInteger M, GenPolynomial A, GenPolynomial B)
            throws NoLiftingException {
        if (C == null || C.isZERO()) {
            return new HenselApprox(C, C, A, B);
        }
        if (A == null || A.isZERO() || B == null || B.isZERO()) {
            throw new IllegalArgumentException("A and B must be nonzero");
        }
        GenPolynomialRing fac = C.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        // one Hensel step on part polynomials
        try {
            GenPolynomial[] gst = A.egcd(B);
            if (!gst[0].isONE()) {
                throw new NoLiftingException("A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B);
            }
            GenPolynomial s = gst[1];
            GenPolynomial t = gst[2];
            HenselApprox ab = HenselUtil. liftHensel(C, M, A, B, s, t);
            return ab;
        } catch (ArithmeticException e) {
            throw new NoLiftingException("coefficient error " + e);
        }
    }


    /**
     * Modular quadratic Hensel lifting algorithm on coefficients. Let p =
     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
     * with gcd(A,B) == 1 mod p and S A + T B == 1 mod p. See algorithm 6.1. in
     * Geddes et.al. and algorithms 3.5.{5,6} in Cohen. Quadratic version, as it
     * also lifts S A + T B == 1 mod p^{e+1}.
     * @param C GenPolynomial
     * @param A GenPolynomial
     * @param B other GenPolynomial
     * @param S GenPolynomial
     * @param T GenPolynomial
     * @param M bound on the coefficients of A1 and B1 as factors of C.
     * @return [A1,B1] = lift(C,A,B), with C = A1 * B1.
     */
    @SuppressWarnings("unchecked")
    public static  & Modular> HenselApprox liftHenselQuadratic(
            GenPolynomial C, BigInteger M, GenPolynomial A, GenPolynomial B,
            GenPolynomial S, GenPolynomial T) throws NoLiftingException {
        if (C == null || C.isZERO()) {
            return new HenselApprox(C, C, A, B);
        }
        if (A == null || A.isZERO() || B == null || B.isZERO()) {
            throw new IllegalArgumentException("A and B must be nonzero");
        }
        GenPolynomialRing fac = C.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        // setup factories
        GenPolynomialRing pfac = A.ring;
        RingFactory p = pfac.coFac;
        RingFactory q = p;
        ModularRingFactory P = (ModularRingFactory) p;
        ModularRingFactory Q = (ModularRingFactory) q;
        BigInteger Qi = Q.getIntegerModul();
        BigInteger M2 = M.multiply(M.fromInteger(2));
        BigInteger Mq = Qi;
        GenPolynomialRing qfac;
        qfac = new GenPolynomialRing(Q, pfac);

        // normalize c and a, b factors, assert p is prime
        GenPolynomial Ai;
        GenPolynomial Bi;
        BigInteger c = C.leadingBaseCoefficient();
        C = C.multiply(c); // sic
        MOD a = A.leadingBaseCoefficient();
        if (!a.isONE()) { // A = A.monic();
            A = A.divide(a);
            S = S.multiply(a);
        }
        MOD b = B.leadingBaseCoefficient();
        if (!b.isONE()) { // B = B.monic();
            B = B.divide(b);
            T = T.multiply(b);
        }
        MOD cm = P.fromInteger(c.getVal());
        A = A.multiply(cm);
        B = B.multiply(cm);
        T = T.divide(cm);
        S = S.divide(cm);
        Ai = PolyUtil.integerFromModularCoefficients(fac, A);
        Bi = PolyUtil.integerFromModularCoefficients(fac, B);
        // replace leading base coefficients
        ExpVector ea = Ai.leadingExpVector();
        ExpVector eb = Bi.leadingExpVector();
        Ai.doPutToMap(ea, c);
        Bi.doPutToMap(eb, c);

        // polynomials mod p
        GenPolynomial Ap;
        GenPolynomial Bp;
        GenPolynomial A1p = A;
        GenPolynomial B1p = B;
        GenPolynomial Ep;
        GenPolynomial Sp = S;
        GenPolynomial Tp = T;

        // polynomials mod q
        GenPolynomial Aq;
        GenPolynomial Bq;
        GenPolynomial Eq;

        // polynomials over the integers
        GenPolynomial E;
        GenPolynomial Ea;
        GenPolynomial Eb;
        GenPolynomial Ea1;
        GenPolynomial Eb1;
        GenPolynomial Si;
        GenPolynomial Ti;

        Si = PolyUtil.integerFromModularCoefficients(fac, S);
        Ti = PolyUtil.integerFromModularCoefficients(fac, T);

        Aq = PolyUtil. fromIntegerCoefficients(qfac, Ai);
        Bq = PolyUtil. fromIntegerCoefficients(qfac, Bi);

        while (Mq.compareTo(M2) < 0) {
            // compute E=(C-AB)/q over the integers
            E = C.subtract(Ai.multiply(Bi));
            if (E.isZERO()) {
                logger.info("leaving on zero E");
                break;
            }
            E = E.divide(Qi);
            // E mod p
            Ep = PolyUtil. fromIntegerCoefficients(qfac, E);
            //logger.info("Ep = " + Ep + ", qfac = " + qfac);
            //if (Ep.isZERO()) {
                //System.out.println("leaving on zero error");
                //??logger.info("leaving on zero Ep");
                //??break;
            //}

            // construct approximation mod p
            Ap = Sp.multiply(Ep); // S,T ++ T,S
            Bp = Tp.multiply(Ep);
            GenPolynomial[] QR;
            //logger.info("Ap = " + Ap + ", Bp = " + Bp + ", fac(Ap) = " + Ap.ring);
            QR = Ap.quotientRemainder(Bq);
            GenPolynomial Qp;
            GenPolynomial Rp;
            Qp = QR[0];
            Rp = QR[1];
            //logger.info("Qp = " + Qp + ", Rp = " + Rp);
            A1p = Rp;
            B1p = Bp.sum(Aq.multiply(Qp));

            // construct q-adic approximation, convert to integer
            Ea = PolyUtil.integerFromModularCoefficients(fac, A1p);
            Eb = PolyUtil.integerFromModularCoefficients(fac, B1p);
            Ea1 = Ea.multiply(Qi);
            Eb1 = Eb.multiply(Qi);
            Ea = Ai.sum(Eb1); // Eb1 and Ea1 are required
            Eb = Bi.sum(Ea1); //--------------------------
            assert (Ea.degree(0) + Eb.degree(0) <= C.degree(0));
            //if ( Ea.degree(0)+Eb.degree(0) > C.degree(0) ) { // debug
            //   throw new RuntimeException("deg(A)+deg(B) > deg(C)");
            //}
            Ai = Ea;
            Bi = Eb;

            // gcd representation factors error --------------------------------
            // compute E=(1-SA-TB)/q over the integers
            E = fac.getONE();
            E = E.subtract(Si.multiply(Ai)).subtract(Ti.multiply(Bi));
            E = E.divide(Qi);
            // E mod q
            Ep = PolyUtil. fromIntegerCoefficients(qfac, E);
            //logger.info("Ep2 = " + Ep);

            // construct approximation mod q
            Ap = Sp.multiply(Ep); // S,T ++ T,S
            Bp = Tp.multiply(Ep);
            QR = Bp.quotientRemainder(Aq); // Ai == A mod p ?
            Qp = QR[0];
            Rp = QR[1];
            B1p = Rp;
            A1p = Ap.sum(Bq.multiply(Qp));

            if (false && debug) {
                Eq = A1p.multiply(Aq).sum(B1p.multiply(Bq)).subtract(Ep);
                if (!Eq.isZERO()) {
                    System.out.println("A*A1p+B*B1p-Ep2 != 0 " + Eq);
                    throw new RuntimeException("A*A1p+B*B1p-Ep2 != 0 mod " + Q.getIntegerModul());
                }
            }

            // construct q-adic approximation, convert to integer
            Ea = PolyUtil.integerFromModularCoefficients(fac, A1p);
            Eb = PolyUtil.integerFromModularCoefficients(fac, B1p);
            Ea1 = Ea.multiply(Qi);
            Eb1 = Eb.multiply(Qi);
            Ea = Si.sum(Ea1); // Eb1 and Ea1 are required
            Eb = Ti.sum(Eb1); //--------------------------
            Si = Ea;
            Ti = Eb;

            // prepare for next iteration
            Mq = Qi;
            Qi = Q.getIntegerModul().multiply(Q.getIntegerModul());
            if (ModLongRing.MAX_LONG.compareTo(Qi.getVal()) > 0) {
                Q = (ModularRingFactory) new ModLongRing(Qi.getVal());
            } else {
                Q = (ModularRingFactory) new ModIntegerRing(Qi.getVal());
            }
            //Q = new ModIntegerRing(Qi.getVal());
            //System.out.println("Q = " + Q + ", from Q = " + Mq);

            qfac = new GenPolynomialRing(Q, pfac);

            Aq = PolyUtil. fromIntegerCoefficients(qfac, Ai);
            Bq = PolyUtil. fromIntegerCoefficients(qfac, Bi);
            Sp = PolyUtil. fromIntegerCoefficients(qfac, Si);
            Tp = PolyUtil. fromIntegerCoefficients(qfac, Ti);
            if (false && debug) {
                E = Ai.multiply(Si).sum(Bi.multiply(Ti));
                Eq = PolyUtil. fromIntegerCoefficients(qfac, E);
                if (!Eq.isONE()) {
                    System.out.println("Ai*Si+Bi*Ti=1 " + Eq);
                    throw new RuntimeException("Ai*Si+Bi*Ti != 1 mod " + Q.getIntegerModul());
                }
            }
        }
        GreatestCommonDivisorAbstract ufd = new GreatestCommonDivisorPrimitive();

        // remove normalization if possible
        BigInteger ai = ufd.baseContent(Ai);
        Ai = Ai.divide(ai);
        BigInteger bi = null;
        try {
            bi = c.divide(ai);
            Bi = Bi.divide(bi); // divide( c/a )
        } catch (RuntimeException e) {
            //System.out.println("C  = " + C );
            //System.out.println("Ai = " + Ai );
            //System.out.println("Bi = " + Bi );
            //System.out.println("c  = " + c );
            //System.out.println("ai = " + ai );
            //System.out.println("bi = " + bi );
            //System.out.println("no exact lifting possible " + e);
            throw new NoLiftingException("no exact lifting possible " + e);
        }
        return new HenselApprox(Ai, Bi, A1p, B1p);
    }


    /**
     * Modular quadratic Hensel lifting algorithm on coefficients. Let p =
     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
     * with gcd(A,B) == 1 mod p. See algorithm 6.1. in Geddes et.al. and
     * algorithms 3.5.{5,6} in Cohen. Quadratic version.
     * @param C GenPolynomial
     * @param A GenPolynomial
     * @param B other GenPolynomial
     * @param M bound on the coefficients of A1 and B1 as factors of C.
     * @return [A1,B1] = lift(C,A,B), with C = A1 * B1.
     */
    @SuppressWarnings("unchecked")
    public static  & Modular> HenselApprox liftHenselQuadratic(
            GenPolynomial C, BigInteger M, GenPolynomial A, GenPolynomial B)
            throws NoLiftingException {
        if (C == null || C.isZERO()) {
            return new HenselApprox(C, C, A, B);
        }
        if (A == null || A.isZERO() || B == null || B.isZERO()) {
            throw new IllegalArgumentException("A and B must be nonzero");
        }
        GenPolynomialRing fac = C.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        // one Hensel step on part polynomials
        try {
            GenPolynomial[] gst = A.egcd(B);
            if (!gst[0].isONE()) {
                throw new NoLiftingException("A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B);
            }
            GenPolynomial s = gst[1];
            GenPolynomial t = gst[2];
            HenselApprox ab = HenselUtil. liftHenselQuadratic(C, M, A, B, s, t);
            return ab;
        } catch (ArithmeticException e) {
            throw new NoLiftingException("coefficient error " + e);
        }
    }


    /**
     * Modular Hensel lifting algorithm on coefficients. Let p =
     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
     * with gcd(A,B) == 1 mod p. See algorithm 6.1. in Geddes et.al. and
     * algorithms 3.5.{5,6} in Cohen. Quadratic version.
     * @param C GenPolynomial
     * @param A GenPolynomial
     * @param B other GenPolynomial
     * @param M bound on the coefficients of A1 and B1 as factors of C.
     * @return [A1,B1] = lift(C,A,B), with C = A1 * B1.
     */
    @SuppressWarnings("unchecked")
    public static  & Modular> HenselApprox liftHenselQuadraticFac(
            GenPolynomial C, BigInteger M, GenPolynomial A, GenPolynomial B)
            throws NoLiftingException {
        if (C == null || C.isZERO()) {
            throw new IllegalArgumentException("C must be nonzero");
        }
        if (A == null || A.isZERO() || B == null || B.isZERO()) {
            throw new IllegalArgumentException("A and B must be nonzero");
        }
        GenPolynomialRing fac = C.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        // one Hensel step on part polynomials
        try {
            GenPolynomial[] gst = A.egcd(B);
            if (!gst[0].isONE()) {
                throw new NoLiftingException("A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B);
            }
            GenPolynomial s = gst[1];
            GenPolynomial t = gst[2];
            HenselApprox ab = HenselUtil. liftHenselQuadraticFac(C, M, A, B, s, t);
            return ab;
        } catch (ArithmeticException e) {
            throw new NoLiftingException("coefficient error " + e);
        }
    }


    /**
     * Modular Hensel lifting algorithm on coefficients. Let p =
     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
     * with gcd(A,B) == 1 mod p and S A + T B == 1 mod p. See algorithm 6.1. in
     * Geddes et.al. and algorithms 3.5.{5,6} in Cohen. Quadratic version, as it
     * also lifts S A + T B == 1 mod p^{e+1}.
     * @param C primitive GenPolynomial
     * @param A GenPolynomial
     * @param B other GenPolynomial
     * @param S GenPolynomial
     * @param T GenPolynomial
     * @param M bound on the coefficients of A1 and B1 as factors of C.
     * @return [A1,B1] = lift(C,A,B), with C = A1 * B1.
     */
    @SuppressWarnings("unchecked")
    public static  & Modular> HenselApprox liftHenselQuadraticFac(
            GenPolynomial C, BigInteger M, GenPolynomial A, GenPolynomial B,
            GenPolynomial S, GenPolynomial T) throws NoLiftingException {
        //System.out.println("*** version for factorization *** ");
        //GenPolynomial[] AB = new GenPolynomial[2];
        if (C == null || C.isZERO()) {
            throw new IllegalArgumentException("C must be nonzero");
        }
        if (A == null || A.isZERO() || B == null || B.isZERO()) {
            throw new IllegalArgumentException("A and B must be nonzero");
        }
        GenPolynomialRing fac = C.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        // setup factories
        GenPolynomialRing pfac = A.ring;
        RingFactory p = pfac.coFac;
        RingFactory q = p;
        ModularRingFactory P = (ModularRingFactory) p;
        ModularRingFactory Q = (ModularRingFactory) q;
        BigInteger PP = Q.getIntegerModul();
        BigInteger Qi = PP;
        BigInteger M2 = M.multiply(M.fromInteger(2));
        if (debug) {
            //System.out.println("M2 =  " + M2);
            logger.debug("M2 =  " + M2);
        }
        BigInteger Mq = Qi;
        GenPolynomialRing qfac;
        qfac = new GenPolynomialRing(Q, pfac); // mod p
        GenPolynomialRing mfac;
        BigInteger Mi = Q.getIntegerModul().multiply(Q.getIntegerModul());
        ModularRingFactory Qmm;
        // = new ModIntegerRing(Mi.getVal());
        if (ModLongRing.MAX_LONG.compareTo(Mi.getVal()) > 0) {
            Qmm = (ModularRingFactory) new ModLongRing(Mi.getVal());
        } else {
            Qmm = (ModularRingFactory) new ModIntegerRing(Mi.getVal());
        }
        mfac = new GenPolynomialRing(Qmm, qfac); // mod p^e
        MOD Qm = Qmm.fromInteger(Qi.getVal());

        // partly normalize c and a, b factors, assert p is prime
        GenPolynomial Ai;
        GenPolynomial Bi;
        BigInteger c = C.leadingBaseCoefficient();
        C = C.multiply(c); // sic
        MOD a = A.leadingBaseCoefficient();
        if (!a.isONE()) { // A = A.monic();
            A = A.divide(a);
            S = S.multiply(a);
        }
        MOD b = B.leadingBaseCoefficient();
        if (!b.isONE()) { // B = B.monic();
            B = B.divide(b);
            T = T.multiply(b);
        }
        MOD cm = P.fromInteger(c.getVal());
        if (cm.isZERO()) {
            System.out.println("c =  " + c);
            System.out.println("P =  " + P);
            throw new ArithmeticException("c mod p == 0 not meaningful");
        }
        // mod p
        A = A.multiply(cm);
        S = S.divide(cm);
        B = B.multiply(cm);
        T = T.divide(cm);
        Ai = PolyUtil.integerFromModularCoefficients(fac, A);
        Bi = PolyUtil.integerFromModularCoefficients(fac, B);
        // replace leading base coefficients
        ExpVector ea = Ai.leadingExpVector();
        ExpVector eb = Bi.leadingExpVector();
        Ai.doPutToMap(ea, c);
        Bi.doPutToMap(eb, c);

        // polynomials mod p
        GenPolynomial Ap;
        GenPolynomial Bp;
        GenPolynomial A1p = A;
        GenPolynomial B1p = B;
        GenPolynomial Sp = S;
        GenPolynomial Tp = T;

        // polynomials mod q
        GenPolynomial Aq;
        GenPolynomial Bq;

        // polynomials mod p^e
        GenPolynomial Cm;
        GenPolynomial Am;
        GenPolynomial Bm;
        GenPolynomial Em;
        GenPolynomial Emp;
        GenPolynomial Sm;
        GenPolynomial Tm;
        GenPolynomial Ema;
        GenPolynomial Emb;
        GenPolynomial Ema1;
        GenPolynomial Emb1;

        // polynomials over the integers
        GenPolynomial Ei;
        GenPolynomial Si;
        GenPolynomial Ti;

        Si = PolyUtil.integerFromModularCoefficients(fac, S);
        Ti = PolyUtil.integerFromModularCoefficients(fac, T);

        Aq = PolyUtil. fromIntegerCoefficients(qfac, Ai);
        Bq = PolyUtil. fromIntegerCoefficients(qfac, Bi);

        // polynomials mod p^e
        Cm = PolyUtil. fromIntegerCoefficients(mfac, C);
        Am = PolyUtil. fromIntegerCoefficients(mfac, Ai);
        Bm = PolyUtil. fromIntegerCoefficients(mfac, Bi);
        Sm = PolyUtil. fromIntegerCoefficients(mfac, Si);
        Tm = PolyUtil. fromIntegerCoefficients(mfac, Ti);
        //System.out.println("Cm =  " + Cm);
        //System.out.println("Am =  " + Am);
        //System.out.println("Bm =  " + Bm);
        //System.out.println("Ai =  " + Ai);
        //System.out.println("Bi =  " + Bi);
        //System.out.println("mfac =  " + mfac);

        while (Mq.compareTo(M2) < 0) {
            // compute E=(C-AB)/p mod p^e
            if (debug) {
                //System.out.println("mfac =  " + Cm.ring);
                logger.debug("mfac =  " + Cm.ring);
            }
            Em = Cm.subtract(Am.multiply(Bm));
            //System.out.println("Em =  " + Em);
            if (Em.isZERO()) {
                if (C.subtract(Ai.multiply(Bi)).isZERO()) {
                    logger.info("leaving on zero E");
                    break;
                }
            }
            // Em = Em.divide( Qm );
            Ei = PolyUtil.integerFromModularCoefficients(fac, Em);
            Ei = Ei.divide(Qi);
            //System.out.println("Ei =  " + Ei);

            // Ei mod p
            Emp = PolyUtil. fromIntegerCoefficients(qfac, Ei);
            //            Emp = PolyUtil.fromIntegerCoefficients(qfac,
            //               PolyUtil.integerFromModularCoefficients(fac,Em) ); 
            //System.out.println("Emp =  " + Emp);
            //logger.info("Emp = " + Emp);
            if (Emp.isZERO()) {
                if (C.subtract(Ai.multiply(Bi)).isZERO()) {
                    logger.info("leaving on zero Emp");
                    break;
                }
            }

            // construct approximation mod p
            Ap = Sp.multiply(Emp); // S,T ++ T,S 
            Bp = Tp.multiply(Emp);
            GenPolynomial[] QR = null;
            QR = Ap.quotientRemainder(Bq); // Bq !
            GenPolynomial Qp = QR[0];
            GenPolynomial Rp = QR[1];
            A1p = Rp;
            B1p = Bp.sum(Aq.multiply(Qp)); // Aq !
            //System.out.println("A1p =  " + A1p);
            //System.out.println("B1p =  " + B1p);

            // construct q-adic approximation
            Ema = PolyUtil. fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac,
                    A1p));
            Emb = PolyUtil. fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac,
                    B1p));
            //System.out.println("Ema =  " + Ema);
            //System.out.println("Emb =  " + Emb);
            Ema1 = Ema.multiply(Qm);
            Emb1 = Emb.multiply(Qm);
            Ema = Am.sum(Emb1); // Eb1 and Ea1 are required
            Emb = Bm.sum(Ema1); //--------------------------
            assert (Ema.degree(0) + Emb.degree(0) <= Cm.degree(0));
            //if ( Ema.degree(0)+Emb.degree(0) > Cm.degree(0) ) { // debug
            //   throw new RuntimeException("deg(A)+deg(B) > deg(C)");
            //}
            Am = Ema;
            Bm = Emb;
            Ai = PolyUtil.integerFromModularCoefficients(fac, Am);
            Bi = PolyUtil.integerFromModularCoefficients(fac, Bm);
            //System.out.println("Am =  " + Am);
            //System.out.println("Bm =  " + Bm);
            //System.out.println("Ai =  " + Ai);
            //System.out.println("Bi =  " + Bi + "\n");

            // gcd representation factors error --------------------------------
            // compute E=(1-SA-TB)/p mod p^e
            Em = mfac.getONE();
            Em = Em.subtract(Sm.multiply(Am)).subtract(Tm.multiply(Bm));
            //System.out.println("Em  =  " + Em);
            // Em = Em.divide( Pm );

            Ei = PolyUtil.integerFromModularCoefficients(fac, Em);
            Ei = Ei.divide(Qi);
            //System.out.println("Ei =  " + Ei);
            // Ei mod p
            Emp = PolyUtil. fromIntegerCoefficients(qfac, Ei);
            //Emp = PolyUtil.fromIntegerCoefficients(qfac,
            //               PolyUtil.integerFromModularCoefficients(fac,Em) );
            //System.out.println("Emp =  " + Emp);

            // construct approximation mod p
            Ap = Sp.multiply(Emp); // S,T ++ T,S // Ep Eqp
            Bp = Tp.multiply(Emp); // Ep Eqp
            QR = Bp.quotientRemainder(Aq); // Ap Aq ! // Ai == A mod p ?
            Qp = QR[0];
            Rp = QR[1];
            B1p = Rp;
            A1p = Ap.sum(Bq.multiply(Qp));
            //System.out.println("A1p =  " + A1p);
            //System.out.println("B1p =  " + B1p);

            // construct p^e-adic approximation
            Ema = PolyUtil. fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac,
                    A1p));
            Emb = PolyUtil. fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac,
                    B1p));
            Ema1 = Ema.multiply(Qm);
            Emb1 = Emb.multiply(Qm);
            Ema = Sm.sum(Ema1); // Emb1 and Ema1 are required
            Emb = Tm.sum(Emb1); //--------------------------
            Sm = Ema;
            Tm = Emb;
            Si = PolyUtil.integerFromModularCoefficients(fac, Sm);
            Ti = PolyUtil.integerFromModularCoefficients(fac, Tm);
            //System.out.println("Sm =  " + Sm);
            //System.out.println("Tm =  " + Tm);
            //System.out.println("Si =  " + Si);
            //System.out.println("Ti =  " + Ti + "\n");

            // prepare for next iteration
            Qi = Q.getIntegerModul().multiply(Q.getIntegerModul());
            Mq = Qi;
            //lmfac = mfac;
            // Q = new ModIntegerRing(Qi.getVal());
            if (ModLongRing.MAX_LONG.compareTo(Qi.getVal()) > 0) {
                Q = (ModularRingFactory) new ModLongRing(Qi.getVal());
            } else {
                Q = (ModularRingFactory) new ModIntegerRing(Qi.getVal());
            }
            qfac = new GenPolynomialRing(Q, pfac);
            BigInteger Qmmi = Qmm.getIntegerModul().multiply(Qmm.getIntegerModul());
            //Qmm = new ModIntegerRing(Qmmi.getVal());
            if (ModLongRing.MAX_LONG.compareTo(Qmmi.getVal()) > 0) {
                Qmm = (ModularRingFactory) new ModLongRing(Qmmi.getVal());
            } else {
                Qmm = (ModularRingFactory) new ModIntegerRing(Qmmi.getVal());
            }
            mfac = new GenPolynomialRing(Qmm, qfac);
            Qm = Qmm.fromInteger(Qi.getVal());

            Cm = PolyUtil. fromIntegerCoefficients(mfac, C);
            Am = PolyUtil. fromIntegerCoefficients(mfac, Ai);
            Bm = PolyUtil. fromIntegerCoefficients(mfac, Bi);
            Sm = PolyUtil. fromIntegerCoefficients(mfac, Si);
            Tm = PolyUtil. fromIntegerCoefficients(mfac, Ti);

            assert (isHenselLift(C, Mi, PP, Ai, Bi));
            Mi = Mi.fromInteger(Qmm.getIntegerModul().getVal());

            Aq = PolyUtil. fromIntegerCoefficients(qfac, Ai);
            Bq = PolyUtil. fromIntegerCoefficients(qfac, Bi);
            Sp = PolyUtil. fromIntegerCoefficients(qfac, Si);
            Tp = PolyUtil. fromIntegerCoefficients(qfac, Ti);

            //System.out.println("Am =  " + Am);
            //System.out.println("Bm =  " + Bm);
            //System.out.println("Sm =  " + Sm);
            //System.out.println("Tm =  " + Tm);
            //System.out.println("mfac =  " + mfac);
            //System.out.println("Qmm = " + Qmm + ", M2   =  " + M2 + ", Mq   =  " + Mq + "\n");
        }
        //System.out.println("*Ai =  " + Ai);
        //System.out.println("*Bi =  " + Bi);

        Em = Cm.subtract(Am.multiply(Bm));
        if (!Em.isZERO()) {
            System.out.println("Em =  " + Em);
            //throw new NoLiftingException("no exact lifting possible");
        }
        // remove normalization not possible when not exact factorization
        GreatestCommonDivisorAbstract ufd = new GreatestCommonDivisorPrimitive();
        // remove normalization if possible
        BigInteger ai = ufd.baseContent(Ai);
        Ai = Ai.divide(ai); // Ai=pp(Ai)
        BigInteger[] qr = c.quotientRemainder(ai);
        BigInteger bi = null;
        boolean exact = true;
        if (qr[1].isZERO()) {
            bi = qr[0];
            try {
                Bi = Bi.divide(bi); // divide( c/a )
            } catch (RuntimeException e) {
                System.out.println("*catch: no exact factorization: " + bi + ", e = " + e);
                exact = false;
            }
        } else {
            System.out.println("*remainder: no exact factorization: q = " + qr[0] + ", r = " + qr[1]);
            exact = false;
        }
        if (!exact) {
            System.out.println("*Ai =  " + Ai);
            System.out.println("*ai =  " + ai);
            System.out.println("*Bi =  " + Bi);
            System.out.println("*bi =  " + bi);
            System.out.println("*c  =  " + c);
            throw new NoLiftingException("no exact lifting possible");
        }
        return new HenselApprox(Ai, Bi, Aq, Bq);
    }


    /**
     * Modular Hensel lifting test. Let p be a prime number and assume C ==
     * prod_{0,...,n-1} g_i mod p with gcd(g_i,g_j) == 1 mod p for i != j.
     * @param C GenPolynomial
     * @param G = [g_0,...,g_{n-1}] list of GenPolynomial
     * @param M bound on the coefficients of g_i as factors of C.
     * @param p prime number.
     * @return true if C = prod_{0,...,n-1} g_i mod p^e, else false.
     */
    public static// & Modular> 
    boolean isHenselLift(GenPolynomial C, BigInteger M, BigInteger p,
            List> G) {
        if (C == null || C.isZERO()) {
            return false;
        }
        GenPolynomialRing pfac = C.ring;
        ModIntegerRing pm = new ModIntegerRing(p.getVal(), true);
        GenPolynomialRing mfac = new GenPolynomialRing(pm, pfac);

        // check mod p
        GenPolynomial cl = mfac.getONE();
        GenPolynomial hlp;
        for (GenPolynomial hl : G) {
            //System.out.println("hl       = " + hl);
            hlp = PolyUtil. fromIntegerCoefficients(mfac, hl);
            //System.out.println("hl mod p = " + hlp);
            cl = cl.multiply(hlp);
        }
        GenPolynomial cp = PolyUtil. fromIntegerCoefficients(mfac, C);
        if (!cp.equals(cl)) {
            System.out.println("Hensel precondition wrong!");
            if (debug) {
                System.out.println("cl      = " + cl);
                System.out.println("cp      = " + cp);
                System.out.println("mon(cl) = " + cl.monic());
                System.out.println("mon(cp) = " + cp.monic());
                System.out.println("cp-cl   = " + cp.subtract(cl));
                System.out.println("M       = " + M + ", p = " + p);
            }
            return false;
        }

        // check mod p^e 
        BigInteger mip = p;
        while (mip.compareTo(M) < 0) {
            mip = mip.multiply(mip); // p
        }
        // mip = mip.multiply(p);
        pm = new ModIntegerRing(mip.getVal(), false);
        mfac = new GenPolynomialRing(pm, pfac);
        cl = mfac.getONE();
        for (GenPolynomial hl : G) {
            //System.out.println("hl         = " + hl);
            hlp = PolyUtil. fromIntegerCoefficients(mfac, hl);
            //System.out.println("hl mod p^e = " + hlp);
            cl = cl.multiply(hlp);
        }
        cp = PolyUtil. fromIntegerCoefficients(mfac, C);
        if (!cp.equals(cl)) {
            System.out.println("Hensel post condition wrong!");
            System.out.println("cl    = " + cl);
            System.out.println("cp    = " + cp);
            System.out.println("cp-cl = " + cp.subtract(cl));
            System.out.println("M = " + M + ", p = " + p + ", p^e = " + mip);
            return false;
        }
        return true;
    }


    /**
     * Modular Hensel lifting test. Let p be a prime number and assume C == A *
     * B mod p with gcd(A,B) == 1 mod p.
     * @param C GenPolynomial
     * @param A GenPolynomial
     * @param B GenPolynomial
     * @param M bound on the coefficients of A and B as factors of C.
     * @param p prime number.
     * @return true if C = A * B mod p**e, else false.
     */
    public static// & Modular>
    boolean isHenselLift(GenPolynomial C, BigInteger M, BigInteger p,
            GenPolynomial A, GenPolynomial B) {
        List> G = new ArrayList>(2);
        G.add(A);
        G.add(B);
        return isHenselLift(C, M, p, G);
    }


    /**
     * Modular Hensel lifting test. Let p be a prime number and assume C == A *
     * B mod p with gcd(A,B) == 1 mod p.
     * @param C GenPolynomial
     * @param Ha Hensel approximation.
     * @param M bound on the coefficients of A and B as factors of C.
     * @param p prime number.
     * @return true if C = A * B mod p^e, else false.
     */
    public static  & Modular> boolean isHenselLift(GenPolynomial C,
            BigInteger M, BigInteger p, HenselApprox Ha) {
        List> G = new ArrayList>(2);
        G.add(Ha.A);
        G.add(Ha.B);
        return isHenselLift(C, M, p, G);
    }


    /**
     * Constructing and lifting algorithm for extended Euclidean relation. Let p
     * = A.ring.coFac.modul() and assume gcd(A,B) == 1 mod p.
     * @param A modular GenPolynomial
     * @param B modular GenPolynomial
     * @param k desired approximation exponent p^k.
     * @return [s,t] with s A + t B = 1 mod p^k.
     */
    @SuppressWarnings("unchecked")
    public static  & Modular> GenPolynomial[] liftExtendedEuclidean(
            GenPolynomial A, GenPolynomial B, long k) throws NoLiftingException {
        if (A == null || A.isZERO() || B == null || B.isZERO()) {
            throw new IllegalArgumentException("A and B must be nonzero, A = " + A + ", B = " + B);
        }
        GenPolynomialRing fac = A.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        // start with extended Euclidean relation mod p
        GenPolynomial[] gst = null;
        try {
            gst = A.egcd(B);
            if (!gst[0].isONE()) {
                throw new NoLiftingException("A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B);
            }
        } catch (ArithmeticException e) {
            throw new NoLiftingException("coefficient error " + e);
        }
        GenPolynomial S = gst[1];
        GenPolynomial T = gst[2];
        //System.out.println("eeS = " + S + ": " + S.ring.coFac);
        //System.out.println("eeT = " + T + ": " + T.ring.coFac);

        // setup integer polynomial ring
        GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(), fac);
        GenPolynomial one = ifac.getONE();
        GenPolynomial Ai = PolyUtil.integerFromModularCoefficients(ifac, A);
        GenPolynomial Bi = PolyUtil.integerFromModularCoefficients(ifac, B);
        GenPolynomial Si = PolyUtil.integerFromModularCoefficients(ifac, S);
        GenPolynomial Ti = PolyUtil.integerFromModularCoefficients(ifac, T);
        //System.out.println("Ai = " + Ai);
        //System.out.println("Bi = " + Bi);
        //System.out.println("Si = " + Si);
        //System.out.println("Ti = " + Ti);

        // approximate mod p^i
        ModularRingFactory mcfac = (ModularRingFactory) fac.coFac;
        BigInteger p = mcfac.getIntegerModul();
        BigInteger modul = p;
        GenPolynomialRing mfac; // = new GenPolynomialRing(mcfac, fac);
        for (int i = 1; i < k; i++) {
            // e = 1 - s a - t b in Z[x]
            GenPolynomial e = one.subtract(Si.multiply(Ai)).subtract(Ti.multiply(Bi));
            //System.out.println("\ne = " + e);
            if (e.isZERO()) {
                logger.info("leaving on zero e in liftExtendedEuclidean");
                break;
            }
            e = e.divide(modul);
            // move to Z_p[x] and compute next approximation 
            GenPolynomial c = PolyUtil. fromIntegerCoefficients(fac, e);
            //System.out.println("c = " + c + ": " + c.ring.coFac);
            GenPolynomial s = S.multiply(c);
            GenPolynomial t = T.multiply(c);
            //System.out.println("s = " + s + ": " + s.ring.coFac);
            //System.out.println("t = " + t + ": " + t.ring.coFac);

            GenPolynomial[] QR = s.quotientRemainder(B); // watch for ordering 
            GenPolynomial q = QR[0];
            s = QR[1];
            t = t.sum(q.multiply(A));
            //System.out.println("s = " + s + ": " + s.ring.coFac);
            //System.out.println("t = " + t + ": " + t.ring.coFac);

            GenPolynomial si = PolyUtil.integerFromModularCoefficients(ifac, s);
            GenPolynomial ti = PolyUtil.integerFromModularCoefficients(ifac, t);
            //System.out.println("si = " + si);
            //System.out.println("ti = " + si);
            // add approximation to solution
            Si = Si.sum(si.multiply(modul));
            Ti = Ti.sum(ti.multiply(modul));
            //System.out.println("Si = " + Si);
            //System.out.println("Ti = " + Ti);
            modul = modul.multiply(p);
            //System.out.println("modul = " + modul + ", " + p + "^" + k + ", p^k = " + Power.power(new BigInteger(),p,k));
        }
        //System.out.println("Si = " + Si + ", Ti = " + Ti);
        // setup ring mod p^i
        if (ModLongRing.MAX_LONG.compareTo(modul.getVal()) > 0) {
            mcfac = (ModularRingFactory) new ModLongRing(modul.getVal());
        } else {
            mcfac = (ModularRingFactory) new ModIntegerRing(modul.getVal());
        }
        //System.out.println("mcfac = " + mcfac);
        mfac = new GenPolynomialRing(mcfac, fac);
        S = PolyUtil. fromIntegerCoefficients(mfac, Si);
        T = PolyUtil. fromIntegerCoefficients(mfac, Ti);
        //System.out.println("S = " + S + ": " + S.ring.coFac);
        //System.out.println("T = " + T + ": " + T.ring.coFac);
        if ( true || debug ) {
            List> AP = new ArrayList>();
            AP.add(B); AP.add(A);
            List> SP = new ArrayList>();
            SP.add(S); SP.add(T);
            if ( !HenselUtil.isExtendedEuclideanLift(AP,SP) ) {
                System.out.println("isExtendedEuclideanLift: false");
            }
        }
        GenPolynomial[] rel = (GenPolynomial[]) new GenPolynomial[2];
        rel[0] = S;
        rel[1] = T;
        return rel;
    }


    /**
     * Constructing and lifting algorithm for extended Euclidean relation. Let p
     * = A_i.ring.coFac.modul() and assume gcd(A_i,A_j) == 1 mod p, i != j.
     * @param A list of modular GenPolynomials
     * @param k desired approximation exponent p^k.
     * @return [s_0,...,s_n-1] with sum_i s_i * B_i = 1 mod p^k, with B_i =
     *         prod_{i!=j} A_j.
     */
    public static  & Modular> List> liftExtendedEuclidean(
            List> A, long k) throws NoLiftingException {
        if (A == null || A.size() == 0) {
            throw new IllegalArgumentException("A must be non null and non empty");
        }
        GenPolynomialRing fac = A.get(0).ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        GenPolynomial zero = fac.getZERO();
        int r = A.size();
        List> Q = new ArrayList>(r);
        for (int i = 0; i < r; i++) {
            Q.add(zero);
        }
        //System.out.println("A = " + A);
        Q.set(r - 2, A.get(r - 1));
        for (int j = r - 3; j >= 0; j--) {
            GenPolynomial q = A.get(j + 1).multiply(Q.get(j + 1));
            Q.set(j, q);
        }
        //System.out.println("Q = " + Q);
        List> B = new ArrayList>(r + 1);
        List> lift = new ArrayList>(r);
        for (int i = 0; i < r; i++) {
            B.add(zero);
            lift.add(zero);
        }
        GenPolynomial one = fac.getONE();
        GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(), fac);
        B.add(0, one);
        //System.out.println("B(0) = " + B.get(0));
        GenPolynomial b = one;
        for (int j = 0; j < r - 1; j++) {
            //System.out.println("Q("+(j)+") = " + Q.get(j));
            //System.out.println("A("+(j)+") = " + A.get(j));
            //System.out.println("B("+(j)+") = " + B.get(j));
            List> S = liftDiophant(Q.get(j), A.get(j), B.get(j), k);
            //System.out.println("\nSb = " + S);
            b = S.get(0);
            GenPolynomial bb = PolyUtil. fromIntegerCoefficients(fac, PolyUtil
                    .integerFromModularCoefficients(ifac, b));
            B.set(j + 1, bb);
            lift.set(j, S.get(1));
            //System.out.println("B("+(j+1)+") = " + B.get(j+1));
            if (debug) {
                logger.info("lift(" + j + ") = " + lift.get(j));
            }
        }
        //System.out.println("liftb = " + lift);
        lift.set(r - 1, b);
        if (debug) {
            logger.info("lift(" + (r - 1) + ") = " + b);
        }
        //System.out.println("B("+(r-1)+") = " + B.get(r-1) + " : " +  B.get(r-1).ring.coFac + ", b = " +  b + " : " +  b.ring.coFac);
        //System.out.println("B = " + B);
        //System.out.println("liftb = " + lift);
        return lift;
    }


    /**
     * Modular diophantine equation solution and lifting algorithm. Let p =
     * A_i.ring.coFac.modul() and assume gcd(A,B) == 1 mod p.
     * @param A modular GenPolynomial, mod p^k
     * @param B modular GenPolynomial, mod p^k
     * @param C modular GenPolynomial, mod p^k
     * @param k desired approximation exponent p^k.
     * @return [s, t] with s A' + t B' = C mod p^k, with A' = B, B' = A.
     */
    public static  & Modular> List> liftDiophant(
            GenPolynomial A, GenPolynomial B, GenPolynomial C, long k)
            throws NoLiftingException {
        if (A == null || A.isZERO() || B == null || B.isZERO()) {
            throw new IllegalArgumentException("A and B must be nonzero, A = " + A + ", B = " + B  + ", C = " + C);
        }
        List> sol = new ArrayList>();
        GenPolynomialRing fac = C.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        //System.out.println("C = " + C);
        GenPolynomial zero = fac.getZERO();
        for (int i = 0; i < 2; i++) {
            sol.add(zero);
        }
        GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(),fac);
        for (Monomial m : C) {
            //System.out.println("monomial = " + m);
            long e = m.e.getVal(0);
            List> S = liftDiophant(A, B, e, k);
            //System.out.println("Se = " + S);
            MOD a = m.c;
            //System.out.println("C.fac = " + fac.toScript());
            a = fac.coFac.fromInteger(a.getSymmetricInteger().getVal());
            int i = 0;
            for (GenPolynomial d : S) {
                //System.out.println("d = " + d);
                d = PolyUtil. fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, d));
                d = d.multiply(a);
                d = sol.get(i).sum(d);
                //System.out.println("d = " + d);
                sol.set(i++, d);
            }
            //System.out.println("sol = " + sol + ", for " + m);
        }
        if (true || debug) {
            //GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(), fac);
            A = PolyUtil. fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, A));
            B = PolyUtil. fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, B));
            C = PolyUtil. fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, C));
            GenPolynomial y = B.multiply(sol.get(0)).sum(A.multiply(sol.get(1)));
            if (!y.equals(C)) {
                System.out.println("A = " + A + ", B = " + B);
                System.out.println("s1 = " + sol.get(0) + ", s2 = " + sol.get(1));
                System.out.println("Error: A*r1 + B*r2 = " + y + " : " + fac.coFac);
            }
        }
        return sol;
    }


    /**
     * Modular diophantine equation solution and lifting algorithm. Let p =
     * A_i.ring.coFac.modul() and assume gcd(a,b) == 1 mod p, for a, b in A.
     * @param A list of modular GenPolynomials, mod p^k
     * @param C modular GenPolynomial, mod p^k
     * @param k desired approximation exponent p^k.
     * @return [s_1,..., s_n] with sum_i s_i A_i' = C mod p^k, with Ai' = prod_{j!=i} A_j.
     */
    public static  & Modular> List> liftDiophant(
            List> A, GenPolynomial C, long k)
            throws NoLiftingException {
        if ( false && A.size() <= 2 ) {
            return HenselUtil. liftDiophant(A.get(0),A.get(1),C,k);
        }
        List> sol = new ArrayList>();
        GenPolynomialRing fac = C.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        //System.out.println("C = " + C);
        GenPolynomial zero = fac.getZERO();
        for (int i = 0; i < A.size(); i++) {
            sol.add(zero);
        }
        GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(),fac);
        for (Monomial m : C) {
            //System.out.println("monomial = " + m);
            long e = m.e.getVal(0);
            List> S = liftDiophant(A, e, k);
            //System.out.println("Se = " + S);
            MOD a = m.c;
            //System.out.println("C.fac = " + fac.toScript());
            a = fac.coFac.fromInteger(a.getSymmetricInteger().getVal());
            int i = 0;
            for (GenPolynomial d : S) {
                //System.out.println("d = " + d);
                d = PolyUtil. fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, d));
                d = d.multiply(a);
                d = sol.get(i).sum(d);
                //System.out.println("d = " + d);
                sol.set(i++, d);
            }
            //System.out.println("sol = " + sol + ", for " + m);
        }
        /*
        if (true || debug) {
            //GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(), fac);
            A = PolyUtil. fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, A));
            B = PolyUtil. fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, B));
            C = PolyUtil. fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, C));
            GenPolynomial y = B.multiply(sol.get(0)).sum(A.multiply(sol.get(1)));
            if (!y.equals(C)) {
                System.out.println("A = " + A + ", B = " + B);
                System.out.println("s1 = " + sol.get(0) + ", s2 = " + sol.get(1));
                System.out.println("Error: A*r1 + B*r2 = " + y + " : " + fac.coFac);
            }
        }
        */
        return sol;
    }


    /**
     * Modular diophantine equation solution and lifting algorithm. Let p =
     * A_i.ring.coFac.modul() and assume gcd(A,B) == 1 mod p.
     * @param A modular GenPolynomial
     * @param B modular GenPolynomial
     * @param e exponent for x^e
     * @param k desired approximation exponent p^k.
     * @return [s, t] with s A' + t B' = x^e mod p^k, with A' = B, B' = A.
     */
    public static  & Modular> List> liftDiophant(
            GenPolynomial A, GenPolynomial B, long e, long k) throws NoLiftingException {
        if (A == null || A.isZERO() || B == null || B.isZERO()) {
            throw new IllegalArgumentException("A and B must be nonzero, A = " + A + ", B = " + B);
        }
        List> sol = new ArrayList>();
        GenPolynomialRing fac = A.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        // lift EE relation to p^k
        GenPolynomial[] lee = liftExtendedEuclidean(B,A, k);
        GenPolynomial s1 = lee[0];
        GenPolynomial s2 = lee[1];
        if ( e == 0L ) {
            sol.add(s1);
            sol.add(s2);
            //System.out.println("sol@0 = " + sol); 
            return sol;
        }
        fac = s1.ring;
        GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(), fac);
        A = PolyUtil. fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, A));
        B = PolyUtil. fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, B));

        //      this is the wrong sequence:
        //         GenPolynomial xe = fac.univariate(0,e);
        //         GenPolynomial q = s1.multiply(xe);
        //         GenPolynomial[] QR = q.quotientRemainder(B);
        //         q = QR[0];
        //         GenPolynomial r1 = QR[1];
        //         GenPolynomial r2 = s2.multiply(xe).sum( q.multiply(A) );

        GenPolynomial xe = fac.univariate(0, e);
        GenPolynomial q = s1.multiply(xe);
        GenPolynomial[] QR = q.quotientRemainder(A);
        q = QR[0];
        //System.out.println("ee coeff qr = " + Arrays.toString(QR));
        GenPolynomial r1 = QR[1];
        GenPolynomial r2 = s2.multiply(xe).sum(q.multiply(B));
        //System.out.println("r1 = " + r1 + ", r2 = " + r2);
        sol.add(r1);
        sol.add(r2);
        //System.out.println("sol@"+ e + " = " + sol); 
        if (true || debug) {
            GenPolynomial y = B.multiply(r1).sum(A.multiply(r2));
            if (!y.equals(xe)) {
                System.out.println("A = " + A + ", B = " + B);
                System.out.println("r1 = " + r1 + ", r2 = " + r2);
                System.out.println("Error: A*r1 + B*r2 = " + y);
            }
        }
        return sol;
    }


    /**
     * Modular diophantine equation solution and lifting algorithm. Let p =
     * A_i.ring.coFac.modul() and assume gcd(a,b) == 1 mod p, for a, b in A.
     * @param A list of modular GenPolynomials
     * @param e exponent for x^e
     * @param k desired approximation exponent p^k.
     * @return [s_1,..., s_n] with sum_i s_i A_i' = x^e mod p^k, with Ai' = prod_{j!=i} A_j.
     */
    public static  & Modular> List> liftDiophant(
            List> A, long e, long k) throws NoLiftingException {
        if ( false && A.size() <= 2 ) {
            return HenselUtil. liftDiophant(A.get(0),A.get(1),e,k);
        }
        List> sol = new ArrayList>();
        GenPolynomialRing fac = A.get(0).ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        // lift EE relation to p^k
        List> lee = liftExtendedEuclidean(A, k);
        if ( e == 0L ) {
            //System.out.println("sol@0 = " + sol); 
            return lee;
        }
        fac = lee.get(0).ring;
        GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(), fac);
        List> S = new ArrayList>(lee.size());
        for ( GenPolynomial a : lee ) {
             a = PolyUtil. fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, a));
             S.add(a);
        }
        GenPolynomial xe = fac.univariate(0, e);
        //List> Sr = new ArrayList>(lee.size());
        int i = 0;
        for ( GenPolynomial s : S ) {
             GenPolynomial q = s.multiply(xe);
             GenPolynomial r = q.remainder(A.get(i++));
             //System.out.println("r = " + r);
             sol.add(r);
        }
        //System.out.println("sol@"+ e + " = " + sol); 
        /*
        if (true || debug) {
            GenPolynomial y = B.multiply(r1).sum(A.multiply(r2));
            if (!y.equals(xe)) {
                System.out.println("A = " + A + ", B = " + B);
                System.out.println("r1 = " + r1 + ", r2 = " + r2);
                System.out.println("Error: A*r1 + B*r2 = " + y);
            }
        }
        */
        return sol;
    }


    /**
     * Modular Diophant relation lifting test.
     * @param A modular GenPolynomial
     * @param B modular GenPolynomial
     * @param C modular GenPolynomial
     * @param S1 modular GenPolynomial
     * @param S2 modular GenPolynomial
     * @return true if A*S1 + B*S2 = C, else false.
     */
    public static  & Modular> boolean isDiophantLift(
            GenPolynomial A, GenPolynomial B, GenPolynomial S1, GenPolynomial S2, GenPolynomial C) {
        GenPolynomialRing fac = C.ring;
        GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(),fac);
        GenPolynomial a  = PolyUtil. fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,A));
        GenPolynomial b  = PolyUtil. fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,B));
        GenPolynomial s1 = PolyUtil. fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,S1));
        GenPolynomial s2 = PolyUtil. fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,S2));
        GenPolynomial t  = a.multiply(s1).sum( b.multiply(s2) );
        if ( t.equals(C) ) {
            return true;
        }
        if (debug) {
            System.out.println("a  = " + a);
            System.out.println("b  = " + b);
            System.out.println("s1 = " + s1);
            System.out.println("s2 = " + s2);
            System.out.println("t  = " + t);
            System.out.println("C  = " + C);
        }
        return false;
    }


    /**
     * Modular extended Euclidean relation lifting test.
     * @param A list of GenPolynomials
     * @param S = [s_0,...,s_{n-1}] list of GenPolynomial
     * @return true if prod_{0,...,n-1} s_i * B_i = 1 mod p^e, with B_i =
     *         prod_{i!=j} A_j, else false.
     */
    public static  & Modular> boolean isExtendedEuclideanLift(
            List> A, List> S) {
        GenPolynomialRing fac = A.get(0).ring;
        GenPolynomial C = fac.getONE();
        return isDiophantLift(A,S,C);
    }


    /**
     * Modular Diophant relation lifting test.
     * @param A list of GenPolynomials
     * @param S = [s_0,...,s_{n-1}] list of GenPolynomials
     * @param C = GenPolynomial
     * @return true if prod_{0,...,n-1} s_i * B_i = C mod p^k, with B_i =
     *         prod_{i!=j} A_j, else false.
     */
    public static  & Modular> boolean isDiophantLift(
                  List> A, List> S, GenPolynomial C) {
        GenPolynomialRing fac = A.get(0).ring;
        GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(),fac);
        List> B = new ArrayList>(A.size());
        int i = 0;
        for (GenPolynomial ai : A) {
            GenPolynomial b = fac.getONE();
            int j = 0;
            for (GenPolynomial aj : A) {
                if (i != j /*!ai.equals(aj)*/) {
                    b = b.multiply(aj);
                }
                j++;
            }
            //System.out.println("b = " + b);
            b = PolyUtil. fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,b));
            B.add(b);
            i++;
        }
        //System.out.println("B = " + B);
        // check mod p^e 
        GenPolynomial t = fac.getZERO();
        i = 0;
        for (GenPolynomial a : B) {
            GenPolynomial b = S.get(i++);
            b = PolyUtil. fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,b));
            GenPolynomial s = a.multiply(b);
            t = t.sum(s);
        }
        if (!t.equals(C)) {
            if (debug) { 
                System.out.println("no diophant lift!");
                System.out.println("A = " + A);
                System.out.println("B = " + B);
                System.out.println("S = " + S);
                System.out.println("C = " + C);
                System.out.println("t = " + t);
            }
            return false;
        }
        return true;
    }


    /**
     * Modular Hensel lifting algorithm on coefficients. Let p =
     * f_i.ring.coFac.modul() and assume C == prod_{0,...,n-1} f_i mod p with
     * gcd(f_i,f_j) == 1 mod p for i != j
     * @param C monic integer polynomial
     * @param F = [f_0,...,f_{n-1}] list of monic modular polynomials.
     * @param k approximation exponent.
     * @return [g_0,...,g_{n-1}] with C = prod_{0,...,n-1} g_i mod p^k.
     */
    @SuppressWarnings("unchecked")
    public static  & Modular> 
        List> liftHenselMonic(GenPolynomial C, List> F, long k) 
                                                 throws NoLiftingException {
        if (C == null || C.isZERO() || F == null || F.size() == 0) {
            throw new IllegalArgumentException("C must be nonzero and F must be nonempty");
        }
        GenPolynomialRing fac = C.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        List> lift = new ArrayList>(F.size());
        GenPolynomialRing pfac = F.get(0).ring;
        RingFactory pcfac = pfac.coFac;
        ModularRingFactory PF = (ModularRingFactory) pcfac;
        BigInteger P = PF.getIntegerModul();
        int n = F.size();
        if (n == 1) { // lift F_0, this case will probably never be used
            GenPolynomial f = F.get(0);
            ModularRingFactory mcfac;
            if (ModLongRing.MAX_LONG.compareTo(P.getVal()) > 0) {
                mcfac = (ModularRingFactory) new ModLongRing(P.getVal());
            } else {
                mcfac = (ModularRingFactory) new ModIntegerRing(P.getVal());
            }
            GenPolynomialRing mfac = new GenPolynomialRing(mcfac, fac);
            f = PolyUtil.fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac, f));
            lift.add(f);
            return lift;
        }
        //         if (n == 2) { // only one step
        //             HenselApprox ab = HenselUtil. liftHenselQuadratic(C, M, F.get(0), F.get(1));
        //             lift.add(ab.Am);
        //             lift.add(ab.Bm);
        //             return lift;
        //         }

        // setup integer polynomial ring
        GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(), fac);
        List> Fi = PolyUtil.integerFromModularCoefficients(ifac, F);
        //System.out.println("Fi = " + Fi);

        List> S = liftExtendedEuclidean(F, k + 1); // lift works for any k, TODO: use this
        //System.out.println("Sext = " + S);
        if (true || debug) {
            logger.info("EE lift = " + S);
            // adjust coefficients
            List> Sx = PolyUtil.fromIntegerCoefficients(pfac, PolyUtil
                    .integerFromModularCoefficients(ifac, S));
            try {
                boolean il = HenselUtil. isExtendedEuclideanLift(F, Sx);
                //System.out.println("islift = " + il);
            } catch (RuntimeException e) {
                e.printStackTrace();
            }
        }
        List> Si = PolyUtil.integerFromModularCoefficients(ifac, S);
        //System.out.println("Si = " + Si);
        //System.out.println("C = " + C);

        // approximate mod p^i
        ModularRingFactory mcfac = PF;
        BigInteger p = mcfac.getIntegerModul();
        BigInteger modul = p;
        GenPolynomialRing mfac = new GenPolynomialRing(mcfac, fac);
        List> Sp = PolyUtil.fromIntegerCoefficients(mfac, Si);
        //System.out.println("Sp = " + Sp);
        for (int i = 1; i < k; i++) {
            //System.out.println("i = " + i);
            GenPolynomial e = fac.getONE();
            for (GenPolynomial fi : Fi) {
                e = e.multiply(fi);
            }
            e = C.subtract(e);
            //System.out.println("\ne = " + e);
            if (e.isZERO()) {
                logger.info("leaving on zero e");
                break;
            }
            try {
                e = e.divide(modul);
            } catch (RuntimeException ex) {
                ex.printStackTrace();
                throw ex;
            }
            //System.out.println("e = " + e);
            // move to in Z_p[x]
            GenPolynomial c = PolyUtil. fromIntegerCoefficients(mfac, e);
            //System.out.println("c = " + c + ": " + c.ring.coFac);

            List> s = new ArrayList>(S.size());
            int j = 0;
            for (GenPolynomial f : Sp) {
                f = f.multiply(c);
                //System.out.println("f = " + f + " : " + f.ring.coFac);
                //System.out.println("F,i = " + F.get(j) + " : " + F.get(j).ring.coFac);
                f = f.remainder(F.get(j++));
                //System.out.println("f = " + f + " : " + f.ring.coFac);
                s.add(f);
            }
            //System.out.println("s = " + s);
            List> si = PolyUtil.integerFromModularCoefficients(ifac, s);
            //System.out.println("si = " + si);

            List> Fii = new ArrayList>(F.size());
            j = 0;
            for (GenPolynomial f : Fi) {
                f = f.sum(si.get(j++).multiply(modul));
                Fii.add(f);
            }
            //System.out.println("Fii = " + Fii);
            Fi = Fii;
            modul = modul.multiply(p);
            if (i >= k - 1) {
                logger.info("e != 0 for k = " + k);
            }
        }
        // setup ring mod p^k
        modul = Power.positivePower(p,k);
        if (ModLongRing.MAX_LONG.compareTo(modul.getVal()) > 0) {
            mcfac = (ModularRingFactory) new ModLongRing(modul.getVal());
        } else {
            mcfac = (ModularRingFactory) new ModIntegerRing(modul.getVal());
        }
        //System.out.println("mcfac = " + mcfac);
        mfac = new GenPolynomialRing(mcfac, fac);
        lift = PolyUtil. fromIntegerCoefficients(mfac, Fi);
        //System.out.println("lift = " + lift + ": " + lift.get(0).ring.coFac);
        return lift;
    }


    /**
     * Modular Hensel lifting algorithm on coefficients. Let p =
     * f_i.ring.coFac.modul() and assume C == prod_{0,...,n-1} f_i mod p with
     * gcd(f_i,f_j) == 1 mod p for i != j
     * @param C integer polynomial
     * @param F = [f_0,...,f_{n-1}] list of monic modular polynomials.
     * @param k approximation exponent.
     * @param g leading coefficient.
     * @return [g_0,...,g_{n-1}] with C = prod_{0,...,n-1} g_i mod p^k.
     */
    @SuppressWarnings("unchecked")
    public static  & Modular> 
        List> liftHensel(GenPolynomial C, List> F, long k, BigInteger g) 
                                            throws NoLiftingException {
        if (C == null || C.isZERO() || F == null || F.size() == 0) {
            throw new IllegalArgumentException("C must be nonzero and F must be nonempty");
        }
        GenPolynomialRing fac = C.ring;
        if (fac.nvar != 1) { // todo assert
            throw new IllegalArgumentException("polynomial ring not univariate");
        }
        List> lift = new ArrayList>(F.size());
        GenPolynomialRing pfac = F.get(0).ring;
        RingFactory pcfac = pfac.coFac;
        ModularRingFactory PF = (ModularRingFactory) pcfac;
        BigInteger P = PF.getIntegerModul();
        int n = F.size();
        if (n == 1) { // lift F_0, this case will probably never be used
            GenPolynomial f = F.get(0);
            ModularRingFactory mcfac;
            if (ModLongRing.MAX_LONG.compareTo(P.getVal()) > 0) {
                mcfac = (ModularRingFactory) new ModLongRing(P.getVal());
            } else {
                mcfac = (ModularRingFactory) new ModIntegerRing(P.getVal());
            }
            GenPolynomialRing mfac = new GenPolynomialRing(mcfac, fac);
            f = PolyUtil.fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac, f));
            lift.add(f);
            return lift;
        }
        // if (n == 2) { // only one step
        //     HenselApprox ab = HenselUtil. liftHenselQuadratic(C, M, F.get(0), F.get(1));
        //     lift.add(ab.Am);
        //     lift.add(ab.Bm);
        //     return lift;
        // }
        
        // normalize C and F_i factors
        BigInteger cc = g;     //C.leadingBaseCoefficient(); // == g ??
        for ( int i = 1; i < F.size(); i++ ) { // #F-1
             C = C.multiply(cc); // sic
        }
        MOD cm = PF.fromInteger(cc.getVal());
        List> Fp = new ArrayList>(F.size());
        for ( GenPolynomial fm : F ) {
            GenPolynomial am = fm.monic();
            am = am.multiply(cm);
            Fp.add(am);
        }
        F = Fp;

        // setup integer polynomial ring
        GenPolynomialRing ifac = new GenPolynomialRing(new BigInteger(), fac);
        List> Fi = PolyUtil.integerFromModularCoefficients(ifac, F);
        //System.out.println("Fi = " + Fi);

        // inplace modify polynomials, replace leading coefficient
        for ( GenPolynomial ai : Fi ) {
            if ( ai.isZERO() ) {
                continue;
            }
            ExpVector ea = ai.leadingExpVector();
            ai.doPutToMap(ea, cc);
        }
        //System.out.println("Fi = " + Fi);

        List> S = liftExtendedEuclidean(F, k + 1); // lift works for any k, TODO: use this
        //System.out.println("Sext = " + S);
        if (true || debug) {
            logger.info("EE lift = " + S);
            // adjust coefficients
            List> Sx = PolyUtil.fromIntegerCoefficients(pfac, PolyUtil
                    .integerFromModularCoefficients(ifac, S));
            try {
                boolean il = HenselUtil. isExtendedEuclideanLift(F, Sx);
                //System.out.println("islift = " + il);
            } catch (RuntimeException e) {
                e.printStackTrace();
            }
        }
        List> Si = PolyUtil.integerFromModularCoefficients(ifac, S);
        //System.out.println("Si = " + Si);
        //System.out.println("C = " + C);

        // approximate mod p^i
        ModularRingFactory mcfac = PF;
        BigInteger p = mcfac.getIntegerModul();
        BigInteger modul = p;
        GenPolynomialRing mfac = new GenPolynomialRing(mcfac, fac);
        List> Sp = PolyUtil.fromIntegerCoefficients(mfac, Si);
        //System.out.println("Sp = " + Sp);
        for (int i = 1; i < k; i++) {
            //System.out.println("i = " + i);
            GenPolynomial e = fac.getONE();
            for (GenPolynomial fi : Fi) {
                e = e.multiply(fi);
            }
            e = C.subtract(e);
            //System.out.println("\ne = " + e);
            if (e.isZERO()) {
                logger.info("leaving on zero e");
                break;
            }
            try {
                e = e.divide(modul);
            } catch (RuntimeException ex) {
                ex.printStackTrace();
                throw ex;
            }
            //System.out.println("e = " + e);
            // move to in Z_p[x]
            GenPolynomial c = PolyUtil. fromIntegerCoefficients(mfac, e);
            //System.out.println("c = " + c + ": " + c.ring.coFac);

            List> s = new ArrayList>(S.size());
            int j = 0;
            for (GenPolynomial f : Sp) {
                f = f.multiply(c);
                //System.out.println("f = " + f + " : " + f.ring.coFac);
                //System.out.println("F,i = " + F.get(j) + " : " + F.get(j).ring.coFac);
                f = f.remainder(F.get(j++));
                //System.out.println("f = " + f + " : " + f.ring.coFac);
                s.add(f);
            }
            //System.out.println("s = " + s);
            List> si = PolyUtil.integerFromModularCoefficients(ifac, s);
            //System.out.println("si = " + si);

            List> Fii = new ArrayList>(F.size());
            j = 0;
            for (GenPolynomial f : Fi) {
                f = f.sum(si.get(j++).multiply(modul));
                Fii.add(f);
            }
            //System.out.println("Fii = " + Fii);
            Fi = Fii;
            modul = modul.multiply(p);
            if (i >= k - 1) {
                logger.info("e != 0 for k = " + k);
            }
        }
        //System.out.println("Fi = " + Fi);
        // remove normalization
        GreatestCommonDivisorAbstract ufd = GCDFactory.getImplementation(cc);
        //BigInteger ai = ufd.baseContent(Fi.get(0));
        //System.out.println("ai = " + ai + ", cc = " + cc);
        List> Fii = new ArrayList>(F.size());
        //int j = 0;
        for ( GenPolynomial bi : Fi ) {
            GenPolynomial ci = null;
            //if ( j++ == 0 ) {
            //    ci = bi.divide(ai);
            //} else {
            //    BigInteger i = cc.divide(ai);
            //    ci = bi.divide(i);
            //}
            ci = ufd.basePrimitivePart(bi); // ??
            //System.out.println("bi = " + bi + ", ci = " + ci);
            Fii.add(ci);
        }
        Fi = Fii;

        // setup ring mod p^k
        modul = Power.positivePower(p,k);
        if (ModLongRing.MAX_LONG.compareTo(modul.getVal()) > 0) {
            mcfac = (ModularRingFactory) new ModLongRing(modul.getVal());
        } else {
            mcfac = (ModularRingFactory) new ModIntegerRing(modul.getVal());
        }
        //System.out.println("mcfac = " + mcfac);
        mfac = new GenPolynomialRing(mcfac, fac);
        lift = PolyUtil. fromIntegerCoefficients(mfac, Fi);
        //System.out.println("lift = " + lift + ": " + lift.get(0).ring.coFac);
        return lift;
    }

}