Source code

001/*
002 * $Id$
003 */
004
005package edu.jas.ufd;
006
007
008import java.util.ArrayList;
009import java.util.List;
010
011import org.apache.logging.log4j.Logger;
012import org.apache.logging.log4j.LogManager; 
013
014import edu.jas.arith.BigInteger;
015import edu.jas.arith.ModInteger;
016import edu.jas.arith.ModIntegerRing;
017import edu.jas.arith.ModLongRing;
018import edu.jas.arith.Modular;
019import edu.jas.arith.ModularRingFactory;
020import edu.jas.poly.ExpVector;
021import edu.jas.poly.GenPolynomial;
022import edu.jas.poly.GenPolynomialRing;
023import edu.jas.poly.Monomial;
024import edu.jas.poly.PolyUtil;
025import edu.jas.structure.GcdRingElem;
026import edu.jas.structure.RingFactory;
027
028
029/**
030 * Hensel utilities for ufd.
031 * @author Heinz Kredel
032 */
033
034public class HenselUtil {
035
036
037    private static final Logger logger = LogManager.getLogger(HenselUtil.class);
038
039
040    private static final boolean debug = logger.isDebugEnabled();
041
042
043    /**
044     * Modular Hensel lifting algorithm on coefficients. Let p =
045     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
046     * with gcd(A,B) == 1 mod p and S A + T B == 1 mod p. See Algorithm 6.1. in
047     * Geddes et.al.. Linear version, as it does not lift S A + T B == 1 mod
048     * p^{e+1}.
049     * @param C GenPolynomial
050     * @param A GenPolynomial
051     * @param B other GenPolynomial
052     * @param S GenPolynomial
053     * @param T GenPolynomial
054     * @param M bound on the coefficients of A1 and B1 as factors of C.
055     * @return [A1,B1,Am,Bm] = lift(C,A,B), with C = A1 * B1 mod p^e, Am = A1
056     *         mod p^e, Bm = B1 mod p^e .
057     */
058    @SuppressWarnings("unchecked")
059    public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHensel(
060                    GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B,
061                    GenPolynomial<MOD> S, GenPolynomial<MOD> T) throws NoLiftingException {
062        if (C == null || C.isZERO()) {
063            return new HenselApprox<MOD>(C, C, A, B);
064        }
065        if (A == null || A.isZERO() || B == null || B.isZERO()) {
066            throw new IllegalArgumentException("A and B must be nonzero");
067        }
068        GenPolynomialRing<BigInteger> fac = C.ring;
069        if (fac.nvar != 1) { // assert ?
070            throw new IllegalArgumentException("polynomial ring not univariate");
071        }
072        // setup factories
073        GenPolynomialRing<MOD> pfac = A.ring;
074        RingFactory<MOD> p = pfac.coFac;
075        RingFactory<MOD> q = p;
076        ModularRingFactory<MOD> P = (ModularRingFactory<MOD>) p;
077        ModularRingFactory<MOD> Q = (ModularRingFactory<MOD>) q;
078        BigInteger Qi = Q.getIntegerModul();
079        BigInteger M2 = M.multiply(M.fromInteger(2));
080        BigInteger Mq = Qi;
081
082        // normalize c and a, b factors, assert p is prime
083        GenPolynomial<BigInteger> Ai;
084        GenPolynomial<BigInteger> Bi;
085        BigInteger c = C.leadingBaseCoefficient();
086        C = C.multiply(c); // sic
087        MOD a = A.leadingBaseCoefficient();
088        if (!a.isONE()) { // A = A.monic();
089            A = A.divide(a);
090            S = S.multiply(a);
091        }
092        MOD b = B.leadingBaseCoefficient();
093        if (!b.isONE()) { // B = B.monic();
094            B = B.divide(b);
095            T = T.multiply(b);
096        }
097        MOD cm = P.fromInteger(c.getVal());
098        A = A.multiply(cm);
099        B = B.multiply(cm);
100        T = T.divide(cm);
101        S = S.divide(cm);
102        Ai = PolyUtil.integerFromModularCoefficients(fac, A);
103        Bi = PolyUtil.integerFromModularCoefficients(fac, B);
104        // replace leading base coefficients
105        ExpVector ea = Ai.leadingExpVector();
106        ExpVector eb = Bi.leadingExpVector();
107        Ai.doPutToMap(ea, c);
108        Bi.doPutToMap(eb, c);
109
110        // polynomials mod p
111        GenPolynomial<MOD> Ap;
112        GenPolynomial<MOD> Bp;
113        GenPolynomial<MOD> A1p = A;
114        GenPolynomial<MOD> B1p = B;
115        GenPolynomial<MOD> Ep;
116
117        // polynomials over the integers
118        GenPolynomial<BigInteger> E;
119        GenPolynomial<BigInteger> Ea;
120        GenPolynomial<BigInteger> Eb;
121        GenPolynomial<BigInteger> Ea1;
122        GenPolynomial<BigInteger> Eb1;
123
124        while (Mq.compareTo(M2) < 0) {
125            // compute E=(C-AB)/q over the integers
126            E = C.subtract(Ai.multiply(Bi));
127            if (E.isZERO()) {
128                logger.info("leaving on zero E");
129                break;
130            }
131            try {
132                E = E.divide(Qi);
133            } catch (RuntimeException e) {
134                // useful in debugging
135                //System.out.println("C  = " + C );
136                //System.out.println("Ai = " + Ai );
137                //System.out.println("Bi = " + Bi );
138                //System.out.println("E  = " + E );
139                //System.out.println("Qi = " + Qi );
140                throw e;
141            }
142            // E mod p
143            Ep = PolyUtil.<MOD> fromIntegerCoefficients(pfac, E);
144            //logger.info("Ep = {}", Ep);
145
146            // construct approximation mod p
147            Ap = S.multiply(Ep); // S,T ++ T,S
148            Bp = T.multiply(Ep);
149            GenPolynomial<MOD>[] QR;
150            QR = Ap.quotientRemainder(B);
151            GenPolynomial<MOD> Qp;
152            GenPolynomial<MOD> Rp;
153            Qp = QR[0];
154            Rp = QR[1];
155            A1p = Rp;
156            B1p = Bp.sum(A.multiply(Qp));
157
158            // construct q-adic approximation, convert to integer
159            Ea = PolyUtil.integerFromModularCoefficients(fac, A1p);
160            Eb = PolyUtil.integerFromModularCoefficients(fac, B1p);
161            Ea1 = Ea.multiply(Qi);
162            Eb1 = Eb.multiply(Qi);
163
164            Ea = Ai.sum(Eb1); // Eb1 and Ea1 are required
165            Eb = Bi.sum(Ea1); //--------------------------
166            assert (Ea.degree(0) + Eb.degree(0) <= C.degree(0));
167            //if ( Ea.degree(0)+Eb.degree(0) > C.degree(0) ) { // debug
168            //   throw new RuntimeException("deg(A)+deg(B) > deg(C)");
169            //}
170
171            // prepare for next iteration
172            Mq = Qi;
173            Qi = Q.getIntegerModul().multiply(P.getIntegerModul());
174            // Q = new ModIntegerRing(Qi.getVal());
175            if (ModLongRing.MAX_LONG.compareTo(Qi.getVal()) > 0) {
176                Q = (ModularRingFactory) new ModLongRing(Qi.getVal());
177            } else {
178                Q = (ModularRingFactory) new ModIntegerRing(Qi.getVal());
179            }
180            Ai = Ea;
181            Bi = Eb;
182        }
183        GreatestCommonDivisorAbstract<BigInteger> ufd = new GreatestCommonDivisorPrimitive<BigInteger>();
184
185        // remove normalization
186        BigInteger ai = ufd.baseContent(Ai);
187        Ai = Ai.divide(ai);
188        BigInteger bi = null;
189        try {
190            bi = c.divide(ai);
191            Bi = Bi.divide(bi); // divide( c/a )
192        } catch (RuntimeException e) {
193            //System.out.println("C  = " + C );
194            //System.out.println("Ai = " + Ai );
195            //System.out.println("Bi = " + Bi );
196            //System.out.println("c  = " + c );
197            //System.out.println("ai = " + ai );
198            //System.out.println("bi = " + bi );
199            //System.out.println("no exact lifting possible " + e);
200            throw new NoLiftingException("no exact lifting possible " + e);
201        }
202        return new HenselApprox<MOD>(Ai, Bi, A1p, B1p);
203    }
204
205
206    /**
207     * Modular Hensel lifting algorithm on coefficients. Let p =
208     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
209     * with gcd(A,B) == 1 mod p. See algorithm 6.1. in Geddes et.al. and
210     * algorithms 3.5.{5,6} in Cohen.
211     * @param C GenPolynomial
212     * @param A GenPolynomial
213     * @param B other GenPolynomial
214     * @param M bound on the coefficients of A1 and B1 as factors of C.
215     * @return [A1,B1] = lift(C,A,B), with C = A1 * B1.
216     */
217    @SuppressWarnings("unchecked")
218    public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHensel(
219                    GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B)
220                    throws NoLiftingException {
221        if (C == null || C.isZERO()) {
222            return new HenselApprox<MOD>(C, C, A, B);
223        }
224        if (A == null || A.isZERO() || B == null || B.isZERO()) {
225            throw new IllegalArgumentException("A and B must be nonzero");
226        }
227        GenPolynomialRing<BigInteger> fac = C.ring;
228        if (fac.nvar != 1) { // assert ?
229            throw new IllegalArgumentException("polynomial ring not univariate");
230        }
231        // one Hensel step on part polynomials
232        try {
233            GenPolynomial<MOD>[] gst = A.egcd(B);
234            if (!gst[0].isONE()) {
235                throw new NoLiftingException(
236                                "A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B);
237            }
238            GenPolynomial<MOD> s = gst[1];
239            GenPolynomial<MOD> t = gst[2];
240            HenselApprox<MOD> ab = HenselUtil.<MOD> liftHensel(C, M, A, B, s, t);
241            return ab;
242        } catch (ArithmeticException e) {
243            throw new NoLiftingException("coefficient error " + e);
244        }
245    }
246
247
248    /**
249     * Modular quadratic Hensel lifting algorithm on coefficients. Let p =
250     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
251     * with gcd(A,B) == 1 mod p and S A + T B == 1 mod p. See algorithm 6.1. in
252     * Geddes et.al. and algorithms 3.5.{5,6} in Cohen. Quadratic version, as it
253     * also lifts S A + T B == 1 mod p^{e+1}.
254     * @param C GenPolynomial
255     * @param A GenPolynomial
256     * @param B other GenPolynomial
257     * @param S GenPolynomial
258     * @param T GenPolynomial
259     * @param M bound on the coefficients of A1 and B1 as factors of C.
260     * @return [A1,B1] = lift(C,A,B), with C = A1 * B1.
261     */
262    @SuppressWarnings("unchecked")
263    public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHenselQuadratic(
264                    GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B,
265                    GenPolynomial<MOD> S, GenPolynomial<MOD> T) throws NoLiftingException {
266        if (C == null || C.isZERO()) {
267            return new HenselApprox<MOD>(C, C, A, B);
268        }
269        if (A == null || A.isZERO() || B == null || B.isZERO()) {
270            throw new IllegalArgumentException("A and B must be nonzero");
271        }
272        GenPolynomialRing<BigInteger> fac = C.ring;
273        if (fac.nvar != 1) { // assert ?
274            throw new IllegalArgumentException("polynomial ring not univariate");
275        }
276        // setup factories
277        GenPolynomialRing<MOD> pfac = A.ring;
278        RingFactory<MOD> p = pfac.coFac;
279        RingFactory<MOD> q = p;
280        ModularRingFactory<MOD> P = (ModularRingFactory<MOD>) p;
281        ModularRingFactory<MOD> Q = (ModularRingFactory<MOD>) q;
282        BigInteger Qi = Q.getIntegerModul();
283        BigInteger M2 = M.multiply(M.fromInteger(2));
284        BigInteger Mq = Qi;
285        GenPolynomialRing<MOD> qfac;
286        qfac = new GenPolynomialRing<MOD>(Q, pfac);
287
288        // normalize c and a, b factors, assert p is prime
289        GenPolynomial<BigInteger> Ai;
290        GenPolynomial<BigInteger> Bi;
291        BigInteger c = C.leadingBaseCoefficient();
292        C = C.multiply(c); // sic
293        MOD a = A.leadingBaseCoefficient();
294        if (!a.isONE()) { // A = A.monic();
295            A = A.divide(a);
296            S = S.multiply(a);
297        }
298        MOD b = B.leadingBaseCoefficient();
299        if (!b.isONE()) { // B = B.monic();
300            B = B.divide(b);
301            T = T.multiply(b);
302        }
303        MOD cm = P.fromInteger(c.getVal());
304        A = A.multiply(cm);
305        B = B.multiply(cm);
306        T = T.divide(cm);
307        S = S.divide(cm);
308        Ai = PolyUtil.integerFromModularCoefficients(fac, A);
309        Bi = PolyUtil.integerFromModularCoefficients(fac, B);
310        // replace leading base coefficients
311        ExpVector ea = Ai.leadingExpVector();
312        ExpVector eb = Bi.leadingExpVector();
313        Ai.doPutToMap(ea, c);
314        Bi.doPutToMap(eb, c);
315
316        // polynomials mod p
317        GenPolynomial<MOD> Ap;
318        GenPolynomial<MOD> Bp;
319        GenPolynomial<MOD> A1p = A;
320        GenPolynomial<MOD> B1p = B;
321        GenPolynomial<MOD> Ep;
322        GenPolynomial<MOD> Sp = S;
323        GenPolynomial<MOD> Tp = T;
324
325        // polynomials mod q
326        GenPolynomial<MOD> Aq;
327        GenPolynomial<MOD> Bq;
328        //GenPolynomial<MOD> Eq;
329
330        // polynomials over the integers
331        GenPolynomial<BigInteger> E;
332        GenPolynomial<BigInteger> Ea;
333        GenPolynomial<BigInteger> Eb;
334        GenPolynomial<BigInteger> Ea1;
335        GenPolynomial<BigInteger> Eb1;
336        GenPolynomial<BigInteger> Si;
337        GenPolynomial<BigInteger> Ti;
338
339        Si = PolyUtil.integerFromModularCoefficients(fac, S);
340        Ti = PolyUtil.integerFromModularCoefficients(fac, T);
341
342        Aq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ai);
343        Bq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Bi);
344
345        while (Mq.compareTo(M2) < 0) {
346            // compute E=(C-AB)/q over the integers
347            E = C.subtract(Ai.multiply(Bi));
348            if (E.isZERO()) {
349                logger.info("leaving on zero E");
350                break;
351            }
352            E = E.divide(Qi);
353            // E mod p
354            Ep = PolyUtil.<MOD> fromIntegerCoefficients(qfac, E);
355            //logger.info("Ep = {}, qfac = {}", Ep, qfac);
356            //if (Ep.isZERO()) {
357            //System.out.println("leaving on zero error");
358            //??logger.info("leaving on zero Ep");
359            //??break;
360            //}
361
362            // construct approximation mod p
363            Ap = Sp.multiply(Ep); // S,T ++ T,S
364            Bp = Tp.multiply(Ep);
365            GenPolynomial<MOD>[] QR;
366            //logger.info("Ap = {}, Bp = {}, fac(Ap) = {}", Ap, Bp, Ap.ring);
367            QR = Ap.quotientRemainder(Bq);
368            GenPolynomial<MOD> Qp;
369            GenPolynomial<MOD> Rp;
370            Qp = QR[0];
371            Rp = QR[1];
372            //logger.info("Qp = {}, Rp = {}", Qp, Rp);
373            A1p = Rp;
374            B1p = Bp.sum(Aq.multiply(Qp));
375
376            // construct q-adic approximation, convert to integer
377            Ea = PolyUtil.integerFromModularCoefficients(fac, A1p);
378            Eb = PolyUtil.integerFromModularCoefficients(fac, B1p);
379            Ea1 = Ea.multiply(Qi);
380            Eb1 = Eb.multiply(Qi);
381            Ea = Ai.sum(Eb1); // Eb1 and Ea1 are required
382            Eb = Bi.sum(Ea1); //--------------------------
383            assert (Ea.degree(0) + Eb.degree(0) <= C.degree(0));
384            //if ( Ea.degree(0)+Eb.degree(0) > C.degree(0) ) { // debug
385            //   throw new RuntimeException("deg(A)+deg(B) > deg(C)");
386            //}
387            Ai = Ea;
388            Bi = Eb;
389
390            // gcd representation factors error --------------------------------
391            // compute E=(1-SA-TB)/q over the integers
392            E = fac.getONE();
393            E = E.subtract(Si.multiply(Ai)).subtract(Ti.multiply(Bi));
394            E = E.divide(Qi);
395            // E mod q
396            Ep = PolyUtil.<MOD> fromIntegerCoefficients(qfac, E);
397            //logger.info("Ep2 = {}", Ep);
398
399            // construct approximation mod q
400            Ap = Sp.multiply(Ep); // S,T ++ T,S
401            Bp = Tp.multiply(Ep);
402            QR = Bp.quotientRemainder(Aq); // Ai == A mod p ?
403            Qp = QR[0];
404            Rp = QR[1];
405            B1p = Rp;
406            A1p = Ap.sum(Bq.multiply(Qp));
407
408            //if (debug) {
409            //    Eq = A1p.multiply(Aq).sum(B1p.multiply(Bq)).subtract(Ep);
410            //    if (!Eq.isZERO()) {
411            //        System.out.println("A*A1p+B*B1p-Ep2 != 0 " + Eq);
412            //        throw new RuntimeException("A*A1p+B*B1p-Ep2 != 0 mod " + Q.getIntegerModul());
413            //    }
414            //}
415
416            // construct q-adic approximation, convert to integer
417            Ea = PolyUtil.integerFromModularCoefficients(fac, A1p);
418            Eb = PolyUtil.integerFromModularCoefficients(fac, B1p);
419            Ea1 = Ea.multiply(Qi);
420            Eb1 = Eb.multiply(Qi);
421            Ea = Si.sum(Ea1); // Eb1 and Ea1 are required
422            Eb = Ti.sum(Eb1); //--------------------------
423            Si = Ea;
424            Ti = Eb;
425
426            // prepare for next iteration
427            Mq = Qi;
428            Qi = Q.getIntegerModul().multiply(Q.getIntegerModul());
429            if (ModLongRing.MAX_LONG.compareTo(Qi.getVal()) > 0) {
430                Q = (ModularRingFactory) new ModLongRing(Qi.getVal());
431            } else {
432                Q = (ModularRingFactory) new ModIntegerRing(Qi.getVal());
433            }
434            //Q = new ModIntegerRing(Qi.getVal());
435            //System.out.println("Q = " + Q + ", from Q = " + Mq);
436
437            qfac = new GenPolynomialRing<MOD>(Q, pfac);
438
439            Aq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ai);
440            Bq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Bi);
441            Sp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Si);
442            Tp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ti);
443            //if (debug) {
444            //    E = Ai.multiply(Si).sum(Bi.multiply(Ti));
445            //    Eq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, E);
446            //    if (!Eq.isONE()) {
447            //        System.out.println("Ai*Si+Bi*Ti=1 " + Eq);
448            //        throw new RuntimeException("Ai*Si+Bi*Ti != 1 mod " + Q.getIntegerModul());
449            //    }
450            //}
451        }
452        GreatestCommonDivisorAbstract<BigInteger> ufd = new GreatestCommonDivisorPrimitive<BigInteger>();
453
454        // remove normalization if possible
455        BigInteger ai = ufd.baseContent(Ai);
456        Ai = Ai.divide(ai);
457        BigInteger bi = null;
458        try {
459            bi = c.divide(ai);
460            Bi = Bi.divide(bi); // divide( c/a )
461        } catch (RuntimeException e) {
462            //System.out.println("C  = " + C );
463            //System.out.println("Ai = " + Ai );
464            //System.out.println("Bi = " + Bi );
465            //System.out.println("c  = " + c );
466            //System.out.println("ai = " + ai );
467            //System.out.println("bi = " + bi );
468            //System.out.println("no exact lifting possible " + e);
469            throw new NoLiftingException("no exact lifting possible " + e);
470        }
471        return new HenselApprox<MOD>(Ai, Bi, A1p, B1p);
472    }
473
474
475    /**
476     * Modular quadratic Hensel lifting algorithm on coefficients. Let p =
477     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
478     * with gcd(A,B) == 1 mod p. See algorithm 6.1. in Geddes et.al. and
479     * algorithms 3.5.{5,6} in Cohen. Quadratic version.
480     * @param C GenPolynomial
481     * @param A GenPolynomial
482     * @param B other GenPolynomial
483     * @param M bound on the coefficients of A1 and B1 as factors of C.
484     * @return [A1,B1] = lift(C,A,B), with C = A1 * B1.
485     */
486    @SuppressWarnings("unchecked")
487    public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHenselQuadratic(
488                    GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B)
489                    throws NoLiftingException {
490        if (C == null || C.isZERO()) {
491            return new HenselApprox<MOD>(C, C, A, B);
492        }
493        if (A == null || A.isZERO() || B == null || B.isZERO()) {
494            throw new IllegalArgumentException("A and B must be nonzero");
495        }
496        GenPolynomialRing<BigInteger> fac = C.ring;
497        if (fac.nvar != 1) { // assert ?
498            throw new IllegalArgumentException("polynomial ring not univariate");
499        }
500        // one Hensel step on part polynomials
501        try {
502            GenPolynomial<MOD>[] gst = A.egcd(B);
503            if (!gst[0].isONE()) {
504                throw new NoLiftingException(
505                                "A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B);
506            }
507            GenPolynomial<MOD> s = gst[1];
508            GenPolynomial<MOD> t = gst[2];
509            HenselApprox<MOD> ab = HenselUtil.<MOD> liftHenselQuadratic(C, M, A, B, s, t);
510            return ab;
511        } catch (ArithmeticException e) {
512            throw new NoLiftingException("coefficient error " + e);
513        }
514    }
515
516
517    /**
518     * Modular Hensel lifting algorithm on coefficients. Let p =
519     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
520     * with gcd(A,B) == 1 mod p. See algorithm 6.1. in Geddes et.al. and
521     * algorithms 3.5.{5,6} in Cohen. Quadratic version.
522     * @param C GenPolynomial
523     * @param A GenPolynomial
524     * @param B other GenPolynomial
525     * @param M bound on the coefficients of A1 and B1 as factors of C.
526     * @return [A1,B1] = lift(C,A,B), with C = A1 * B1.
527     */
528    @SuppressWarnings("unchecked")
529    public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHenselQuadraticFac(
530                    GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B)
531                    throws NoLiftingException {
532        if (C == null || C.isZERO()) {
533            throw new IllegalArgumentException("C must be nonzero");
534        }
535        if (A == null || A.isZERO() || B == null || B.isZERO()) {
536            throw new IllegalArgumentException("A and B must be nonzero");
537        }
538        GenPolynomialRing<BigInteger> fac = C.ring;
539        if (fac.nvar != 1) { // assert ?
540            throw new IllegalArgumentException("polynomial ring not univariate");
541        }
542        // one Hensel step on part polynomials
543        try {
544            GenPolynomial<MOD>[] gst = A.egcd(B);
545            if (!gst[0].isONE()) {
546                throw new NoLiftingException(
547                                "A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B);
548            }
549            GenPolynomial<MOD> s = gst[1];
550            GenPolynomial<MOD> t = gst[2];
551            HenselApprox<MOD> ab = HenselUtil.<MOD> liftHenselQuadraticFac(C, M, A, B, s, t);
552            return ab;
553        } catch (ArithmeticException e) {
554            throw new NoLiftingException("coefficient error " + e);
555        }
556    }
557
558
559    /**
560     * Modular Hensel lifting algorithm on coefficients. Let p =
561     * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p
562     * with gcd(A,B) == 1 mod p and S A + T B == 1 mod p. See algorithm 6.1. in
563     * Geddes et.al. and algorithms 3.5.{5,6} in Cohen. Quadratic version, as it
564     * also lifts S A + T B == 1 mod p^{e+1}.
565     * @param C primitive GenPolynomial
566     * @param A GenPolynomial
567     * @param B other GenPolynomial
568     * @param S GenPolynomial
569     * @param T GenPolynomial
570     * @param M bound on the coefficients of A1 and B1 as factors of C.
571     * @return [A1,B1] = lift(C,A,B), with C = A1 * B1.
572     */
573    @SuppressWarnings("unchecked")
574    public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHenselQuadraticFac(
575                    GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B,
576                    GenPolynomial<MOD> S, GenPolynomial<MOD> T) throws NoLiftingException {
577        //System.out.println("*** version for factorization *** ");
578        //GenPolynomial<BigInteger>[] AB = new GenPolynomial[2];
579        if (C == null || C.isZERO()) {
580            throw new IllegalArgumentException("C must be nonzero");
581        }
582        if (A == null || A.isZERO() || B == null || B.isZERO()) {
583            throw new IllegalArgumentException("A and B must be nonzero");
584        }
585        GenPolynomialRing<BigInteger> fac = C.ring;
586        if (fac.nvar != 1) { // assert ?
587            throw new IllegalArgumentException("polynomial ring not univariate");
588        }
589        // setup factories
590        GenPolynomialRing<MOD> pfac = A.ring;
591        RingFactory<MOD> p = pfac.coFac;
592        RingFactory<MOD> q = p;
593        ModularRingFactory<MOD> P = (ModularRingFactory<MOD>) p;
594        ModularRingFactory<MOD> Q = (ModularRingFactory<MOD>) q;
595        BigInteger PP = Q.getIntegerModul();
596        BigInteger Qi = PP;
597        BigInteger M2 = M.multiply(M.fromInteger(2));
598        if (debug) {
599            //System.out.println("M2 =  " + M2);
600            logger.debug("M2 =  {}", M2);
601        }
602        BigInteger Mq = Qi;
603        GenPolynomialRing<MOD> qfac;
604        qfac = new GenPolynomialRing<MOD>(Q, pfac); // mod p
605        GenPolynomialRing<MOD> mfac;
606        BigInteger Mi = Q.getIntegerModul().multiply(Q.getIntegerModul());
607        ModularRingFactory<MOD> Qmm;
608        // = new ModIntegerRing(Mi.getVal());
609        if (ModLongRing.MAX_LONG.compareTo(Mi.getVal()) > 0) {
610            Qmm = (ModularRingFactory) new ModLongRing(Mi.getVal());
611        } else {
612            Qmm = (ModularRingFactory) new ModIntegerRing(Mi.getVal());
613        }
614        mfac = new GenPolynomialRing<MOD>(Qmm, qfac); // mod p^e
615        MOD Qm = Qmm.fromInteger(Qi.getVal());
616
617        // partly normalize c and a, b factors, assert p is prime
618        GenPolynomial<BigInteger> Ai;
619        GenPolynomial<BigInteger> Bi;
620        BigInteger c = C.leadingBaseCoefficient();
621        C = C.multiply(c); // sic
622        MOD a = A.leadingBaseCoefficient();
623        if (!a.isONE()) { // A = A.monic();
624            A = A.divide(a);
625            S = S.multiply(a);
626        }
627        MOD b = B.leadingBaseCoefficient();
628        if (!b.isONE()) { // B = B.monic();
629            B = B.divide(b);
630            T = T.multiply(b);
631        }
632        MOD cm = P.fromInteger(c.getVal());
633        if (cm.isZERO()) {
634            System.out.println("c =  " + c);
635            System.out.println("P =  " + P);
636            throw new ArithmeticException("c mod p == 0 not meaningful");
637        }
638        // mod p
639        A = A.multiply(cm);
640        S = S.divide(cm);
641        B = B.multiply(cm);
642        T = T.divide(cm);
643        Ai = PolyUtil.integerFromModularCoefficients(fac, A);
644        Bi = PolyUtil.integerFromModularCoefficients(fac, B);
645        // replace leading base coefficients
646        ExpVector ea = Ai.leadingExpVector();
647        ExpVector eb = Bi.leadingExpVector();
648        Ai.doPutToMap(ea, c);
649        Bi.doPutToMap(eb, c);
650
651        // polynomials mod p
652        GenPolynomial<MOD> Ap;
653        GenPolynomial<MOD> Bp;
654        GenPolynomial<MOD> A1p = A;
655        GenPolynomial<MOD> B1p = B;
656        GenPolynomial<MOD> Sp = S;
657        GenPolynomial<MOD> Tp = T;
658
659        // polynomials mod q
660        GenPolynomial<MOD> Aq;
661        GenPolynomial<MOD> Bq;
662
663        // polynomials mod p^e
664        GenPolynomial<MOD> Cm;
665        GenPolynomial<MOD> Am;
666        GenPolynomial<MOD> Bm;
667        GenPolynomial<MOD> Em;
668        GenPolynomial<MOD> Emp;
669        GenPolynomial<MOD> Sm;
670        GenPolynomial<MOD> Tm;
671        GenPolynomial<MOD> Ema;
672        GenPolynomial<MOD> Emb;
673        GenPolynomial<MOD> Ema1;
674        GenPolynomial<MOD> Emb1;
675
676        // polynomials over the integers
677        GenPolynomial<BigInteger> Ei;
678        GenPolynomial<BigInteger> Si;
679        GenPolynomial<BigInteger> Ti;
680
681        Si = PolyUtil.integerFromModularCoefficients(fac, S);
682        Ti = PolyUtil.integerFromModularCoefficients(fac, T);
683
684        Aq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ai);
685        Bq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Bi);
686
687        // polynomials mod p^e
688        Cm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, C);
689        Am = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Ai);
690        Bm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Bi);
691        Sm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Si);
692        Tm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Ti);
693        //System.out.println("Cm =  " + Cm);
694        //System.out.println("Am =  " + Am);
695        //System.out.println("Bm =  " + Bm);
696        //System.out.println("Ai =  " + Ai);
697        //System.out.println("Bi =  " + Bi);
698        //System.out.println("mfac =  " + mfac);
699
700        while (Mq.compareTo(M2) < 0) {
701            // compute E=(C-AB)/p mod p^e
702            if (debug) {
703                //System.out.println("mfac =  " + Cm.ring);
704                logger.debug("mfac =  {}", Cm.ring);
705            }
706            Em = Cm.subtract(Am.multiply(Bm));
707            //System.out.println("Em =  " + Em);
708            if (Em.isZERO()) {
709                if (C.subtract(Ai.multiply(Bi)).isZERO()) {
710                    logger.info("leaving on zero E");
711                    break;
712                }
713            }
714            // Em = Em.divide( Qm );
715            Ei = PolyUtil.integerFromModularCoefficients(fac, Em);
716            Ei = Ei.divide(Qi);
717            //System.out.println("Ei =  " + Ei);
718
719            // Ei mod p
720            Emp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ei);
721            //            Emp = PolyUtil.<MOD>fromIntegerCoefficients(qfac,
722            //               PolyUtil.integerFromModularCoefficients(fac,Em) ); 
723            //System.out.println("Emp =  " + Emp);
724            //logger.info("Emp = {}", Emp);
725            if (Emp.isZERO()) {
726                if (C.subtract(Ai.multiply(Bi)).isZERO()) {
727                    logger.info("leaving on zero Emp");
728                    break;
729                }
730            }
731
732            // construct approximation mod p
733            Ap = Sp.multiply(Emp); // S,T ++ T,S 
734            Bp = Tp.multiply(Emp);
735            GenPolynomial<MOD>[] QR = null;
736            QR = Ap.quotientRemainder(Bq); // Bq !
737            GenPolynomial<MOD> Qp = QR[0];
738            GenPolynomial<MOD> Rp = QR[1];
739            A1p = Rp;
740            B1p = Bp.sum(Aq.multiply(Qp)); // Aq !
741            //System.out.println("A1p =  " + A1p);
742            //System.out.println("B1p =  " + B1p);
743
744            // construct q-adic approximation
745            Ema = PolyUtil.<MOD> fromIntegerCoefficients(mfac,
746                            PolyUtil.integerFromModularCoefficients(fac, A1p));
747            Emb = PolyUtil.<MOD> fromIntegerCoefficients(mfac,
748                            PolyUtil.integerFromModularCoefficients(fac, B1p));
749            //System.out.println("Ema =  " + Ema);
750            //System.out.println("Emb =  " + Emb);
751            Ema1 = Ema.multiply(Qm);
752            Emb1 = Emb.multiply(Qm);
753            Ema = Am.sum(Emb1); // Eb1 and Ea1 are required
754            Emb = Bm.sum(Ema1); //--------------------------
755            assert (Ema.degree(0) + Emb.degree(0) <= Cm.degree(0));
756            //if ( Ema.degree(0)+Emb.degree(0) > Cm.degree(0) ) { // debug
757            //   throw new RuntimeException("deg(A)+deg(B) > deg(C)");
758            //}
759            Am = Ema;
760            Bm = Emb;
761            Ai = PolyUtil.integerFromModularCoefficients(fac, Am);
762            Bi = PolyUtil.integerFromModularCoefficients(fac, Bm);
763            //System.out.println("Am =  " + Am);
764            //System.out.println("Bm =  " + Bm);
765            //System.out.println("Ai =  " + Ai);
766            //System.out.println("Bi =  " + Bi + "\n");
767
768            // gcd representation factors error --------------------------------
769            // compute E=(1-SA-TB)/p mod p^e
770            Em = mfac.getONE();
771            Em = Em.subtract(Sm.multiply(Am)).subtract(Tm.multiply(Bm));
772            //System.out.println("Em  =  " + Em);
773            // Em = Em.divide( Pm );
774
775            Ei = PolyUtil.integerFromModularCoefficients(fac, Em);
776            Ei = Ei.divide(Qi);
777            //System.out.println("Ei =  " + Ei);
778            // Ei mod p
779            Emp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ei);
780            //Emp = PolyUtil.<MOD>fromIntegerCoefficients(qfac,
781            //               PolyUtil.integerFromModularCoefficients(fac,Em) );
782            //System.out.println("Emp =  " + Emp);
783
784            // construct approximation mod p
785            Ap = Sp.multiply(Emp); // S,T ++ T,S // Ep Eqp
786            Bp = Tp.multiply(Emp); // Ep Eqp
787            QR = Bp.quotientRemainder(Aq); // Ap Aq ! // Ai == A mod p ?
788            Qp = QR[0];
789            Rp = QR[1];
790            B1p = Rp;
791            A1p = Ap.sum(Bq.multiply(Qp));
792            //System.out.println("A1p =  " + A1p);
793            //System.out.println("B1p =  " + B1p);
794
795            // construct p^e-adic approximation
796            Ema = PolyUtil.<MOD> fromIntegerCoefficients(mfac,
797                            PolyUtil.integerFromModularCoefficients(fac, A1p));
798            Emb = PolyUtil.<MOD> fromIntegerCoefficients(mfac,
799                            PolyUtil.integerFromModularCoefficients(fac, B1p));
800            Ema1 = Ema.multiply(Qm);
801            Emb1 = Emb.multiply(Qm);
802            Ema = Sm.sum(Ema1); // Emb1 and Ema1 are required
803            Emb = Tm.sum(Emb1); //--------------------------
804            Sm = Ema;
805            Tm = Emb;
806            Si = PolyUtil.integerFromModularCoefficients(fac, Sm);
807            Ti = PolyUtil.integerFromModularCoefficients(fac, Tm);
808            //System.out.println("Sm =  " + Sm);
809            //System.out.println("Tm =  " + Tm);
810            //System.out.println("Si =  " + Si);
811            //System.out.println("Ti =  " + Ti + "\n");
812
813            // prepare for next iteration
814            Qi = Q.getIntegerModul().multiply(Q.getIntegerModul());
815            Mq = Qi;
816            //lmfac = mfac;
817            // Q = new ModIntegerRing(Qi.getVal());
818            if (ModLongRing.MAX_LONG.compareTo(Qi.getVal()) > 0) {
819                Q = (ModularRingFactory) new ModLongRing(Qi.getVal());
820            } else {
821                Q = (ModularRingFactory) new ModIntegerRing(Qi.getVal());
822            }
823            qfac = new GenPolynomialRing<MOD>(Q, pfac);
824            BigInteger Qmmi = Qmm.getIntegerModul().multiply(Qmm.getIntegerModul());
825            //Qmm = new ModIntegerRing(Qmmi.getVal());
826            if (ModLongRing.MAX_LONG.compareTo(Qmmi.getVal()) > 0) {
827                Qmm = (ModularRingFactory) new ModLongRing(Qmmi.getVal());
828            } else {
829                Qmm = (ModularRingFactory) new ModIntegerRing(Qmmi.getVal());
830            }
831            mfac = new GenPolynomialRing<MOD>(Qmm, qfac);
832            Qm = Qmm.fromInteger(Qi.getVal());
833
834            Cm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, C);
835            Am = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Ai);
836            Bm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Bi);
837            Sm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Si);
838            Tm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Ti);
839
840            assert (isHenselLift(C, Mi, PP, Ai, Bi));
841            Mi = Mi.fromInteger(Qmm.getIntegerModul().getVal());
842
843            Aq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ai);
844            Bq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Bi);
845            Sp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Si);
846            Tp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ti);
847
848            //System.out.println("Am =  " + Am);
849            //System.out.println("Bm =  " + Bm);
850            //System.out.println("Sm =  " + Sm);
851            //System.out.println("Tm =  " + Tm);
852            //System.out.println("mfac =  " + mfac);
853            //System.out.println("Qmm = " + Qmm + ", M2   =  " + M2 + ", Mq   =  " + Mq + "\n");
854        }
855        //System.out.println("*Ai =  " + Ai);
856        //System.out.println("*Bi =  " + Bi);
857
858        Em = Cm.subtract(Am.multiply(Bm));
859        if (!Em.isZERO()) {
860            System.out.println("Em =  " + Em);
861            //throw new NoLiftingException("no exact lifting possible");
862        }
863        // remove normalization not possible when not exact factorization
864        GreatestCommonDivisorAbstract<BigInteger> ufd = new GreatestCommonDivisorPrimitive<BigInteger>();
865        // remove normalization if possible
866        BigInteger ai = ufd.baseContent(Ai);
867        Ai = Ai.divide(ai); // Ai=pp(Ai)
868        BigInteger[] qr = c.quotientRemainder(ai);
869        BigInteger bi = null;
870        boolean exact = true;
871        if (qr[1].isZERO()) {
872            bi = qr[0];
873            try {
874                Bi = Bi.divide(bi); // divide( c/a )
875            } catch (RuntimeException e) {
876                System.out.println("*catch: no exact factorization: " + bi + ", e = " + e);
877                exact = false;
878            }
879        } else {
880            System.out.println("*remainder: no exact factorization: q = " + qr[0] + ", r = " + qr[1]);
881            exact = false;
882        }
883        if (!exact) {
884            System.out.println("*Ai =  " + Ai);
885            System.out.println("*ai =  " + ai);
886            System.out.println("*Bi =  " + Bi);
887            System.out.println("*bi =  " + bi);
888            System.out.println("*c  =  " + c);
889            throw new NoLiftingException("no exact lifting possible");
890        }
891        return new HenselApprox<MOD>(Ai, Bi, Aq, Bq);
892    }
893
894
895    /**
896     * Modular Hensel lifting test. Let p be a prime number and assume C ==
897     * prod_{0,...,n-1} g_i mod p with gcd(g_i,g_j) == 1 mod p for i != j.
898     * @param C GenPolynomial
899     * @param G = [g_0,...,g_{n-1}] list of GenPolynomial
900     * @param M bound on the coefficients of g_i as factors of C.
901     * @param p prime number.
902     * @return true if C = prod_{0,...,n-1} g_i mod p^e, else false.
903     */
904    public static//<MOD extends GcdRingElem<MOD> & Modular> 
905    boolean isHenselLift(GenPolynomial<BigInteger> C, BigInteger M, BigInteger p,
906                    List<GenPolynomial<BigInteger>> G) {
907        if (C == null || C.isZERO()) {
908            return false;
909        }
910        GenPolynomialRing<BigInteger> pfac = C.ring;
911        ModIntegerRing pm = new ModIntegerRing(p.getVal(), true);
912        GenPolynomialRing<ModInteger> mfac = new GenPolynomialRing<ModInteger>(pm, pfac);
913
914        // check mod p
915        GenPolynomial<ModInteger> cl = mfac.getONE();
916        GenPolynomial<ModInteger> hlp;
917        for (GenPolynomial<BigInteger> hl : G) {
918            //System.out.println("hl       = " + hl);
919            hlp = PolyUtil.<ModInteger> fromIntegerCoefficients(mfac, hl);
920            //System.out.println("hl mod p = " + hlp);
921            cl = cl.multiply(hlp);
922        }
923        GenPolynomial<ModInteger> cp = PolyUtil.<ModInteger> fromIntegerCoefficients(mfac, C);
924        if (!cp.equals(cl)) {
925            System.out.println("Hensel precondition wrong!");
926            if (debug) {
927                System.out.println("cl      = " + cl);
928                System.out.println("cp      = " + cp);
929                System.out.println("mon(cl) = " + cl.monic());
930                System.out.println("mon(cp) = " + cp.monic());
931                System.out.println("cp-cl   = " + cp.subtract(cl));
932                System.out.println("M       = " + M + ", p = " + p);
933            }
934            return false;
935        }
936
937        // check mod p^e 
938        BigInteger mip = p;
939        while (mip.compareTo(M) < 0) {
940            mip = mip.multiply(mip); // p
941        }
942        // mip = mip.multiply(p);
943        pm = new ModIntegerRing(mip.getVal(), false);
944        mfac = new GenPolynomialRing<ModInteger>(pm, pfac);
945        cl = mfac.getONE();
946        for (GenPolynomial<BigInteger> hl : G) {
947            //System.out.println("hl         = " + hl);
948            hlp = PolyUtil.<ModInteger> fromIntegerCoefficients(mfac, hl);
949            //System.out.println("hl mod p^e = " + hlp);
950            cl = cl.multiply(hlp);
951        }
952        cp = PolyUtil.<ModInteger> fromIntegerCoefficients(mfac, C);
953        if (!cp.equals(cl)) {
954            System.out.println("Hensel post condition wrong!");
955            System.out.println("cl    = " + cl);
956            System.out.println("cp    = " + cp);
957            System.out.println("cp-cl = " + cp.subtract(cl));
958            System.out.println("M = " + M + ", p = " + p + ", p^e = " + mip);
959            return false;
960        }
961        return true;
962    }
963
964
965    /**
966     * Modular Hensel lifting test. Let p be a prime number and assume C == A *
967     * B mod p with gcd(A,B) == 1 mod p.
968     * @param C GenPolynomial
969     * @param A GenPolynomial
970     * @param B GenPolynomial
971     * @param M bound on the coefficients of A and B as factors of C.
972     * @param p prime number.
973     * @return true if C = A * B mod p**e, else false.
974     */
975    public static//<MOD extends GcdRingElem<MOD> & Modular>
976    boolean isHenselLift(GenPolynomial<BigInteger> C, BigInteger M, BigInteger p, GenPolynomial<BigInteger> A,
977                    GenPolynomial<BigInteger> B) {
978        List<GenPolynomial<BigInteger>> G = new ArrayList<GenPolynomial<BigInteger>>(2);
979        G.add(A);
980        G.add(B);
981        return isHenselLift(C, M, p, G);
982    }
983
984
985    /**
986     * Modular Hensel lifting test. Let p be a prime number and assume C == A *
987     * B mod p with gcd(A,B) == 1 mod p.
988     * @param C GenPolynomial
989     * @param Ha Hensel approximation.
990     * @param M bound on the coefficients of A and B as factors of C.
991     * @param p prime number.
992     * @return true if C = A * B mod p^e, else false.
993     */
994    public static <MOD extends GcdRingElem<MOD> & Modular> boolean isHenselLift(GenPolynomial<BigInteger> C,
995                    BigInteger M, BigInteger p, HenselApprox<MOD> Ha) {
996        List<GenPolynomial<BigInteger>> G = new ArrayList<GenPolynomial<BigInteger>>(2);
997        G.add(Ha.A);
998        G.add(Ha.B);
999        return isHenselLift(C, M, p, G);
1000    }
1001
1002
1003    /**
1004     * Constructing and lifting algorithm for extended Euclidean relation. Let p
1005     * = A.ring.coFac.modul() and assume gcd(A,B) == 1 mod p.
1006     * @param A modular GenPolynomial
1007     * @param B modular GenPolynomial
1008     * @param k desired approximation exponent p^k.
1009     * @return [s,t] with s A + t B = 1 mod p^k.
1010     */
1011    @SuppressWarnings("unchecked")
1012    public static <MOD extends GcdRingElem<MOD> & Modular> GenPolynomial<MOD>[] liftExtendedEuclidean(
1013                    GenPolynomial<MOD> A, GenPolynomial<MOD> B, long k) throws NoLiftingException {
1014        if (A == null || A.isZERO() || B == null || B.isZERO()) {
1015            throw new IllegalArgumentException("A and B must be nonzero, A = " + A + ", B = " + B);
1016        }
1017        GenPolynomialRing<MOD> fac = A.ring;
1018        if (fac.nvar != 1) { // assert ?
1019            throw new IllegalArgumentException("polynomial ring not univariate");
1020        }
1021        // start with extended Euclidean relation mod p
1022        GenPolynomial<MOD>[] gst = null;
1023        try {
1024            gst = A.egcd(B);
1025            if (!gst[0].isONE()) {
1026                throw new NoLiftingException(
1027                                "A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B);
1028            }
1029        } catch (ArithmeticException e) {
1030            throw new NoLiftingException("coefficient error " + e);
1031        }
1032        GenPolynomial<MOD> S = gst[1];
1033        GenPolynomial<MOD> T = gst[2];
1034        //System.out.println("eeS = " + S + ": " + S.ring.coFac);
1035        //System.out.println("eeT = " + T + ": " + T.ring.coFac);
1036
1037        // setup integer polynomial ring
1038        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1039        GenPolynomial<BigInteger> one = ifac.getONE();
1040        GenPolynomial<BigInteger> Ai = PolyUtil.integerFromModularCoefficients(ifac, A);
1041        GenPolynomial<BigInteger> Bi = PolyUtil.integerFromModularCoefficients(ifac, B);
1042        GenPolynomial<BigInteger> Si = PolyUtil.integerFromModularCoefficients(ifac, S);
1043        GenPolynomial<BigInteger> Ti = PolyUtil.integerFromModularCoefficients(ifac, T);
1044        //System.out.println("Ai = " + Ai);
1045        //System.out.println("Bi = " + Bi);
1046        //System.out.println("Si = " + Si);
1047        //System.out.println("Ti = " + Ti);
1048
1049        // approximate mod p^i
1050        ModularRingFactory<MOD> mcfac = (ModularRingFactory<MOD>) fac.coFac;
1051        BigInteger p = mcfac.getIntegerModul();
1052        BigInteger modul = p;
1053        GenPolynomialRing<MOD> mfac; // = new GenPolynomialRing<MOD>(mcfac, fac);
1054        for (int i = 1; i < k; i++) {
1055            // e = 1 - s a - t b in Z[x]
1056            GenPolynomial<BigInteger> e = one.subtract(Si.multiply(Ai)).subtract(Ti.multiply(Bi));
1057            //System.out.println("\ne = " + e);
1058            if (e.isZERO()) {
1059                logger.info("leaving on zero e in liftExtendedEuclidean");
1060                break;
1061            }
1062            e = e.divide(modul);
1063            // move to Z_p[x] and compute next approximation 
1064            GenPolynomial<MOD> c = PolyUtil.<MOD> fromIntegerCoefficients(fac, e);
1065            //System.out.println("c = " + c + ": " + c.ring.coFac);
1066            GenPolynomial<MOD> s = S.multiply(c);
1067            GenPolynomial<MOD> t = T.multiply(c);
1068            //System.out.println("s = " + s + ": " + s.ring.coFac);
1069            //System.out.println("t = " + t + ": " + t.ring.coFac);
1070
1071            GenPolynomial<MOD>[] QR = s.quotientRemainder(B); // watch for ordering 
1072            GenPolynomial<MOD> q = QR[0];
1073            s = QR[1];
1074            t = t.sum(q.multiply(A));
1075            //System.out.println("s = " + s + ": " + s.ring.coFac);
1076            //System.out.println("t = " + t + ": " + t.ring.coFac);
1077
1078            GenPolynomial<BigInteger> si = PolyUtil.integerFromModularCoefficients(ifac, s);
1079            GenPolynomial<BigInteger> ti = PolyUtil.integerFromModularCoefficients(ifac, t);
1080            //System.out.println("si = " + si);
1081            //System.out.println("ti = " + si);
1082            // add approximation to solution
1083            Si = Si.sum(si.multiply(modul));
1084            Ti = Ti.sum(ti.multiply(modul));
1085            //System.out.println("Si = " + Si);
1086            //System.out.println("Ti = " + Ti);
1087            modul = modul.multiply(p);
1088            //System.out.println("modul = " + modul + ", " + p + "^" + i + ", p^i = " + p.power(i));
1089        }
1090        //System.out.println("Si = " + Si + ", Ti = " + Ti);
1091        // setup ring mod p^i
1092        if (ModLongRing.MAX_LONG.compareTo(modul.getVal()) > 0) {
1093            mcfac = (ModularRingFactory) new ModLongRing(modul.getVal());
1094        } else {
1095            mcfac = (ModularRingFactory) new ModIntegerRing(modul.getVal());
1096        }
1097        //System.out.println("mcfac = " + mcfac);
1098        mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1099        S = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Si);
1100        T = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Ti);
1101        //System.out.println("S = " + S + ": " + S.ring.coFac);
1102        //System.out.println("T = " + T + ": " + T.ring.coFac);
1103        if (debug) {
1104            List<GenPolynomial<MOD>> AP = new ArrayList<GenPolynomial<MOD>>();
1105            AP.add(B);
1106            AP.add(A);
1107            List<GenPolynomial<MOD>> SP = new ArrayList<GenPolynomial<MOD>>();
1108            SP.add(S);
1109            SP.add(T);
1110            if (!HenselUtil.<MOD> isExtendedEuclideanLift(AP, SP)) {
1111                System.out.println("isExtendedEuclideanLift: false");
1112            }
1113        }
1114        @SuppressWarnings("cast")
1115        GenPolynomial<MOD>[] rel = (GenPolynomial<MOD>[]) new GenPolynomial[2];
1116        rel[0] = S;
1117        rel[1] = T;
1118        return rel;
1119    }
1120
1121
1122    /**
1123     * Constructing and lifting algorithm for extended Euclidean relation. Let p
1124     * = A_i.ring.coFac.modul() and assume gcd(A_i,A_j) == 1 mod p, i != j.
1125     * @param A list of modular GenPolynomials
1126     * @param k desired approximation exponent p^k.
1127     * @return [s_0,...,s_n-1] with sum_i s_i * B_i = 1 mod p^k, with B_i =
1128     *         prod_{i!=j} A_j.
1129     */
1130    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftExtendedEuclidean(
1131                    List<GenPolynomial<MOD>> A, long k) throws NoLiftingException {
1132        if (A == null || A.size() <= 1) {
1133            throw new IllegalArgumentException("A must be non null and non empty");
1134        }
1135        GenPolynomialRing<MOD> fac = A.get(0).ring;
1136        if (fac.nvar != 1) { // assert ?
1137            throw new IllegalArgumentException("polynomial ring not univariate");
1138        }
1139        GenPolynomial<MOD> zero = fac.getZERO();
1140        int r = A.size();
1141        List<GenPolynomial<MOD>> Q = new ArrayList<GenPolynomial<MOD>>(r);
1142        for (int i = 0; i < r; i++) {
1143            Q.add(zero);
1144        }
1145        //System.out.println("A = " + A);
1146        Q.set(r - 2, A.get(r - 1));
1147        for (int j = r - 3; j >= 0; j--) {
1148            GenPolynomial<MOD> q = A.get(j + 1).multiply(Q.get(j + 1));
1149            Q.set(j, q);
1150        }
1151        //System.out.println("Q = " + Q);
1152        List<GenPolynomial<MOD>> B = new ArrayList<GenPolynomial<MOD>>(r + 1);
1153        List<GenPolynomial<MOD>> lift = new ArrayList<GenPolynomial<MOD>>(r);
1154        for (int i = 0; i < r; i++) {
1155            B.add(zero);
1156            lift.add(zero);
1157        }
1158        GenPolynomial<MOD> one = fac.getONE();
1159        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1160        B.add(0, one);
1161        //System.out.println("B(0) = " + B.get(0));
1162        GenPolynomial<MOD> b = one;
1163        for (int j = 0; j < r - 1; j++) {
1164            //System.out.println("Q("+(j)+") = " + Q.get(j));
1165            //System.out.println("A("+(j)+") = " + A.get(j));
1166            //System.out.println("B("+(j)+") = " + B.get(j));
1167            List<GenPolynomial<MOD>> S = liftDiophant(Q.get(j), A.get(j), B.get(j), k);
1168            //System.out.println("\nSb = " + S);
1169            b = S.get(0);
1170            GenPolynomial<MOD> bb = PolyUtil.<MOD> fromIntegerCoefficients(fac,
1171                            PolyUtil.integerFromModularCoefficients(ifac, b));
1172            B.set(j + 1, bb);
1173            lift.set(j, S.get(1));
1174            //System.out.println("B("+(j+1)+") = " + B.get(j+1));
1175            if (debug) {
1176                logger.info("lift({}) = {}", j, lift.get(j));
1177            }
1178        }
1179        //System.out.println("liftb = " + lift);
1180        lift.set(r - 1, b);
1181        if (debug) {
1182            logger.info("lift({}) = {}", (r - 1), b);
1183        }
1184        //System.out.println("B("+(r-1)+") = " + B.get(r-1) + " : " +  B.get(r-1).ring.coFac + ", b = " +  b + " : " +  b.ring.coFac);
1185        //System.out.println("B = " + B);
1186        //System.out.println("liftb = " + lift);
1187        return lift;
1188    }
1189
1190
1191    /**
1192     * Modular diophantine equation solution and lifting algorithm. Let p =
1193     * A_i.ring.coFac.modul() and assume gcd(A,B) == 1 mod p.
1194     * @param A modular GenPolynomial, mod p^k
1195     * @param B modular GenPolynomial, mod p^k
1196     * @param C modular GenPolynomial, mod p^k
1197     * @param k desired approximation exponent p^k.
1198     * @return [s, t] with s A' + t B' = C mod p^k, with A' = B, B' = A.
1199     */
1200    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant(
1201                    GenPolynomial<MOD> A, GenPolynomial<MOD> B, GenPolynomial<MOD> C, long k)
1202                    throws NoLiftingException {
1203        if (A == null || A.isZERO() || B == null || B.isZERO()) {
1204            throw new IllegalArgumentException(
1205                            "A and B must be nonzero, A = " + A + ", B = " + B + ", C = " + C);
1206        }
1207        List<GenPolynomial<MOD>> sol = new ArrayList<GenPolynomial<MOD>>();
1208        GenPolynomialRing<MOD> fac = C.ring;
1209        if (fac.nvar != 1) { // assert ?
1210            throw new IllegalArgumentException("polynomial ring not univariate");
1211        }
1212        //System.out.println("C = " + C);
1213        GenPolynomial<MOD> zero = fac.getZERO();
1214        for (int i = 0; i < 2; i++) {
1215            sol.add(zero);
1216        }
1217        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1218        for (Monomial<MOD> m : C) {
1219            //System.out.println("monomial = " + m);
1220            long e = m.e.getVal(0);
1221            List<GenPolynomial<MOD>> S = liftDiophant(A, B, e, k);
1222            //System.out.println("Se = " + S);
1223            MOD a = m.c;
1224            //System.out.println("C.fac = " + fac.toScript());
1225            a = fac.coFac.fromInteger(a.getSymmetricInteger().getVal());
1226            int i = 0;
1227            for (GenPolynomial<MOD> d : S) {
1228                //System.out.println("d = " + d);
1229                d = PolyUtil.<MOD> fromIntegerCoefficients(fac,
1230                                PolyUtil.integerFromModularCoefficients(ifac, d));
1231                d = d.multiply(a);
1232                d = sol.get(i).sum(d);
1233                //System.out.println("d = " + d);
1234                sol.set(i++, d);
1235            }
1236            //System.out.println("sol = " + sol + ", for " + m);
1237        }
1238        if (debug) {
1239            //GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1240            A = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, A));
1241            B = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, B));
1242            C = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, C));
1243            GenPolynomial<MOD> y = B.multiply(sol.get(0)).sum(A.multiply(sol.get(1)));
1244            if (!y.equals(C)) {
1245                System.out.println("A = " + A + ", B = " + B);
1246                System.out.println("s1 = " + sol.get(0) + ", s2 = " + sol.get(1));
1247                System.out.println("Error: A*r1 + B*r2 = " + y + " : " + fac.coFac);
1248            }
1249        }
1250        return sol;
1251    }
1252
1253
1254    /**
1255     * Modular diophantine equation solution and lifting algorithm. Let p =
1256     * A_i.ring.coFac.modul() and assume gcd(a,b) == 1 mod p, for a, b in A.
1257     * @param A list of modular GenPolynomials, mod p^k
1258     * @param C modular GenPolynomial, mod p^k
1259     * @param k desired approximation exponent p^k.
1260     * @return [s_1,..., s_n] with sum_i s_i A_i' = C mod p^k, with Ai' =
1261     *         prod_{j!=i} A_j.
1262     */
1263    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant(
1264                    List<GenPolynomial<MOD>> A, GenPolynomial<MOD> C, long k) throws NoLiftingException {
1265        if (false && A.size() <= 2) {
1266            return HenselUtil.<MOD> liftDiophant(A.get(0), A.get(1), C, k);
1267        }
1268        List<GenPolynomial<MOD>> sol = new ArrayList<GenPolynomial<MOD>>();
1269        GenPolynomialRing<MOD> fac = C.ring;
1270        if (fac.nvar != 1) { // assert ?
1271            throw new IllegalArgumentException("polynomial ring not univariate");
1272        }
1273        //System.out.println("C = " + C);
1274        GenPolynomial<MOD> zero = fac.getZERO();
1275        for (int i = 0; i < A.size(); i++) {
1276            sol.add(zero);
1277        }
1278        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1279        for (Monomial<MOD> m : C) {
1280            //System.out.println("monomial = " + m);
1281            long e = m.e.getVal(0);
1282            List<GenPolynomial<MOD>> S = liftDiophant(A, e, k);
1283            //System.out.println("Se = " + S);
1284            MOD a = m.c;
1285            //System.out.println("C.fac = " + fac.toScript());
1286            a = fac.coFac.fromInteger(a.getSymmetricInteger().getVal());
1287            int i = 0;
1288            for (GenPolynomial<MOD> d : S) {
1289                //System.out.println("d = " + d);
1290                d = PolyUtil.<MOD> fromIntegerCoefficients(fac,
1291                                PolyUtil.integerFromModularCoefficients(ifac, d));
1292                d = d.multiply(a);
1293                d = sol.get(i).sum(d);
1294                //System.out.println("d = " + d);
1295                sol.set(i++, d);
1296            }
1297            //System.out.println("sol = " + sol + ", for " + m);
1298        }
1299        /*
1300        if (true || debug) {
1301            //GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1302            A = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, A));
1303            B = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, B));
1304            C = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, C));
1305            GenPolynomial<MOD> y = B.multiply(sol.get(0)).sum(A.multiply(sol.get(1)));
1306            if (!y.equals(C)) {
1307                System.out.println("A = " + A + ", B = " + B);
1308                System.out.println("s1 = " + sol.get(0) + ", s2 = " + sol.get(1));
1309                System.out.println("Error: A*r1 + B*r2 = " + y + " : " + fac.coFac);
1310            }
1311        }
1312        */
1313        return sol;
1314    }
1315
1316
1317    /**
1318     * Modular diophantine equation solution and lifting algorithm. Let p =
1319     * A_i.ring.coFac.modul() and assume gcd(A,B) == 1 mod p.
1320     * @param A modular GenPolynomial
1321     * @param B modular GenPolynomial
1322     * @param e exponent for x^e
1323     * @param k desired approximation exponent p^k.
1324     * @return [s, t] with s A' + t B' = x^e mod p^k, with A' = B, B' = A.
1325     */
1326    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant(
1327                    GenPolynomial<MOD> A, GenPolynomial<MOD> B, long e, long k) throws NoLiftingException {
1328        if (A == null || A.isZERO() || B == null || B.isZERO()) {
1329            throw new IllegalArgumentException("A and B must be nonzero, A = " + A + ", B = " + B);
1330        }
1331        List<GenPolynomial<MOD>> sol = new ArrayList<GenPolynomial<MOD>>();
1332        GenPolynomialRing<MOD> fac = A.ring;
1333        if (fac.nvar != 1) { // assert ?
1334            throw new IllegalArgumentException("polynomial ring not univariate");
1335        }
1336        // lift EE relation to p^k
1337        GenPolynomial<MOD>[] lee = liftExtendedEuclidean(B, A, k);
1338        GenPolynomial<MOD> s1 = lee[0];
1339        GenPolynomial<MOD> s2 = lee[1];
1340        if (e == 0L) {
1341            sol.add(s1);
1342            sol.add(s2);
1343            //System.out.println("sol@0 = " + sol); 
1344            return sol;
1345        }
1346        fac = s1.ring;
1347        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1348        A = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, A));
1349        B = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, B));
1350
1351        //      this is the wrong sequence:
1352        //         GenPolynomial<MOD> xe = fac.univariate(0,e);
1353        //         GenPolynomial<MOD> q = s1.multiply(xe);
1354        //         GenPolynomial<MOD>[] QR = q.quotientRemainder(B);
1355        //         q = QR[0];
1356        //         GenPolynomial<MOD> r1 = QR[1];
1357        //         GenPolynomial<MOD> r2 = s2.multiply(xe).sum( q.multiply(A) );
1358
1359        GenPolynomial<MOD> xe = fac.univariate(0, e);
1360        GenPolynomial<MOD> q = s1.multiply(xe);
1361        GenPolynomial<MOD>[] QR = q.quotientRemainder(A);
1362        q = QR[0];
1363        //System.out.println("ee coeff qr = " + Arrays.toString(QR));
1364        GenPolynomial<MOD> r1 = QR[1];
1365        GenPolynomial<MOD> r2 = s2.multiply(xe).sum(q.multiply(B));
1366        //System.out.println("r1 = " + r1 + ", r2 = " + r2);
1367        sol.add(r1);
1368        sol.add(r2);
1369        //System.out.println("sol@"+ e + " = " + sol); 
1370        if (debug) {
1371            GenPolynomial<MOD> y = B.multiply(r1).sum(A.multiply(r2));
1372            if (!y.equals(xe)) {
1373                System.out.println("A = " + A + ", B = " + B);
1374                System.out.println("r1 = " + r1 + ", r2 = " + r2);
1375                System.out.println("Error: A*r1 + B*r2 = " + y);
1376            }
1377        }
1378        return sol;
1379    }
1380
1381
1382    /**
1383     * Modular diophantine equation solution and lifting algorithm. Let p =
1384     * A_i.ring.coFac.modul() and assume gcd(a,b) == 1 mod p, for a, b in A.
1385     * @param A list of modular GenPolynomials
1386     * @param e exponent for x^e
1387     * @param k desired approximation exponent p^k.
1388     * @return [s_1,..., s_n] with sum_i s_i A_i' = x^e mod p^k, with Ai' =
1389     *         prod_{j!=i} A_j.
1390     */
1391    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant(
1392                    List<GenPolynomial<MOD>> A, long e, long k) throws NoLiftingException {
1393        if (false && A.size() <= 2) {
1394            return HenselUtil.<MOD> liftDiophant(A.get(0), A.get(1), e, k);
1395        }
1396        List<GenPolynomial<MOD>> sol = new ArrayList<GenPolynomial<MOD>>();
1397        GenPolynomialRing<MOD> fac = A.get(0).ring;
1398        if (fac.nvar != 1) { // assert ?
1399            throw new IllegalArgumentException("polynomial ring not univariate");
1400        }
1401        // lift EE relation to p^k
1402        List<GenPolynomial<MOD>> lee = liftExtendedEuclidean(A, k);
1403        if (e == 0L) {
1404            //System.out.println("sol@0 = " + sol); 
1405            return lee;
1406        }
1407        fac = lee.get(0).ring;
1408        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1409        List<GenPolynomial<MOD>> S = new ArrayList<GenPolynomial<MOD>>(lee.size());
1410        for (GenPolynomial<MOD> a : lee) {
1411            a = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, a));
1412            S.add(a);
1413        }
1414        GenPolynomial<MOD> xe = fac.univariate(0, e);
1415        //List<GenPolynomial<MOD>> Sr = new ArrayList<GenPolynomial<MOD>>(lee.size());
1416        int i = 0;
1417        for (GenPolynomial<MOD> s : S) {
1418            GenPolynomial<MOD> q = s.multiply(xe);
1419            GenPolynomial<MOD> r = q.remainder(A.get(i++));
1420            //System.out.println("r = " + r);
1421            sol.add(r);
1422        }
1423        //System.out.println("sol@"+ e + " = " + sol); 
1424        /*
1425        if (true || debug) {
1426            GenPolynomial<MOD> y = B.multiply(r1).sum(A.multiply(r2));
1427            if (!y.equals(xe)) {
1428                System.out.println("A = " + A + ", B = " + B);
1429                System.out.println("r1 = " + r1 + ", r2 = " + r2);
1430                System.out.println("Error: A*r1 + B*r2 = " + y);
1431            }
1432        }
1433        */
1434        return sol;
1435    }
1436
1437
1438    /**
1439     * Modular Diophant relation lifting test.
1440     * @param A modular GenPolynomial
1441     * @param B modular GenPolynomial
1442     * @param C modular GenPolynomial
1443     * @param S1 modular GenPolynomial
1444     * @param S2 modular GenPolynomial
1445     * @return true if A*S1 + B*S2 = C, else false.
1446     */
1447    public static <MOD extends GcdRingElem<MOD> & Modular> boolean isDiophantLift(GenPolynomial<MOD> A,
1448                    GenPolynomial<MOD> B, GenPolynomial<MOD> S1, GenPolynomial<MOD> S2,
1449                    GenPolynomial<MOD> C) {
1450        GenPolynomialRing<MOD> fac = C.ring;
1451        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1452        GenPolynomial<MOD> a = PolyUtil.<MOD> fromIntegerCoefficients(fac,
1453                        PolyUtil.integerFromModularCoefficients(ifac, A));
1454        GenPolynomial<MOD> b = PolyUtil.<MOD> fromIntegerCoefficients(fac,
1455                        PolyUtil.integerFromModularCoefficients(ifac, B));
1456        GenPolynomial<MOD> s1 = PolyUtil.<MOD> fromIntegerCoefficients(fac,
1457                        PolyUtil.integerFromModularCoefficients(ifac, S1));
1458        GenPolynomial<MOD> s2 = PolyUtil.<MOD> fromIntegerCoefficients(fac,
1459                        PolyUtil.integerFromModularCoefficients(ifac, S2));
1460        GenPolynomial<MOD> t = a.multiply(s1).sum(b.multiply(s2));
1461        if (t.equals(C)) {
1462            return true;
1463        }
1464        if (debug) {
1465            System.out.println("a  = " + a);
1466            System.out.println("b  = " + b);
1467            System.out.println("s1 = " + s1);
1468            System.out.println("s2 = " + s2);
1469            System.out.println("t  = " + t);
1470            System.out.println("C  = " + C);
1471        }
1472        return false;
1473    }
1474
1475
1476    /**
1477     * Modular extended Euclidean relation lifting test.
1478     * @param A list of GenPolynomials
1479     * @param S = [s_0,...,s_{n-1}] list of GenPolynomial
1480     * @return true if prod_{0,...,n-1} s_i * B_i = 1 mod p^e, with B_i =
1481     *         prod_{i!=j} A_j, else false.
1482     */
1483    public static <MOD extends GcdRingElem<MOD> & Modular> boolean isExtendedEuclideanLift(
1484                    List<GenPolynomial<MOD>> A, List<GenPolynomial<MOD>> S) {
1485        GenPolynomialRing<MOD> fac = A.get(0).ring;
1486        GenPolynomial<MOD> C = fac.getONE();
1487        return isDiophantLift(A, S, C);
1488    }
1489
1490
1491    /**
1492     * Modular Diophant relation lifting test.
1493     * @param A list of GenPolynomials
1494     * @param S = [s_0,...,s_{n-1}] list of GenPolynomials
1495     * @param C = GenPolynomial
1496     * @return true if prod_{0,...,n-1} s_i * B_i = C mod p^k, with B_i =
1497     *         prod_{i!=j} A_j, else false.
1498     */
1499    public static <MOD extends GcdRingElem<MOD> & Modular> boolean isDiophantLift(List<GenPolynomial<MOD>> A,
1500                    List<GenPolynomial<MOD>> S, GenPolynomial<MOD> C) {
1501        GenPolynomialRing<MOD> fac = A.get(0).ring;
1502        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1503        List<GenPolynomial<MOD>> B = new ArrayList<GenPolynomial<MOD>>(A.size());
1504        int i = 0;
1505        for (GenPolynomial<MOD> ai : A) {
1506            GenPolynomial<MOD> b = fac.getONE();
1507            int j = 0;
1508            for (GenPolynomial<MOD> aj : A) {
1509                if (i != j /*!ai.equals(aj)*/) {
1510                    b = b.multiply(aj);
1511                }
1512                j++;
1513            }
1514            //System.out.println("b = " + b);
1515            b = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, b));
1516            B.add(b);
1517            i++;
1518        }
1519        //System.out.println("B = " + B);
1520        // check mod p^e 
1521        GenPolynomial<MOD> t = fac.getZERO();
1522        i = 0;
1523        for (GenPolynomial<MOD> a : B) {
1524            GenPolynomial<MOD> b = S.get(i++);
1525            b = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, b));
1526            GenPolynomial<MOD> s = a.multiply(b);
1527            t = t.sum(s);
1528        }
1529        if (!t.equals(C)) {
1530            if (debug) {
1531                System.out.println("no diophant lift!");
1532                System.out.println("A = " + A);
1533                System.out.println("B = " + B);
1534                System.out.println("S = " + S);
1535                System.out.println("C = " + C);
1536                System.out.println("t = " + t);
1537            }
1538            return false;
1539        }
1540        return true;
1541    }
1542
1543
1544    /**
1545     * Modular Hensel lifting algorithm on coefficients. Let p =
1546     * f_i.ring.coFac.modul() and assume C == prod_{0,...,n-1} f_i mod p with
1547     * gcd(f_i,f_j) == 1 mod p for i != j
1548     * @param C monic integer polynomial
1549     * @param F = [f_0,...,f_{n-1}] list of monic modular polynomials.
1550     * @param k approximation exponent.
1551     * @return [g_0,...,g_{n-1}] with C = prod_{0,...,n-1} g_i mod p^k.
1552     */
1553    @SuppressWarnings("unchecked")
1554    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHenselMonic(
1555                    GenPolynomial<BigInteger> C, List<GenPolynomial<MOD>> F, long k)
1556                    throws NoLiftingException {
1557        if (C == null || C.isZERO() || F == null || F.size() == 0) {
1558            throw new IllegalArgumentException("C must be nonzero and F must be nonempty");
1559        }
1560        GenPolynomialRing<BigInteger> fac = C.ring;
1561        if (fac.nvar != 1) { // assert ?
1562            throw new IllegalArgumentException("polynomial ring not univariate");
1563        }
1564        List<GenPolynomial<MOD>> lift = new ArrayList<GenPolynomial<MOD>>(F.size());
1565        GenPolynomialRing<MOD> pfac = F.get(0).ring;
1566        RingFactory<MOD> pcfac = pfac.coFac;
1567        ModularRingFactory<MOD> PF = (ModularRingFactory<MOD>) pcfac;
1568        BigInteger P = PF.getIntegerModul();
1569        int n = F.size();
1570        if (n == 1) { // lift F_0, this case will probably never be used
1571            GenPolynomial<MOD> f = F.get(0);
1572            ModularRingFactory<MOD> mcfac;
1573            if (ModLongRing.MAX_LONG.compareTo(P.getVal()) > 0) {
1574                mcfac = (ModularRingFactory) new ModLongRing(P.getVal());
1575            } else {
1576                mcfac = (ModularRingFactory) new ModIntegerRing(P.getVal());
1577            }
1578            GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1579            f = PolyUtil.fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac, f));
1580            lift.add(f);
1581            return lift;
1582        }
1583        //         if (n == 2) { // only one step
1584        //             HenselApprox<MOD> ab = HenselUtil.<MOD> liftHenselQuadratic(C, M, F.get(0), F.get(1));
1585        //             lift.add(ab.Am);
1586        //             lift.add(ab.Bm);
1587        //             return lift;
1588        //         }
1589
1590        // setup integer polynomial ring
1591        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1592        List<GenPolynomial<BigInteger>> Fi = PolyUtil.integerFromModularCoefficients(ifac, F);
1593        //System.out.println("Fi = " + Fi);
1594
1595        List<GenPolynomial<MOD>> S = liftExtendedEuclidean(F, k + 1); // lift works for any k, use this
1596        //System.out.println("Sext = " + S);
1597        if (debug) {
1598            logger.info("EE lift = {}", S);
1599            // adjust coefficients
1600            List<GenPolynomial<MOD>> Sx = PolyUtil.fromIntegerCoefficients(pfac,
1601                            PolyUtil.integerFromModularCoefficients(ifac, S));
1602            try {
1603                boolean il = HenselUtil.<MOD> isExtendedEuclideanLift(F, Sx);
1604                //System.out.println("islift = " + il);
1605            } catch (RuntimeException e) {
1606                e.printStackTrace();
1607            }
1608        }
1609        List<GenPolynomial<BigInteger>> Si = PolyUtil.integerFromModularCoefficients(ifac, S);
1610        //System.out.println("Si = " + Si);
1611        //System.out.println("C = " + C);
1612
1613        // approximate mod p^i
1614        ModularRingFactory<MOD> mcfac = PF;
1615        BigInteger p = mcfac.getIntegerModul();
1616        BigInteger modul = p;
1617        GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1618        List<GenPolynomial<MOD>> Sp = PolyUtil.fromIntegerCoefficients(mfac, Si);
1619        //System.out.println("Sp = " + Sp);
1620        for (int i = 1; i < k; i++) {
1621            //System.out.println("i = " + i);
1622            GenPolynomial<BigInteger> e = fac.getONE();
1623            for (GenPolynomial<BigInteger> fi : Fi) {
1624                e = e.multiply(fi);
1625            }
1626            e = C.subtract(e);
1627            //System.out.println("\ne = " + e);
1628            if (e.isZERO()) {
1629                logger.info("leaving on zero e");
1630                break;
1631            }
1632            try {
1633                e = e.divide(modul);
1634            } catch (RuntimeException ex) {
1635                ex.printStackTrace();
1636                throw ex;
1637            }
1638            //System.out.println("e = " + e);
1639            // move to in Z_p[x]
1640            GenPolynomial<MOD> c = PolyUtil.<MOD> fromIntegerCoefficients(mfac, e);
1641            //System.out.println("c = " + c + ": " + c.ring.coFac);
1642
1643            List<GenPolynomial<MOD>> s = new ArrayList<GenPolynomial<MOD>>(S.size());
1644            int j = 0;
1645            for (GenPolynomial<MOD> f : Sp) {
1646                f = f.multiply(c);
1647                //System.out.println("f = " + f + " : " + f.ring.coFac);
1648                //System.out.println("F,i = " + F.get(j) + " : " + F.get(j).ring.coFac);
1649                f = f.remainder(F.get(j++));
1650                //System.out.println("f = " + f + " : " + f.ring.coFac);
1651                s.add(f);
1652            }
1653            //System.out.println("s = " + s);
1654            List<GenPolynomial<BigInteger>> si = PolyUtil.integerFromModularCoefficients(ifac, s);
1655            //System.out.println("si = " + si);
1656
1657            List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(F.size());
1658            j = 0;
1659            for (GenPolynomial<BigInteger> f : Fi) {
1660                f = f.sum(si.get(j++).multiply(modul));
1661                Fii.add(f);
1662            }
1663            //System.out.println("Fii = " + Fii);
1664            Fi = Fii;
1665            modul = modul.multiply(p);
1666            if (i >= k - 1) {
1667                logger.info("e != 0 for k = {}", k);
1668            }
1669        }
1670        // setup ring mod p^k
1671        modul = p.power(k);
1672        if (ModLongRing.MAX_LONG.compareTo(modul.getVal()) > 0) {
1673            mcfac = (ModularRingFactory) new ModLongRing(modul.getVal());
1674        } else {
1675            mcfac = (ModularRingFactory) new ModIntegerRing(modul.getVal());
1676        }
1677        //System.out.println("mcfac = " + mcfac);
1678        mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1679        lift = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Fi);
1680        //System.out.println("lift = " + lift + ": " + lift.get(0).ring.coFac);
1681        return lift;
1682    }
1683
1684
1685    /**
1686     * Modular Hensel lifting algorithm on coefficients. Let p =
1687     * f_i.ring.coFac.modul() and assume C == prod_{0,...,n-1} f_i mod p with
1688     * gcd(f_i,f_j) == 1 mod p for i != j
1689     * @param C integer polynomial
1690     * @param F = [f_0,...,f_{n-1}] list of monic modular polynomials.
1691     * @param k approximation exponent.
1692     * @param g leading coefficient.
1693     * @return [g_0,...,g_{n-1}] with C = prod_{0,...,n-1} g_i mod p^k.
1694     */
1695    @SuppressWarnings("unchecked")
1696    public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHensel(
1697                    GenPolynomial<BigInteger> C, List<GenPolynomial<MOD>> F, long k, BigInteger g)
1698                    throws NoLiftingException {
1699        if (C == null || C.isZERO() || F == null || F.size() == 0) {
1700            throw new IllegalArgumentException("C must be nonzero and F must be nonempty");
1701        }
1702        GenPolynomialRing<BigInteger> fac = C.ring;
1703        if (fac.nvar != 1) { // assert ?
1704            throw new IllegalArgumentException("polynomial ring not univariate");
1705        }
1706        List<GenPolynomial<MOD>> lift = new ArrayList<GenPolynomial<MOD>>(F.size());
1707        GenPolynomialRing<MOD> pfac = F.get(0).ring;
1708        RingFactory<MOD> pcfac = pfac.coFac;
1709        ModularRingFactory<MOD> PF = (ModularRingFactory<MOD>) pcfac;
1710        BigInteger P = PF.getIntegerModul();
1711        int n = F.size();
1712        if (n == 1) { // lift F_0, this case will probably never be used
1713            GenPolynomial<MOD> f = F.get(0);
1714            ModularRingFactory<MOD> mcfac;
1715            if (ModLongRing.MAX_LONG.compareTo(P.getVal()) > 0) {
1716                mcfac = (ModularRingFactory) new ModLongRing(P.getVal());
1717            } else {
1718                mcfac = (ModularRingFactory) new ModIntegerRing(P.getVal());
1719            }
1720            GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1721            f = PolyUtil.fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac, f));
1722            lift.add(f);
1723            return lift;
1724        }
1725        // if (n == 2) { // only one step
1726        //     HenselApprox<MOD> ab = HenselUtil.<MOD> liftHenselQuadratic(C, M, F.get(0), F.get(1));
1727        //     lift.add(ab.Am);
1728        //     lift.add(ab.Bm);
1729        //     return lift;
1730        // }
1731
1732        // normalize C and F_i factors
1733        BigInteger cc = g; //C.leadingBaseCoefficient(); // == g ??
1734        for (int i = 1; i < F.size(); i++) { // #F-1
1735            C = C.multiply(cc); // sic
1736        }
1737        MOD cm = PF.fromInteger(cc.getVal());
1738        List<GenPolynomial<MOD>> Fp = new ArrayList<GenPolynomial<MOD>>(F.size());
1739        for (GenPolynomial<MOD> fm : F) {
1740            GenPolynomial<MOD> am = fm.monic();
1741            am = am.multiply(cm);
1742            Fp.add(am);
1743        }
1744        F = Fp;
1745
1746        // setup integer polynomial ring
1747        GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1748        List<GenPolynomial<BigInteger>> Fi = PolyUtil.integerFromModularCoefficients(ifac, F);
1749        //System.out.println("Fi = " + Fi);
1750
1751        // inplace modify polynomials, replace leading coefficient
1752        for (GenPolynomial<BigInteger> ai : Fi) {
1753            if (ai.isZERO()) {
1754                continue;
1755            }
1756            ExpVector ea = ai.leadingExpVector();
1757            ai.doPutToMap(ea, cc);
1758        }
1759        //System.out.println("Fi = " + Fi);
1760
1761        List<GenPolynomial<MOD>> S = liftExtendedEuclidean(F, k + 1); // lift works for any k, use this
1762        //System.out.println("Sext = " + S);
1763        if (debug) {
1764            logger.info("EE lift = {}", S);
1765            // adjust coefficients
1766            List<GenPolynomial<MOD>> Sx = PolyUtil.fromIntegerCoefficients(pfac,
1767                            PolyUtil.integerFromModularCoefficients(ifac, S));
1768            try {
1769                boolean il = HenselUtil.<MOD> isExtendedEuclideanLift(F, Sx);
1770                //System.out.println("islift = " + il);
1771            } catch (RuntimeException e) {
1772                e.printStackTrace();
1773            }
1774        }
1775        List<GenPolynomial<BigInteger>> Si = PolyUtil.integerFromModularCoefficients(ifac, S);
1776        //System.out.println("Si = " + Si);
1777        //System.out.println("C = " + C);
1778
1779        // approximate mod p^i
1780        ModularRingFactory<MOD> mcfac = PF;
1781        BigInteger p = mcfac.getIntegerModul();
1782        BigInteger modul = p;
1783        GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1784        List<GenPolynomial<MOD>> Sp = PolyUtil.fromIntegerCoefficients(mfac, Si);
1785        //System.out.println("Sp = " + Sp);
1786        for (int i = 1; i < k; i++) {
1787            //System.out.println("i = " + i);
1788            GenPolynomial<BigInteger> e = fac.getONE();
1789            for (GenPolynomial<BigInteger> fi : Fi) {
1790                e = e.multiply(fi);
1791            }
1792            e = C.subtract(e);
1793            //System.out.println("\ne = " + e);
1794            if (e.isZERO()) {
1795                logger.info("leaving on zero e");
1796                break;
1797            }
1798            try {
1799                e = e.divide(modul);
1800            } catch (RuntimeException ex) {
1801                ex.printStackTrace();
1802                throw ex;
1803            }
1804            //System.out.println("e = " + e);
1805            // move to in Z_p[x]
1806            GenPolynomial<MOD> c = PolyUtil.<MOD> fromIntegerCoefficients(mfac, e);
1807            //System.out.println("c = " + c + ": " + c.ring.coFac);
1808
1809            List<GenPolynomial<MOD>> s = new ArrayList<GenPolynomial<MOD>>(S.size());
1810            int j = 0;
1811            for (GenPolynomial<MOD> f : Sp) {
1812                f = f.multiply(c);
1813                //System.out.println("f = " + f + " : " + f.ring.coFac);
1814                //System.out.println("F,i = " + F.get(j) + " : " + F.get(j).ring.coFac);
1815                f = f.remainder(F.get(j++));
1816                //System.out.println("f = " + f + " : " + f.ring.coFac);
1817                s.add(f);
1818            }
1819            //System.out.println("s = " + s);
1820            List<GenPolynomial<BigInteger>> si = PolyUtil.integerFromModularCoefficients(ifac, s);
1821            //System.out.println("si = " + si);
1822
1823            List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(F.size());
1824            j = 0;
1825            for (GenPolynomial<BigInteger> f : Fi) {
1826                f = f.sum(si.get(j++).multiply(modul));
1827                Fii.add(f);
1828            }
1829            //System.out.println("Fii = " + Fii);
1830            Fi = Fii;
1831            modul = modul.multiply(p);
1832            if (i >= k - 1) {
1833                logger.info("e != 0 for k = {}", k);
1834            }
1835        }
1836        //System.out.println("Fi = " + Fi);
1837        // remove normalization
1838        GreatestCommonDivisorAbstract<BigInteger> ufd = GCDFactory.getImplementation(cc);
1839        //BigInteger ai = ufd.baseContent(Fi.get(0));
1840        //System.out.println("ai = " + ai + ", cc = " + cc);
1841        List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(F.size());
1842        //int j = 0;
1843        for (GenPolynomial<BigInteger> bi : Fi) {
1844            GenPolynomial<BigInteger> ci = null;
1845            //if ( j++ == 0 ) {
1846            //    ci = bi.divide(ai);
1847            //} else {
1848            //    BigInteger i = cc.divide(ai);
1849            //    ci = bi.divide(i);
1850            //}
1851            ci = ufd.basePrimitivePart(bi); // ??
1852            //System.out.println("bi = " + bi + ", ci = " + ci);
1853            Fii.add(ci);
1854        }
1855        Fi = Fii;
1856
1857        // setup ring mod p^k
1858        modul = p.power(k);
1859        if (ModLongRing.MAX_LONG.compareTo(modul.getVal()) > 0) {
1860            mcfac = (ModularRingFactory) new ModLongRing(modul.getVal());
1861        } else {
1862            mcfac = (ModularRingFactory) new ModIntegerRing(modul.getVal());
1863        }
1864        //System.out.println("mcfac = " + mcfac);
1865        mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1866        lift = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Fi);
1867        //System.out.println("lift = " + lift + ": " + lift.get(0).ring.coFac);
1868        return lift;
1869    }
1870
1871}