Source code

001/*
002 * $Id: GreatestCommonDivisorModular.java 4965 2014-10-17 20:07:51Z kredel $
003 */
004
005package edu.jas.ufd;
006
007
008import org.apache.log4j.Logger;
009
010import edu.jas.arith.BigInteger;
011import edu.jas.arith.Combinatoric;
012import edu.jas.arith.ModIntegerRing;
013import edu.jas.arith.ModLongRing;
014import edu.jas.arith.Modular;
015import edu.jas.arith.ModularRingFactory;
016import edu.jas.arith.PrimeList;
017import edu.jas.poly.ExpVector;
018import edu.jas.poly.GenPolynomial;
019import edu.jas.poly.GenPolynomialRing;
020import edu.jas.poly.PolyUtil;
021import edu.jas.structure.GcdRingElem;
022import edu.jas.structure.Power;
023
024
025/**
026 * Greatest common divisor algorithms with modular computation and chinese
027 * remainder algorithm.
028 * @author Heinz Kredel
029 */
030
031public class GreatestCommonDivisorModular<MOD extends GcdRingElem<MOD> & Modular> extends
032                GreatestCommonDivisorAbstract<BigInteger> {
033
034
035    private static final Logger logger = Logger.getLogger(GreatestCommonDivisorModular.class);
036
037
038    private final boolean debug = logger.isDebugEnabled(); //logger.isInfoEnabled();
039
040
041    /*
042     * Modular gcd algorithm to use.
043     */
044    protected final GreatestCommonDivisorAbstract<MOD> mufd;
045
046
047    /*
048     * Integer gcd algorithm for fall back.
049     */
050    protected final GreatestCommonDivisorAbstract<BigInteger> iufd = new GreatestCommonDivisorSubres<BigInteger>();
051
052
053    /**
054     * Constructor to set recursive algorithm. Use modular evaluation GCD
055     * algorithm.
056     */
057    public GreatestCommonDivisorModular() {
058        this(false);
059    }
060
061
062    /**
063     * Constructor to set recursive algorithm.
064     * @param simple , true if the simple PRS should be used.
065     */
066    public GreatestCommonDivisorModular(boolean simple) {
067        if (simple) {
068            mufd = new GreatestCommonDivisorSimple<MOD>();
069        } else {
070            mufd = new GreatestCommonDivisorModEval<MOD>();
071        }
072    }
073
074
075    /**
076     * Univariate GenPolynomial greatest comon divisor. Delegate to subresultant
077     * baseGcd, should not be needed.
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        return iufd.baseGcd(P, S);
085    }
086
087
088    /**
089     * Univariate GenPolynomial recursive greatest comon divisor. Delegate to
090     * subresultant recursiveGcd, should not be needed.
091     * @param P univariate recursive GenPolynomial.
092     * @param S univariate recursive GenPolynomial.
093     * @return gcd(P,S).
094     */
095    @Override
096    public GenPolynomial<GenPolynomial<BigInteger>> recursiveUnivariateGcd(
097                    GenPolynomial<GenPolynomial<BigInteger>> P, GenPolynomial<GenPolynomial<BigInteger>> S) {
098        return iufd.recursiveUnivariateGcd(P, S);
099    }
100
101
102    /**
103     * GenPolynomial greatest comon divisor, modular algorithm.
104     * @param P GenPolynomial.
105     * @param S GenPolynomial.
106     * @return gcd(P,S).
107     */
108    @Override
109    @SuppressWarnings("unchecked")
110    public GenPolynomial<BigInteger> gcd(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) {
111        if (S == null || S.isZERO()) {
112            return P;
113        }
114        if (P == null || P.isZERO()) {
115            return S;
116        }
117        GenPolynomialRing<BigInteger> fac = P.ring;
118        // special case for univariate polynomials
119        if (fac.nvar <= 1) {
120            GenPolynomial<BigInteger> T = baseGcd(P, S);
121            return T;
122        }
123        long e = P.degree(0);
124        long f = S.degree(0);
125        GenPolynomial<BigInteger> q;
126        GenPolynomial<BigInteger> r;
127        if (f > e) {
128            r = P;
129            q = S;
130            long g = f;
131            f = e;
132            e = g;
133        } else {
134            q = P;
135            r = S;
136        }
137        if (debug) {
138            logger.debug("degrees: e = " + e + ", f = " + f);
139        }
140        r = r.abs();
141        q = q.abs();
142        // compute contents and primitive parts
143        BigInteger a = baseContent(r);
144        BigInteger b = baseContent(q);
145        // gcd of coefficient contents
146        BigInteger c = gcd(a, b); // indirection
147        r = divide(r, a); // indirection
148        q = divide(q, b); // indirection
149        if (r.isONE()) {
150            return r.multiply(c);
151        }
152        if (q.isONE()) {
153            return q.multiply(c);
154        }
155        // compute normalization factor
156        BigInteger ac = r.leadingBaseCoefficient();
157        BigInteger bc = q.leadingBaseCoefficient();
158        BigInteger cc = gcd(ac, bc); // indirection
159        // compute norms
160        BigInteger an = r.maxNorm();
161        BigInteger bn = q.maxNorm();
162        BigInteger n = (an.compareTo(bn) < 0 ? bn : an);
163        n = n.multiply(cc).multiply(n.fromInteger(2));
164        // compute degree vectors
165        ExpVector rdegv = r.degreeVector();
166        ExpVector qdegv = q.degreeVector();
167        //compute factor coefficient bounds
168        BigInteger af = an.multiply(PolyUtil.factorBound(rdegv));
169        BigInteger bf = bn.multiply(PolyUtil.factorBound(qdegv));
170        BigInteger cf = (af.compareTo(bf) < 0 ? bf : af);
171        cf = cf.multiply(cc.multiply(cc.fromInteger(8)));
172        //initialize prime list and degree vector
173        PrimeList primes = new PrimeList();
174        int pn = 10; //primes.size();
175        ExpVector wdegv = rdegv.subst(0, rdegv.getVal(0) + 1);
176        // +1 seems to be a hack for the unlucky prime test
177        ModularRingFactory<MOD> cofac;
178        ModularRingFactory<MOD> cofacM = null;
179        GenPolynomial<MOD> qm;
180        GenPolynomial<MOD> rm;
181        GenPolynomialRing<MOD> mfac;
182        GenPolynomialRing<MOD> rfac = null;
183        int i = 0;
184        BigInteger M = null;
185        BigInteger cfe = null;
186        GenPolynomial<MOD> cp = null;
187        GenPolynomial<MOD> cm = null;
188        GenPolynomial<BigInteger> cpi = null;
189        if (debug) {
190            logger.debug("c = " + c);
191            logger.debug("cc = " + cc);
192            logger.debug("n  = " + n);
193            logger.debug("cf = " + cf);
194            logger.info("wdegv = " + wdegv);
195        }
196        for (java.math.BigInteger p : primes) {
197            //System.out.println("next run ++++++++++++++++++++++++++++++++++");
198            if (p.longValue() == 2L) { // skip 2
199                continue;
200            }
201            if (++i >= pn) {
202                logger.warn("prime list exhausted, pn = " + pn);
203                return iufd.gcd(P, S);
204                //throw new ArithmeticException("prime list exhausted");
205            }
206            // initialize coefficient factory and map normalization factor
207            if (ModLongRing.MAX_LONG.compareTo(p) > 0) {
208                cofac = (ModularRingFactory) new ModLongRing(p, true);
209            } else {
210                cofac = (ModularRingFactory) new ModIntegerRing(p, true);
211            }
212            MOD nf = cofac.fromInteger(cc.getVal());
213            if (nf.isZERO()) {
214                continue;
215            }
216            // initialize polynomial factory and map polynomials
217            mfac = new GenPolynomialRing<MOD>(cofac, fac.nvar, fac.tord, fac.getVars());
218            qm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, q);
219            if (qm.isZERO() || !qm.degreeVector().equals(qdegv)) {
220                continue;
221            }
222            rm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, r);
223            if (rm.isZERO() || !rm.degreeVector().equals(rdegv)) {
224                continue;
225            }
226            if (debug) {
227                logger.info("cofac = " + cofac.getIntegerModul());
228            }
229            // compute modular gcd
230            cm = mufd.gcd(rm, qm);
231            // test for constant g.c.d
232            if (cm.isConstant()) {
233                logger.debug("cm, constant = " + cm);
234                return fac.getONE().multiply(c);
235                //return cm.abs().multiply( c );
236            }
237            // test for unlucky prime
238            ExpVector mdegv = cm.degreeVector();
239            if (wdegv.equals(mdegv)) { // TL = 0
240                // prime ok, next round
241                if (M != null) {
242                    if (M.compareTo(cfe) > 0) {
243                        System.out.println("M > cfe: " + M + " > " + cfe);
244                        // continue; // why should this be required?
245                    }
246                }
247            } else { // TL = 3
248                boolean ok = false;
249                if (wdegv.multipleOf(mdegv)) { // TL = 2 // EVMT(wdegv,mdegv)
250                    M = null; // init chinese remainder
251                    ok = true; // prime ok
252                }
253                if (mdegv.multipleOf(wdegv)) { // TL = 1 // EVMT(mdegv,wdegv)
254                    continue; // skip this prime
255                }
256                if (!ok) {
257                    M = null; // discard chinese remainder and previous work
258                    continue; // prime not ok
259                }
260            }
261            //--wdegv = mdegv;
262            // prepare chinese remainder algorithm
263            cm = cm.multiply(nf);
264            if (M == null) {
265                // initialize chinese remainder algorithm
266                M = new BigInteger(p);
267                cofacM = cofac;
268                rfac = mfac;
269                cp = cm;
270                wdegv = wdegv.gcd(mdegv); //EVGCD(wdegv,mdegv);
271                cfe = cf;
272                for (int k = 0; k < wdegv.length(); k++) {
273                    cfe = cfe.multiply(new BigInteger(wdegv.getVal(k) + 1));
274                }
275            } else {
276                // apply chinese remainder algorithm
277                BigInteger Mp = M;
278                MOD mi = cofac.fromInteger(Mp.getVal());
279                mi = mi.inverse(); // mod p
280                M = M.multiply(new BigInteger(p));
281                if (ModLongRing.MAX_LONG.compareTo(M.getVal()) > 0) {
282                    cofacM = (ModularRingFactory) new ModLongRing(M.getVal());
283                } else {
284                    cofacM = (ModularRingFactory) new ModIntegerRing(M.getVal());
285                }
286                rfac = new GenPolynomialRing<MOD>(cofacM, fac);
287                if (!cofac.getClass().equals(cofacM.getClass())) {
288                    logger.info("adjusting coefficents: cofacM = " + cofacM.getClass() + ", cofacP = "
289                                    + cofac.getClass());
290                    cofac = (ModularRingFactory) new ModIntegerRing(p);
291                    mfac = new GenPolynomialRing<MOD>(cofac, fac);
292                    GenPolynomial<BigInteger> mm = PolyUtil.<MOD> integerFromModularCoefficients(fac, cm);
293                    cm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, mm);
294                    mi = cofac.fromInteger(Mp.getVal());
295                    mi = mi.inverse(); // mod p
296                }
297                if (!cp.ring.coFac.getClass().equals(cofacM.getClass())) {
298                    logger.info("adjusting coefficents: cofacM = " + cofacM.getClass() + ", cofacM' = "
299                                    + cp.ring.coFac.getClass());
300                    ModularRingFactory cop = (ModularRingFactory) cp.ring.coFac;
301                    cofac = (ModularRingFactory) new ModIntegerRing(cop.getIntegerModul().getVal());
302                    mfac = new GenPolynomialRing<MOD>(cofac, fac);
303                    GenPolynomial<BigInteger> mm = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp);
304                    cp = PolyUtil.<MOD> fromIntegerCoefficients(mfac, mm);
305                }
306                cp = PolyUtil.<MOD> chineseRemainder(rfac, cp, mi, cm);
307            }
308            // test for completion
309            if (n.compareTo(M) <= 0) {
310                break;
311            }
312            // must use integer.sumNorm
313            cpi = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp);
314            BigInteger cmn = cpi.sumNorm();
315            cmn = cmn.multiply(cmn.fromInteger(4));
316            //if ( cmn.compareTo( M ) <= 0 ) {
317            // does not work: only if cofactors are also considered?
318            // break;
319            //}
320            if (i % 2 != 0 && !cp.isZERO()) {
321                // check if done on every second prime
322                GenPolynomial<BigInteger> x;
323                x = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp);
324                x = basePrimitivePart(x);
325                if (!PolyUtil.<BigInteger> baseSparsePseudoRemainder(q, x).isZERO()) {
326                    continue;
327                }
328                if (!PolyUtil.<BigInteger> baseSparsePseudoRemainder(r, x).isZERO()) {
329                    continue;
330                }
331                logger.info("done on exact division, #primes = " + i);
332                break;
333            }
334        }
335        if (debug) {
336            logger.info("done on M = " + M + ", #primes = " + i);
337        }
338        // remove normalization
339        q = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp);
340        q = basePrimitivePart(q);
341        return q.abs().multiply(c);
342    }
343
344
345    /**
346     * Univariate GenPolynomial resultant.
347     * @param P univariate GenPolynomial.
348     * @param S univariate GenPolynomial.
349     * @return res(P,S).
350     */
351    @Override
352    public GenPolynomial<BigInteger> baseResultant(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) {
353        // not a special case here
354        return resultant(P, S);
355    }
356
357
358    /**
359     * Univariate GenPolynomial recursive resultant.
360     * @param P univariate recursive GenPolynomial.
361     * @param S univariate recursive GenPolynomial.
362     * @return res(P,S).
363     */
364    @Override
365    public GenPolynomial<GenPolynomial<BigInteger>> recursiveUnivariateResultant(
366                    GenPolynomial<GenPolynomial<BigInteger>> P, GenPolynomial<GenPolynomial<BigInteger>> S) {
367        // only in this class
368        return recursiveResultant(P, S);
369    }
370
371
372    /**
373     * GenPolynomial resultant, modular algorithm.
374     * @param P GenPolynomial.
375     * @param S GenPolynomial.
376     * @return res(P,S).
377     */
378    @Override
379    @SuppressWarnings("unchecked")
380    public GenPolynomial<BigInteger> resultant(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) {
381        if (S == null || S.isZERO()) {
382            return S;
383        }
384        if (P == null || P.isZERO()) {
385            return P;
386        }
387        GenPolynomialRing<BigInteger> fac = P.ring;
388        // no special case for univariate polynomials in this class !
389        //if (fac.nvar <= 1) {
390        //    GenPolynomial<BigInteger> T = iufd.baseResultant(P, S);
391        //    return T;
392        //}
393        long e = P.degree(0);
394        long f = S.degree(0);
395        GenPolynomial<BigInteger> q;
396        GenPolynomial<BigInteger> r;
397        if (f > e) {
398            r = P;
399            q = S;
400            long g = f;
401            f = e;
402            e = g;
403        } else {
404            q = P;
405            r = S;
406        }
407        // compute norms
408        BigInteger an = r.maxNorm();
409        BigInteger bn = q.maxNorm();
410        an = Power.<BigInteger> power(fac.coFac, an, f);
411        bn = Power.<BigInteger> power(fac.coFac, bn, e);
412        BigInteger cn = Combinatoric.factorial(e + f);
413        BigInteger n = cn.multiply(an).multiply(bn);
414
415        // compute degree vectors
416        ExpVector rdegv = r.leadingExpVector(); //degreeVector();
417        ExpVector qdegv = q.leadingExpVector(); //degreeVector();
418
419        //initialize prime list and degree vector
420        PrimeList primes = new PrimeList();
421        int pn = 30; //primes.size();
422        ModularRingFactory<MOD> cofac;
423        ModularRingFactory<MOD> cofacM = null;
424        GenPolynomial<MOD> qm;
425        GenPolynomial<MOD> rm;
426        GenPolynomialRing<MOD> mfac;
427        GenPolynomialRing<MOD> rfac = null;
428        int i = 0;
429        BigInteger M = null;
430        GenPolynomial<MOD> cp = null;
431        GenPolynomial<MOD> cm = null;
432        //GenPolynomial<BigInteger> cpi = null;
433        if (debug) {
434            logger.debug("an  = " + an);
435            logger.debug("bn  = " + bn);
436            logger.debug("e+f = " + (e + f));
437            logger.debug("cn  = " + cn);
438            logger.info("n     = " + n);
439            //logger.info("q     = " + q);
440            //logger.info("r     = " + r);
441            //logger.info("rdegv = " + rdegv.toString(fac.getVars()));
442            //logger.info("qdegv = " + qdegv.toString(fac.getVars()));
443        }
444        for (java.math.BigInteger p : primes) {
445            //System.out.println("next run ++++++++++++++++++++++++++++++++++");
446            if (p.longValue() == 2L) { // skip 2
447                continue;
448            }
449            if (++i >= pn) {
450                logger.warn("prime list exhausted, pn = " + pn);
451                return iufd.resultant(P, S);
452                //throw new ArithmeticException("prime list exhausted");
453            }
454            // initialize coefficient factory and map normalization factor
455            if (ModLongRing.MAX_LONG.compareTo(p) > 0) {
456                cofac = (ModularRingFactory) new ModLongRing(p, true);
457            } else {
458                cofac = (ModularRingFactory) new ModIntegerRing(p, true);
459            }
460            // initialize polynomial factory and map polynomials
461            mfac = new GenPolynomialRing<MOD>(cofac, fac);
462            qm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, q);
463            if (qm.isZERO() || !qm.leadingExpVector().equals(qdegv)) { //degreeVector()
464                //logger.info("qm = " + qm);
465                if (debug) {
466                    logger.info("unlucky prime = " + cofac.getIntegerModul() + ", degv = "
467                                    + qm.leadingExpVector());
468                }
469                continue;
470            }
471            rm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, r);
472            if (rm.isZERO() || !rm.leadingExpVector().equals(rdegv)) { //degreeVector()
473                //logger.info("rm = " + rm);
474                if (debug) {
475                    logger.info("unlucky prime = " + cofac.getIntegerModul() + ", degv = "
476                                    + rm.leadingExpVector());
477                }
478                continue;
479            }
480            logger.info("lucky prime = " + cofac.getIntegerModul());
481
482            // compute modular resultant
483            cm = mufd.resultant(qm, rm);
484            if (debug) {
485                logger.info("res_p = " + cm);
486            }
487
488            // prepare chinese remainder algorithm
489            if (M == null) {
490                // initialize chinese remainder algorithm
491                M = new BigInteger(p);
492                cofacM = cofac;
493                //rfac = mfac;
494                cp = cm;
495            } else {
496                // apply chinese remainder algorithm
497                BigInteger Mp = M;
498                MOD mi = cofac.fromInteger(Mp.getVal());
499                mi = mi.inverse(); // mod p
500                M = M.multiply(new BigInteger(p));
501                if (ModLongRing.MAX_LONG.compareTo(M.getVal()) > 0) {
502                    cofacM = (ModularRingFactory) new ModLongRing(M.getVal());
503                } else {
504                    cofacM = (ModularRingFactory) new ModIntegerRing(M.getVal());
505                }
506                rfac = new GenPolynomialRing<MOD>(cofacM, fac);
507                if (!cofac.getClass().equals(cofacM.getClass())) {
508                    logger.info("adjusting coefficents: cofacM = " + cofacM.getClass() + ", cofacP = "
509                                    + cofac.getClass());
510                    cofac = (ModularRingFactory) new ModIntegerRing(p);
511                    mfac = new GenPolynomialRing<MOD>(cofac, fac);
512                    GenPolynomial<BigInteger> mm = PolyUtil.<MOD> integerFromModularCoefficients(fac, cm);
513                    cm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, mm);
514                    mi = cofac.fromInteger(Mp.getVal());
515                    mi = mi.inverse(); // mod p
516                }
517                if (!cp.ring.coFac.getClass().equals(cofacM.getClass())) {
518                    logger.info("adjusting coefficents: cofacM = " + cofacM.getClass() + ", cofacM' = "
519                                    + cp.ring.coFac.getClass());
520                    ModularRingFactory cop = (ModularRingFactory) cp.ring.coFac;
521                    cofac = (ModularRingFactory) new ModIntegerRing(cop.getIntegerModul().getVal());
522                    mfac = new GenPolynomialRing<MOD>(cofac, fac);
523                    GenPolynomial<BigInteger> mm = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp);
524                    cp = PolyUtil.<MOD> fromIntegerCoefficients(mfac, mm);
525                }
526                cp = PolyUtil.<MOD> chineseRemainder(rfac, cp, mi, cm);
527            }
528            // test for completion
529            if (n.compareTo(M) <= 0) {
530                break;
531            }
532        }
533        if (debug) {
534            logger.info("done on M = " + M + ", #primes = " + i);
535        }
536        // convert to integer polynomial
537        q = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp);
538        return q;
539    }
540
541}