Source code

001/*
002 * $Id: GreatestCommonDivisorHensel.java 4067 2012-07-27 16:17:35Z 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.Power;
026import edu.jas.structure.RingFactory;
027
028
029/**
030 * Greatest common divisor algorithms with subresultant polynomial remainder
031 * sequence and univariate Hensel lifting.
032 * @author Heinz Kredel
033 */
034
035public class GreatestCommonDivisorHensel<MOD extends GcdRingElem<MOD> & Modular> extends
036                GreatestCommonDivisorAbstract<BigInteger> {
037
038
039    private static final Logger logger = Logger.getLogger(GreatestCommonDivisorHensel.class);
040
041
042    private final boolean debug = logger.isDebugEnabled();
043
044
045    /**
046     * Flag for linear or quadratic Hensel lift.
047     */
048    public final boolean quadratic;
049
050
051    /**
052     * Fall back gcd algorithm.
053     */
054    public final GreatestCommonDivisorAbstract<BigInteger> iufd;
055
056
057    /**
058     * Constructor.
059     */
060    public GreatestCommonDivisorHensel() {
061        this(true);
062    }
063
064
065    /**
066     * Constructor.
067     * @param quadratic use quadratic Hensel lift.
068     */
069    public GreatestCommonDivisorHensel(boolean quadratic) {
070        this.quadratic = quadratic;
071        iufd = new GreatestCommonDivisorSubres<BigInteger>();
072    }
073
074
075    /**
076     * Univariate GenPolynomial greatest comon divisor. Uses univariate Hensel
077     * lifting.
078     * @param P univariate GenPolynomial.
079     * @param S univariate GenPolynomial.
080     * @return gcd(P,S).
081     */
082    @Override
083    public GenPolynomial<BigInteger> baseGcd(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) {
084        if (S == null || S.isZERO()) {
085            return P;
086        }
087        if (P == null || P.isZERO()) {
088            return S;
089        }
090        if (P.ring.nvar > 1) {
091            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
092        }
093        GenPolynomialRing<BigInteger> fac = P.ring;
094        long e = P.degree(0);
095        long f = S.degree(0);
096        GenPolynomial<BigInteger> q;
097        GenPolynomial<BigInteger> r;
098        if (f > e) {
099            r = P;
100            q = S;
101            long g = f;
102            f = e;
103            e = g;
104        } else {
105            q = P;
106            r = S;
107        }
108        if (debug) {
109            logger.debug("degrees: e = " + e + ", f = " + f);
110        }
111        r = r.abs();
112        q = q.abs();
113        // compute contents and primitive parts
114        BigInteger a = baseContent(r);
115        BigInteger b = baseContent(q);
116        // gcd of coefficient contents
117        BigInteger c = gcd(a, b); // indirection
118        r = divide(r, a); // indirection
119        q = divide(q, b); // indirection
120        if (r.isONE()) {
121            return r.multiply(c);
122        }
123        if (q.isONE()) {
124            return q.multiply(c);
125        }
126        // compute normalization factor
127        BigInteger ac = r.leadingBaseCoefficient();
128        BigInteger bc = q.leadingBaseCoefficient();
129        BigInteger cc = gcd(ac, bc); // indirection
130        // compute degree vectors, only univeriate
131        ExpVector rdegv = r.degreeVector();
132        ExpVector qdegv = q.degreeVector();
133        //initialize prime list and degree vector
134        PrimeList primes = new PrimeList(PrimeList.Range.medium);
135        int pn = 50; //primes.size();
136
137        ModularRingFactory<MOD> cofac;
138        GenPolynomial<MOD> qm;
139        GenPolynomial<MOD> qmf;
140        GenPolynomial<MOD> rm;
141        GenPolynomial<MOD> rmf;
142        GenPolynomial<MOD> cmf;
143        GenPolynomialRing<MOD> mfac;
144        GenPolynomial<MOD> cm = null;
145        GenPolynomial<MOD>[] ecm = null;
146        GenPolynomial<MOD> sm = null;
147        GenPolynomial<MOD> tm = null;
148        HenselApprox<MOD> lift = null;
149        if (debug) {
150            logger.debug("c = " + c);
151            logger.debug("cc = " + cc);
152            logger.debug("primes = " + primes);
153        }
154
155        int i = 0;
156        for (java.math.BigInteger p : primes) {
157            //System.out.println("next run ++++++++++++++++++++++++++++++++++");
158            if (++i >= pn) {
159                logger.error("prime list exhausted, pn = " + pn);
160                logger.info("primes = " + primes);
161                return iufd.baseGcd(P, S);
162                //throw new ArithmeticException("prime list exhausted");
163            }
164            // initialize coefficient factory and map normalization factor
165            //cofac = new ModIntegerRing(p, true);
166            if (ModLongRing.MAX_LONG.compareTo(p) > 0) {
167                cofac = (ModularRingFactory) new ModLongRing(p, true);
168            } else {
169                cofac = (ModularRingFactory) new ModIntegerRing(p, true);
170            }
171            MOD nf = cofac.fromInteger(cc.getVal());
172            if (nf.isZERO()) {
173                continue;
174            }
175            nf = cofac.fromInteger(q.leadingBaseCoefficient().getVal());
176            if (nf.isZERO()) {
177                continue;
178            }
179            nf = cofac.fromInteger(r.leadingBaseCoefficient().getVal());
180            if (nf.isZERO()) {
181                continue;
182            }
183            // initialize polynomial factory and map polynomials
184            mfac = new GenPolynomialRing<MOD>(cofac, fac.nvar, fac.tord, fac.getVars());
185            qm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, q);
186            if (!qm.degreeVector().equals(qdegv)) {
187                continue;
188            }
189            rm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, r);
190            if (!rm.degreeVector().equals(rdegv)) {
191                continue;
192            }
193            if (debug) {
194                logger.info("cofac = " + cofac.getIntegerModul());
195            }
196
197            // compute univariate modular gcd
198            cm = qm.gcd(rm);
199
200            // test for constant g.c.d
201            if (cm.isConstant()) {
202                logger.debug("cm, constant = " + cm);
203                return fac.getONE().multiply(c);
204            }
205
206            // compute factors and gcd with factor
207            GenPolynomial<BigInteger> crq;
208            rmf = rm.divide(cm); // rm = cm * rmf
209            ecm = cm.egcd(rmf);
210            if (ecm[0].isONE()) {
211                //logger.debug("gcd() first factor " + rmf);
212                crq = r;
213                cmf = rmf;
214                sm = ecm[1];
215                tm = ecm[2];
216            } else {
217                qmf = qm.divide(cm); // qm = cm * qmf
218                ecm = cm.egcd(qmf);
219                if (ecm[0].isONE()) {
220                    //logger.debug("gcd() second factor " + qmf);
221                    crq = q;
222                    cmf = qmf;
223                    sm = ecm[1];
224                    tm = ecm[2];
225                } else {
226                    logger.info("giving up on Hensel gcd reverting to Subres gcd");
227                    return iufd.baseGcd(P, S);
228                }
229            }
230            BigInteger cn = crq.maxNorm();
231            cn = cn.multiply(crq.leadingBaseCoefficient().abs());
232            cn = cn.multiply(cn.fromInteger(2));
233            if (debug) {
234                System.out.println("crq = " + crq);
235                System.out.println("cm  = " + cm);
236                System.out.println("cmf = " + cmf);
237                System.out.println("sm  = " + sm);
238                System.out.println("tm  = " + tm);
239                System.out.println("cn  = " + cn);
240            }
241            try {
242                if (quadratic) {
243                    lift = HenselUtil.liftHenselQuadratic(crq, cn, cm, cmf, sm, tm);
244                } else {
245                    lift = HenselUtil.liftHensel(crq, cn, cm, cmf, sm, tm);
246                }
247            } catch (NoLiftingException nle) {
248                logger.info("giving up on Hensel gcd reverting to Subres gcd " + nle);
249                return iufd.baseGcd(P, S);
250            }
251            q = lift.A;
252            if (debug) {
253                System.out.println("q   = " + q);
254                System.out.println("qf  = " + lift.B);
255            }
256            q = basePrimitivePart(q);
257            q = q.multiply(c).abs();
258            if (PolyUtil.<BigInteger> baseSparsePseudoRemainder(P, q).isZERO()
259                            && PolyUtil.<BigInteger> baseSparsePseudoRemainder(S, q).isZERO()) {
260                break;
261            }
262            logger.info("final devision not successfull");
263            //System.out.println("P rem q = " + PolyUtil.<BigInteger>basePseudoRemainder(P,q));
264            //System.out.println("S rem q = " + PolyUtil.<BigInteger>basePseudoRemainder(S,q));
265            //break;
266        }
267        return q;
268    }
269
270
271    /**
272     * Univariate GenPolynomial recursive greatest comon divisor. Uses
273     * multivariate Hensel list.
274     * @param P univariate recursive GenPolynomial.
275     * @param S univariate recursive GenPolynomial.
276     * @return gcd(P,S).
277     */
278    @Override
279    public GenPolynomial<GenPolynomial<BigInteger>> recursiveUnivariateGcd(
280                    GenPolynomial<GenPolynomial<BigInteger>> P, GenPolynomial<GenPolynomial<BigInteger>> S) {
281        if (S == null || S.isZERO()) {
282            return P;
283        }
284        if (P == null || P.isZERO()) {
285            return S;
286        }
287        if (P.ring.nvar > 1) {
288            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
289        }
290        long e = P.degree(0);
291        long f = S.degree(0);
292        GenPolynomial<GenPolynomial<BigInteger>> q;
293        GenPolynomial<GenPolynomial<BigInteger>> r;
294        if (f > e) {
295            r = P;
296            q = S;
297            long g = f;
298            f = e;
299            e = g;
300        } else {
301            q = P;
302            r = S;
303        }
304        if (debug) {
305            logger.debug("degrees: e = " + e + ", f = " + f);
306        }
307        r = r.abs();
308        q = q.abs();
309        GenPolynomial<BigInteger> a = recursiveContent(r);
310        GenPolynomial<BigInteger> b = recursiveContent(q);
311
312        GenPolynomial<BigInteger> c = gcd(a, b); // go to recursion
313        //System.out.println("rgcd c = " + c);
314        r = PolyUtil.<BigInteger> recursiveDivide(r, a);
315        q = PolyUtil.<BigInteger> recursiveDivide(q, b);
316        a = PolyUtil.<BigInteger> basePseudoDivide(a, c);
317        b = PolyUtil.<BigInteger> basePseudoDivide(b, c);
318        if (r.isONE()) {
319            return r.multiply(c);
320        }
321        if (q.isONE()) {
322            return q.multiply(c);
323        }
324        // convert from Z[y1,...,yr][x] to Z[x][y1,...,yr] to Z[x,y1,...,yr]
325        GenPolynomial<GenPolynomial<BigInteger>> qs = PolyUtil.<BigInteger> switchVariables(q);
326        GenPolynomial<GenPolynomial<BigInteger>> rs = PolyUtil.<BigInteger> switchVariables(r);
327
328        GenPolynomialRing<GenPolynomial<BigInteger>> rfac = qs.ring;
329        RingFactory<GenPolynomial<BigInteger>> rrfac = rfac.coFac;
330        GenPolynomialRing<BigInteger> cfac = (GenPolynomialRing<BigInteger>) rrfac;
331        GenPolynomialRing<BigInteger> dfac = cfac.extend(rfac.getVars());
332        //System.out.println("pfac = " + P.ring.toScript());
333        //System.out.println("rfac = " + rfac.toScript());
334        //System.out.println("dfac = " + dfac.toScript());
335
336        GenPolynomial<BigInteger> qd = PolyUtil.<BigInteger> distribute(dfac, qs);
337        GenPolynomial<BigInteger> rd = PolyUtil.<BigInteger> distribute(dfac, rs);
338        //System.out.println("qd = " + qd);
339        //System.out.println("rd = " + rd);
340        //System.out.println("a  = " + a);
341        //System.out.println("b  = " + b);
342        //System.out.println("c  = " + c);
343
344        // compute normalization factor
345        BigInteger ac = rd.leadingBaseCoefficient();
346        BigInteger bc = qd.leadingBaseCoefficient();
347        BigInteger cc = gcd(ac, bc); // indirection
348
349        //initialize prime list
350        PrimeList primes = new PrimeList(PrimeList.Range.medium); // PrimeList.Range.medium);
351        Iterator<java.math.BigInteger> primeIter = primes.iterator();
352        int pn = 50; //primes.size();
353
354        // double check variables
355        // need qe,re,qd,rd,a,b
356        GenPolynomial<MOD> qe0;
357        GenPolynomial<MOD> re0;
358        GenPolynomial<MOD> ce0 = null;
359
360        for (int i = 0; i < 11; i++) { // meta loop
361            //System.out.println("======================================================= run " 
362            //                   + dfac.nvar + ", " + i);
363            java.math.BigInteger p = null; //new java.math.BigInteger("19"); //primes.next();
364            // 5 small, 5 medium and 1 large size primes
365            if (i == 0) { // medium size
366                primes = new PrimeList(PrimeList.Range.medium);
367                primeIter = primes.iterator();
368            }
369            if (i == 5) { // smal size
370                primes = new PrimeList(PrimeList.Range.small);
371                primeIter = primes.iterator();
372                p = primeIter.next(); // 2
373                p = primeIter.next(); // 3
374                p = primeIter.next(); // 5
375                p = primeIter.next(); // 7
376            }
377            if (i == 10) { // large size
378                primes = new PrimeList(PrimeList.Range.large);
379                primeIter = primes.iterator();
380            }
381            ModularRingFactory<MOD> cofac = null;
382            int pi = 0;
383            while (pi < pn && primeIter.hasNext()) {
384                p = primeIter.next();
385                //p = new java.math.BigInteger("19");
386                logger.info("prime = " + p);
387                // initialize coefficient factory and map normalization factor and polynomials
388                ModularRingFactory<MOD> cf = null;
389                if (ModLongRing.MAX_LONG.compareTo(p) > 0) {
390                    cf = (ModularRingFactory) new ModLongRing(p, true);
391                } else {
392                    cf = (ModularRingFactory) new ModIntegerRing(p, true);
393                }
394                MOD nf = cf.fromInteger(cc.getVal());
395                if (nf.isZERO()) {
396                    continue;
397                }
398                nf = cf.fromInteger(q.leadingBaseCoefficient().leadingBaseCoefficient().getVal());
399                if (nf.isZERO()) {
400                    continue;
401                }
402                nf = cf.fromInteger(r.leadingBaseCoefficient().leadingBaseCoefficient().getVal());
403                if (nf.isZERO()) {
404                    continue;
405                }
406                cofac = cf;
407                break;
408            }
409            if (cofac == null) { // no lucky prime found
410                System.out.println("giving up on Hensel gcd reverting to Subres gcd");
411                GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(P, S);
412                return T.abs().multiply(c); //.abs();
413            }
414            //System.out.println("cofac = " + cofac);
415
416            List<MOD> V = new ArrayList<MOD>(P.ring.nvar);
417            GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(cofac, dfac);
418            //System.out.println("mfac = " + mfac.toScript());
419            GenPolynomial<MOD> qm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, qd);
420            GenPolynomial<MOD> rm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, rd);
421            //System.out.println("qm = " + qm);
422            //System.out.println("rm = " + rm);
423
424            // search evaluation point and evaluate
425            GenPolynomialRing<MOD> ckfac = mfac;
426            GenPolynomial<MOD> qe = qm;
427            GenPolynomial<MOD> re = rm;
428            GenPolynomial<MOD> qep;
429            GenPolynomial<MOD> rep;
430            for (int j = dfac.nvar; j > 1; j--) {
431                // evaluation to univariate case
432                long degq = qe.degree(ckfac.nvar - 2);
433                long degr = re.degree(ckfac.nvar - 2);
434                ckfac = ckfac.contract(1);
435                long vi = 1L; //(long)(dfac.nvar-j); // 1L; 0 not so good for small p
436                if (p.longValue() > 1000L) {
437                    //vi = (long)j+1L;
438                    vi = 0L;
439                }
440                // search small evaluation point
441                while (true) {
442                    MOD vp = cofac.fromInteger(vi++);
443                    //System.out.println("vp = " + vp);
444                    if (vp.isZERO() && vi != 1L) { // all elements of Z_p exhausted
445                        qe = null;
446                        re = null;
447                        break;
448                        //GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(P, S);
449                        //return T.abs().multiply(c); //.abs();
450                    }
451                    qep = PolyUtil.<MOD> evaluateMain(ckfac, qe, vp);
452                    rep = PolyUtil.<MOD> evaluateMain(ckfac, re, vp);
453                    //System.out.println("qep = " + qep);
454                    //System.out.println("rep = " + rep);
455
456                    // check lucky evaluation point 
457                    MOD ql = qep.leadingBaseCoefficient();
458                    MOD rl = rep.leadingBaseCoefficient();
459                    if (ql.isZERO() || rl.isZERO()) { // nearly non sense
460                        continue;
461                    }
462                    if (degq != qep.degree(ckfac.nvar - 1)) {
463                        //System.out.println("degv(qe) = " + qe.degreeVector());
464                        //System.out.println("deg(qe) = " + degq + ", deg(qe) = " + qep.degree(ckfac.nvar-1));
465                        continue;
466                    }
467                    if (degr != rep.degree(ckfac.nvar - 1)) {
468                        //System.out.println("degv(re) = " + re.degreeVector());
469                        //System.out.println("deg(re) = " + degr + ", deg(re) = " + rep.degree(ckfac.nvar-1));
470                        continue;
471                    }
472                    V.add(vp);
473                    qe = qep;
474                    re = rep;
475                    break;
476                }
477                if (qe == null && re == null) {
478                    break;
479                }
480            }
481            if (qe == null && re == null) {
482                continue;
483            }
484            logger.info("evaluation points  = " + V);
485            //System.out.println("qe = " + qe);
486            //System.out.println("re = " + re);
487
488            // recursion base:
489            GreatestCommonDivisorAbstract<MOD> mufd = GCDFactory.getImplementation(cofac);
490            GenPolynomial<MOD> ce = mufd.baseGcd(qe, re);
491            if (ce.isConstant()) {
492                return P.ring.getONE().multiply(c);
493            }
494            //System.out.println("ce = " + ce);
495            logger.info("base gcd = " + ce);
496            //System.out.println("c = " + c);
497
498            // double check 
499            // need qe,re,qd,rd,a,b
500            if (i == 0) {
501                qe0 = qe;
502                re0 = re;
503                ce0 = ce;
504                continue;
505            }
506            long d0 = ce0.degree(0);
507            long d1 = ce.degree(0);
508            //System.out.println("d0, d1 = " + d0 + ", " + d1);
509            if (d1 < d0) {
510                qe0 = qe;
511                re0 = re;
512                ce0 = ce;
513                continue;
514            } else if (d1 > d0) {
515                continue;
516            }
517            // d0 == d1 is ok
518            long dx = r.degree(0);
519            //System.out.println("d0, dx = " + d0 + ", " + dx);
520            if (d0 == dx) { // gcd == r ??
521                if (PolyUtil.<BigInteger> recursivePseudoRemainder(q, r).isZERO()) {
522                    r = r.abs().multiply(c); //.abs();
523                    logger.info("exit with r | q : " + r);
524                    return r;
525                }
526                continue;
527            }
528
529            // norm
530            BigInteger mn = null;
531            //mn = mn.multiply(cc).multiply(mn.fromInteger(2));
532
533            // prepare lifting
534            GenPolynomial<MOD> re1 = PolyUtil.<MOD> basePseudoDivide(re, ce);
535            GenPolynomial<MOD> qe1 = PolyUtil.<MOD> basePseudoDivide(qe, ce);
536            GenPolynomial<BigInteger> ui;
537            GenPolynomial<MOD> he;
538            GenPolynomial<BigInteger> g;
539            if (mufd.baseGcd(re1, ce).isONE()) {
540                ui = rd;
541                he = re1;
542                BigInteger an = rd.maxNorm();
543                mn = an.multiply(cc).multiply(new BigInteger(2L));
544                g = a;
545            } else if (mufd.baseGcd(qe1, ce).isONE()) {
546                ui = qd;
547                he = qe1;
548                BigInteger bn = qd.maxNorm();
549                mn = bn.multiply(cc).multiply(new BigInteger(2L));
550                g = b;
551            } else {
552                break;
553                //System.out.println("giving up on Hensel gcd reverting to Subres gcd");
554            }
555            //System.out.println("ui = " + ui);
556
557            long k = Power.logarithm(new BigInteger(p), mn);
558            //System.out.println("mn = " + mn);
559            //System.out.println("k = " + k);
560
561            List<GenPolynomial<MOD>> F = new ArrayList<GenPolynomial<MOD>>(2);
562            F.add(ce);
563            F.add(he);
564            List<GenPolynomial<BigInteger>> G = new ArrayList<GenPolynomial<BigInteger>>(2);
565            G.add(g.ring.getONE()); // TODO
566            G.add(g.ring.getONE());
567            List<GenPolynomial<MOD>> lift;
568            try {
569                lift = HenselMultUtil.<MOD> liftHenselFull(ui, F, V, k, G);
570                logger.info("lift = " + lift);
571            } catch (NoLiftingException nle) {
572                //System.out.println("exception : " + nle);
573                continue;
574            } catch (ArithmeticException ae) {
575                //System.out.println("exception : " + ae);
576                continue;
577            }
578
579            // convert Ci from Z_{p^k}[x,y1,...,yr] to Z[x,y1,...,yr] to Z[x][y1,...,yr] to Z[y1,...,yr][x]
580            GenPolynomial<BigInteger> ci = PolyUtil.integerFromModularCoefficients(dfac, lift.get(0));
581            GenPolynomial<GenPolynomial<BigInteger>> Cr = PolyUtil.<BigInteger> recursive(rfac, ci);
582            GenPolynomial<GenPolynomial<BigInteger>> Cs = PolyUtil.<BigInteger> switchVariables(Cr);
583            if (!Cs.ring.equals(P.ring)) {
584                System.out.println("Cs.ring = " + Cs.ring + ", P.ring = " + P.ring);
585            }
586            q = recursivePrimitivePart(Cs);
587            q = q.abs().multiply(c); //.abs();
588            if (PolyUtil.<BigInteger> recursivePseudoRemainder(P, q).isZERO()
589                            && PolyUtil.<BigInteger> recursivePseudoRemainder(S, q).isZERO()) {
590                return q;
591            }
592            //System.out.println("bad q = " + q);
593            //System.out.println("bad P/q = " + PolyUtil.<BigInteger> recursivePseudoRemainder(P,q));
594            //System.out.println("bad S/q = " + PolyUtil.<BigInteger> recursivePseudoRemainder(S,q));
595        } // end for meta loop
596          // Hensel gcd failed
597        GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(P, S);
598        T = T.abs().multiply(c); //.abs();
599        logger.info("giving up on Hensel gcd reverting to Subres gcd " + T);
600        //System.out.println("giving up on Hensel gcd reverted to Subres gcd: " + T);
601        return T;
602    }
603
604}