Source code

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