Source code

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