Source code

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