Source code

001/*
002 * $Id: GreatestCommonDivisorHensel.java 5047 2014-12-30 17:44:11Z kredel $
003 */
004
005package edu.jas.ufd;
006
007
008import java.util.ArrayList;
009import java.util.Iterator;
010import java.util.List;
011
012import org.apache.log4j.Logger;
013
014import edu.jas.arith.BigInteger;
015import edu.jas.arith.ModIntegerRing;
016import edu.jas.arith.ModLongRing;
017import edu.jas.arith.Modular;
018import edu.jas.arith.ModularRingFactory;
019import edu.jas.arith.PrimeList;
020import edu.jas.poly.ExpVector;
021import edu.jas.poly.GenPolynomial;
022import edu.jas.poly.GenPolynomialRing;
023import edu.jas.poly.PolyUtil;
024import edu.jas.structure.GcdRingElem;
025import edu.jas.structure.NotInvertibleException;
026import edu.jas.structure.Power;
027import edu.jas.structure.RingFactory;
028
029
030// import edu.jas.application.Ideal;
031
032
033/**
034 * Greatest common divisor algorithms with subresultant polynomial remainder
035 * sequence and univariate Hensel lifting.
036 * @author Heinz Kredel
037 */
038
039public class GreatestCommonDivisorHensel<MOD extends GcdRingElem<MOD> & Modular> extends
040                GreatestCommonDivisorAbstract<BigInteger> {
041
042
043    private static final Logger logger = Logger.getLogger(GreatestCommonDivisorHensel.class);
044
045
046    private final boolean debug = logger.isDebugEnabled();
047
048
049    /**
050     * Flag for linear or quadratic Hensel lift.
051     */
052    public final boolean quadratic;
053
054
055    /**
056     * Fall back gcd algorithm.
057     */
058    public final GreatestCommonDivisorAbstract<BigInteger> iufd;
059
060
061    /*
062     * Internal dispatcher.
063     */
064    private final GreatestCommonDivisorAbstract<BigInteger> ufd;
065
066
067    /**
068     * Constructor.
069     */
070    public GreatestCommonDivisorHensel() {
071        this(true);
072    }
073
074
075    /**
076     * Constructor.
077     * @param quadratic use quadratic Hensel lift.
078     */
079    public GreatestCommonDivisorHensel(boolean quadratic) {
080        this.quadratic = quadratic;
081        iufd = new GreatestCommonDivisorSubres<BigInteger>();
082        ufd = this; //iufd;
083    }
084
085
086    /**
087     * Univariate GenPolynomial greatest comon divisor. Uses univariate Hensel
088     * lifting.
089     * @param P univariate GenPolynomial.
090     * @param S univariate GenPolynomial.
091     * @return gcd(P,S).
092     */
093    @Override
094    @SuppressWarnings("unchecked")
095    public GenPolynomial<BigInteger> baseGcd(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) {
096        if (S == null || S.isZERO()) {
097            return P;
098        }
099        if (P == null || P.isZERO()) {
100            return S;
101        }
102        if (P.ring.nvar > 1) {
103            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
104        }
105        GenPolynomialRing<BigInteger> fac = P.ring;
106        long e = P.degree(0);
107        long f = S.degree(0);
108        GenPolynomial<BigInteger> q;
109        GenPolynomial<BigInteger> r;
110        if (f > e) {
111            r = P;
112            q = S;
113            long g = f;
114            f = e;
115            e = g;
116        } else {
117            q = P;
118            r = S;
119        }
120        if (debug) {
121            logger.debug("degrees: e = " + e + ", f = " + f);
122        }
123        r = r.abs();
124        q = q.abs();
125        // compute contents and primitive parts
126        BigInteger a = baseContent(r);
127        BigInteger b = baseContent(q);
128        // gcd of coefficient contents
129        BigInteger c = gcd(a, b); // indirection
130        r = divide(r, a); // indirection
131        q = divide(q, b); // indirection
132        if (r.isONE()) {
133            return r.multiply(c);
134        }
135        if (q.isONE()) {
136            return q.multiply(c);
137        }
138        // compute normalization factor
139        BigInteger ac = r.leadingBaseCoefficient();
140        BigInteger bc = q.leadingBaseCoefficient();
141        BigInteger cc = gcd(ac, bc); // indirection
142        // compute degree vectors, only univeriate
143        ExpVector rdegv = r.degreeVector();
144        ExpVector qdegv = q.degreeVector();
145        //initialize prime list and degree vector
146        PrimeList primes = new PrimeList(PrimeList.Range.medium);
147        int pn = 50; //primes.size();
148
149        ModularRingFactory<MOD> cofac;
150        GenPolynomial<MOD> qm;
151        GenPolynomial<MOD> qmf;
152        GenPolynomial<MOD> rm;
153        GenPolynomial<MOD> rmf;
154        GenPolynomial<MOD> cmf;
155        GenPolynomialRing<MOD> mfac;
156        GenPolynomial<MOD> cm = null;
157        GenPolynomial<MOD>[] ecm = null;
158        GenPolynomial<MOD> sm = null;
159        GenPolynomial<MOD> tm = null;
160        HenselApprox<MOD> lift = null;
161        if (debug) {
162            logger.debug("c = " + c);
163            logger.debug("cc = " + cc);
164            logger.debug("primes = " + primes);
165        }
166
167        int i = 0;
168        for (java.math.BigInteger p : primes) {
169            //System.out.println("next run ++++++++++++++++++++++++++++++++++");
170            if (++i >= pn) {
171                logger.error("prime list exhausted, pn = " + pn);
172                //logger.info("primes = " + primes);
173                return iufd.baseGcd(P, S);
174                //throw new ArithmeticException("prime list exhausted");
175            }
176            // initialize coefficient factory and map normalization factor
177            //cofac = new ModIntegerRing(p, true);
178            if (ModLongRing.MAX_LONG.compareTo(p) > 0) {
179                cofac = (ModularRingFactory) new ModLongRing(p, true);
180            } else {
181                cofac = (ModularRingFactory) new ModIntegerRing(p, true);
182            }
183            MOD nf = cofac.fromInteger(cc.getVal());
184            if (nf.isZERO()) {
185                continue;
186            }
187            nf = cofac.fromInteger(q.leadingBaseCoefficient().getVal());
188            if (nf.isZERO()) {
189                continue;
190            }
191            nf = cofac.fromInteger(r.leadingBaseCoefficient().getVal());
192            if (nf.isZERO()) {
193                continue;
194            }
195            // initialize polynomial factory and map polynomials
196            mfac = new GenPolynomialRing<MOD>(cofac, fac.nvar, fac.tord, fac.getVars());
197            qm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, q);
198            if (!qm.degreeVector().equals(qdegv)) {
199                continue;
200            }
201            rm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, r);
202            if (!rm.degreeVector().equals(rdegv)) {
203                continue;
204            }
205            if (debug) {
206                logger.info("cofac = " + cofac.getIntegerModul());
207            }
208
209            // compute univariate modular gcd
210            cm = qm.gcd(rm);
211
212            // test for constant g.c.d
213            if (cm.isConstant()) {
214                logger.debug("cm, constant = " + cm);
215                return fac.getONE().multiply(c);
216            }
217
218            // compute factors and gcd with factor
219            GenPolynomial<BigInteger> crq;
220            rmf = rm.divide(cm); // rm = cm * rmf
221            ecm = cm.egcd(rmf);
222            if (ecm[0].isONE()) {
223                //logger.debug("gcd() first factor " + rmf);
224                crq = r;
225                cmf = rmf;
226                sm = ecm[1];
227                tm = ecm[2];
228            } else {
229                qmf = qm.divide(cm); // qm = cm * qmf
230                ecm = cm.egcd(qmf);
231                if (ecm[0].isONE()) {
232                    //logger.debug("gcd() second factor " + qmf);
233                    crq = q;
234                    cmf = qmf;
235                    sm = ecm[1];
236                    tm = ecm[2];
237                } else {
238                    logger.info("both gcd != 1: Hensel not applicable");
239                    return iufd.baseGcd(P, S);
240                }
241            }
242            BigInteger cn = crq.maxNorm();
243            cn = cn.multiply(crq.leadingBaseCoefficient().abs());
244            cn = cn.multiply(cn.fromInteger(2));
245            if (debug) {
246                System.out.println("crq = " + crq);
247                System.out.println("cm  = " + cm);
248                System.out.println("cmf = " + cmf);
249                System.out.println("sm  = " + sm);
250                System.out.println("tm  = " + tm);
251                System.out.println("cn  = " + cn);
252            }
253            try {
254                if (quadratic) {
255                    lift = HenselUtil.liftHenselQuadratic(crq, cn, cm, cmf, sm, tm);
256                } else {
257                    lift = HenselUtil.liftHensel(crq, cn, cm, cmf, sm, tm);
258                }
259            } catch (NoLiftingException nle) {
260                logger.info("giving up on Hensel gcd reverting to Subres gcd " + nle);
261                return iufd.baseGcd(P, S);
262            }
263            q = lift.A;
264            if (debug) {
265                System.out.println("q   = " + q);
266                System.out.println("qf  = " + lift.B);
267            }
268            q = basePrimitivePart(q);
269            q = q.multiply(c).abs();
270            if (PolyUtil.<BigInteger> baseSparsePseudoRemainder(P, q).isZERO()
271                            && PolyUtil.<BigInteger> baseSparsePseudoRemainder(S, q).isZERO()) {
272                break;
273            }
274            logger.info("final devision not successfull");
275            //System.out.println("P rem q = " + PolyUtil.<BigInteger>basePseudoRemainder(P,q));
276            //System.out.println("S rem q = " + PolyUtil.<BigInteger>basePseudoRemainder(S,q));
277            //break;
278        }
279        return q;
280    }
281
282
283    /**
284     * Univariate GenPolynomial recursive greatest comon divisor. Uses
285     * multivariate Hensel list.
286     * @param P univariate recursive GenPolynomial.
287     * @param S univariate recursive GenPolynomial.
288     * @return gcd(P,S).
289     */
290    @Override
291    @SuppressWarnings("unchecked")
292    public GenPolynomial<GenPolynomial<BigInteger>> recursiveUnivariateGcd(
293                    GenPolynomial<GenPolynomial<BigInteger>> P, GenPolynomial<GenPolynomial<BigInteger>> S) {
294        if (S == null || S.isZERO()) {
295            return P;
296        }
297        if (P == null || P.isZERO()) {
298            return S;
299        }
300        if (P.ring.nvar > 1) {
301            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
302        }
303        long e = P.degree(0);
304        long f = S.degree(0);
305        GenPolynomial<GenPolynomial<BigInteger>> q, r, s;
306        if (f > e) {
307            r = P;
308            q = S;
309            long g = f;
310            f = e;
311            e = g;
312        } else {
313            q = P;
314            r = S;
315        }
316        if (debug) {
317            logger.debug("degrees: e = " + e + ", f = " + f);
318        }
319        r = r.abs();
320        q = q.abs();
321        //logger.info("r: " + r + ", q: " + q);
322
323        GenPolynomial<BigInteger> a = ufd.recursiveContent(r);
324        GenPolynomial<BigInteger> b = ufd.recursiveContent(q);
325
326        GenPolynomial<BigInteger> c = ufd.gcd(a, b); // go to recursion
327        //System.out.println("rgcd c = " + c);
328        r = PolyUtil.<BigInteger> recursiveDivide(r, a);
329        q = PolyUtil.<BigInteger> recursiveDivide(q, b);
330        a = PolyUtil.<BigInteger> basePseudoDivide(a, c); // unused ?
331        b = PolyUtil.<BigInteger> basePseudoDivide(b, c); // unused ?
332        if (r.isONE()) {
333            return r.multiply(c);
334        }
335        if (q.isONE()) {
336            return q.multiply(c);
337        }
338        // check constant ldcf, TODO general case
339        GenPolynomial<BigInteger> la, lb, lc, lh;
340        la = r.leadingBaseCoefficient();
341        lb = q.leadingBaseCoefficient();
342        lc = ufd.gcd(la, lb);
343        //logger.info("la = " + la + ", lb = " + lb + ", lc = " + lc);
344        if (!lc.isConstant()) {
345            //continue; // easy way out
346            GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(r, q);
347            T = T.abs().multiply(c);
348            logger.info("non monic ldcf (" + lc + ") not implemented: " + T + "= gcd(" + r + "," + q + ") * "
349                            + c);
350            return T;
351        }
352
353        // convert from Z[y1,...,yr][x] to Z[x][y1,...,yr] to Z[x,y1,...,yr]
354        GenPolynomial<GenPolynomial<BigInteger>> qs = PolyUtil.<BigInteger> switchVariables(q);
355        GenPolynomial<GenPolynomial<BigInteger>> rs = PolyUtil.<BigInteger> switchVariables(r);
356
357        GenPolynomialRing<GenPolynomial<BigInteger>> rfac = qs.ring;
358        RingFactory<GenPolynomial<BigInteger>> rrfac = rfac.coFac;
359        GenPolynomialRing<BigInteger> cfac = (GenPolynomialRing<BigInteger>) rrfac;
360        GenPolynomialRing<BigInteger> dfac = cfac.extend(rfac.getVars());
361        //System.out.println("pfac = " + P.ring.toScript());
362        //System.out.println("rfac = " + rfac.toScript());
363        //System.out.println("dfac = " + dfac.toScript());
364        GenPolynomial<BigInteger> qd = PolyUtil.<BigInteger> distribute(dfac, qs);
365        GenPolynomial<BigInteger> rd = PolyUtil.<BigInteger> distribute(dfac, rs);
366
367        // compute normalization factor
368        BigInteger ac = rd.leadingBaseCoefficient();
369        BigInteger bc = qd.leadingBaseCoefficient();
370        BigInteger cc = gcd(ac, bc); // indirection
371
372        //initialize prime list
373        PrimeList primes = new PrimeList(PrimeList.Range.medium);
374        Iterator<java.math.BigInteger> primeIter = primes.iterator();
375        int pn = 50; //primes.size();
376
377        // double check variables
378        // need qe,re,qd,rd,a,b
379        GenPolynomial<BigInteger> qe0, re0, ce0 = null;
380
381        for (int i = 0; i < 11; i++) { // meta loop
382            //System.out.println("======== run " + dfac.nvar + ", " + i);
383            java.math.BigInteger p = null; //new java.math.BigInteger("19"); //primes.next();
384            // 5 small, 4 medium and 2 large size primes
385            if (i == 0) { // medium size
386                primes = new PrimeList(PrimeList.Range.medium);
387                primeIter = primes.iterator();
388            }
389            if (i == 4) { // small size
390                primes = new PrimeList(PrimeList.Range.small);
391                primeIter = primes.iterator();
392                p = primeIter.next(); // 2
393                p = primeIter.next(); // 3
394                p = primeIter.next(); // 5
395                p = primeIter.next(); // 7
396            }
397            if (i == 9) { // large size
398                primes = new PrimeList(PrimeList.Range.large);
399                primeIter = primes.iterator();
400            }
401            ModularRingFactory<MOD> cofac = null;
402            int pi = 0;
403            while (pi < pn && primeIter.hasNext()) {
404                p = primeIter.next();
405                //p = new java.math.BigInteger("19");
406                logger.info("prime = " + p);
407                // initialize coefficient factory and map normalization factor and polynomials
408                ModularRingFactory<MOD> cf = null;
409                if (ModLongRing.MAX_LONG.compareTo(p) > 0) {
410                    cf = (ModularRingFactory) new ModLongRing(p, true);
411                } else {
412                    cf = (ModularRingFactory) new ModIntegerRing(p, true);
413                }
414                MOD nf = cf.fromInteger(cc.getVal());
415                if (nf.isZERO()) {
416                    continue;
417                }
418                nf = cf.fromInteger(q.leadingBaseCoefficient().leadingBaseCoefficient().getVal());
419                if (nf.isZERO()) {
420                    continue;
421                }
422                nf = cf.fromInteger(r.leadingBaseCoefficient().leadingBaseCoefficient().getVal());
423                if (nf.isZERO()) {
424                    continue;
425                }
426                cofac = cf;
427                break;
428            }
429            if (cofac == null) { // no lucky prime found
430                GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(q, r);
431                logger.info("no lucky prime, gave up on Hensel: " + T + "= gcd(" + r + "," + q + ")");
432                return T.abs().multiply(c); //.abs();
433            }
434            //System.out.println("cofac = " + cofac);
435
436            // search evaluation points and evaluate
437            List<BigInteger> V = new ArrayList<BigInteger>(P.ring.nvar);
438            GenPolynomialRing<BigInteger> ckfac = dfac;
439            GenPolynomial<BigInteger> qe = qd;
440            GenPolynomial<BigInteger> re = rd;
441            GenPolynomial<BigInteger> qei;
442            GenPolynomial<BigInteger> rei;
443            for (int j = dfac.nvar; j > 1; j--) {
444                // evaluation to univariate case
445                long degq = qe.degree(ckfac.nvar - 2);
446                long degr = re.degree(ckfac.nvar - 2);
447                ckfac = ckfac.contract(1);
448                long vi = 1L; //(long)(dfac.nvar-j); // 1L; 0 not so good for small p
449                if (p.longValue() > 1000L) {
450                    //vi = (long)j+1L;
451                    vi = 0L;
452                }
453                // search small evaluation point
454                while (true) {
455                    MOD vp = cofac.fromInteger(vi++);
456                    //System.out.println("vp = " + vp);
457                    if (vp.isZERO() && vi != 1L) { // all elements of Z_p exhausted
458                        qe = null;
459                        re = null;
460                        break;
461                    }
462                    BigInteger vii = new BigInteger(vi - 1);
463                    qei = PolyUtil.<BigInteger> evaluateMain(ckfac, qe, vii);
464                    rei = PolyUtil.<BigInteger> evaluateMain(ckfac, re, vii);
465                    //System.out.println("qei = " + qei);
466                    //System.out.println("rei = " + rei);
467
468                    // check lucky evaluation point 
469                    if (degq != qei.degree(ckfac.nvar - 1)) {
470                        //System.out.println("degv(qe) = " + qe.degreeVector());
471                        //System.out.println("deg(qe) = " + degq + ", deg(qe) = " + qei.degree(ckfac.nvar-1));
472                        continue;
473                    }
474                    if (degr != rei.degree(ckfac.nvar - 1)) {
475                        //System.out.println("degv(re) = " + re.degreeVector());
476                        //System.out.println("deg(re) = " + degr + ", deg(re) = " + rei.degree(ckfac.nvar-1));
477                        continue;
478                    }
479                    V.add(vii);
480                    qe = qei;
481                    re = rei;
482                    break;
483                }
484                if (qe == null && re == null) {
485                    break;
486                }
487            }
488            if (qe == null && re == null) {
489                continue;
490            }
491            logger.info("evaluation points  = " + V);
492
493            // recursion base:
494            GenPolynomial<BigInteger> ce = ufd.baseGcd(qe, re);
495            if (ce.isConstant()) {
496                return P.ring.getONE().multiply(c);
497            }
498            logger.info("base gcd = " + ce);
499
500            // double check 
501            // need qe,re,qd,rd,a,b
502            if (i == 0) {
503                qe0 = qe;
504                re0 = re;
505                ce0 = ce;
506                continue;
507            }
508            long d0 = ce0.degree(0);
509            long d1 = ce.degree(0);
510            //System.out.println("d0, d1 = " + d0 + ", " + d1);
511            if (d1 < d0) {
512                qe0 = qe;
513                re0 = re;
514                ce0 = ce;
515                continue;
516            } else if (d1 > d0) {
517                continue;
518            }
519            // d0 == d1 is ok
520            long dx = r.degree(0);
521            //System.out.println("d0, dx = " + d0 + ", " + dx);
522            if (d0 == dx) { // gcd == r ?
523                if (PolyUtil.<BigInteger> recursivePseudoRemainder(q, r).isZERO()) {
524                    r = r.abs().multiply(c); //.abs();
525                    logger.info("exit with r | q : " + r);
526                    return r;
527                }
528                continue;
529            }
530            // norm
531            BigInteger mn = null; //mn = mn.multiply(cc).multiply(mn.fromInteger(2));
532            // prepare lifting, chose factor polynomials
533            GenPolynomial<BigInteger> re1 = PolyUtil.<BigInteger> basePseudoDivide(re, ce);
534            GenPolynomial<BigInteger> qe1 = PolyUtil.<BigInteger> basePseudoDivide(qe, ce);
535            GenPolynomial<BigInteger> ui, he, pe;
536            GenPolynomial<BigInteger> g, gi, lui;
537            GenPolynomial<BigInteger> gcr, gcq;
538            gcr = ufd.baseGcd(re1, ce);
539            gcq = ufd.baseGcd(qe1, ce);
540            if (gcr.isONE() && gcq.isONE()) { // both gcds == 1: chose smaller ldcf
541                if (la.totalDegree() > lb.totalDegree()) {
542                    ui = qd;
543                    s = q;
544                    he = qe1;
545                    pe = qe;
546                    BigInteger bn = qd.maxNorm();
547                    mn = bn.multiply(cc).multiply(new BigInteger(2L));
548                    g = lb;
549                    logger.debug("select deg: ui = qd, g = b"); //, qe1 = " + qe1); // + ", qe = " + qe);
550                } else {
551                    ui = rd;
552                    s = r;
553                    he = re1;
554                    pe = re;
555                    BigInteger an = rd.maxNorm();
556                    mn = an.multiply(cc).multiply(new BigInteger(2L));
557                    g = la;
558                    logger.debug("select deg: ui = rd, g = a"); //, re1 = " + re1); // + ", re = " + re);
559                }
560            } else if (gcr.isONE()) {
561                ui = rd;
562                s = r;
563                he = re1;
564                pe = re;
565                BigInteger an = rd.maxNorm();
566                mn = an.multiply(cc).multiply(new BigInteger(2L));
567                g = la;
568                logger.debug("select: ui = rd, g = a"); //, re1 = " + re1); // + ", re = " + re);
569            } else if (gcq.isONE()) {
570                ui = qd;
571                s = q;
572                he = qe1;
573                pe = qe;
574                BigInteger bn = qd.maxNorm();
575                mn = bn.multiply(cc).multiply(new BigInteger(2L));
576                g = lb;
577                logger.debug("select: ui = qd, g = b"); //, qe1 = " + qe1); // + ", qe = " + qe);
578            } else { // both gcds != 1: method not applicable
579                logger.info("both gcds != 1: method not applicable");
580                break;
581            }
582            lui = lc; //s.leadingBaseCoefficient();
583            lh = PolyUtil.<BigInteger> basePseudoDivide(g, lui);
584            BigInteger ge = PolyUtil.<BigInteger> evaluateAll(g.ring.coFac, lui, V);
585            if (ge.isZERO()) {
586                continue;
587            }
588            BigInteger geh = PolyUtil.<BigInteger> evaluateAll(g.ring.coFac, lh, V);
589            if (geh.isZERO()) {
590                continue;
591            }
592            BigInteger gg = PolyUtil.<BigInteger> evaluateAll(g.ring.coFac, g, V);
593            if (gg.isZERO()) {
594                continue;
595            }
596            //System.out.println("ge = " + ge + ", geh = " + geh + ", gg = " + gg + ", pe = " + pe); 
597            // 
598            //ce = ce.multiply(geh); //ge);
599            // 
600            he = he.multiply(ge); //gg); //ge); //geh);
601            //
602            gi = lui.extendLower(dfac, 0, 0L); //lui. // g.
603            //
604            ui = ui.multiply(gi); // gi !.multiply(gi) 
605            //System.out.println("ui = " + ui + ", deg(ui) = " + ui.degreeVector());
606            //System.out.println("ce = " + ce + ", he = " + he + ", ge = " + ge);
607            logger.info("gcd(ldcf): " + lui + ", ldcf cofactor: " + lh + ", base cofactor: " + he);
608
609            long k = Power.logarithm(new BigInteger(p), mn);
610            //System.out.println("mn = " + mn);
611            //System.out.println("k = " + k);
612
613            BigInteger qp = Power.positivePower(cofac.getIntegerModul(), k); // == p
614            ModularRingFactory<MOD> muqfac;
615            if (ModLongRing.MAX_LONG.compareTo(qp.getVal()) > 0) {
616                muqfac = (ModularRingFactory) new ModLongRing(qp.getVal(), true); // nearly a field
617            } else {
618                muqfac = (ModularRingFactory) new ModIntegerRing(qp.getVal(), true); // nearly a field
619            }
620            GenPolynomialRing<MOD> mucpfac = new GenPolynomialRing<MOD>(muqfac, ckfac);
621            //System.out.println("mucpfac = " + mucpfac.toScript());
622            if (muqfac.fromInteger(ge.getVal()).isZERO()) {
623                continue;
624            }
625            //GenPolynomial<BigInteger> xxx = invertPoly(muqfac,lui,V);
626            //System.out.println("inv(lui) = " + xxx + ", muqfac = " + muqfac + ", lui = " + lui);
627            //ce = ce.multiply(xxx); //.leadingBaseCoefficient()); 
628            //xxx = invertPoly(muqfac,lh,V);
629            //System.out.println("inv(lh) = " + xxx + ", muqfac = " + muqfac + ", lh = " + lh);
630            //he = he.multiply(xxx); //.leadingBaseCoefficient()); 
631
632            GenPolynomial<MOD> cm = PolyUtil.<MOD> fromIntegerCoefficients(mucpfac, ce);
633            GenPolynomial<MOD> hm = PolyUtil.<MOD> fromIntegerCoefficients(mucpfac, he);
634            if (cm.degree(0) != ce.degree(0) || hm.degree(0) != he.degree(0)) {
635                continue;
636            }
637            if (cm.isZERO() || hm.isZERO()) {
638                continue;
639            }
640            logger.info("univariate modulo p^k: " + cm + ", " + hm);
641
642            // convert C from Z[...] to Z_q[...]
643            GenPolynomialRing<MOD> qcfac = new GenPolynomialRing<MOD>(muqfac, dfac);
644            GenPolynomial<MOD> uq = PolyUtil.<MOD> fromIntegerCoefficients(qcfac, ui);
645            if (!ui.leadingExpVector().equals(uq.leadingExpVector())) {
646                logger.info("ev(ui) = " + ui.leadingExpVector() + ", ev(uq) = " + uq.leadingExpVector());
647                continue;
648            }
649            logger.info("multivariate modulo p^k: " + uq);
650
651            List<GenPolynomial<MOD>> F = new ArrayList<GenPolynomial<MOD>>(2);
652            F.add(cm);
653            F.add(hm);
654            List<GenPolynomial<BigInteger>> G = new ArrayList<GenPolynomial<BigInteger>>(2);
655            G.add(lui.ring.getONE()); //lui: lui.ring.getONE()); // TODO 
656            G.add(lui.ring.getONE()); //lh: lui);
657            List<GenPolynomial<MOD>> lift;
658            try {
659                //lift = HenselMultUtil.<MOD> liftHenselFull(ui, F, V, k, G);
660                lift = HenselMultUtil.<MOD> liftHensel(ui, uq, F, V, k, G);
661                logger.info("lift = " + lift);
662            } catch (NoLiftingException nle) {
663                logger.info("NoLiftingException");
664                //System.out.println("exception : " + nle);
665                continue;
666            } catch (ArithmeticException ae) {
667                logger.info("ArithmeticException");
668                //System.out.println("exception : " + ae);
669                continue;
670            } catch (NotInvertibleException ni) {
671                logger.info("NotInvertibleException");
672                //System.out.println("exception : " + ni);
673                continue;
674            }
675            //if (!HenselMultUtil.<MOD> isHenselLift(ui, uq, F, k, lift)) { // not meaningfull test
676            //    logger.info("isHenselLift: false");
677            //    //continue;
678            //}
679
680            // convert Ci from Z_{p^k}[x,y1,...,yr] to Z[x,y1,...,yr] to Z[x][y1,...,yr] to Z[y1,...,yr][x]
681            GenPolynomial<BigInteger> ci = PolyUtil.integerFromModularCoefficients(dfac, lift.get(0));
682            ci = basePrimitivePart(ci);
683            GenPolynomial<GenPolynomial<BigInteger>> Cr = PolyUtil.<BigInteger> recursive(rfac, ci);
684            GenPolynomial<GenPolynomial<BigInteger>> Cs = PolyUtil.<BigInteger> switchVariables(Cr);
685            if (!Cs.ring.equals(P.ring)) {
686                System.out.println("Cs.ring = " + Cs.ring + ", P.ring = " + P.ring);
687            }
688            GenPolynomial<GenPolynomial<BigInteger>> Q = ufd.recursivePrimitivePart(Cs);
689            Q = ufd.baseRecursivePrimitivePart(Q);
690            Q = Q.abs().multiply(c); //.abs();
691            GenPolynomial<GenPolynomial<BigInteger>> Pq, Sq;
692            Pq = PolyUtil.<BigInteger> recursivePseudoRemainder(P, Q);
693            Sq = PolyUtil.<BigInteger> recursivePseudoRemainder(S, Q);
694            if (Pq.isZERO() && Sq.isZERO()) {
695                logger.info("gcd normal exit: " + Q);
696                return Q;
697            }
698            logger.info("bad Q = " + Q); // + ", Pq = " + Pq + ", Sq = " + Sq);
699        } // end for meta loop
700          // Hensel gcd failed
701        GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(r, q);
702        T = T.abs().multiply(c);
703        logger.info("no lucky prime or evaluation points, gave up on Hensel: " + T + "= gcd(" + r + "," + q
704                        + ")");
705        return T;
706    }
707
708
709    GenPolynomial<BigInteger> invertPoly(ModularRingFactory<MOD> mfac, GenPolynomial<BigInteger> li,
710                    List<BigInteger> V) {
711        if (li == null || li.isZERO()) {
712            throw new RuntimeException("li not invertible: " + li);
713        }
714        if (li.isONE()) {
715            return li;
716        }
717        //System.out.println("mfac = " + mfac + ", V = " + V +", li = " + li);
718        GenPolynomialRing<BigInteger> pfac = li.ring;
719        GenPolynomialRing<MOD> mpfac = new GenPolynomialRing<MOD>(mfac, pfac);
720        GenPolynomial<MOD> lm = PolyUtil.<MOD> fromIntegerCoefficients(mpfac, li);
721        //System.out.println("pfac = " + pfac + ", lm = " + lm);
722        List<GenPolynomial<MOD>> lid = new ArrayList<GenPolynomial<MOD>>(V.size());
723        int i = 0;
724        for (BigInteger bi : V) {
725            MOD m = mfac.fromInteger(bi.getVal());
726            GenPolynomial<MOD> mp = mpfac.univariate(i);
727            mp = mp.subtract(m); // X_i - v_i
728            lid.add(mp);
729            i++;
730        }
731        //System.out.println("lid = " + lid);
732        //Ideal<MOD> id = new Ideal<MOD>(mpfac,lid,true); // is a GB
733        //System.out.println("id = " + id);
734        GenPolynomial<MOD> mi = lm; //id.inverse(lm);
735        //System.out.println("mi = " + mi);
736        GenPolynomial<BigInteger> inv = PolyUtil.integerFromModularCoefficients(pfac, mi);
737        return inv;
738    }
739
740}