Source code

001/*
002 * $Id: GreatestCommonDivisorModEval.java 4025 2012-07-23 16:41:43Z kredel $
003 */
004
005package edu.jas.ufd;
006
007
008import org.apache.log4j.Logger;
009
010import edu.jas.arith.Modular;
011import edu.jas.arith.ModularRingFactory;
012import edu.jas.poly.ExpVector;
013import edu.jas.poly.GenPolynomial;
014import edu.jas.poly.GenPolynomialRing;
015import edu.jas.poly.PolyUtil;
016import edu.jas.structure.GcdRingElem;
017
018
019/**
020 * Greatest common divisor algorithms with modular evaluation algorithm for
021 * recursion.
022 * @author Heinz Kredel
023 */
024
025public class GreatestCommonDivisorModEval <MOD extends GcdRingElem<MOD> & Modular> 
026        extends GreatestCommonDivisorAbstract<MOD> {
027
028
029    private static final Logger logger = Logger.getLogger(GreatestCommonDivisorModEval.class);
030
031
032    private final boolean debug = logger.isDebugEnabled();
033
034
035    /**
036     * Modular gcd algorithm to use.
037     */
038    protected final GreatestCommonDivisorAbstract<MOD> mufd 
039       = new GreatestCommonDivisorSimple<MOD>();
040    // = new GreatestCommonDivisorPrimitive<MOD>();
041    // not okay: = new GreatestCommonDivisorSubres<MOD>();
042
043
044    /**
045     * Univariate GenPolynomial greatest common divisor. 
046     * @param P univariate GenPolynomial.
047     * @param S univariate GenPolynomial.
048     * @return gcd(P,S).
049     */
050    @Override
051    public GenPolynomial<MOD> baseGcd(GenPolynomial<MOD> P, GenPolynomial<MOD> S) {
052        // required as recursion base
053        return mufd.baseGcd(P, S);
054    }
055
056
057    /**
058     * Recursive univariate GenPolynomial greatest common divisor. 
059     * @param P univariate recursive GenPolynomial.
060     * @param S univariate recursive GenPolynomial.
061     * @return gcd(P,S).
062     */
063    @Override
064    public GenPolynomial<GenPolynomial<MOD>> recursiveUnivariateGcd(
065            GenPolynomial<GenPolynomial<MOD>> P, GenPolynomial<GenPolynomial<MOD>> S) {
066        return mufd.recursiveUnivariateGcd(P, S);
067    }
068
069
070    /**
071     * GenPolynomial greatest common divisor, modular evaluation algorithm.
072     * @param P GenPolynomial.
073     * @param S GenPolynomial.
074     * @return gcd(P,S).
075     */
076    @Override
077    public GenPolynomial<MOD> gcd(GenPolynomial<MOD> P, GenPolynomial<MOD> S) {
078        if (S == null || S.isZERO()) {
079            return P;
080        }
081        if (P == null || P.isZERO()) {
082            return S;
083        }
084        GenPolynomialRing<MOD> fac = P.ring;
085        // recusion base case for univariate polynomials
086        if (fac.nvar <= 1) {
087            GenPolynomial<MOD> T = baseGcd(P, S);
088            return T;
089        }
090        long e = P.degree(fac.nvar-1);
091        long f = S.degree(fac.nvar-1);
092        if ( e == 0 && f == 0 ) {
093            GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(1);
094            GenPolynomial<MOD> Pc = PolyUtil.<MOD> recursive(rfac, P).leadingBaseCoefficient();
095            GenPolynomial<MOD> Sc = PolyUtil.<MOD> recursive(rfac, S).leadingBaseCoefficient();
096            GenPolynomial<MOD> r = gcd(Pc,Sc);
097            return r.extend(fac,0,0L);
098        }
099        GenPolynomial<MOD> q;
100        GenPolynomial<MOD> r;
101        if (f > e) {
102            r = P;
103            q = S;
104            long g = f;
105            f = e;
106            e = g;
107        } else {
108            q = P;
109            r = S;
110        }
111        if (debug) {
112            logger.debug("degrees: e = " + e + ", f = " + f);
113        }
114        r = r.abs();
115        q = q.abs();
116        // setup factories
117        ModularRingFactory<MOD> cofac = (ModularRingFactory<MOD>) P.ring.coFac;
118        if (!cofac.isField()) {
119            logger.warn("cofac is not a field: " + cofac);
120        }
121        GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(fac.nvar - 1);
122        GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(cofac, rfac);
123        GenPolynomialRing<MOD> ufac = (GenPolynomialRing<MOD>) rfac.coFac;
124        //GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(cofac, fac.nvar - 1, fac.tord);
125        //GenPolynomialRing<MOD> ufac = new GenPolynomialRing<MOD>(cofac, 1, fac.tord);
126        //GenPolynomialRing<GenPolynomial<MOD>> rfac = new GenPolynomialRing<GenPolynomial<MOD>>(ufac, fac.nvar - 1, fac.tord);
127        // convert polynomials
128        GenPolynomial<GenPolynomial<MOD>> qr;
129        GenPolynomial<GenPolynomial<MOD>> rr;
130        qr = PolyUtil.<MOD> recursive(rfac, q);
131        rr = PolyUtil.<MOD> recursive(rfac, r);
132
133        // compute univariate contents and primitive parts
134        GenPolynomial<MOD> a = recursiveContent(rr);
135        GenPolynomial<MOD> b = recursiveContent(qr);
136        // gcd of univariate coefficient contents
137        GenPolynomial<MOD> c = gcd(a, b);
138        rr = PolyUtil.<MOD> recursiveDivide(rr, a);
139        qr = PolyUtil.<MOD> recursiveDivide(qr, b);
140        if (rr.isONE()) {
141            rr = rr.multiply(c);
142            r = PolyUtil.<MOD> distribute(fac, rr);
143            return r;
144        }
145        if (qr.isONE()) {
146            qr = qr.multiply(c);
147            q = PolyUtil.<MOD> distribute(fac, qr);
148            return q;
149        }
150        // compute normalization factor
151        GenPolynomial<MOD> ac = rr.leadingBaseCoefficient();
152        GenPolynomial<MOD> bc = qr.leadingBaseCoefficient();
153        GenPolynomial<MOD> cc = gcd(ac, bc);
154        // compute degrees and degree vectors
155        ExpVector rdegv = rr.degreeVector();
156        ExpVector qdegv = qr.degreeVector();
157        long rd0 = PolyUtil.<MOD> coeffMaxDegree(rr);
158        long qd0 = PolyUtil.<MOD> coeffMaxDegree(qr);
159        long cd0 = cc.degree(0);
160        long G = (rd0 >= qd0 ? rd0 : qd0) + cd0;
161
162        // initialize element and degree vector
163        ExpVector wdegv = rdegv.subst(0, rdegv.getVal(0) + 1);
164        // +1 seems to be a hack for the unlucky prime test
165        MOD inc = cofac.getONE();
166        long i = 0;
167        long en = cofac.getIntegerModul().longValue() - 1; // just a stopper
168        MOD end = cofac.fromInteger(en);
169        MOD mi;
170        GenPolynomial<MOD> M = null;
171        GenPolynomial<MOD> mn;
172        GenPolynomial<MOD> qm;
173        GenPolynomial<MOD> rm;
174        GenPolynomial<MOD> cm;
175        GenPolynomial<GenPolynomial<MOD>> cp = null;
176        if (debug) {
177            logger.debug("c = " + c);
178            logger.debug("cc = " + cc);
179            logger.debug("G = " + G);
180            logger.info("wdegv = " + wdegv);
181        }
182        for (MOD d = cofac.getZERO(); d.compareTo(end) <= 0; d = d.sum(inc)) {
183            if (++i >= en) {
184                logger.warn("elements of Z_p exhausted, en = " + en);
185                return mufd.gcd(P, S);
186                //throw new ArithmeticException("prime list exhausted");
187            }
188            // map normalization factor
189            MOD nf = PolyUtil.<MOD> evaluateMain(cofac, cc, d);
190            if (nf.isZERO()) {
191                continue;
192            }
193            // map polynomials
194            qm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, qr, d);
195            if (qm.isZERO() || !qm.degreeVector().equals(qdegv)) {
196                continue;
197            }
198            rm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, rr, d);
199            if (rm.isZERO() || !rm.degreeVector().equals(rdegv)) {
200                continue;
201            }
202            if (debug) {
203                logger.debug("eval d = " + d);
204            }
205            // compute modular gcd in recursion
206            cm = gcd(rm, qm);
207            //System.out.println("cm = " + cm);
208            // test for constant g.c.d
209            if (cm.isConstant()) {
210                logger.debug("cm.isConstant = " + cm + ", c = " + c);
211                if (c.ring.nvar < cm.ring.nvar) {
212                    c = c.extend(mfac, 0, 0);
213                }
214                cm = cm.abs().multiply(c);
215                q = cm.extend(fac, 0, 0);
216                logger.debug("q             = " + q + ", c = " + c);
217                return q;
218            }
219            // test for unlucky prime
220            ExpVector mdegv = cm.degreeVector();
221            if (wdegv.equals(mdegv)) { // TL = 0
222                // prime ok, next round
223                if (M != null) {
224                    if (M.degree(0) > G) {
225                        logger.info("deg(M) > G: " + M.degree(0) + " > " + G);
226                        // continue; // why should this be required?
227                    }
228                }
229            } else { // TL = 3
230                boolean ok = false;
231                if (wdegv.multipleOf(mdegv)) { // TL = 2
232                    M = null; // init chinese remainder
233                    ok = true; // prime ok
234                }
235                if (mdegv.multipleOf(wdegv)) { // TL = 1
236                    continue; // skip this prime
237                }
238                if (!ok) {
239                    M = null; // discard chinese remainder and previous work
240                    continue; // prime not ok
241                }
242            }
243            // prepare interpolation algorithm
244            cm = cm.multiply(nf);
245            if (M == null) {
246                // initialize interpolation
247                M = ufac.getONE();
248                cp = rfac.getZERO();
249                wdegv = wdegv.gcd(mdegv); //EVGCD(wdegv,mdegv);
250            }
251            // interpolate
252            mi = PolyUtil.<MOD> evaluateMain(cofac, M, d);
253            mi = mi.inverse(); // mod p
254            cp = PolyUtil.interpolate(rfac, cp, M, mi, cm, d);
255            mn = ufac.getONE().multiply(d);
256            mn = ufac.univariate(0).subtract(mn);
257            M = M.multiply(mn);
258            // test for completion
259            if (M.degree(0) > G) {
260                break;
261            }
262            //long cmn = PolyUtil.<MOD>coeffMaxDegree(cp);
263            //if ( M.degree(0) > cmn ) {
264            // does not work: only if cofactors are also considered?
265            // break;
266            //}
267        }
268        // remove normalization
269        cp = recursivePrimitivePart(cp).abs();
270        cp = cp.multiply(c);
271        q = PolyUtil.<MOD> distribute(fac, cp);
272        return q;
273    }
274
275
276    /**
277     * Univariate GenPolynomial resultant. 
278     * @param P univariate GenPolynomial.
279     * @param S univariate GenPolynomial.
280     * @return res(P,S).
281     */
282    @Override
283    public GenPolynomial<MOD> baseResultant(GenPolynomial<MOD> P, GenPolynomial<MOD> S) { 
284        // required as recursion base
285        return mufd.baseResultant(P, S);
286    }
287
288
289    /**
290     * Univariate GenPolynomial recursive resultant. 
291     * @param P univariate recursive GenPolynomial.
292     * @param S univariate recursive GenPolynomial.
293     * @return res(P,S).
294     */
295    @Override
296    public GenPolynomial<GenPolynomial<MOD>> recursiveUnivariateResultant(GenPolynomial<GenPolynomial<MOD>> P,
297            GenPolynomial<GenPolynomial<MOD>> S) { 
298        // only in this class
299        return recursiveResultant(P,S);
300    }
301
302
303    /**
304     * GenPolynomial resultant, modular evaluation algorithm.
305     * @param P GenPolynomial.
306     * @param S GenPolynomial.
307     * @return res(P,S).
308     */
309    @Override
310    public GenPolynomial<MOD> resultant(GenPolynomial<MOD> P, GenPolynomial<MOD> S) {
311        if (S == null || S.isZERO()) {
312            return S;
313        }
314        if (P == null || P.isZERO()) {
315            return P;
316        }
317        GenPolynomialRing<MOD> fac = P.ring;
318        // recusion base case for univariate polynomials
319        if (fac.nvar <= 1) {
320            GenPolynomial<MOD> T = mufd.baseResultant(P, S);
321            return T;
322        }
323        long e = P.degree(fac.nvar-1);
324        long f = S.degree(fac.nvar-1);
325        if ( e == 0 && f == 0 ) {
326            GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(1);
327            GenPolynomial<MOD> Pc = PolyUtil.<MOD> recursive(rfac, P).leadingBaseCoefficient();
328            GenPolynomial<MOD> Sc = PolyUtil.<MOD> recursive(rfac, S).leadingBaseCoefficient();
329            GenPolynomial<MOD> r = resultant(Pc,Sc);
330            return r.extend(fac,0,0L);
331        }
332        GenPolynomial<MOD> q;
333        GenPolynomial<MOD> r;
334        if (f > e) {
335            r = P;
336            q = S;
337            long g = f;
338            f = e;
339            e = g;
340        } else {
341            q = P;
342            r = S;
343        }
344        if (debug) {
345            logger.debug("degrees: e = " + e + ", f = " + f);
346        }
347        // setup factories
348        ModularRingFactory<MOD> cofac = (ModularRingFactory<MOD>) P.ring.coFac;
349        if (!cofac.isField()) {
350            logger.warn("cofac is not a field: " + cofac);
351        }
352        GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(fac.nvar - 1);
353        GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(cofac, rfac);
354        GenPolynomialRing<MOD> ufac = (GenPolynomialRing<MOD>) rfac.coFac;
355
356        // convert polynomials
357        GenPolynomial<GenPolynomial<MOD>> qr = PolyUtil.<MOD> recursive(rfac, q);
358        GenPolynomial<GenPolynomial<MOD>> rr = PolyUtil.<MOD> recursive(rfac, r);
359
360        // compute degrees and degree vectors
361        ExpVector qdegv = qr.degreeVector();
362        ExpVector rdegv = rr.degreeVector();
363
364        long qd0 = PolyUtil.<MOD> coeffMaxDegree(qr);
365        long rd0 = PolyUtil.<MOD> coeffMaxDegree(rr);
366        qd0 = ( qd0 == 0 ? 1 : qd0 );
367        rd0 = ( rd0 == 0 ? 1 : rd0 );
368        long qd1 = qr.degree(); 
369        long rd1 = rr.degree();
370        qd1 = ( qd1 == 0 ? 1 : qd1 );
371        rd1 = ( rd1 == 0 ? 1 : rd1 );
372        long G = qd0 * rd1 + rd0 * qd1 + 1;
373
374        // initialize element
375        MOD inc = cofac.getONE();
376        long i = 0;
377        long en = cofac.getIntegerModul().longValue() - 1; // just a stopper
378        MOD end = cofac.fromInteger(en);
379        MOD mi;
380        GenPolynomial<MOD> M = null;
381        GenPolynomial<MOD> mn;
382        GenPolynomial<MOD> qm;
383        GenPolynomial<MOD> rm;
384        GenPolynomial<MOD> cm;
385        GenPolynomial<GenPolynomial<MOD>> cp = null;
386        if (debug) {
387            //logger.info("qr    = " + qr + ", q = " + q);
388            //logger.info("rr    = " + rr + ", r = " + r);
389            //logger.info("qd0   = " + qd0);
390            //logger.info("rd0   = " + rd0);
391            logger.info("G     = " + G);
392            //logger.info("rdegv = " + rdegv); // + ", rr.degree(0) = " + rr.degree(0));
393            //logger.info("qdegv = " + qdegv); // + ", qr.degree(0) = " + qr.degree(0));
394        }
395        for (MOD d = cofac.getZERO(); d.compareTo(end) <= 0; d = d.sum(inc)) {
396            if (++i >= en) {
397                logger.warn("elements of Z_p exhausted, en = " + en + ", p = " + cofac.getIntegerModul());
398                return mufd.resultant(P, S);
399                //throw new ArithmeticException("prime list exhausted");
400            }
401            // map polynomials
402            qm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, qr, d);
403            //logger.info("qr(" + d + ") = " + qm + ", qr = " + qr);
404            if (qm.isZERO() || !qm.degreeVector().equals(qdegv)) {
405                if (debug) {
406                   logger.info("un-lucky evaluation point " + d + ", qm = " + qm.degreeVector() + " < " + qdegv);
407                }
408                continue;
409            }
410            rm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, rr, d);
411            //logger.info("rr(" + d + ") = " + rm + ", rr = " + rr);
412            if (rm.isZERO() || !rm.degreeVector().equals(rdegv)) {
413                if (debug) {
414                    logger.info("un-lucky evaluation point " + d + ", rm = " + rm.degreeVector() + " < " + rdegv);
415                }
416                continue;
417            }
418            // compute modular resultant in recursion
419            cm = resultant(rm, qm);
420            //System.out.println("cm = " + cm);
421
422            // prepare interpolation algorithm
423            if (M == null) {
424                // initialize interpolation
425                M = ufac.getONE();
426                cp = rfac.getZERO();
427            }
428            // interpolate
429            mi = PolyUtil.<MOD> evaluateMain(cofac, M, d);
430            mi = mi.inverse(); // mod p
431            cp = PolyUtil.interpolate(rfac, cp, M, mi, cm, d);
432            //logger.info("cp = " + cp);
433            mn = ufac.getONE().multiply(d);
434            mn = ufac.univariate(0).subtract(mn);
435            M = M.multiply(mn);
436            // test for completion
437            if (M.degree(0) > G) {
438                if (debug) {
439                    logger.info("last lucky evaluation point " + d);
440                }
441                break;
442            }
443            //logger.info("M  = " + M);
444        }
445        // distribute
446        q = PolyUtil.<MOD> distribute(fac, cp);
447        return q;
448    }
449
450}