Source code

001/*
002 * $Id$
003 */
004
005package edu.jas.ufd;
006
007
008import java.util.ArrayList;
009import java.util.List;
010
011import org.apache.logging.log4j.Logger;
012import org.apache.logging.log4j.LogManager; 
013
014import edu.jas.arith.BigInteger;
015import edu.jas.arith.ModIntegerRing;
016import edu.jas.arith.ModLongRing;
017import edu.jas.arith.Modular;
018import edu.jas.arith.ModularRingFactory;
019import edu.jas.poly.GenPolynomial;
020import edu.jas.poly.GenPolynomialRing;
021import edu.jas.poly.PolyUtil;
022import edu.jas.ps.PolynomialTaylorFunction;
023import edu.jas.ps.TaylorFunction;
024import edu.jas.ps.UnivPowerSeries;
025import edu.jas.ps.UnivPowerSeriesRing;
026import edu.jas.structure.GcdRingElem;
027import edu.jas.structure.RingFactory;
028
029
030/**
031 * Hensel multivariate lifting utilities.
032 * @author Heinz Kredel
033 */
034
035public class HenselMultUtil {
036
037
038    private static final Logger logger = LogManager.getLogger(HenselMultUtil.class);
039
040
041    private static final boolean debug = logger.isInfoEnabled();
042
043
044    /**
045     * Modular diophantine equation solution and lifting algorithm. Let p =
046     * A_i.ring.coFac.modul() and assume ggt(A,B) == 1 mod p.
047     * @param A modular GenPolynomial, mod p^k
048     * @param B modular GenPolynomial, mod p^k
049     * @param C modular GenPolynomial, mod p^k
050     * @param V list of substitution values, mod p^k
051     * @param d desired approximation exponent (x_i-v_i)^d.
052     * @param k desired approximation exponent p^k.
053     * @return [s, t] with s A' + t B' = C mod p^k, with A' = B, B' = A.
054     */
055    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant(
056                    GenPolynomial<MOD> A, GenPolynomial<MOD> B, GenPolynomial<MOD> C, List<MOD> V, long d,
057                    long k) throws NoLiftingException {
058        GenPolynomialRing<MOD> pkfac = C.ring;
059        if (pkfac.nvar == 1) { // V, d ignored
060            return HenselUtil.<MOD> liftDiophant(A, B, C, k);
061        }
062        if (!pkfac.equals(A.ring)) {
063            throw new IllegalArgumentException("A.ring != pkfac: " + A.ring + " != " + pkfac);
064        }
065
066        // evaluate at v_n:
067        List<MOD> Vp = new ArrayList<MOD>(V);
068        MOD v = Vp.remove(Vp.size() - 1);
069        //GenPolynomial<MOD> zero = pkfac.getZERO();
070        // (x_n - v)
071        GenPolynomial<MOD> mon = pkfac.getONE();
072        GenPolynomial<MOD> xv = pkfac.univariate(0, 1);
073        xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal()));
074        //System.out.println("xv = " + xv);
075        // A(v), B(v), C(v)
076        ModularRingFactory<MOD> cf = (ModularRingFactory<MOD>) pkfac.coFac;
077        MOD vp = cf.fromInteger(v.getSymmetricInteger().getVal());
078        //System.out.println("v = " + v + ", vp = " + vp);
079        GenPolynomialRing<MOD> ckfac = pkfac.contract(1);
080        GenPolynomial<MOD> Ap = PolyUtil.<MOD> evaluateMain(ckfac, A, vp);
081        GenPolynomial<MOD> Bp = PolyUtil.<MOD> evaluateMain(ckfac, B, vp);
082        GenPolynomial<MOD> Cp = PolyUtil.<MOD> evaluateMain(ckfac, C, vp);
083        //System.out.println("Ap = " + Ap);
084        //System.out.println("Bp = " + Bp);
085        //System.out.println("Cp = " + Cp);
086
087        // recursion:
088        List<GenPolynomial<MOD>> su = HenselMultUtil.<MOD> liftDiophant(Ap, Bp, Cp, Vp, d, k);
089        //System.out.println("su@p^" + k + " = " + su);
090        //System.out.println("coFac = " + su.get(0).ring.coFac.toScript());
091        if (pkfac.nvar == 2 && !HenselUtil.<MOD> isDiophantLift(Bp, Ap, su.get(0), su.get(1), Cp)) {
092            //System.out.println("isDiophantLift: false");
093            throw new NoLiftingException("isDiophantLift: false");
094        }
095        if (!ckfac.equals(su.get(0).ring)) {
096            throw new IllegalArgumentException("qfac != ckfac: " + su.get(0).ring + " != " + ckfac);
097        }
098        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac);
099        //System.out.println("ifac = " + ifac.toScript());
100        String[] mn = new String[] { pkfac.getVars()[pkfac.nvar - 1] };
101        GenPolynomialRing<GenPolynomial<MOD>> qrfac = new GenPolynomialRing<GenPolynomial<MOD>>(ckfac, 1, mn);
102        //System.out.println("qrfac = " + qrfac);
103
104        List<GenPolynomial<MOD>> sup = new ArrayList<GenPolynomial<MOD>>(su.size());
105        List<GenPolynomial<BigInteger>> supi = new ArrayList<GenPolynomial<BigInteger>>(su.size());
106        for (GenPolynomial<MOD> s : su) {
107            GenPolynomial<MOD> sp = s.extend(pkfac, 0, 0L);
108            sup.add(sp);
109            GenPolynomial<BigInteger> spi = PolyUtil.integerFromModularCoefficients(ifac, sp);
110            supi.add(spi);
111        }
112        //System.out.println("sup  = " + sup);
113        //System.out.println("supi = " + supi);
114        GenPolynomial<BigInteger> Ai = PolyUtil.integerFromModularCoefficients(ifac, A);
115        GenPolynomial<BigInteger> Bi = PolyUtil.integerFromModularCoefficients(ifac, B);
116        GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, C);
117        //System.out.println("Ai = " + Ai);
118        //System.out.println("Bi = " + Bi);
119        //System.out.println("Ci = " + Ci);
120        //GenPolynomial<MOD> aq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, Ai);
121        //GenPolynomial<MOD> bq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, Bi);
122        //System.out.println("aq = " + aq);
123        //System.out.println("bq = " + bq);
124
125        // compute error:
126        GenPolynomial<BigInteger> E = Ci; // - sum_i s_i b_i
127        E = E.subtract(Bi.multiply(supi.get(0)));
128        E = E.subtract(Ai.multiply(supi.get(1)));
129        //System.out.println("E     = " + E);
130        if (E.isZERO()) {
131            logger.info("liftDiophant leaving on zero E");
132            return sup;
133        }
134        GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
135        logger.info("Ep(0,{}) = {}", pkfac.nvar, Ep);
136        if (Ep.isZERO()) {
137            logger.info("liftDiophant leaving on zero Ep mod p^k");
138            return sup;
139        }
140        for (int e = 1; e <= d; e++) {
141            //System.out.println("\ne = " + e + " -------------------------------------- " + pkfac.nvar);
142            GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(qrfac, Ep);
143            //System.out.println("Epr   = " + Epr);
144            UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>(
145                            qrfac);
146            //System.out.println("psfac = " + psfac);
147            TaylorFunction<GenPolynomial<MOD>> F = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr);
148            //System.out.println("F     = " + F);
149            //List<GenPolynomial<MOD>> Vs = new ArrayList<GenPolynomial<MOD>>(1);
150            GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal());
151            //Vs.add(vq);
152            //System.out.println("Vs    = " + Vs);
153            UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(F, vq);
154            //System.out.println("Epst  = " + Epst);
155            GenPolynomial<MOD> cm = Epst.coefficient(e);
156            //System.out.println("cm   = " + cm + ", cm.ring   = " + cm.ring.toScript());
157
158            // recursion:
159            List<GenPolynomial<MOD>> S = HenselMultUtil.<MOD> liftDiophant(Ap, Bp, cm, Vp, d, k);
160            //System.out.println("S    = " + S);
161            if (!ckfac.coFac.equals(S.get(0).ring.coFac)) {
162                throw new IllegalArgumentException(
163                                "ckfac != pkfac: " + ckfac.coFac + " != " + S.get(0).ring.coFac);
164            }
165            if (pkfac.nvar == 2 && !HenselUtil.<MOD> isDiophantLift(Ap, Bp, S.get(1), S.get(0), cm)) {
166                //System.out.println("isDiophantLift: false");
167                throw new NoLiftingException("isDiophantLift: false");
168            }
169            mon = mon.multiply(xv); // Power.<GenPolynomial<MOD>> power(pkfac,xv,e);
170            //System.out.println("mon  = " + mon);
171            //List<GenPolynomial<MOD>> Sp = new ArrayList<GenPolynomial<MOD>>(S.size());
172            int i = 0;
173            supi = new ArrayList<GenPolynomial<BigInteger>>(su.size());
174            for (GenPolynomial<MOD> dd : S) {
175                //System.out.println("dd = " + dd);
176                GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L);
177                GenPolynomial<MOD> dm = de.multiply(mon);
178                //Sp.add(dm);
179                de = sup.get(i).sum(dm);
180                //System.out.println("dd = " + dd);
181                sup.set(i++, de);
182                GenPolynomial<BigInteger> spi = PolyUtil.integerFromModularCoefficients(ifac, dm);
183                supi.add(spi);
184            }
185            //System.out.println("Sp   = " + Sp);
186            //System.out.println("sup  = " + sup);
187            //System.out.println("supi = " + supi);
188            // compute new error
189            //E = E; // - sum_i s_i b_i
190            E = E.subtract(Bi.multiply(supi.get(0)));
191            E = E.subtract(Ai.multiply(supi.get(1)));
192            //System.out.println("E     = " + E);
193            if (E.isZERO()) {
194                logger.info("liftDiophant leaving on zero E");
195                return sup;
196            }
197            Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
198            //System.out.println("Ep(" + e + "," + pkfac.nvar + ") = " + Ep); 
199            logger.info("Ep({},{}) = {}", e, pkfac.nvar, Ep);
200            if (Ep.isZERO()) {
201                logger.info("liftDiophant leaving on zero Ep mod p^k");
202                return sup;
203            }
204        }
205        //System.out.println("*** done: " + pkfac.nvar);
206        return sup;
207    }
208
209
210    /**
211     * Modular diophantine equation solution and lifting algorithm. Let p =
212     * A_i.ring.coFac.modul() and assume ggt(a,b) == 1 mod p, for a, b in A.
213     * @param A list of modular GenPolynomials, mod p^k
214     * @param C modular GenPolynomial, mod p^k
215     * @param V list of substitution values, mod p^k
216     * @param d desired approximation exponent (x_i-v_i)^d.
217     * @param k desired approximation exponent p^k.
218     * @return [s_1,..., s_n] with sum_i s_i A_i' = C mod p^k, with Ai' =
219     *         prod_{j!=i} A_j.
220     */
221    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant(
222                    List<GenPolynomial<MOD>> A, GenPolynomial<MOD> C, List<MOD> V, long d, long k)
223                    throws NoLiftingException {
224        GenPolynomialRing<MOD> pkfac = C.ring;
225        if (pkfac.nvar == 1) { // V, d ignored
226            return HenselUtil.<MOD> liftDiophant(A, C, k);
227        }
228        if (!pkfac.equals(A.get(0).ring)) {
229            throw new IllegalArgumentException("A.ring != pkfac: " + A.get(0).ring + " != " + pkfac);
230        }
231        // co-products
232        GenPolynomial<MOD> As = pkfac.getONE();
233        for (GenPolynomial<MOD> a : A) {
234            As = As.multiply(a);
235        }
236        List<GenPolynomial<MOD>> Bp = new ArrayList<GenPolynomial<MOD>>(A.size());
237        for (GenPolynomial<MOD> a : A) {
238            GenPolynomial<MOD> b = PolyUtil.<MOD> basePseudoDivide(As, a);
239            Bp.add(b);
240        }
241
242        // evaluate at v_n:
243        List<MOD> Vp = new ArrayList<MOD>(V);
244        MOD v = Vp.remove(Vp.size() - 1);
245        // (x_n - v)
246        GenPolynomial<MOD> mon = pkfac.getONE();
247        GenPolynomial<MOD> xv = pkfac.univariate(0, 1);
248        xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal()));
249        //System.out.println("xv = " + xv);
250        // A(v), B(v), C(v)
251        ModularRingFactory<MOD> cf = (ModularRingFactory<MOD>) pkfac.coFac;
252        MOD vp = cf.fromInteger(v.getSymmetricInteger().getVal());
253        //System.out.println("v = " + v + ", vp = " + vp);
254        GenPolynomialRing<MOD> ckfac = pkfac.contract(1);
255        List<GenPolynomial<MOD>> Ap = new ArrayList<GenPolynomial<MOD>>(A.size());
256        for (GenPolynomial<MOD> a : A) {
257            GenPolynomial<MOD> ap = PolyUtil.<MOD> evaluateMain(ckfac, a, vp);
258            Ap.add(ap);
259        }
260        GenPolynomial<MOD> Cp = PolyUtil.<MOD> evaluateMain(ckfac, C, vp);
261        //System.out.println("Ap = " + Ap);
262        //System.out.println("Cp = " + Cp);
263
264        // recursion:
265        List<GenPolynomial<MOD>> su = HenselMultUtil.<MOD> liftDiophant(Ap, Cp, Vp, d, k);
266        //System.out.println("su@p^" + k + " = " + su);
267        //System.out.println("coFac = " + su.get(0).ring.coFac.toScript());
268        if (pkfac.nvar == 2 && !HenselUtil.<MOD> isDiophantLift(Ap, su, Cp)) {
269            //System.out.println("isDiophantLift: false");
270            throw new NoLiftingException("isDiophantLift: false");
271        }
272        if (!ckfac.equals(su.get(0).ring)) {
273            throw new IllegalArgumentException("qfac != ckfac: " + su.get(0).ring + " != " + ckfac);
274        }
275        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac);
276        //System.out.println("ifac = " + ifac.toScript());
277        String[] mn = new String[] { pkfac.getVars()[pkfac.nvar - 1] };
278        GenPolynomialRing<GenPolynomial<MOD>> qrfac = new GenPolynomialRing<GenPolynomial<MOD>>(ckfac, 1, mn);
279        //System.out.println("qrfac = " + qrfac);
280
281        List<GenPolynomial<MOD>> sup = new ArrayList<GenPolynomial<MOD>>(su.size());
282        List<GenPolynomial<BigInteger>> supi = new ArrayList<GenPolynomial<BigInteger>>(su.size());
283        for (GenPolynomial<MOD> s : su) {
284            GenPolynomial<MOD> sp = s.extend(pkfac, 0, 0L);
285            sup.add(sp);
286            GenPolynomial<BigInteger> spi = PolyUtil.integerFromModularCoefficients(ifac, sp);
287            supi.add(spi);
288        }
289        //System.out.println("sup  = " + sup);
290        //System.out.println("supi = " + supi);
291        //List<GenPolynomial<BigInteger>> Ai = new ArrayList<GenPolynomial<BigInteger>>(A.size());
292        //for (GenPolynomial<MOD> a : A) {
293        //    GenPolynomial<BigInteger> ai = PolyUtil.integerFromModularCoefficients(ifac, a);
294        //    Ai.add(ai);
295        //}
296        List<GenPolynomial<BigInteger>> Bi = new ArrayList<GenPolynomial<BigInteger>>(A.size());
297        for (GenPolynomial<MOD> b : Bp) {
298            GenPolynomial<BigInteger> bi = PolyUtil.integerFromModularCoefficients(ifac, b);
299            Bi.add(bi);
300        }
301        GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, C);
302        //System.out.println("Ai = " + Ai);
303        //System.out.println("Ci = " + Ci);
304
305        //List<GenPolynomial<MOD>> Aq = new ArrayList<GenPolynomial<MOD>>(A.size());
306        //for (GenPolynomial<BigInteger> ai : Ai) {
307        //    GenPolynomial<MOD> aq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, ai);
308        //    Aq.add(aq);
309        //}
310        //System.out.println("Aq = " + Aq);
311
312        // compute error:
313        GenPolynomial<BigInteger> E = Ci; // - sum_i s_i b_i
314        int i = 0;
315        for (GenPolynomial<BigInteger> bi : Bi) {
316            E = E.subtract(bi.multiply(supi.get(i++)));
317        }
318        //System.out.println("E     = " + E);
319        if (E.isZERO()) {
320            logger.info("liftDiophant leaving on zero E");
321            return sup;
322        }
323        GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
324        logger.info("Ep(0,{}) = {}", pkfac.nvar, Ep);
325        if (Ep.isZERO()) {
326            logger.info("liftDiophant leaving on zero Ep mod p^k");
327            return sup;
328        }
329        for (int e = 1; e <= d; e++) {
330            //System.out.println("\ne = " + e + " -------------------------------------- " + pkfac.nvar);
331            GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(qrfac, Ep);
332            //System.out.println("Epr   = " + Epr);
333            UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>(
334                            qrfac);
335            //System.out.println("psfac = " + psfac);
336            TaylorFunction<GenPolynomial<MOD>> F = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr);
337            //System.out.println("F     = " + F);
338            //List<GenPolynomial<MOD>> Vs = new ArrayList<GenPolynomial<MOD>>(1);
339            GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal());
340            //Vs.add(vq);
341            //System.out.println("Vs    = " + Vs);
342            UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(F, vq);
343            //System.out.println("Epst  = " + Epst);
344            GenPolynomial<MOD> cm = Epst.coefficient(e);
345            //System.out.println("cm   = " + cm + ", cm.ring   = " + cm.ring.toScript());
346            if (cm.isZERO()) {
347                continue;
348            }
349            // recursion:
350            List<GenPolynomial<MOD>> S = HenselMultUtil.<MOD> liftDiophant(Ap, cm, Vp, d, k);
351            //System.out.println("S    = " + S);
352            if (!ckfac.coFac.equals(S.get(0).ring.coFac)) {
353                throw new IllegalArgumentException(
354                                "ckfac != pkfac: " + ckfac.coFac + " != " + S.get(0).ring.coFac);
355            }
356            if (pkfac.nvar == 2 && !HenselUtil.<MOD> isDiophantLift(Ap, S, cm)) {
357                //System.out.println("isDiophantLift: false");
358                throw new NoLiftingException("isDiophantLift: false");
359            }
360            mon = mon.multiply(xv); // Power.<GenPolynomial<MOD>> power(pkfac,xv,e);
361            //System.out.println("mon  = " + mon);
362            //List<GenPolynomial<MOD>> Sp = new ArrayList<GenPolynomial<MOD>>(S.size());
363            i = 0;
364            supi = new ArrayList<GenPolynomial<BigInteger>>(su.size());
365            for (GenPolynomial<MOD> dd : S) {
366                //System.out.println("dd = " + dd);
367                GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L);
368                GenPolynomial<MOD> dm = de.multiply(mon);
369                //Sp.add(dm);
370                de = sup.get(i).sum(dm);
371                //System.out.println("dd = " + dd);
372                sup.set(i++, de);
373                GenPolynomial<BigInteger> spi = PolyUtil.integerFromModularCoefficients(ifac, dm);
374                supi.add(spi);
375            }
376            //System.out.println("Sp   = " + Sp);
377            //System.out.println("sup  = " + sup);
378            //System.out.println("supi = " + supi);
379            // compute new error
380            //E = E; // - sum_i s_i b_i
381            i = 0;
382            for (GenPolynomial<BigInteger> bi : Bi) {
383                E = E.subtract(bi.multiply(supi.get(i++)));
384            }
385            //System.out.println("E     = " + E);
386            if (E.isZERO()) {
387                logger.info("liftDiophant leaving on zero E");
388                return sup;
389            }
390            Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
391            logger.info("Ep({},{}) = {}", e, pkfac.nvar, Ep);
392            if (Ep.isZERO()) {
393                logger.info("liftDiophant leaving on zero Ep mod p^k");
394                return sup;
395            }
396        }
397        //System.out.println("*** done: " + pkfac.nvar);
398        return sup;
399    }
400
401
402    /**
403     * Modular Hensel lifting algorithm on coefficients test. Let p =
404     * f_i.ring.coFac.modul() and assume C == prod_{0,...,n-1} f_i mod p with
405     * gcd(f_i,f_j) == 1 mod p for i != j
406     * @param C integer polynomial
407     * @param Cp GenPolynomial mod p^k
408     * @param F = [f_0,...,f_{n-1}] list of monic modular polynomials.
409     * @param L = [g_0,...,g_{n-1}] list of lifted modular polynomials.
410     * @return true if C = prod_{0,...,n-1} g_i mod p^k, else false.
411     */
412    @SuppressWarnings("unused")
413    public static <MOD extends GcdRingElem<MOD> & Modular> boolean isHenselLift(GenPolynomial<BigInteger> C,
414                    GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F, List<GenPolynomial<MOD>> L) {
415        boolean t = true;
416        GenPolynomialRing<MOD> qfac = L.get(0).ring;
417        GenPolynomial<MOD> q = qfac.getONE();
418        for (GenPolynomial<MOD> fi : L) {
419            q = q.multiply(fi);
420        }
421        t = Cp.equals(q);
422        if (!t) {
423            System.out.println("Cp     = " + Cp);
424            System.out.println("q      = " + q);
425            System.out.println("Cp != q: " + Cp.subtract(q));
426            return t;
427        }
428        GenPolynomialRing<BigInteger> dfac = C.ring;
429        GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(dfac, q);
430        t = C.equals(Ci);
431        if (!t) {
432            System.out.println("C      = " + C);
433            System.out.println("Ci     = " + Ci);
434            System.out.println("C != Ci: " + C.subtract(Ci));
435            return t;
436        }
437        // test L mod id(V) == F
438        return t;
439    }
440
441
442    /**
443     * Modular Hensel lifting algorithm, monic case. Let p =
444     * A_i.ring.coFac.modul() and assume ggt(a,b) == 1 mod p, for a, b in A.
445     * @param C monic GenPolynomial with integer coefficients
446     * @param Cp GenPolynomial mod p^k
447     * @param F list of modular GenPolynomials, mod (I_v, p^k )
448     * @param V list of integer substitution values
449     * @param k desired approximation exponent p^k.
450     * @return [g'_1,..., g'_n] with prod_i g'_i = Cp mod p^k.
451     */
452    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHenselMonic(
453                    GenPolynomial<BigInteger> C, GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F,
454                    List<BigInteger> V, long k) throws NoLiftingException {
455        GenPolynomialRing<MOD> pkfac = Cp.ring;
456        //if (pkfac.nvar == 1) { // V ignored
457        //    return HenselUtil.<MOD> liftHenselMonic(C,F,k);
458        //}
459        long d = C.degree();
460        //System.out.println("d = " + d);
461        // prepare stack of polynomial rings and polynomials
462        List<GenPolynomialRing<MOD>> Pfac = new ArrayList<GenPolynomialRing<MOD>>();
463        List<GenPolynomial<MOD>> Ap = new ArrayList<GenPolynomial<MOD>>();
464        List<MOD> Vb = new ArrayList<MOD>();
465        MOD v = pkfac.coFac.fromInteger(V.get(0).getVal());
466        Pfac.add(pkfac);
467        Ap.add(Cp);
468        Vb.add(v);
469        GenPolynomialRing<MOD> pf = pkfac;
470        GenPolynomial<MOD> ap = Cp;
471        for (int j = pkfac.nvar; j > 2; j--) {
472            pf = pf.contract(1);
473            Pfac.add(0, pf);
474            //MOD vp = pkfac.coFac.fromInteger(V.get(j - 2).getSymmetricInteger().getVal());
475            MOD vp = pkfac.coFac.fromInteger(V.get(j - 2).getVal());
476            //System.out.println("vp     = " + vp);
477            Vb.add(1, vp);
478            ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp);
479            Ap.add(0, ap);
480        }
481        logger.debug("Pfac   = {}", Pfac);
482        //System.out.println("Ap     = " + Ap);
483        //System.out.println("V      = " + V);
484        //System.out.println("Vb     = " + Vb);
485        // setup bi-variate base case
486        GenPolynomialRing<MOD> pk1fac = F.get(0).ring;
487        if (!pkfac.coFac.equals(pk1fac.coFac)) {
488            throw new IllegalArgumentException("F.ring != pkfac: " + pk1fac + " != " + pkfac);
489        }
490        // TODO: adjust leading coefficients
491        pkfac = Pfac.get(0);
492        //Cp = Ap.get(0);
493        //System.out.println("pkfac  = " + pkfac.toScript());
494        //System.out.println("pk1fac = " + pk1fac.toScript());
495        GenPolynomialRing<BigInteger> i1fac = new GenPolynomialRing<BigInteger>(new BigInteger(), pk1fac);
496        //System.out.println("i1fac = " + i1fac.toScript());
497        List<GenPolynomial<BigInteger>> Bi = new ArrayList<GenPolynomial<BigInteger>>(F.size());
498        for (GenPolynomial<MOD> b : F) {
499            GenPolynomial<BigInteger> bi = PolyUtil.integerFromModularCoefficients(i1fac, b);
500            Bi.add(bi);
501        }
502        //System.out.println("Bi = " + Bi);
503        // evaluate Cp at v_n:
504        //ModularRingFactory<MOD> cf = (ModularRingFactory<MOD>) pkfac.coFac;
505        //MOD vp = cf.fromInteger(v.getSymmetricInteger().getVal());
506        //System.out.println("v = " + v + ", vp = " + vp);
507        GenPolynomialRing<MOD> ckfac; // = pkfac.contract(1);
508        //GenPolynomial<MOD> Cs = PolyUtil.<MOD> evaluateMain(ckfac, Cp, vp);
509        //System.out.println("Cp = " + Cp);
510        //System.out.println("Cs = " + Cs);
511
512        List<GenPolynomial<MOD>> U = new ArrayList<GenPolynomial<MOD>>(F.size());
513        for (GenPolynomial<MOD> b : F) {
514            GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L);
515            U.add(bi);
516        }
517        //System.out.println("U  = " + U);
518        List<GenPolynomial<MOD>> U1 = F;
519        //System.out.println("U1 = " + U1);
520
521        GenPolynomial<BigInteger> E = C.ring.getZERO();
522        List<MOD> Vh = new ArrayList<MOD>();
523
524        while (Pfac.size() > 0) { // loop through stack of polynomial rings
525            pkfac = Pfac.remove(0);
526            Cp = Ap.remove(0);
527            v = Vb.remove(0);
528            //Vh.add(0,v);
529            //System.out.println("\npkfac = " + pkfac.toScript() + " ================================== " + Vh);
530
531            // (x_n - v)
532            GenPolynomial<MOD> mon = pkfac.getONE();
533            GenPolynomial<MOD> xv = pkfac.univariate(0, 1);
534            xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal()));
535            //System.out.println("xv = " + xv);
536
537            long deg = Cp.degree(pkfac.nvar - 1);
538            //System.out.println("deg = " + deg);
539
540            GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac);
541            //System.out.println("ifac = " + ifac.toScript());
542            List<GenPolynomial<BigInteger>> Bip = new ArrayList<GenPolynomial<BigInteger>>(F.size());
543            for (GenPolynomial<BigInteger> b : Bi) {
544                GenPolynomial<BigInteger> bi = b.extend(ifac, 0, 0L);
545                Bip.add(bi);
546            }
547            Bi = Bip;
548            //System.out.println("Bi = " + Bi);
549            GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cp);
550            //System.out.println("Ci = " + Ci);
551
552            // compute error:
553            E = ifac.getONE();
554            for (GenPolynomial<BigInteger> bi : Bi) {
555                E = E.multiply(bi);
556            }
557            E = Ci.subtract(E);
558            //System.out.println("E     = " + E);
559            GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
560            logger.info("Ep(0,{},{}) = {}", deg, pkfac.nvar, Ep);
561
562            String[] mn = new String[] { pkfac.getVars()[pkfac.nvar - 1] };
563            ckfac = pkfac.contract(1);
564            GenPolynomialRing<GenPolynomial<MOD>> pkrfac = new GenPolynomialRing<GenPolynomial<MOD>>(ckfac, 1,
565                            mn);
566            //System.out.println("pkrfac = " + pkrfac.toScript());
567
568            for (int e = 1; e <= deg && !Ep.isZERO(); e++) {
569                //System.out.println("\ne = " + e + " -------------------------------------- " + pkfac.nvar);
570                GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(pkrfac, Ep);
571                //System.out.println("Epr   = " + Epr);
572                UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>(
573                                pkrfac);
574                //System.out.println("psfac = " + psfac);
575                TaylorFunction<GenPolynomial<MOD>> T = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr);
576                //System.out.println("T     = " + T);
577                //List<GenPolynomial<MOD>> Vs = new ArrayList<GenPolynomial<MOD>>(1);
578                GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal());
579                //Vs.add(vq);
580                //System.out.println("Vs    = " + Vs + ", Vh = " + Vh);
581                UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(T, vq);
582                //System.out.println("Epst  = " + Epst);
583                logger.info("Epst({},{},{}) = {}", e, deg, pkfac.nvar, Epst);
584                GenPolynomial<MOD> cm = Epst.coefficient(e);
585                //System.out.println("cm   = " + cm);
586                if (cm.isZERO()) {
587                    continue;
588                }
589                List<GenPolynomial<MOD>> Ud = HenselMultUtil.<MOD> liftDiophant(U1, cm, Vh, d, k);
590                //System.out.println("Ud = " + Ud);
591
592                mon = mon.multiply(xv);
593                //System.out.println("mon  = " + mon);
594                //List<GenPolynomial<MOD>> Sd = new ArrayList<GenPolynomial<MOD>>(Ud.size());
595                int i = 0;
596                List<GenPolynomial<BigInteger>> Si = new ArrayList<GenPolynomial<BigInteger>>(Ud.size());
597                for (GenPolynomial<MOD> dd : Ud) {
598                    //System.out.println("dd = " + dd);
599                    GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L);
600                    GenPolynomial<MOD> dm = de.multiply(mon);
601                    //Sd.add(dm);
602                    de = U.get(i).sum(dm);
603                    //System.out.println("de = " + de);
604                    U.set(i++, de);
605                    GenPolynomial<BigInteger> si = PolyUtil.integerFromModularCoefficients(ifac, de);
606                    Si.add(si);
607                }
608                //System.out.println("Sd   = " + Sd);
609                //System.out.println("U    = " + U);
610                //System.out.println("Si   = " + Si);
611
612                // compute new error:
613                E = ifac.getONE();
614                for (GenPolynomial<BigInteger> bi : Si) {
615                    E = E.multiply(bi);
616                }
617                E = Ci.subtract(E);
618                //System.out.println("E     = " + E);
619                Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
620                //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep);
621                logger.info("Ep({},{},{}) = {}", e, deg, pkfac.nvar, Ep);
622            }
623            Vh.add(v);
624            U1 = U;
625            if (Pfac.size() > 0) {
626                List<GenPolynomial<MOD>> U2 = new ArrayList<GenPolynomial<MOD>>(U.size());
627                pkfac = Pfac.get(0);
628                for (GenPolynomial<MOD> b : U) {
629                    GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L);
630                    U2.add(bi);
631                }
632                U = U2;
633                //System.out.println("U  = " + U);
634            }
635        }
636        if (E.isZERO()) {
637            logger.info("liftHensel leaving with zero E");
638        }
639        return U;
640    }
641
642
643    /**
644     * Modular Hensel lifting algorithm. Let p = A_i.ring.coFac.modul() and
645     * assume ggt(a,b) == 1 mod p, for a, b in A.
646     * @param C GenPolynomial with integer coefficients
647     * @param Cp GenPolynomial C mod p^k
648     * @param F list of modular GenPolynomials, mod (I_v, p^k )
649     * @param V list of integral substitution values
650     * @param k desired approximation exponent p^k.
651     * @param G list of leading coefficients of the factors of C.
652     * @return [g'_1,..., g'_n] with prod_i g'_i = Cp mod p^k.
653     */
654    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHensel(
655                    GenPolynomial<BigInteger> C, GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F,
656                    List<BigInteger> V, long k, List<GenPolynomial<BigInteger>> G) throws NoLiftingException {
657        GenPolynomialRing<MOD> pkfac = Cp.ring;
658        long d = C.degree();
659        //System.out.println("C = " + C);
660        //System.out.println("Cp = " + Cp);
661        //System.out.println("G = " + G);
662
663        //GenPolynomial<BigInteger> cd = G.get(0); // 1
664        //System.out.println("cd = " + cd + ", ring = " + C.ring);
665        //if ( cd.equals(C.ring.univariate(0)) ) {
666        //    System.out.println("cd == G[1]");
667        //}
668        // G mod p^k, in all variables
669        GenPolynomialRing<MOD> pkfac1 = new GenPolynomialRing<MOD>(pkfac.coFac, G.get(0).ring);
670        List<GenPolynomial<MOD>> Lp = new ArrayList<GenPolynomial<MOD>>(G.size());
671        for (GenPolynomial<BigInteger> cd1 : G) {
672            GenPolynomial<MOD> cdq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac1, cd1);
673            cdq = cdq.extendLower(pkfac, 0, 0L); // reintroduce lower variable
674            Lp.add(cdq);
675        }
676        logger.info("G modulo p^k: {}", Lp); // + ", ring = {}", pkfac1);
677
678        // prepare stack of polynomial rings, polynomials and evaluated leading coefficients
679        List<GenPolynomialRing<MOD>> Pfac = new ArrayList<GenPolynomialRing<MOD>>();
680        List<GenPolynomial<MOD>> Ap = new ArrayList<GenPolynomial<MOD>>();
681        List<List<GenPolynomial<MOD>>> Gp = new ArrayList<List<GenPolynomial<MOD>>>();
682        List<MOD> Vb = new ArrayList<MOD>();
683        //MOD v = V.get(0); // fromInteger
684        Pfac.add(pkfac);
685        Ap.add(Cp);
686        Gp.add(Lp);
687        GenPolynomialRing<MOD> pf = pkfac;
688        //GenPolynomialRing<MOD> pf1 = pkfac1;
689        GenPolynomial<MOD> ap = Cp;
690        List<GenPolynomial<MOD>> Lpp = Lp;
691        for (int j = pkfac.nvar; j > 2; j--) {
692            pf = pf.contract(1);
693            Pfac.add(0, pf);
694            //MOD vp = pkfac.coFac.fromInteger(V.get(pkfac.nvar - j).getSymmetricInteger().getVal());
695            MOD vp = pkfac.coFac.fromInteger(V.get(pkfac.nvar - j).getVal());
696            //System.out.println("vp     = " + vp);
697            Vb.add(vp);
698            ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp);
699            Ap.add(0, ap);
700            List<GenPolynomial<MOD>> Lps = new ArrayList<GenPolynomial<MOD>>(Lpp.size());
701            for (GenPolynomial<MOD> qp : Lpp) {
702                GenPolynomial<MOD> qpe = PolyUtil.<MOD> evaluateMain(pf, qp, vp);
703                Lps.add(qpe);
704            }
705            //System.out.println("Lps = " + Lps);
706            Lpp = Lps;
707            Gp.add(0, Lpp);
708        }
709        Vb.add(pkfac.coFac.fromInteger(V.get(pkfac.nvar - 2).getVal()));
710        //System.out.println("Pfac   = " + Pfac);
711        logger.debug("Pfac   = {}", Pfac);
712        //System.out.println("Ap     = " + Ap);
713        //System.out.println("Gp     = " + Gp);
714        //System.out.println("Gp[0]  = " + Gp.get(0) + ", Gp[0].ring = " + Gp.get(0).get(0).ring);
715        //System.out.println("V      = " + V);
716        //System.out.println("Vb     = " + Vb + ", V == Vb: " + V.equals(Vb));
717
718        // check bi-variate base case
719        GenPolynomialRing<MOD> pk1fac = F.get(0).ring;
720        if (!pkfac.coFac.equals(pk1fac.coFac)) {
721            throw new IllegalArgumentException("F.ring != pkfac: " + pk1fac + " != " + pkfac);
722        }
723
724        // init recursion
725        List<GenPolynomial<MOD>> U = F;
726        //logger.info("to lift U = {}", U); // + ", U1.ring = {}", U1.get(0).ring);
727        GenPolynomial<BigInteger> E = C.ring.getZERO();
728        List<MOD> Vh = new ArrayList<MOD>();
729        List<GenPolynomial<BigInteger>> Si; // = new ArrayList<GenPolynomial<BigInteger>>(F.size());
730        MOD v = null;
731
732        while (Pfac.size() > 0) { // loop through stack of polynomial rings
733            pkfac = Pfac.remove(0);
734            Cp = Ap.remove(0);
735            Lpp = Gp.remove(0);
736            v = Vb.remove(Vb.size() - 1); // last in stack
737            //System.out.println("\npkfac = " + pkfac.toScript() + " ================================== " + v);
738            logger.info("stack loop: pkfac = {} v = {}", pkfac.toScript(), v);
739
740            List<GenPolynomial<MOD>> U1 = U;
741            logger.info("to lift U1 = {}", U1); // + ", U1.ring = {}", U1.get(0).ring);
742            U = new ArrayList<GenPolynomial<MOD>>(U1.size());
743
744            // update U, replace leading coefficient if required
745            int j = 0;
746            for (GenPolynomial<MOD> b : U1) {
747                //System.out.println("b = " + b + ", b.ring = " + b.ring);
748                GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L);
749                GenPolynomial<MOD> li = Lpp.get(j);
750                if (!li.isONE()) {
751                    //System.out.println("li = " + li + ", li.ring = " + li.ring);
752                    //System.out.println("bi = " + bi);
753                    GenPolynomialRing<GenPolynomial<MOD>> pkrfac = pkfac.recursive(pkfac.nvar - 1);
754                    //System.out.println("pkrfac = " + pkrfac);
755                    GenPolynomial<GenPolynomial<MOD>> br = PolyUtil.<MOD> recursive(pkrfac, bi);
756                    //System.out.println("br = " + br);
757                    GenPolynomial<GenPolynomial<MOD>> bs = PolyUtil.<MOD> switchVariables(br);
758                    //System.out.println("bs = " + bs + ", bs.ring = " + bs.ring);
759
760                    GenPolynomial<GenPolynomial<MOD>> lr = PolyUtil.<MOD> recursive(pkrfac, li);
761                    //System.out.println("lr = " + lr);
762                    GenPolynomial<GenPolynomial<MOD>> ls = PolyUtil.<MOD> switchVariables(lr);
763                    //System.out.println("ls = " + ls + ", ls.ring = " + ls.ring);
764                    if (!ls.isConstant() && !ls.isZERO()) {
765                        throw new RuntimeException("ls not constant " + ls + ", li = " + li);
766                    }
767                    bs.doPutToMap(bs.leadingExpVector(), ls.leadingBaseCoefficient());
768                    //System.out.println("bs = " + bs + ", bs.ring = " + bs.ring);
769                    br = PolyUtil.<MOD> switchVariables(bs);
770                    //System.out.println("br = " + br);
771                    bi = PolyUtil.<MOD> distribute(pkfac, br);
772                    //System.out.println("bi = " + bi);
773                }
774                U.add(bi);
775                j++;
776            }
777            logger.info("U with leading coefficient replaced = {}", U); // + ", U.ring = {}", U.get(0).ring);
778
779            // (x_n - v)
780            GenPolynomial<MOD> mon = pkfac.getONE();
781            GenPolynomial<MOD> xv = pkfac.univariate(0, 1);
782            xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal()));
783            //System.out.println("xv = " + xv);
784
785            long deg = Cp.degree(pkfac.nvar - 1);
786            //System.out.println("deg = " + deg + ", degv = " + Cp.degreeVector());
787
788            // convert to integer polynomials
789            GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac);
790            //System.out.println("ifac = " + ifac.toScript());
791            List<GenPolynomial<BigInteger>> Bi = PolyUtil.integerFromModularCoefficients(ifac, U);
792            //System.out.println("Bi = " + Bi);
793            GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cp);
794            //System.out.println("Ci = " + Ci);
795
796            // compute error:
797            E = ifac.getONE();
798            for (GenPolynomial<BigInteger> bi : Bi) {
799                E = E.multiply(bi);
800            }
801            //System.out.println("E  = " + E);
802            E = Ci.subtract(E);
803            //System.out.println("E  = " + E);
804            GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
805            logger.info("Ep(0,{},{}) = {}", deg, pkfac.nvar, Ep);
806
807            GenPolynomialRing<GenPolynomial<MOD>> pkrfac = pkfac.recursive(1);
808            GenPolynomialRing<MOD> ckfac = (GenPolynomialRing<MOD>) pkrfac.coFac;
809            //System.out.println("pkrfac = " + pkrfac.toScript());
810
811            for (int e = 1; e <= deg && !Ep.isZERO(); e++) {
812                //System.out.println("\ne = " + e + " -------------------------------------- " + deg);
813                logger.info("approximation loop: e = {} of deg = {}", e, deg);
814                GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(pkrfac, Ep);
815                //System.out.println("Epr   = " + Epr);
816                UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>(
817                                pkrfac);
818                //System.out.println("psfac = " + psfac);
819                TaylorFunction<GenPolynomial<MOD>> T = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr);
820                //System.out.println("T     = " + T);
821                GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal());
822                //System.out.println("vq    = " + vq + ", Vh = " + Vh);
823                UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(T, vq);
824                //System.out.println("Epst  = " + Epst);
825                logger.info("Epst({},{},{}) = {}", e, deg, pkfac.nvar, Epst);
826                GenPolynomial<MOD> cm = Epst.coefficient(e);
827                if (cm.isZERO()) {
828                    //System.out.println("cm   = " + cm);
829                    continue;
830                }
831                List<GenPolynomial<MOD>> Ud = HenselMultUtil.<MOD> liftDiophant(U1, cm, Vh, d, k);
832                //System.out.println("Ud = " + Ud);
833
834                mon = mon.multiply(xv);
835                //System.out.println("mon  = " + mon);
836                //List<GenPolynomial<MOD>> Sd = new ArrayList<GenPolynomial<MOD>>(Ud.size());
837                int i = 0;
838                Si = new ArrayList<GenPolynomial<BigInteger>>(Ud.size());
839                for (GenPolynomial<MOD> dd : Ud) {
840                    //System.out.println("dd = " + dd);
841                    GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L);
842                    GenPolynomial<MOD> dm = de.multiply(mon);
843                    //Sd.add(dm);
844                    de = U.get(i).sum(dm);
845                    //System.out.println("de = " + de);
846                    U.set(i++, de);
847                    GenPolynomial<BigInteger> si = PolyUtil.integerFromModularCoefficients(ifac, de);
848                    Si.add(si);
849                }
850                //System.out.println("Sd   = " + Sd);
851                //System.out.println("U    = " + U + ", U.ring = " + U.get(0).ring);
852                //System.out.println("Si   = " + Si);
853
854                // compute new error:
855                E = ifac.getONE();
856                for (GenPolynomial<BigInteger> bi : Si) {
857                    E = E.multiply(bi);
858                }
859                E = Ci.subtract(E);
860                //System.out.println("E = " + E);
861                Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
862                //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep);
863                logger.info("Ep({},{},{}) = {}", e, deg, pkfac.nvar, Ep);
864            }
865            Vh.add(v);
866            GenPolynomial<MOD> Uf = U.get(0).ring.getONE();
867            for (GenPolynomial<MOD> Upp : U) {
868                Uf = Uf.multiply(Upp);
869            }
870            if (false && !Cp.leadingExpVector().equals(Uf.leadingExpVector())) { // not meaningful test
871                System.out.println("\nU    = " + U);
872                System.out.println("Cp   = " + Cp);
873                System.out.println("Uf   = " + Uf);
874                //System.out.println("Cp.ring = " + Cp.ring.toScript() + ", Uf.ring = " + Uf.ring.toScript() + "\n");
875                System.out.println("");
876                //throw new NoLiftingException("no factorization, Cp != Uf");
877            }
878        }
879        if (E.isZERO()) {
880            logger.info("liftHensel leaving with zero E, Ep");
881        }
882        if (false && debug) {
883            // remove normalization required ??
884            GreatestCommonDivisorAbstract<BigInteger> ufd = GCDFactory.getImplementation(new BigInteger());
885            List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(U.size());
886            for (GenPolynomial<BigInteger> bi : Si) {
887                GenPolynomial<BigInteger> ci = ufd.content(bi); //ufd.primitivePart(bi); // ??
888                if (!ci.isONE()) {
889                    System.out.println("bi = " + bi + ", cont(bi) = " + ci);
890                }
891                //Fii.add(ci);
892            }
893            //Si = Fii;
894            //System.out.println("Si  = " + Si);
895        }
896        logger.info("multivariate lift: U = {}, of {}", U, F);
897        return U;
898    }
899
900
901    /**
902     * Modular Hensel full lifting algorithm. Let p = A_i.ring.coFac.modul() and
903     * assume ggt(a,b) == 1 mod p, for a, b in A.
904     * @param C GenPolynomial with integer coefficients
905     * @param F list of modular GenPolynomials, mod (I_v, p )
906     * @param V list of integer substitution values
907     * @param k desired approximation exponent p^k.
908     * @param G = [g_1,...,g_n] list of factors of leading coefficients.
909     * @return [c_1,..., c_n] with prod_i c_i = C mod p^k.
910     */
911    @SuppressWarnings("unchecked")
912    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHenselFull(
913                    GenPolynomial<BigInteger> C, List<GenPolynomial<MOD>> F, List<BigInteger> V, long k,
914                    List<GenPolynomial<BigInteger>> G) throws NoLiftingException {
915        if (F == null || F.size() == 0) {
916            return new ArrayList<GenPolynomial<MOD>>();
917        }
918        GenPolynomialRing<MOD> pkfac = F.get(0).ring;
919        //long d = C.degree();
920        // setup q = p^k
921        RingFactory<MOD> cfac = pkfac.coFac;
922        ModularRingFactory<MOD> pcfac = (ModularRingFactory<MOD>) cfac;
923        //System.out.println("pcfac = " + pcfac);
924        BigInteger p = pcfac.getIntegerModul();
925        BigInteger q = p.power(k);
926        ModularRingFactory<MOD> mcfac;
927        if (ModLongRing.MAX_LONG.compareTo(q.getVal()) > 0) {
928            mcfac = (ModularRingFactory) new ModLongRing(q.getVal());
929        } else {
930            mcfac = (ModularRingFactory) new ModIntegerRing(q.getVal());
931        }
932        //System.out.println("mcfac = " + mcfac);
933
934        // convert C from Z[...] to Z_q[...]
935        GenPolynomialRing<MOD> qcfac = new GenPolynomialRing<MOD>(mcfac, C.ring);
936        GenPolynomial<MOD> Cq = PolyUtil.<MOD> fromIntegerCoefficients(qcfac, C);
937        //System.out.println("C  = " + C);
938        //System.out.println("Cq = " + Cq);
939
940        // convert g_i from Z[...] to Z_q[...]
941        GenPolynomialRing<MOD> gcfac = new GenPolynomialRing<MOD>(mcfac, G.get(0).ring);
942        List<GenPolynomial<MOD>> GQ = new ArrayList<GenPolynomial<MOD>>();
943        boolean allOnes = true;
944        for (GenPolynomial<BigInteger> g : G) {
945            if (!g.isONE()) {
946                allOnes = false;
947            }
948            GenPolynomial<MOD> gq = PolyUtil.<MOD> fromIntegerCoefficients(gcfac, g);
949            GQ.add(gq);
950        }
951        //System.out.println("G  = " + G);
952        //System.out.println("GQ = " + GQ);
953
954        // evaluate C to Z_q[x]
955        GenPolynomialRing<MOD> pf = qcfac;
956        GenPolynomial<MOD> ap = Cq;
957        for (int j = C.ring.nvar; j > 1; j--) {
958            pf = pf.contract(1);
959            //MOD vp = mcfac.fromInteger(V.get(C.ring.nvar - j).getSymmetricInteger().getVal());
960            MOD vp = mcfac.fromInteger(V.get(C.ring.nvar - j).getVal());
961            //System.out.println("vp     = " + vp);
962            ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp);
963            //System.out.println("ap     = " + ap);
964        }
965        GenPolynomial<MOD> Cq1 = ap;
966        //System.out.println("Cq1 = " + Cq1);
967        if (Cq1.isZERO()) {
968            throw new NoLiftingException("C mod (I, p^k) == 0: " + C);
969        }
970        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pf);
971        GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cq1);
972        //System.out.println("Ci  = " + Ci);
973        GreatestCommonDivisorAbstract<BigInteger> ufd = GCDFactory.getImplementation(new BigInteger());
974        Ci = Ci.abs();
975        BigInteger cCi = ufd.baseContent(Ci);
976        Ci = Ci.divide(cCi);
977        //System.out.println("cCi = " + cCi);
978        //System.out.println("Ci  = " + Ci);
979        ////System.out.println("F.fac = " + F.get(0).ring);
980
981        // evaluate G to Z_q
982        //List<GenPolynomial<MOD>> GP = new ArrayList<GenPolynomial<MOD>>();
983        for (GenPolynomial<MOD> gq : GQ) {
984            GenPolynomialRing<MOD> gf = gcfac;
985            GenPolynomial<MOD> gp = gq;
986            for (int j = gcfac.nvar; j > 1; j--) {
987                gf = gf.contract(1);
988                //MOD vp = mcfac.fromInteger(V.get(gcfac.nvar - j).getSymmetricInteger().getVal());
989                MOD vp = mcfac.fromInteger(V.get(gcfac.nvar - j).getVal());
990                //System.out.println("vp     = " + vp);
991                gp = PolyUtil.<MOD> evaluateMain(gf, gp, vp);
992                //System.out.println("gp     = " + gp);
993            }
994            //GP.add(gp);
995        }
996        //System.out.println("GP = " + GP); // + ", GP.ring = " + GP.get(0).ring);
997
998        // leading coefficient for recursion base, for Cq1 and list GP 
999        BigInteger gi0 = Ci.leadingBaseCoefficient(); // gq0.getSymmetricInteger();
1000        //System.out.println("gi0 = " + gi0);
1001
1002        // lift F to Z_{p^k}[x]
1003        //System.out.println("Ci = " + Ci + ", F = " + F + ", k = " + k + ", p = " + F.get(0).ring + ", gi0 = " + gi0);
1004        List<GenPolynomial<MOD>> U1 = null;
1005        if (gi0.isONE()) {
1006            U1 = HenselUtil.<MOD> liftHenselMonic(Ci, F, k);
1007        } else {
1008            U1 = HenselUtil.<MOD> liftHensel(Ci, F, k, gi0); // gi0 TODO ??
1009        }
1010        logger.info("univariate lift: Ci = {}, F = {}, U1 = {}", Ci, F, U1);
1011        //System.out.println("U1.fac = " + U1.get(0).ring);
1012
1013        // adjust leading coefficients of U1 with F
1014        List<GenPolynomial<BigInteger>> U1i = PolyUtil.<MOD> integerFromModularCoefficients(Ci.ring, U1);
1015        //System.out.println("U1i = " + U1i);
1016        boolean t = HenselUtil.isHenselLift(Ci, q, p, U1i);
1017        //System.out.println("isLift(U1) = " + t);
1018        if (!t) {
1019            //System.out.println("NoLiftingException, Ci = " + Ci + ", U1i = " + U1i);
1020            throw new NoLiftingException("Ci = " + Ci + ", U1i = " + U1i);
1021        }
1022        MOD cC = mcfac.fromInteger(cCi.getVal());
1023        List<GenPolynomial<MOD>> U1f = PolyUtil.<MOD> fromIntegerCoefficients(F.get(0).ring, U1i);
1024        //System.out.println("U1f = " + U1f);
1025        List<GenPolynomial<MOD>> U1s = new ArrayList<GenPolynomial<MOD>>(U1.size());
1026        int j = 0;
1027        int s = 0;
1028        for (GenPolynomial<MOD> u : U1) {
1029            GenPolynomial<MOD> uf = U1f.get(j);
1030            GenPolynomial<MOD> f = F.get(j);
1031            GenPolynomial<BigInteger> ui = U1i.get(j);
1032            GenPolynomial<BigInteger> gi = G.get(j);
1033            if (ui.signum() != gi.signum()) {
1034                //System.out.println("ui = " + ui + ", gi = " + gi);
1035                u = u.negate();
1036                uf = uf.negate();
1037                s++;
1038            }
1039            j++;
1040            if (uf.isConstant()) {
1041                //System.out.println("u   = " + u);
1042                u = u.monic();
1043                //System.out.println("u  = " + u);
1044                u = u.multiply(cC);
1045                cC = cC.divide(cC);
1046                //System.out.println("u   = " + u);
1047            } else {
1048                MOD x = f.leadingBaseCoefficient().divide(uf.leadingBaseCoefficient());
1049                //System.out.println("x   = " + x + ", xi = " + x.getSymmetricInteger());
1050                if (!x.isONE()) {
1051                    MOD xq = mcfac.fromInteger(x.getSymmetricInteger().getVal());
1052                    //System.out.println("xq  = " + xq);
1053                    u = u.multiply(xq);
1054                    cC = cC.divide(xq);
1055                    //System.out.println("cC  = " + cC);
1056                }
1057            }
1058            U1s.add(u);
1059        }
1060        //if ( s % 2 != 0 || !cC.isONE()) {
1061        if (!cC.isONE()) {
1062            throw new NoLiftingException("s = " + s + ", Ci = " + Ci + ", U1i = " + U1i + ", cC = " + cC);
1063        }
1064        U1 = U1s;
1065        U1i = PolyUtil.<MOD> integerFromModularCoefficients(Ci.ring, U1);
1066        //System.out.println("U1i = " + U1i);
1067        U1f = PolyUtil.<MOD> fromIntegerCoefficients(F.get(0).ring, U1i);
1068        if (!F.equals(U1f)) { // evtl loop until reached
1069            System.out.println("F   = " + F);
1070            System.out.println("U1f = " + U1f);
1071            throw new NoLiftingException("F = " + F + ", U1f = " + U1f);
1072        }
1073        logger.info("multivariate lift: U1 = {}", U1);
1074
1075        // lift U to Z_{p^k}[x,...]
1076        //System.out.println("C = " + C + ", U1 = " + U1 + ", V = " + V + ", k = " + k + ", q = " + U1.get(0).ring + ", G = " + G);
1077        List<GenPolynomial<MOD>> U = null;
1078        if (allOnes) {
1079            U = HenselMultUtil.<MOD> liftHenselMonic(C, Cq, U1, V, k);
1080        } else {
1081            U = HenselMultUtil.<MOD> liftHensel(C, Cq, U1, V, k, G);
1082        }
1083        logger.info("multivariate lift: C = {}, U1 = {}, U = {}", C, U1, U);
1084        //System.out.println("U  = " + U);
1085        //System.out.println("U.fac = " + U.get(0).ring);
1086        return U;
1087    }
1088
1089}