Source code

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