Source code

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