Source code

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