Source code

001/*
002 * $Id: GreatestCommonDivisorModEval.java 6007 2020-03-29 13:34:49Z 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.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 
041       = new GreatestCommonDivisorSimple<MOD>();
042    // not okay = new GreatestCommonDivisorPrimitive<MOD>();
043    // not okay = new GreatestCommonDivisorSubres<MOD>();
044
045
046    /**
047     * Univariate GenPolynomial greatest common divisor. 
048     * @param P univariate GenPolynomial.
049     * @param S univariate GenPolynomial.
050     * @return gcd(P,S).
051     */
052    @Override
053    public GenPolynomial<MOD> baseGcd(GenPolynomial<MOD> P, GenPolynomial<MOD> S) {
054        // required as recursion base
055        return mufd.baseGcd(P, S);
056    }
057
058
059    /**
060     * Recursive univariate GenPolynomial greatest common divisor. 
061     * @param P univariate recursive GenPolynomial.
062     * @param S univariate recursive GenPolynomial.
063     * @return gcd(P,S).
064     */
065    @Override
066    public GenPolynomial<GenPolynomial<MOD>> recursiveUnivariateGcd(
067            GenPolynomial<GenPolynomial<MOD>> P, GenPolynomial<GenPolynomial<MOD>> S) {
068        //return mufd.recursiveUnivariateGcd(P, S);
069        // distributed polynomials gcd
070        GenPolynomialRing<GenPolynomial<MOD>> rfac = P.ring;
071        RingFactory<GenPolynomial<MOD>> rrfac = rfac.coFac;
072        GenPolynomialRing<MOD> cfac = (GenPolynomialRing<MOD>) rrfac;
073        GenPolynomialRing<MOD> dfac = cfac.extend(rfac.nvar);
074        GenPolynomial<MOD> Pd = PolyUtil.<MOD> distribute(dfac, P);
075        GenPolynomial<MOD> Sd = PolyUtil.<MOD> distribute(dfac, S);
076        GenPolynomial<MOD> Dd = gcd(Pd, Sd);
077        // convert to recursive
078        GenPolynomial<GenPolynomial<MOD>> C = PolyUtil.<MOD> recursive(rfac, Dd);
079        return C;
080    }
081
082
083    /**
084     * GenPolynomial greatest common divisor, modular evaluation algorithm.
085     * @param P GenPolynomial.
086     * @param S GenPolynomial.
087     * @return gcd(P,S).
088     */
089    @Override
090    public GenPolynomial<MOD> gcd(GenPolynomial<MOD> P, GenPolynomial<MOD> S) {
091        if (S == null || S.isZERO()) {
092            return P;
093        }
094        if (P == null || P.isZERO()) {
095            return S;
096        }
097        GenPolynomialRing<MOD> fac = P.ring;
098        // recusion base case for univariate polynomials
099        if (fac.nvar <= 1) {
100            GenPolynomial<MOD> T = baseGcd(P, S);
101            return T;
102        }
103        long e = P.degree(fac.nvar-1);
104        long f = S.degree(fac.nvar-1);
105        if ( e == 0 && f == 0 ) {
106            GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(1);
107            GenPolynomial<MOD> Pc = PolyUtil.<MOD> recursive(rfac, P).leadingBaseCoefficient();
108            GenPolynomial<MOD> Sc = PolyUtil.<MOD> recursive(rfac, S).leadingBaseCoefficient();
109            GenPolynomial<MOD> r = gcd(Pc,Sc);
110            return r.extend(fac,0,0L);
111        }
112        GenPolynomial<MOD> q;
113        GenPolynomial<MOD> r;
114        if (f > e) {
115            r = P;
116            q = S;
117            long g = f;
118            f = e;
119            e = g;
120        } else {
121            q = P;
122            r = S;
123        }
124        if (debug) {
125            logger.debug("degrees: e = " + e + ", f = " + f);
126        }
127        r = r.abs();
128        q = q.abs();
129        // setup factories
130        ModularRingFactory<MOD> cofac = (ModularRingFactory<MOD>) P.ring.coFac;
131        if (!cofac.isField()) {
132            logger.warn("cofac is not a field: " + cofac);
133        }
134        GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(fac.nvar - 1);
135        GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(cofac, rfac);
136        GenPolynomialRing<MOD> ufac = (GenPolynomialRing<MOD>) rfac.coFac;
137        // convert polynomials
138        GenPolynomial<GenPolynomial<MOD>> qr;
139        GenPolynomial<GenPolynomial<MOD>> rr;
140        qr = PolyUtil.<MOD> recursive(rfac, q);
141        rr = PolyUtil.<MOD> recursive(rfac, r);
142
143        // compute univariate contents and primitive parts
144        GenPolynomial<MOD> a = recursiveContent(rr);
145        GenPolynomial<MOD> b = recursiveContent(qr);
146        // gcd of univariate coefficient contents
147        GenPolynomial<MOD> c = gcd(a, b);
148        rr = PolyUtil.<MOD> recursiveDivide(rr, a);
149        qr = PolyUtil.<MOD> recursiveDivide(qr, b);
150        if (rr.isONE()) {
151            rr = rr.multiply(c);
152            r = PolyUtil.<MOD> distribute(fac, rr);
153            return r;
154        }
155        if (qr.isONE()) {
156            qr = qr.multiply(c);
157            q = PolyUtil.<MOD> distribute(fac, qr);
158            return q;
159        }
160        // compute normalization factor
161        GenPolynomial<MOD> ac = rr.leadingBaseCoefficient();
162        GenPolynomial<MOD> bc = qr.leadingBaseCoefficient();
163        GenPolynomial<MOD> cc = gcd(ac, bc);
164        // compute degrees and degree vectors
165        ExpVector rdegv = rr.degreeVector();
166        ExpVector qdegv = qr.degreeVector();
167        long rd0 = PolyUtil.<MOD> coeffMaxDegree(rr);
168        long qd0 = PolyUtil.<MOD> coeffMaxDegree(qr);
169        long cd0 = cc.degree(0);
170        long G = (rd0 >= qd0 ? rd0 : qd0) + cd0;
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().longValue() - 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);
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("prime list exhausted");
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            //System.out.println("cm = " + cm);
218            // test for constant g.c.d
219            if (cm.isConstant()) {
220                logger.debug("cm.isConstant = " + cm + ", c = " + c);
221                if (c.ring.nvar < cm.ring.nvar) {
222                    c = c.extend(mfac, 0, 0);
223                }
224                cm = cm.abs().multiply(c);
225                q = cm.extend(fac, 0, 0);
226                logger.debug("q             = " + q + ", c = " + c);
227                return q;
228            }
229            // test for unlucky evaluation point
230            ExpVector mdegv = cm.degreeVector();
231            if (wdegv.equals(mdegv)) { // TL = 0
232                // evaluation point ok, next round
233                if (M != null) {
234                    if (M.degree(0) > G) {
235                        logger.info("deg(M) > G: " + M.degree(0) + " > " + G);
236                        // continue; // why should this be required?
237                    }
238                }
239            } else { // TL = 3
240                boolean ok = false;
241                if (wdegv.multipleOf(mdegv)) { // TL = 2
242                    M = null; // init chinese remainder
243                    ok = true; // evaluation point ok
244                }
245                if (mdegv.multipleOf(wdegv)) { // TL = 1
246                    continue; // skip this evaluation point
247                }
248                if (!ok) {
249                    M = null; // discard chinese remainder and previous work
250                    continue; // evaluation point not ok
251                }
252            }
253            // prepare interpolation algorithm
254            cm = cm.multiply(nf);
255            if (M == null) {
256                // initialize interpolation
257                M = ufac.getONE();
258                cp = rfac.getZERO();
259                wdegv = wdegv.gcd(mdegv); //EVGCD(wdegv,mdegv);
260            }
261            // interpolate
262            mi = PolyUtil.<MOD> evaluateMain(cofac, M, d);
263            mi = mi.inverse(); // mod p
264            cp = PolyUtil.interpolate(rfac, cp, M, mi, cm, d);
265            mn = ufac.getONE().multiply(d);
266            mn = ufac.univariate(0).subtract(mn);
267            M = M.multiply(mn);
268            // test for completion
269            if (M.degree(0) > G) {
270                break;
271            }
272            //long cmn = PolyUtil.<MOD>coeffMaxDegree(cp);
273            //if ( M.degree(0) > cmn ) {
274            // does not work: only if cofactors are also considered?
275            // break;
276            //}
277        }
278        // remove normalization
279        cp = recursivePrimitivePart(cp).abs();
280        cp = cp.multiply(c);
281        q = PolyUtil.<MOD> distribute(fac, cp);
282        return q;
283    }
284
285
286    /**
287     * Univariate GenPolynomial resultant. 
288     * @param P univariate GenPolynomial.
289     * @param S univariate GenPolynomial.
290     * @return res(P,S).
291     */
292    @Override
293    public GenPolynomial<MOD> baseResultant(GenPolynomial<MOD> P, GenPolynomial<MOD> S) { 
294        // required as recursion base
295        return mufd.baseResultant(P, S);
296    }
297
298
299    /**
300     * Univariate GenPolynomial recursive resultant. 
301     * @param P univariate recursive GenPolynomial.
302     * @param S univariate recursive GenPolynomial.
303     * @return res(P,S).
304     */
305    @Override
306    public GenPolynomial<GenPolynomial<MOD>> recursiveUnivariateResultant(GenPolynomial<GenPolynomial<MOD>> P,
307            GenPolynomial<GenPolynomial<MOD>> S) { 
308        // only in this class
309        return recursiveResultant(P,S);
310    }
311
312
313    /**
314     * GenPolynomial resultant, modular evaluation algorithm.
315     * @param P GenPolynomial.
316     * @param S GenPolynomial.
317     * @return res(P,S).
318     */
319    @Override
320    public GenPolynomial<MOD> resultant(GenPolynomial<MOD> P, GenPolynomial<MOD> S) {
321        if (S == null || S.isZERO()) {
322            return S;
323        }
324        if (P == null || P.isZERO()) {
325            return P;
326        }
327        GenPolynomialRing<MOD> fac = P.ring;
328        // recusion base case for univariate polynomials
329        if (fac.nvar <= 1) {
330            GenPolynomial<MOD> T = baseResultant(P, S);
331            return T;
332        }
333        long e = P.degree(fac.nvar-1);
334        long f = S.degree(fac.nvar-1);
335        if ( e == 0 && f == 0 ) {
336            GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(1);
337            GenPolynomial<MOD> Pc = PolyUtil.<MOD> recursive(rfac, P).leadingBaseCoefficient();
338            GenPolynomial<MOD> Sc = PolyUtil.<MOD> recursive(rfac, S).leadingBaseCoefficient();
339            GenPolynomial<MOD> r = resultant(Pc,Sc);
340            return r.extend(fac,0,0L);
341        }
342        GenPolynomial<MOD> q;
343        GenPolynomial<MOD> r;
344        if (f > e) {
345            r = P;
346            q = S;
347            long g = f;
348            f = e;
349            e = g;
350        } else {
351            q = P;
352            r = S;
353        }
354        if (debug) {
355            logger.debug("degrees: e = " + e + ", f = " + f);
356        }
357        // setup factories
358        ModularRingFactory<MOD> cofac = (ModularRingFactory<MOD>) P.ring.coFac;
359        if (!cofac.isField()) {
360            logger.warn("cofac is not a field: " + cofac);
361        }
362        GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(fac.nvar - 1);
363        GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(cofac, rfac);
364        GenPolynomialRing<MOD> ufac = (GenPolynomialRing<MOD>) rfac.coFac;
365
366        // convert polynomials
367        GenPolynomial<GenPolynomial<MOD>> qr = PolyUtil.<MOD> recursive(rfac, q);
368        GenPolynomial<GenPolynomial<MOD>> rr = PolyUtil.<MOD> recursive(rfac, r);
369
370        // compute degrees and degree vectors
371        ExpVector qdegv = qr.degreeVector();
372        ExpVector rdegv = rr.degreeVector();
373
374        long qd0 = PolyUtil.<MOD> coeffMaxDegree(qr);
375        long rd0 = PolyUtil.<MOD> coeffMaxDegree(rr);
376        qd0 = ( qd0 == 0 ? 1 : qd0 );
377        rd0 = ( rd0 == 0 ? 1 : rd0 );
378        long qd1 = qr.degree(); 
379        long rd1 = rr.degree();
380        qd1 = ( qd1 == 0 ? 1 : qd1 );
381        rd1 = ( rd1 == 0 ? 1 : rd1 );
382        long G = qd0 * rd1 + rd0 * qd1 + 1;
383
384        // initialize element
385        MOD inc = cofac.getONE();
386        long i = 0;
387        long en = cofac.getIntegerModul().longValue() - 1; // just a stopper
388        MOD end = cofac.fromInteger(en);
389        MOD mi;
390        GenPolynomial<MOD> M = null;
391        GenPolynomial<MOD> mn;
392        GenPolynomial<MOD> qm;
393        GenPolynomial<MOD> rm;
394        GenPolynomial<MOD> cm;
395        GenPolynomial<GenPolynomial<MOD>> cp = null;
396        if (debug) {
397            //logger.info("qr    = " + qr + ", q = " + q);
398            //logger.info("rr    = " + rr + ", r = " + r);
399            //logger.info("qd0   = " + qd0);
400            //logger.info("rd0   = " + rd0);
401            logger.info("G     = " + G);
402            //logger.info("rdegv = " + rdegv); // + ", rr.degree(0) = " + rr.degree(0));
403            //logger.info("qdegv = " + qdegv); // + ", qr.degree(0) = " + qr.degree(0));
404        }
405        for (MOD d = cofac.getZERO(); d.compareTo(end) <= 0; d = d.sum(inc)) {
406            if (++i >= en) {
407                logger.warn("elements of Z_p exhausted, en = " + en + ", p = " + cofac.getIntegerModul());
408                return mufd.resultant(P, S);
409                //throw new ArithmeticException("prime list exhausted");
410            }
411            // map polynomials
412            qm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, qr, d);
413            //logger.info("qr(" + d + ") = " + qm + ", qr = " + qr);
414            if (qm.isZERO() || !qm.degreeVector().equals(qdegv)) {
415                if (debug) {
416                   logger.info("un-lucky evaluation point " + d + ", qm = " + qm.degreeVector() + " < " + qdegv);
417                }
418                continue;
419            }
420            rm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, rr, d);
421            //logger.info("rr(" + d + ") = " + rm + ", rr = " + rr);
422            if (rm.isZERO() || !rm.degreeVector().equals(rdegv)) {
423                if (debug) {
424                    logger.info("un-lucky evaluation point " + d + ", rm = " + rm.degreeVector() + " < " + rdegv);
425                }
426                continue;
427            }
428            // compute modular resultant in recursion
429            cm = resultant(rm, qm);
430            //System.out.println("cm = " + cm);
431
432            // prepare interpolation algorithm
433            if (M == null) {
434                // initialize interpolation
435                M = ufac.getONE();
436                cp = rfac.getZERO();
437            }
438            // interpolate
439            mi = PolyUtil.<MOD> evaluateMain(cofac, M, d);
440            mi = mi.inverse(); // mod p
441            cp = PolyUtil.interpolate(rfac, cp, M, mi, cm, d);
442            //logger.info("cp = " + cp);
443            mn = ufac.getONE().multiply(d);
444            mn = ufac.univariate(0).subtract(mn);
445            M = M.multiply(mn);
446            // test for completion
447            if (M.degree(0) > G) {
448                if (debug) {
449                    logger.info("last lucky evaluation point " + d);
450                }
451                break;
452            }
453            //logger.info("M  = " + M);
454        }
455        // distribute
456        q = PolyUtil.<MOD> distribute(fac, cp);
457        return q;
458    }
459
460}