Source code

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