Source code

001/*
002 * $Id: PolyGBUtil.java 5935 2018-09-30 11:47:20Z kredel $
003 */
004
005package edu.jas.gbufd;
006
007
008import java.util.ArrayList;
009import java.util.List;
010
011import org.apache.logging.log4j.Logger;
012import org.apache.logging.log4j.LogManager; 
013
014import edu.jas.gb.GroebnerBaseAbstract;
015import edu.jas.gb.SolvableGroebnerBaseAbstract;
016import edu.jas.gb.SolvableGroebnerBaseSeq;
017import edu.jas.gb.SolvableReductionAbstract;
018import edu.jas.gb.SolvableReductionSeq;
019import edu.jas.gb.WordGroebnerBaseAbstract;
020import edu.jas.gb.WordGroebnerBaseSeq;
021import edu.jas.poly.ExpVector;
022import edu.jas.poly.GenPolynomial;
023import edu.jas.poly.GenPolynomialRing;
024import edu.jas.poly.GenSolvablePolynomial;
025import edu.jas.poly.GenSolvablePolynomialRing;
026import edu.jas.poly.WordFactory;
027import edu.jas.poly.GenWordPolynomial;
028import edu.jas.poly.GenWordPolynomialRing;
029import edu.jas.poly.PolyUtil;
030import edu.jas.structure.GcdRingElem;
031import edu.jas.structure.RingElem;
032
033
034/**
035 * Package gbufd utilities.
036 * @author Heinz Kredel
037 */
038
039public class PolyGBUtil {
040
041
042    private static final Logger logger = LogManager.getLogger(PolyGBUtil.class);
043
044
045    private static final boolean debug = logger.isDebugEnabled();
046
047
048    /**
049     * Test for resultant.
050     * @param A generic polynomial.
051     * @param B generic polynomial.
052     * @param r generic polynomial.
053     * @return true if res(A,B) isContained in ideal(A,B), else false.
054     */
055    public static <C extends GcdRingElem<C>> boolean isResultant(GenPolynomial<C> A, GenPolynomial<C> B,
056                    GenPolynomial<C> r) {
057        if (r == null || r.isZERO()) {
058            return true;
059        }
060        GroebnerBaseAbstract<C> bb = GBFactory.<C> getImplementation(r.ring.coFac);
061        List<GenPolynomial<C>> F = new ArrayList<GenPolynomial<C>>(2);
062        F.add(A);
063        F.add(B);
064        List<GenPolynomial<C>> G = bb.GB(F);
065        //System.out.println("G = " + G);
066        GenPolynomial<C> n = bb.red.normalform(G, r);
067        //System.out.println("n = " + n);
068        return n.isZERO();
069    }
070
071
072    /**
073     * Top pseudo reduction wrt the main variables.
074     * @param P generic polynomial.
075     * @param A list of generic polynomials sorted according to appearing main
076     *            variables.
077     * @return top pseudo remainder of P wrt. A for the appearing variables.
078     */
079    public static <C extends RingElem<C>> GenPolynomial<C> topPseudoRemainder(List<GenPolynomial<C>> A,
080                    GenPolynomial<C> P) {
081        if (A == null || A.isEmpty()) {
082            return P.monic();
083        }
084        if (P.isZERO()) {
085            return P;
086        }
087        //System.out.println("remainder, P = " + P);
088        GenPolynomialRing<C> pfac = A.get(0).ring;
089        if (pfac.nvar <= 1) { // recursion base 
090            GenPolynomial<C> R = PolyUtil.<C> baseSparsePseudoRemainder(P, A.get(0));
091            return R.monic();
092        }
093        // select polynomials according to the main variable
094        GenPolynomialRing<GenPolynomial<C>> rfac = pfac.recursive(1);
095        GenPolynomial<C> Q = A.get(0); // wrong, must eventually search polynomial
096        GenPolynomial<GenPolynomial<C>> qr = PolyUtil.<C> recursive(rfac, Q);
097        GenPolynomial<GenPolynomial<C>> pr = PolyUtil.<C> recursive(rfac, P);
098        GenPolynomial<GenPolynomial<C>> rr;
099        if (qr.isONE()) {
100            return P.ring.getZERO();
101        }
102        if (qr.degree(0) > 0) {
103            rr = PolyUtil.<C> recursiveSparsePseudoRemainder(pr, qr);
104            //System.out.println("remainder, pr = " + pr);
105            //System.out.println("remainder, qr = " + qr);
106            //System.out.println("remainder, rr = " + rr);
107        } else {
108            rr = pr;
109        }
110        if (rr.degree(0) > 0) {
111            GenPolynomial<C> R = PolyUtil.<C> distribute(pfac, rr);
112            return R.monic();
113            // not further reduced wrt. other variables = top-reduction only
114        }
115        List<GenPolynomial<C>> zeroDeg = zeroDegrees(A);
116        GenPolynomial<C> R = topPseudoRemainder(zeroDeg, rr.leadingBaseCoefficient());
117        R = R.extend(pfac, 0, 0L);
118        return R.monic();
119    }
120
121
122    /**
123     * Top coefficient pseudo remainder of the leading coefficient of P wrt A in
124     * the main variables.
125     * @param P generic polynomial in n+1 variables.
126     * @param A list of generic polynomials in n variables sorted according to
127     *            appearing main variables.
128     * @return pseudo remainder of the leading coefficient of P wrt A.
129     */
130    public static <C extends RingElem<C>> GenPolynomial<C> topCoefficientPseudoRemainder(
131                    List<GenPolynomial<C>> A, GenPolynomial<C> P) {
132        if (A == null || A.isEmpty()) {
133            return P.monic();
134        }
135        if (P.isZERO()) {
136            return P;
137        }
138        GenPolynomialRing<C> pfac = P.ring;
139        GenPolynomialRing<C> pfac1 = A.get(0).ring;
140        if (pfac1.nvar <= 1) { // recursion base 
141            GenPolynomial<C> a = A.get(0);
142            GenPolynomialRing<GenPolynomial<C>> rfac = pfac.recursive(pfac.nvar - 1);
143            GenPolynomial<GenPolynomial<C>> pr = PolyUtil.<C> recursive(rfac, P);
144            // ldcf(P,x_m) = q a + r 
145            GenPolynomial<GenPolynomial<C>> rr = PolyGBUtil.<C> coefficientPseudoRemainderBase(pr, a);
146            GenPolynomial<C> R = PolyUtil.<C> distribute(pfac, rr);
147            return R.monic();
148        }
149        // select polynomials according to the main variable
150        GenPolynomialRing<GenPolynomial<C>> rfac1 = pfac1.recursive(1);
151        int nv = pfac.nvar - pfac1.nvar;
152        GenPolynomialRing<GenPolynomial<C>> rfac = pfac.recursive(1 + nv);
153        GenPolynomialRing<GenPolynomial<GenPolynomial<C>>> rfac2 = rfac.recursive(nv);
154        if (debug) {
155            logger.info("rfac =" + rfac);
156        }
157        GenPolynomial<GenPolynomial<C>> pr = PolyUtil.<C> recursive(rfac, P);
158        GenPolynomial<GenPolynomial<GenPolynomial<C>>> pr2 = PolyUtil.<GenPolynomial<C>> recursive(rfac2, pr);
159        //System.out.println("recursion, pr2 = " + pr2);
160        GenPolynomial<C> Q = A.get(0);
161        GenPolynomial<GenPolynomial<C>> qr = PolyUtil.<C> recursive(rfac1, Q);
162        GenPolynomial<GenPolynomial<GenPolynomial<C>>> rr;
163        if (qr.isONE()) {
164            return P.ring.getZERO();
165        }
166        if (qr.degree(0) > 0) {
167            // pseudo remainder:  ldcf(P,x_m) = a q + r 
168            rr = PolyGBUtil.<C> coefficientPseudoRemainder(pr2, qr);
169            //System.out.println("recursion, qr  = " + qr);
170            //System.out.println("recursion, pr  = " + pr2);
171            //System.out.println("recursion, rr  = " + rr);
172        } else {
173            rr = pr2;
174        }
175        // reduction wrt. the other variables
176        List<GenPolynomial<C>> zeroDeg = zeroDegrees(A);
177        GenPolynomial<GenPolynomial<C>> Rr = PolyUtil.<GenPolynomial<C>> distribute(rfac, rr);
178        GenPolynomial<C> R = PolyUtil.<C> distribute(pfac, Rr);
179        R = topCoefficientPseudoRemainder(zeroDeg, R);
180        return R.monic();
181    }
182
183
184    /**
185     * Polynomial leading coefficient pseudo remainder.
186     * @param P generic polynomial in n+1 variables.
187     * @param A generic polynomial in n variables.
188     * @return pseudo remainder of the leading coefficient of P wrt A, with
189     *         ldcf(A)<sup>m'</sup> P = quotient * A + remainder.
190     */
191    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<GenPolynomial<C>>> coefficientPseudoRemainder(
192                    GenPolynomial<GenPolynomial<GenPolynomial<C>>> P, GenPolynomial<GenPolynomial<C>> A) {
193        if (A == null || A.isZERO()) { // findbugs
194            throw new ArithmeticException(P + " division by zero " + A);
195        }
196        if (A.isONE()) {
197            return P.ring.getZERO();
198        }
199        if (P.isZERO() || P.isONE()) {
200            return P;
201        }
202        GenPolynomialRing<GenPolynomial<GenPolynomial<C>>> pfac = P.ring;
203        GenPolynomialRing<GenPolynomial<C>> afac = A.ring; // == pfac.coFac
204        GenPolynomial<GenPolynomial<GenPolynomial<C>>> r = P;
205        GenPolynomial<GenPolynomial<C>> h;
206        GenPolynomial<GenPolynomial<GenPolynomial<C>>> hr;
207        GenPolynomial<GenPolynomial<C>> ldcf = P.leadingBaseCoefficient();
208        long m = ldcf.degree(0);
209        long n = A.degree(0);
210        GenPolynomial<C> c = A.leadingBaseCoefficient();
211        GenPolynomial<GenPolynomial<C>> cc = afac.getZERO().sum(c);
212        //System.out.println("cc = " + cc);
213        ExpVector e = A.leadingExpVector();
214        for (long i = m; i >= n; i--) {
215            if (r.isZERO()) {
216                return r;
217            }
218            GenPolynomial<GenPolynomial<C>> p = r.leadingBaseCoefficient();
219            ExpVector g = r.leadingExpVector();
220            long k = p.degree(0);
221            if (i == k) {
222                GenPolynomial<C> pl = p.leadingBaseCoefficient();
223                ExpVector f = p.leadingExpVector();
224                f = f.subtract(e);
225                r = r.multiply(cc); // coeff cc
226                h = A.multiply(pl, f); // coeff ac
227                hr = new GenPolynomial<GenPolynomial<GenPolynomial<C>>>(pfac, h, g);
228                r = r.subtract(hr);
229            } else {
230                r = r.multiply(cc);
231            }
232            //System.out.println("r = " + r);
233        }
234        if (r.degree(0) < P.degree(0)) { // recursion for degree
235            r = coefficientPseudoRemainder(r, A);
236        }
237        return r;
238    }
239
240
241    /**
242     * Polynomial leading coefficient pseudo remainder, base case.
243     * @param P generic polynomial in 1+1 variables.
244     * @param A generic polynomial in 1 variable.
245     * @return pseudo remainder of the leading coefficient of P wrt. A, with
246     *         ldcf(A)<sup>m'</sup> P = quotient * A + remainder.
247     */
248    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> coefficientPseudoRemainderBase(
249                    GenPolynomial<GenPolynomial<C>> P, GenPolynomial<C> A) {
250        if (A == null || A.isZERO()) { // findbugs
251            throw new ArithmeticException(P + " division by zero " + A);
252        }
253        if (A.isONE()) {
254            return P.ring.getZERO();
255        }
256        if (P.isZERO() || P.isONE()) {
257            return P;
258        }
259        GenPolynomialRing<GenPolynomial<C>> pfac = P.ring;
260        GenPolynomialRing<C> afac = A.ring; // == pfac.coFac
261        GenPolynomial<GenPolynomial<C>> r = P;
262        GenPolynomial<C> h;
263        GenPolynomial<GenPolynomial<C>> hr;
264        GenPolynomial<C> ldcf = P.leadingBaseCoefficient();
265        long m = ldcf.degree(0);
266        long n = A.degree(0);
267        C c = A.leadingBaseCoefficient();
268        GenPolynomial<C> cc = afac.getZERO().sum(c);
269        //System.out.println("cc = " + cc);
270        ExpVector e = A.leadingExpVector();
271        for (long i = m; i >= n; i--) {
272            if (r.isZERO()) {
273                return r;
274            }
275            GenPolynomial<C> p = r.leadingBaseCoefficient();
276            ExpVector g = r.leadingExpVector();
277            long k = p.degree(0);
278            if (i == k) {
279                C pl = p.leadingBaseCoefficient();
280                ExpVector f = p.leadingExpVector();
281                f = f.subtract(e);
282                r = r.multiply(cc); // coeff cc
283                h = A.multiply(pl, f); // coeff ac
284                hr = new GenPolynomial<GenPolynomial<C>>(pfac, h, g);
285                r = r.subtract(hr);
286            } else {
287                r = r.multiply(cc);
288            }
289            //System.out.println("r = " + r);
290        }
291        if (r.degree(0) < P.degree(0)) { // recursion for degree
292            r = coefficientPseudoRemainderBase(r, A);
293        }
294        return r;
295    }
296
297
298    /**
299     * Extract polynomials with degree zero in the main variable.
300     * @param A list of generic polynomials in n variables.
301     * @return Z = [a_i] with deg(a_i,x_n) = 0 and in n-1 variables.
302     */
303    public static <C extends RingElem<C>> List<GenPolynomial<C>> zeroDegrees(List<GenPolynomial<C>> A) {
304        if (A == null || A.isEmpty()) {
305            return A;
306        }
307        GenPolynomialRing<C> pfac = A.get(0).ring;
308        GenPolynomialRing<GenPolynomial<C>> rfac = pfac.recursive(1);
309        List<GenPolynomial<C>> zeroDeg = new ArrayList<GenPolynomial<C>>(A.size());
310        for (int i = 0; i < A.size(); i++) {
311            GenPolynomial<C> q = A.get(i);
312            GenPolynomial<GenPolynomial<C>> fr = PolyUtil.<C> recursive(rfac, q);
313            if (fr.degree(0) == 0) {
314                zeroDeg.add(fr.leadingBaseCoefficient());
315            }
316        }
317        return zeroDeg;
318    }
319
320
321    /**
322     * Intersection. Generators for the intersection of ideals.
323     * @param pfac polynomial ring
324     * @param A list of polynomials
325     * @param B list of polynomials
326     * @return generators for (A \cap B)
327     */
328    public static <C extends GcdRingElem<C>> List<GenPolynomial<C>> intersect(GenPolynomialRing<C> pfac,
329                    List<GenPolynomial<C>> A, List<GenPolynomial<C>> B) {
330        if (A == null || A.isEmpty()) { // (0)
331            return B;
332        }
333        if (B == null || B.isEmpty()) { // (0)
334            return A;
335        }
336        int s = A.size() + B.size();
337        List<GenPolynomial<C>> c = new ArrayList<GenPolynomial<C>>(s);
338        GenPolynomialRing<C> tfac = pfac.extend(1);
339        // term order is also adjusted
340        for (GenPolynomial<C> p : A) {
341            p = p.extend(tfac, 0, 1L); // t*p
342            c.add(p);
343        }
344        for (GenPolynomial<C> p : B) {
345            GenPolynomial<C> q = p.extend(tfac, 0, 1L);
346            GenPolynomial<C> r = p.extend(tfac, 0, 0L);
347            p = r.subtract(q); // (1-t)*p
348            c.add(p);
349        }
350        GroebnerBaseAbstract<C> bb = GBFactory.<C> getImplementation(tfac.coFac);
351        logger.warn("intersect computing GB");
352        List<GenPolynomial<C>> G = bb.GB(c);
353        if (debug) {
354            logger.debug("intersect GB = " + G);
355        }
356        List<GenPolynomial<C>> I = PolyUtil.<C> intersect(pfac, G);
357        return I;
358    }
359
360
361    /**
362     * Intersection. Generators for the intersection of ideals.
363     * @param pfac solvable polynomial ring
364     * @param A list of polynomials
365     * @param B list of polynomials
366     * @return generators for (A \cap B)
367     */
368    public static <C extends GcdRingElem<C>> List<GenSolvablePolynomial<C>> intersect(
369                    GenSolvablePolynomialRing<C> pfac, List<GenSolvablePolynomial<C>> A,
370                    List<GenSolvablePolynomial<C>> B) {
371        if (A == null || A.isEmpty()) { // (0)
372            return B;
373        }
374        if (B == null || B.isEmpty()) { // (0)
375            return A;
376        }
377        int s = A.size() + B.size();
378        List<GenSolvablePolynomial<C>> c = new ArrayList<GenSolvablePolynomial<C>>(s);
379        GenSolvablePolynomialRing<C> tfac = pfac.extend(1);
380        // term order is also adjusted
381        for (GenSolvablePolynomial<C> p : A) {
382            p = (GenSolvablePolynomial<C>) p.extend(tfac, 0, 1L); // t*p
383            c.add(p);
384        }
385        for (GenSolvablePolynomial<C> p : B) {
386            GenSolvablePolynomial<C> q = (GenSolvablePolynomial<C>) p.extend(tfac, 0, 1L);
387            GenSolvablePolynomial<C> r = (GenSolvablePolynomial<C>) p.extend(tfac, 0, 0L);
388            p = (GenSolvablePolynomial<C>) r.subtract(q); // (1-t)*p
389            c.add(p);
390        }
391        SolvableGroebnerBaseAbstract<C> sbb = SGBFactory.<C> getImplementation(tfac.coFac);
392        //new SolvableGroebnerBaseSeq<C>();
393        logger.warn("intersect computing GB");
394        List<GenSolvablePolynomial<C>> g = sbb.leftGB(c);
395        //List<GenSolvablePolynomial<C>> g = sbb.twosidedGB(c);
396        if (debug) {
397            logger.debug("intersect GB = " + g);
398        }
399        List<GenSolvablePolynomial<C>> I = PolyUtil.<C> intersect(pfac, g);
400        return I;
401    }
402
403
404    /**
405     * Intersection. Generators for the intersection of word ideals.
406     * @param pfac word polynomial ring
407     * @param A list of word polynomials
408     * @param B list of word polynomials
409     * @return generators for (A \cap B) if it exists
410     */
411    public static <C extends GcdRingElem<C>> List<GenWordPolynomial<C>> intersect(GenWordPolynomialRing<C> pfac,
412                    List<GenWordPolynomial<C>> A, List<GenWordPolynomial<C>> B) {
413        if (A == null || A.isEmpty()) { // (0)
414            return B;
415        }
416        if (B == null || B.isEmpty()) { // (0)
417            return A;
418        }
419        int s = A.size() + B.size();
420        List<GenWordPolynomial<C>> L = new ArrayList<GenWordPolynomial<C>>(s);
421        GenWordPolynomialRing<C> tfac = pfac.extend(1);
422        List<GenWordPolynomial<C>> gens = tfac.univariateList();
423        //System.out.println("gens = " + gens);
424        GenWordPolynomial<C> t = gens.get(gens.size()-1);
425        //System.out.println("t = " + t);
426        // make t commute with other variables
427        for (GenWordPolynomial<C> p : gens) {
428            if (t == p) {
429                continue;
430            }
431            GenWordPolynomial<C> c = t.multiply(p).subtract(p.multiply(t)); // t p - p t
432            L.add(c);
433        }       
434        for (GenWordPolynomial<C> p : A) {
435            p = tfac.valueOf(p).multiply(t); // t p
436            L.add(p);
437        }
438        for (GenWordPolynomial<C> p : B) {
439            GenWordPolynomial<C> q = tfac.valueOf(p).multiply(t); 
440            GenWordPolynomial<C> r = tfac.valueOf(p);
441            p = r.subtract(q); // (1-t) p
442            L.add(p);
443        }
444        //System.out.println("L = " + L);
445        WordGroebnerBaseAbstract<C> bb = new WordGroebnerBaseSeq<C>();
446        logger.warn("intersect computing GB");
447        List<GenWordPolynomial<C>> G = bb.GB(L);
448        //System.out.println("G = " + G);
449        if (debug) {
450            logger.debug("intersect GB = " + G);
451        }
452        List<GenWordPolynomial<C>> I = PolyUtil.<C> intersect(pfac, G);
453        return I;
454    }
455
456    
457    /**
458     * Solvable quotient and remainder via reduction.
459     * @param n first solvable polynomial.
460     * @param d second solvable polynomial.
461     * @return [ n/d, n - (n/d)*d ]
462     */
463    @SuppressWarnings("unchecked")
464    public static <C extends GcdRingElem<C>> GenSolvablePolynomial<C>[] quotientRemainder(
465                    GenSolvablePolynomial<C> n, GenSolvablePolynomial<C> d) {
466        GenSolvablePolynomial<C>[] res = (GenSolvablePolynomial<C>[]) new GenSolvablePolynomial[2];
467        if (d.isZERO()) {
468            throw new RuntimeException("division by zero: " + n + "/" + d);
469        }
470        if (n.isZERO()) {
471            res[0] = n;
472            res[1] = n;
473            return res;
474        }
475        GenSolvablePolynomialRing<C> r = n.ring;
476        if (d.isONE()) {
477            res[0] = n;
478            res[1] = r.getZERO();
479            return res;
480        }
481        // divide
482        List<GenSolvablePolynomial<C>> Q = new ArrayList<GenSolvablePolynomial<C>>(1);
483        Q.add(r.getZERO());
484        List<GenSolvablePolynomial<C>> D = new ArrayList<GenSolvablePolynomial<C>>(1);
485        D.add(d);
486        SolvableReductionAbstract<C> sred = new SolvableReductionSeq<C>();
487        res[1] = sred.rightNormalform(Q, D, n); // left
488        res[0] = Q.get(0);
489        return res;
490    }
491
492}