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