Source code

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