Source code

001/*
002 * $Id: GreatestCommonDivisorAbstract.java 4148 2012-08-31 19:49:27Z kredel $
003 */
004
005package edu.jas.ufd;
006
007
008import java.util.ArrayList;
009import java.util.List;
010
011import org.apache.log4j.Logger;
012
013import edu.jas.poly.GenPolynomial;
014import edu.jas.poly.GenPolynomialRing;
015import edu.jas.poly.PolyUtil;
016import edu.jas.structure.GcdRingElem;
017import edu.jas.structure.RingFactory;
018
019
020/**
021 * Greatest common divisor algorithms.
022 * @author Heinz Kredel
023 */
024
025public abstract class GreatestCommonDivisorAbstract<C extends GcdRingElem<C>> implements
026                GreatestCommonDivisor<C> {
027
028
029    private static final Logger logger = Logger.getLogger(GreatestCommonDivisorAbstract.class);
030
031
032    private final boolean debug = logger.isDebugEnabled();
033
034
035    /**
036     * Get the String representation.
037     * @see java.lang.Object#toString()
038     */
039    @Override
040    public String toString() {
041        return getClass().getName();
042    }
043
044
045    /**
046     * GenPolynomial base coefficient content.
047     * @param P GenPolynomial.
048     * @return cont(P).
049     */
050    public C baseContent(GenPolynomial<C> P) {
051        if (P == null) {
052            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
053        }
054        if (P.isZERO()) {
055            return P.ring.getZEROCoefficient();
056        }
057        C d = null;
058        for (C c : P.getMap().values()) {
059            if (d == null) {
060                d = c;
061            } else {
062                d = d.gcd(c);
063            }
064            if (d.isONE()) {
065                return d;
066            }
067        }
068        if (d.signum() < 0) {
069            d = d.negate();
070        }
071        return d;
072    }
073
074
075    /**
076     * GenPolynomial base coefficient primitive part.
077     * @param P GenPolynomial.
078     * @return pp(P).
079     */
080    public GenPolynomial<C> basePrimitivePart(GenPolynomial<C> P) {
081        if (P == null) {
082            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
083        }
084        if (P.isZERO()) {
085            return P;
086        }
087        C d = baseContent(P);
088        if (d.isONE()) {
089            return P;
090        }
091        GenPolynomial<C> pp = P.divide(d);
092        if (debug) {
093            GenPolynomial<C> p = pp.multiply(d);
094            if (!p.equals(P)) {
095                throw new ArithmeticException("pp(p)*cont(p) != p: ");
096            }
097        }
098        return pp;
099    }
100
101
102    /**
103     * Univariate GenPolynomial greatest common divisor. Uses sparse
104     * pseudoRemainder for remainder.
105     * @param P univariate GenPolynomial.
106     * @param S univariate GenPolynomial.
107     * @return gcd(P,S).
108     */
109    public abstract GenPolynomial<C> baseGcd(GenPolynomial<C> P, GenPolynomial<C> S);
110
111
112    /**
113     * GenPolynomial recursive content.
114     * @param P recursive GenPolynomial.
115     * @return cont(P).
116     */
117    public GenPolynomial<C> recursiveContent(GenPolynomial<GenPolynomial<C>> P) {
118        if (P == null) {
119            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
120        }
121        if (P.isZERO()) {
122            return P.ring.getZEROCoefficient();
123        }
124        GenPolynomial<C> d = null;
125        for (GenPolynomial<C> c : P.getMap().values()) {
126            if (d == null) {
127                d = c;
128            } else {
129                d = gcd(d, c); // go to recursion
130            }
131            if (d.isONE()) {
132                return d;
133            }
134        }
135        return d.abs();
136    }
137
138
139    /**
140     * GenPolynomial recursive primitive part.
141     * @param P recursive GenPolynomial.
142     * @return pp(P).
143     */
144    public GenPolynomial<GenPolynomial<C>> recursivePrimitivePart(GenPolynomial<GenPolynomial<C>> P) {
145        if (P == null) {
146            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
147        }
148        if (P.isZERO()) {
149            return P;
150        }
151        GenPolynomial<C> d = recursiveContent(P);
152        if (d.isONE()) {
153            return P;
154        }
155        GenPolynomial<GenPolynomial<C>> pp = PolyUtil.<C> recursiveDivide(P, d);
156        return pp;
157    }
158
159
160    /**
161     * GenPolynomial base recursive content.
162     * @param P recursive GenPolynomial.
163     * @return baseCont(P).
164     */
165    public C baseRecursiveContent(GenPolynomial<GenPolynomial<C>> P) {
166        if (P == null) {
167            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
168        }
169        if (P.isZERO()) {
170            GenPolynomialRing<C> cf = (GenPolynomialRing<C>) P.ring.coFac;
171            return cf.coFac.getZERO();
172        }
173        C d = null;
174        for (GenPolynomial<C> c : P.getMap().values()) {
175            C cc = baseContent(c);
176            if (d == null) {
177                d = cc;
178            } else {
179                d = gcd(d, cc);
180            }
181            if (d.isONE()) {
182                return d;
183            }
184        }
185        if (d.signum() < 0) {
186            d = d.negate();
187        }
188        return d;
189    }
190
191
192    /**
193     * GenPolynomial base recursive primitive part.
194     * @param P recursive GenPolynomial.
195     * @return basePP(P).
196     */
197    public GenPolynomial<GenPolynomial<C>> baseRecursivePrimitivePart(GenPolynomial<GenPolynomial<C>> P) {
198        if (P == null) {
199            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
200        }
201        if (P.isZERO()) {
202            return P;
203        }
204        C d = baseRecursiveContent(P);
205        if (d.isONE()) {
206            return P;
207        }
208        GenPolynomial<GenPolynomial<C>> pp = PolyUtil.<C> baseRecursiveDivide(P, d);
209        return pp;
210    }
211
212
213    /**
214     * GenPolynomial recursive greatest common divisor. Uses pseudoRemainder for
215     * remainder.
216     * @param P recursive GenPolynomial.
217     * @param S recursive GenPolynomial.
218     * @return gcd(P,S).
219     */
220    public GenPolynomial<GenPolynomial<C>> recursiveGcd(GenPolynomial<GenPolynomial<C>> P,
221                    GenPolynomial<GenPolynomial<C>> S) {
222        if (S == null || S.isZERO()) {
223            return P;
224        }
225        if (P == null || P.isZERO()) {
226            return S;
227        }
228        if (P.ring.nvar <= 1) {
229            return recursiveUnivariateGcd(P, S);
230        }
231        // distributed polynomials gcd
232        GenPolynomialRing<GenPolynomial<C>> rfac = P.ring;
233        RingFactory<GenPolynomial<C>> rrfac = rfac.coFac;
234        GenPolynomialRing<C> cfac = (GenPolynomialRing<C>) rrfac;
235        GenPolynomialRing<C> dfac = cfac.extend(rfac.nvar);
236        GenPolynomial<C> Pd = PolyUtil.<C> distribute(dfac, P);
237        GenPolynomial<C> Sd = PolyUtil.<C> distribute(dfac, S);
238        GenPolynomial<C> Dd = gcd(Pd, Sd);
239        // convert to recursive
240        GenPolynomial<GenPolynomial<C>> C = PolyUtil.<C> recursive(rfac, Dd);
241        return C;
242    }
243
244
245    /**
246     * Univariate GenPolynomial recursive greatest common divisor. Uses
247     * pseudoRemainder for remainder.
248     * @param P univariate recursive GenPolynomial.
249     * @param S univariate recursive GenPolynomial.
250     * @return gcd(P,S).
251     */
252    public abstract GenPolynomial<GenPolynomial<C>> recursiveUnivariateGcd(GenPolynomial<GenPolynomial<C>> P,
253                    GenPolynomial<GenPolynomial<C>> S);
254
255
256    /**
257     * GenPolynomial content.
258     * @param P GenPolynomial.
259     * @return cont(P).
260     */
261    public GenPolynomial<C> content(GenPolynomial<C> P) {
262        if (P == null) {
263            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
264        }
265        GenPolynomialRing<C> pfac = P.ring;
266        if (pfac.nvar <= 1) {
267            // baseContent not possible by return type
268            throw new IllegalArgumentException(this.getClass().getName()
269                            + " use baseContent for univariate polynomials");
270
271        }
272        GenPolynomialRing<C> cfac = pfac.contract(1);
273        GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(cfac, 1);
274
275        GenPolynomial<GenPolynomial<C>> Pr = PolyUtil.<C> recursive(rfac, P);
276        GenPolynomial<C> D = recursiveContent(Pr);
277        return D;
278    }
279
280
281    /**
282     * GenPolynomial primitive part.
283     * @param P GenPolynomial.
284     * @return pp(P).
285     */
286    public GenPolynomial<C> primitivePart(GenPolynomial<C> P) {
287        if (P == null) {
288            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
289        }
290        if (P.isZERO()) {
291            return P;
292        }
293        GenPolynomialRing<C> pfac = P.ring;
294        if (pfac.nvar <= 1) {
295            return basePrimitivePart(P);
296        }
297        GenPolynomialRing<C> cfac = pfac.contract(1);
298        GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(cfac, 1);
299
300        GenPolynomial<GenPolynomial<C>> Pr = PolyUtil.<C> recursive(rfac, P);
301        GenPolynomial<GenPolynomial<C>> PP = recursivePrimitivePart(Pr);
302
303        GenPolynomial<C> D = PolyUtil.<C> distribute(pfac, PP);
304        return D;
305    }
306
307
308    /**
309     * GenPolynomial division. Indirection to GenPolynomial method.
310     * @param a GenPolynomial.
311     * @param b coefficient.
312     * @return a/b.
313     */
314    public GenPolynomial<C> divide(GenPolynomial<C> a, C b) {
315        if (b == null || b.isZERO()) {
316            throw new IllegalArgumentException("division by zero");
317
318        }
319        if (a == null || a.isZERO()) {
320            return a;
321        }
322        return a.divide(b);
323    }
324
325
326    /**
327     * Coefficient greatest common divisor. Indirection to coefficient method.
328     * @param a coefficient.
329     * @param b coefficient.
330     * @return gcd(a,b).
331     */
332    public C gcd(C a, C b) {
333        if (b == null || b.isZERO()) {
334            return a;
335        }
336        if (a == null || a.isZERO()) {
337            return b;
338        }
339        return a.gcd(b);
340    }
341
342
343    /**
344     * GenPolynomial greatest common divisor.
345     * @param P GenPolynomial.
346     * @param S GenPolynomial.
347     * @return gcd(P,S).
348     */
349    public GenPolynomial<C> gcd(GenPolynomial<C> P, GenPolynomial<C> S) {
350        if (S == null || S.isZERO()) {
351            return P;
352        }
353        if (P == null || P.isZERO()) {
354            return S;
355        }
356        GenPolynomialRing<C> pfac = P.ring;
357        if (pfac.nvar <= 1) {
358            GenPolynomial<C> T = baseGcd(P, S);
359            return T;
360        }
361        GenPolynomialRing<C> cfac = pfac.contract(1);
362        GenPolynomialRing<GenPolynomial<C>> rfac;
363        if (pfac.getVars() != null && pfac.getVars().length > 0) {
364            String[] v = new String[] { pfac.getVars()[pfac.nvar - 1] };
365            rfac = new GenPolynomialRing<GenPolynomial<C>>(cfac, 1, v);
366        } else {
367            rfac = new GenPolynomialRing<GenPolynomial<C>>(cfac, 1);
368        }
369        GenPolynomial<GenPolynomial<C>> Pr = PolyUtil.<C> recursive(rfac, P);
370        GenPolynomial<GenPolynomial<C>> Sr = PolyUtil.<C> recursive(rfac, S);
371        GenPolynomial<GenPolynomial<C>> Dr = recursiveUnivariateGcd(Pr, Sr);
372        GenPolynomial<C> D = PolyUtil.<C> distribute(pfac, Dr);
373        return D;
374    }
375
376
377    /**
378     * GenPolynomial least common multiple.
379     * @param P GenPolynomial.
380     * @param S GenPolynomial.
381     * @return lcm(P,S).
382     */
383    public GenPolynomial<C> lcm(GenPolynomial<C> P, GenPolynomial<C> S) {
384        if (S == null || S.isZERO()) {
385            return S;
386        }
387        if (P == null || P.isZERO()) {
388            return P;
389        }
390        GenPolynomial<C> C = gcd(P, S);
391        GenPolynomial<C> A = P.multiply(S);
392        return PolyUtil.<C> basePseudoDivide(A, C);
393    }
394
395
396    /**
397     * List of GenPolynomials greatest common divisor.
398     * @param A non empty list of GenPolynomials.
399     * @return gcd(A_i).
400     */
401    public GenPolynomial<C> gcd(List<GenPolynomial<C>> A) {
402        if (A == null || A.isEmpty()) {
403            throw new IllegalArgumentException("A may not be empty");
404        }
405        GenPolynomial<C> g = A.get(0);
406        for (int i = 1; i < A.size(); i++) {
407            GenPolynomial<C> f = A.get(i);
408            g = gcd(g, f);
409        }
410        return g;
411    }
412
413
414    /**
415     * Univariate GenPolynomial resultant.
416     * @param P univariate GenPolynomial.
417     * @param S univariate GenPolynomial.
418     * @return res(P,S).
419     * @throws UnsupportedOperationException if there is no implementation in
420     *             the sub-class.
421     */
422    @SuppressWarnings("unused")
423    public GenPolynomial<C> baseResultant(GenPolynomial<C> P, GenPolynomial<C> S) {
424        // can not be abstract
425        throw new UnsupportedOperationException("not implmented");
426    }
427
428
429    /**
430     * Univariate GenPolynomial recursive resultant.
431     * @param P univariate recursive GenPolynomial.
432     * @param S univariate recursive GenPolynomial.
433     * @return res(P,S).
434     * @throws UnsupportedOperationException if there is no implementation in
435     *             the sub-class.
436     */
437    @SuppressWarnings("unused")
438    public GenPolynomial<GenPolynomial<C>> recursiveUnivariateResultant(GenPolynomial<GenPolynomial<C>> P,
439                    GenPolynomial<GenPolynomial<C>> S) {
440        // can not be abstract
441        throw new UnsupportedOperationException("not implmented");
442    }
443
444
445    /**
446     * GenPolynomial recursive resultant.
447     * @param P univariate recursive GenPolynomial.
448     * @param S univariate recursive GenPolynomial.
449     * @return res(P,S).
450     * @throws UnsupportedOperationException if there is no implementation in
451     *             the sub-class.
452     */
453    public GenPolynomial<GenPolynomial<C>> recursiveResultant(GenPolynomial<GenPolynomial<C>> P,
454                    GenPolynomial<GenPolynomial<C>> S) {
455        if (S == null || S.isZERO()) {
456            return S;
457        }
458        if (P == null || P.isZERO()) {
459            return P;
460        }
461        GenPolynomialRing<GenPolynomial<C>> rfac = P.ring;
462        GenPolynomialRing<C> cfac = (GenPolynomialRing<C>) rfac.coFac;
463        GenPolynomialRing<C> dfac = cfac.extend(rfac.getVars());
464        GenPolynomial<C> Pp = PolyUtil.<C> distribute(dfac, P);
465        GenPolynomial<C> Sp = PolyUtil.<C> distribute(dfac, S);
466        GenPolynomial<C> res = resultant(Pp, Sp);
467        GenPolynomial<GenPolynomial<C>> Rr = PolyUtil.<C> recursive(rfac, res);
468        return Rr;
469    }
470
471
472    /**
473     * GenPolynomial resultant. The input polynomials are considered as
474     * univariate polynomials in the main variable.
475     * @param P GenPolynomial.
476     * @param S GenPolynomial.
477     * @return res(P,S).
478     * @see edu.jas.ufd.GreatestCommonDivisorSubres#recursiveResultant
479     * @throws UnsupportedOperationException if there is no implementation in
480     *             the sub-class.
481     */
482    public GenPolynomial<C> resultant(GenPolynomial<C> P, GenPolynomial<C> S) {
483        if (S == null || S.isZERO()) {
484            return S;
485        }
486        if (P == null || P.isZERO()) {
487            return P;
488        }
489        // no more hacked: GreatestCommonDivisorSubres<C> ufd_sr = new GreatestCommonDivisorSubres<C>();
490        GenPolynomialRing<C> pfac = P.ring;
491        if (pfac.nvar <= 1) {
492            return baseResultant(P, S);
493        }
494        GenPolynomialRing<C> cfac = pfac.contract(1);
495        GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(cfac, 1);
496
497        GenPolynomial<GenPolynomial<C>> Pr = PolyUtil.<C> recursive(rfac, P);
498        GenPolynomial<GenPolynomial<C>> Sr = PolyUtil.<C> recursive(rfac, S);
499
500        GenPolynomial<GenPolynomial<C>> Dr = recursiveUnivariateResultant(Pr, Sr);
501        GenPolynomial<C> D = PolyUtil.<C> distribute(pfac, Dr);
502        return D;
503    }
504
505
506    /**
507     * GenPolynomial co-prime list.
508     * @param A list of GenPolynomials.
509     * @return B with gcd(b,c) = 1 for all b != c in B and for all non-constant
510     *         a in A there exists b in B with b|a. B does not contain zero or
511     *         constant polynomials.
512     */
513    public List<GenPolynomial<C>> coPrime(List<GenPolynomial<C>> A) {
514        if (A == null || A.isEmpty()) {
515            return A;
516        }
517        List<GenPolynomial<C>> B = new ArrayList<GenPolynomial<C>>(A.size());
518        // make a coprime to rest of list
519        GenPolynomial<C> a = A.get(0);
520        //System.out.println("a = " + a);
521        if (!a.isZERO() && !a.isConstant()) {
522            for (int i = 1; i < A.size(); i++) {
523                GenPolynomial<C> b = A.get(i);
524                GenPolynomial<C> g = gcd(a, b).abs();
525                if (!g.isONE()) {
526                    a = PolyUtil.<C> basePseudoDivide(a, g);
527                    b = PolyUtil.<C> basePseudoDivide(b, g);
528                    GenPolynomial<C> gp = gcd(a, g).abs();
529                    while (!gp.isONE()) {
530                        a = PolyUtil.<C> basePseudoDivide(a, gp);
531                        g = PolyUtil.<C> basePseudoDivide(g, gp);
532                        B.add(g); // gcd(a,g) == 1
533                        g = gp;
534                        gp = gcd(a, gp).abs();
535                    }
536                    if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
537                        B.add(g); // gcd(a,g) == 1
538                    }
539                }
540                if (!b.isZERO() && !b.isConstant()) {
541                    B.add(b); // gcd(a,b) == 1
542                }
543            }
544        } else {
545            B.addAll(A.subList(1, A.size()));
546        }
547        // make rest coprime
548        B = coPrime(B);
549        //System.out.println("B = " + B);
550        if (!a.isZERO() && !a.isConstant() /*&& !B.contains(a)*/) {
551            a = a.abs();
552            B.add(a);
553        }
554        return B;
555    }
556
557
558    /**
559     * GenPolynomial co-prime list.
560     * @param A list of GenPolynomials.
561     * @return B with gcd(b,c) = 1 for all b != c in B and for all non-constant
562     *         a in A there exists b in B with b|a. B does not contain zero or
563     *         constant polynomials.
564     */
565    public List<GenPolynomial<C>> coPrimeRec(List<GenPolynomial<C>> A) {
566        if (A == null || A.isEmpty()) {
567            return A;
568        }
569        List<GenPolynomial<C>> B = new ArrayList<GenPolynomial<C>>();
570        // make a co-prime to rest of list
571        for (GenPolynomial<C> a : A) {
572            //System.out.println("a = " + a);
573            B = coPrime(a, B);
574            //System.out.println("B = " + B);
575        }
576        return B;
577    }
578
579
580    /**
581     * GenPolynomial co-prime list.
582     * @param a GenPolynomial.
583     * @param P co-prime list of GenPolynomials.
584     * @return B with gcd(b,c) = 1 for all b != c in B and for non-constant a
585     *         there exists b in P with b|a. B does not contain zero or constant
586     *         polynomials.
587     */
588    public List<GenPolynomial<C>> coPrime(GenPolynomial<C> a, List<GenPolynomial<C>> P) {
589        if (a == null || a.isZERO() || a.isConstant()) {
590            return P;
591        }
592        List<GenPolynomial<C>> B = new ArrayList<GenPolynomial<C>>(P.size() + 1);
593        // make a coprime to elements of the list P
594        for (int i = 0; i < P.size(); i++) {
595            GenPolynomial<C> b = P.get(i);
596            GenPolynomial<C> g = gcd(a, b).abs();
597            if (!g.isONE()) {
598                a = PolyUtil.<C> basePseudoDivide(a, g);
599                b = PolyUtil.<C> basePseudoDivide(b, g);
600                // make g co-prime to new a, g is co-prime to c != b in P, B
601                GenPolynomial<C> gp = gcd(a, g).abs();
602                while (!gp.isONE()) {
603                    a = PolyUtil.<C> basePseudoDivide(a, gp);
604                    g = PolyUtil.<C> basePseudoDivide(g, gp);
605                    if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
606                        B.add(g); // gcd(a,g) == 1 and gcd(g,c) == 1 for c != b in P, B
607                    }
608                    g = gp;
609                    gp = gcd(a, gp).abs();
610                }
611                // make new g co-prime to new b
612                gp = gcd(b, g).abs();
613                while (!gp.isONE()) {
614                    b = PolyUtil.<C> basePseudoDivide(b, gp);
615                    g = PolyUtil.<C> basePseudoDivide(g, gp);
616                    if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
617                        B.add(g); // gcd(a,g) == 1 and gcd(g,c) == 1 for c != b in P, B
618                    }
619                    g = gp;
620                    gp = gcd(b, gp).abs();
621                }
622                if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
623                    B.add(g); // gcd(a,g) == 1 and gcd(g,c) == 1 for c != b in P, B
624                }
625            }
626            if (!b.isZERO() && !b.isConstant() /*&& !B.contains(b)*/) {
627                B.add(b); // gcd(a,b) == 1 and gcd(b,c) == 1 for c != b in P, B
628            }
629        }
630        if (!a.isZERO() && !a.isConstant() /*&& !B.contains(a)*/) {
631            B.add(a);
632        }
633        return B;
634    }
635
636
637    /**
638     * GenPolynomial test for co-prime list.
639     * @param A list of GenPolynomials.
640     * @return true if gcd(b,c) = 1 for all b != c in B, else false.
641     */
642    public boolean isCoPrime(List<GenPolynomial<C>> A) {
643        if (A == null || A.isEmpty()) {
644            return true;
645        }
646        if (A.size() == 1) {
647            return true;
648        }
649        for (int i = 0; i < A.size(); i++) {
650            GenPolynomial<C> a = A.get(i);
651            for (int j = i + 1; j < A.size(); j++) {
652                GenPolynomial<C> b = A.get(j);
653                GenPolynomial<C> g = gcd(a, b);
654                if (!g.isONE()) {
655                    System.out.println("not co-prime, a: " + a);
656                    System.out.println("not co-prime, b: " + b);
657                    System.out.println("not co-prime, g: " + g);
658                    return false;
659                }
660            }
661        }
662        return true;
663    }
664
665
666    /**
667     * GenPolynomial test for co-prime list of given list.
668     * @param A list of GenPolynomials.
669     * @param P list of co-prime GenPolynomials.
670     * @return true if isCoPrime(P) and for all a in A exists p in P with p | a,
671     *         else false.
672     */
673    public boolean isCoPrime(List<GenPolynomial<C>> P, List<GenPolynomial<C>> A) {
674        if (!isCoPrime(P)) {
675            return false;
676        }
677        if (A == null || A.isEmpty()) {
678            return true;
679        }
680        for (GenPolynomial<C> q : A) {
681            if (q.isZERO() || q.isConstant()) {
682                continue;
683            }
684            boolean divides = false;
685            for (GenPolynomial<C> p : P) {
686                GenPolynomial<C> a = PolyUtil.<C> baseSparsePseudoRemainder(q, p);
687                if (a.isZERO()) { // p divides q
688                    divides = true;
689                    break;
690                }
691            }
692            if (!divides) {
693                System.out.println("no divisor for: " + q);
694                return false;
695            }
696        }
697        return true;
698    }
699
700
701    /**
702     * Univariate GenPolynomial extended greatest common divisor. Uses sparse
703     * pseudoRemainder for remainder.
704     * @param P univariate GenPolynomial.
705     * @param S univariate GenPolynomial.
706     * @return [ gcd(P,S), a, b ] with a*P + b*S = gcd(P,S).
707     */
708    public GenPolynomial<C>[] baseExtendedGcd(GenPolynomial<C> P, GenPolynomial<C> S) {
709        //return P.egcd(S);
710        GenPolynomial<C>[] hegcd = baseHalfExtendedGcd(P, S);
711        GenPolynomial<C>[] ret = (GenPolynomial<C>[]) new GenPolynomial[3];
712        ret[0] = hegcd[0];
713        ret[1] = hegcd[1];
714        GenPolynomial<C> x = hegcd[0].subtract(hegcd[1].multiply(P));
715        GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(x, S);
716        // assert qr[1].isZERO() 
717        ret[2] = qr[0];
718        return ret;
719    }
720
721
722    /**
723     * Univariate GenPolynomial half extended greatest comon divisor. Uses
724     * sparse pseudoRemainder for remainder.
725     * @param S GenPolynomial.
726     * @return [ gcd(P,S), a ] with a*P + b*S = gcd(P,S).
727     */
728    public GenPolynomial<C>[] baseHalfExtendedGcd(GenPolynomial<C> P, GenPolynomial<C> S) {
729        //if ( P == null ) {
730        //    throw new IllegalArgumentException("null P not allowed");
731        //}
732        GenPolynomial<C>[] ret = (GenPolynomial<C>[]) new GenPolynomial[2];
733        ret[0] = null;
734        ret[1] = null;
735        if (S == null || S.isZERO()) {
736            ret[0] = P;
737            ret[1] = P.ring.getONE();
738            return ret;
739        }
740        if (P == null || P.isZERO()) {
741            ret[0] = S;
742            ret[1] = S.ring.getZERO();
743            return ret;
744        }
745        if (P.ring.nvar != 1) {
746            throw new IllegalArgumentException(this.getClass().getName() + " not univariate polynomials "
747                            + P.ring);
748        }
749        GenPolynomial<C> q = P;
750        GenPolynomial<C> r = S;
751        GenPolynomial<C> c1 = P.ring.getONE().copy();
752        GenPolynomial<C> d1 = P.ring.getZERO().copy();
753        while (!r.isZERO()) {
754            GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(q, r);
755            //q.divideAndRemainder(r);
756            q = qr[0];
757            GenPolynomial<C> x = c1.subtract(q.multiply(d1));
758            c1 = d1;
759            d1 = x;
760            q = r;
761            r = qr[1];
762        }
763        // normalize ldcf(q) to 1, i.e. make monic
764        C g = q.leadingBaseCoefficient();
765        if (g.isUnit()) {
766            C h = g.inverse();
767            q = q.multiply(h);
768            c1 = c1.multiply(h);
769        }
770        //assert ( ((c1.multiply(P)).remainder(S).equals(q) )); 
771        ret[0] = q;
772        ret[1] = c1;
773        return ret;
774    }
775
776
777    /**
778     * Univariate GenPolynomial greatest common divisor diophantine version.
779     * @param P univariate GenPolynomial.
780     * @param S univariate GenPolynomial.
781     * @param c univariate GenPolynomial.
782     * @return [ a, b ] with a*P + b*S = c and deg(a) < deg(S).
783     */
784    public GenPolynomial<C>[] baseGcdDiophant(GenPolynomial<C> P, GenPolynomial<C> S, GenPolynomial<C> c) {
785        GenPolynomial<C>[] egcd = baseExtendedGcd(P, S);
786        GenPolynomial<C> g = egcd[0];
787        GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(c, g);
788        if (!qr[1].isZERO()) {
789            throw new ArithmeticException("not solvable, r = " + qr[1] + ", c = " + c + ", g = " + g);
790        }
791        GenPolynomial<C> q = qr[0];
792        GenPolynomial<C> a = egcd[1].multiply(q);
793        GenPolynomial<C> b = egcd[2].multiply(q);
794        if (!a.isZERO() && a.degree(0) >= S.degree(0)) {
795            qr = PolyUtil.<C> basePseudoQuotientRemainder(a, S);
796            a = qr[1];
797            b = b.sum(P.multiply(qr[0]));
798        }
799        GenPolynomial<C>[] ret = (GenPolynomial<C>[]) new GenPolynomial[2];
800        ret[0] = a;
801        ret[1] = b;
802
803        if (debug) {
804            GenPolynomial<C> y = ret[0].multiply(P).sum(ret[1].multiply(S));
805            if (!y.equals(c)) {
806                System.out.println("P  = " + P);
807                System.out.println("S  = " + S);
808                System.out.println("c  = " + c);
809                System.out.println("a  = " + a);
810                System.out.println("b  = " + b);
811                System.out.println("y  = " + y);
812                throw new ArithmeticException("not diophant, x = " + y.subtract(c));
813            }
814        }
815        return ret;
816    }
817
818
819    /**
820     * Univariate GenPolynomial partial fraction decomposition.
821     * @param A univariate GenPolynomial.
822     * @param P univariate GenPolynomial.
823     * @param S univariate GenPolynomial.
824     * @return [ A0, Ap, As ] with A/(P*S) = A0 + Ap/P + As/S with deg(Ap) <
825     *         deg(P) and deg(As) < deg(S).
826     */
827    public GenPolynomial<C>[] basePartialFraction(GenPolynomial<C> A, GenPolynomial<C> P, GenPolynomial<C> S) {
828        GenPolynomial<C>[] ret = (GenPolynomial<C>[]) new GenPolynomial[3];
829        ret[0] = null;
830        ret[1] = null;
831        ret[2] = null;
832        GenPolynomial<C> ps = P.multiply(S);
833        GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(A, ps);
834        ret[0] = qr[0];
835        GenPolynomial<C> r = qr[1];
836        GenPolynomial<C>[] diop = baseGcdDiophant(S, P, r); // switch arguments
837
838        //         GenPolynomial<C> x = diop[0].multiply(S).sum( diop[1].multiply(P) );
839        //         if ( !x.equals(r) ) {
840        //             System.out.println("r  = " + r);
841        //             System.out.println("x  = " + x);
842        //             throw new RuntimeException("not partial fraction, x = " + x);
843        //         }
844
845        ret[1] = diop[0];
846        ret[2] = diop[1];
847        if (ret[1].degree(0) >= P.degree(0)) {
848            qr = PolyUtil.<C> basePseudoQuotientRemainder(ret[1], P);
849            ret[0] = ret[0].sum(qr[0]);
850            ret[1] = qr[1];
851        }
852        if (ret[2].degree(0) >= S.degree(0)) {
853            qr = PolyUtil.<C> basePseudoQuotientRemainder(ret[2], S);
854            ret[0] = ret[0].sum(qr[0]);
855            ret[2] = qr[1];
856        }
857        return ret;
858    }
859
860
861    /**
862     * Univariate GenPolynomial partial fraction decomposition.
863     * @param A univariate GenPolynomial.
864     * @param P univariate GenPolynomial.
865     * @param e exponent for P.
866     * @return [ F0, F1, ..., Fe ] with A/(P^e) = sum( Fi / P^i ) with deg(Fi) <
867     *         deg(P).
868     */
869    public List<GenPolynomial<C>> basePartialFraction(GenPolynomial<C> A, GenPolynomial<C> P, int e) {
870        if (A == null || P == null || e == 0) {
871            throw new IllegalArgumentException("null A, P or e = 0 not allowed");
872        }
873        List<GenPolynomial<C>> pf = new ArrayList<GenPolynomial<C>>(e);
874        if (A.isZERO()) {
875            for (int i = 0; i < e; i++) {
876                pf.add(A);
877            }
878            return pf;
879        }
880        if (e == 1) {
881            GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(A, P);
882            pf.add(qr[0]);
883            pf.add(qr[1]);
884            return pf;
885        }
886        GenPolynomial<C> a = A;
887        for (int j = e; j > 0; j--) {
888            GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(a, P);
889            a = qr[0];
890            pf.add(0, qr[1]);
891        }
892        pf.add(0, a);
893        return pf;
894    }
895
896
897    /**
898     * Univariate GenPolynomial partial fraction decomposition.
899     * @param A univariate GenPolynomial.
900     * @param D list of co-prime univariate GenPolynomials.
901     * @return [ A0, A1,..., An ] with A/prod(D) = A0 + sum( Ai/Di ) with
902     *         deg(Ai) < deg(Di).
903     */
904    public List<GenPolynomial<C>> basePartialFraction(GenPolynomial<C> A, List<GenPolynomial<C>> D) {
905        if (D == null || A == null) {
906            throw new IllegalArgumentException("null A or D not allowed");
907        }
908        List<GenPolynomial<C>> pf = new ArrayList<GenPolynomial<C>>(D.size() + 1);
909        if (A.isZERO() || D.size() == 0) {
910            pf.add(A);
911            for (int i = 0; i < D.size(); i++) {
912                pf.add(A);
913            }
914            return pf;
915        }
916        List<GenPolynomial<C>> Dp = new ArrayList<GenPolynomial<C>>(D.size() - 1);
917        GenPolynomial<C> P = A.ring.getONE();
918        GenPolynomial<C> d1 = null;
919        for (GenPolynomial<C> d : D) {
920            if (d1 == null) {
921                d1 = d;
922            } else {
923                P = P.multiply(d);
924                Dp.add(d);
925            }
926        }
927        GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(A, P.multiply(d1));
928        GenPolynomial<C> A0 = qr[0];
929        GenPolynomial<C> r = qr[1];
930        if (D.size() == 1) {
931            pf.add(A0);
932            pf.add(r);
933            return pf;
934        }
935        GenPolynomial<C>[] diop = baseGcdDiophant(P, d1, r); // switch arguments
936        GenPolynomial<C> A1 = diop[0];
937        GenPolynomial<C> S = diop[1];
938        List<GenPolynomial<C>> Fr = basePartialFraction(S, Dp);
939        A0 = A0.sum(Fr.remove(0));
940        pf.add(A0);
941        pf.add(A1);
942        pf.addAll(Fr);
943        return pf;
944    }
945
946
947    /**
948     * Test for Univariate GenPolynomial partial fraction decomposition.
949     * @param A univariate GenPolynomial.
950     * @param D list of (co-prime) univariate GenPolynomials.
951     * @param F list of univariate GenPolynomials from a partial fraction
952     *            computation.
953     * @return true if A/prod(D) = F0 + sum( Fi/Di ) with deg(Fi) < deg(Di), Fi
954     *         in F, else false.
955     */
956    public boolean isBasePartialFraction(GenPolynomial<C> A, List<GenPolynomial<C>> D,
957                    List<GenPolynomial<C>> F) {
958        if (D == null || A == null || F == null) {
959            throw new IllegalArgumentException("null A, F or D not allowed");
960        }
961        if (D.size() != F.size() - 1) {
962            return false;
963        }
964        // A0*prod(D) + sum( Ai * Dip ), Dip = prod(D,j!=i)
965        GenPolynomial<C> P = A.ring.getONE();
966        for (GenPolynomial<C> d : D) {
967            P = P.multiply(d);
968        }
969        List<GenPolynomial<C>> Fp = new ArrayList<GenPolynomial<C>>(F);
970        GenPolynomial<C> A0 = Fp.remove(0).multiply(P);
971        //System.out.println("A0 = " + A0);
972        int j = 0;
973        for (GenPolynomial<C> Fi : Fp) {
974            P = A.ring.getONE();
975            int i = 0;
976            for (GenPolynomial<C> d : D) {
977                if (i != j) {
978                    P = P.multiply(d);
979                }
980                i++;
981            }
982            //System.out.println("Fi = " + Fi);
983            //System.out.println("P  = " + P);
984            A0 = A0.sum(Fi.multiply(P));
985            //System.out.println("A0 = " + A0);
986            j++;
987        }
988        boolean t = A.equals(A0);
989        if (!t) {
990            System.out.println("not isPartFrac = " + A0);
991        }
992        return t;
993    }
994
995
996    /**
997     * Test for Univariate GenPolynomial partial fraction decomposition.
998     * @param A univariate GenPolynomial.
999     * @param P univariate GenPolynomial.
1000     * @param e exponent for P.
1001     * @param F list of univariate GenPolynomials from a partial fraction
1002     *            computation.
1003     * @return true if A/(P^e) = F0 + sum( Fi / P^i ) with deg(Fi) < deg(P), Fi
1004     *         in F, else false.
1005     */
1006    public boolean isBasePartialFraction(GenPolynomial<C> A, GenPolynomial<C> P, int e,
1007                    List<GenPolynomial<C>> F) {
1008        if (A == null || P == null || F == null || e == 0) {
1009            throw new IllegalArgumentException("null A, P, F or e = 0 not allowed");
1010        }
1011        GenPolynomial<C> A0 = basePartialFractionValue(P, e, F);
1012        //         A.ring.getZERO();
1013        //         for ( GenPolynomial<C> Fi : F ) {
1014        //             A0 = A0.multiply(P);
1015        //             A0 = A0.sum(Fi);
1016        //             //System.out.println("A0 = " + A0);
1017        //         }
1018        boolean t = A.equals(A0);
1019        if (!t) {
1020            System.out.println("not isPartFrac = " + A0);
1021        }
1022        return t;
1023    }
1024
1025
1026    /**
1027     * Test for Univariate GenPolynomial partial fraction decomposition.
1028     * @param P univariate GenPolynomial.
1029     * @param e exponent for P.
1030     * @param F list of univariate GenPolynomials from a partial fraction
1031     *            computation.
1032     * @return (F0 + sum( Fi / P^i )) * P^e.
1033     */
1034    public GenPolynomial<C> basePartialFractionValue(GenPolynomial<C> P, int e, List<GenPolynomial<C>> F) {
1035        if (P == null || F == null || e == 0) {
1036            throw new IllegalArgumentException("null P, F or e = 0 not allowed");
1037        }
1038        GenPolynomial<C> A0 = P.ring.getZERO();
1039        for (GenPolynomial<C> Fi : F) {
1040            A0 = A0.multiply(P);
1041            A0 = A0.sum(Fi);
1042            //System.out.println("A0 = " + A0);
1043        }
1044        return A0;
1045    }
1046
1047}