001 /*
002 * $Id: GreatestCommonDivisorAbstract.java 3652 2011-06-02 18:17:04Z kredel $
003 */
004
005 package edu.jas.ufd;
006
007
008 import java.util.ArrayList;
009 import java.util.List;
010
011 import org.apache.log4j.Logger;
012
013 import edu.jas.poly.GenPolynomial;
014 import edu.jas.poly.GenPolynomialRing;
015 import edu.jas.poly.PolyUtil;
016 import edu.jas.structure.GcdRingElem;
017 import edu.jas.structure.RingFactory;
018
019
020 /**
021 * Greatest common divisor algorithms.
022 * @author Heinz Kredel
023 */
024
025 public 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(this.getClass().getName() + " 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 * GenPolynomial resultant.
398 * The input polynomials are considered as univariate polynomials in the main variable.
399 * @param P GenPolynomial.
400 * @param S GenPolynomial.
401 * @return res(P,S).
402 * @see edu.jas.ufd.GreatestCommonDivisorSubres#recursiveResultant
403 */
404 public GenPolynomial<C> resultant(GenPolynomial<C> P, GenPolynomial<C> S) {
405 if (S == null || S.isZERO()) {
406 return S;
407 }
408 if (P == null || P.isZERO()) {
409 return P;
410 }
411 // hack for now:
412 GreatestCommonDivisorSubres<C> ufd_sr = new GreatestCommonDivisorSubres<C>();
413 GenPolynomialRing<C> pfac = P.ring;
414 if (pfac.nvar <= 1) {
415 GenPolynomial<C> T = ufd_sr.baseResultant(P, S);
416 return T;
417 }
418 GenPolynomialRing<C> cfac = pfac.contract(1);
419 GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(cfac, 1);
420
421 GenPolynomial<GenPolynomial<C>> Pr = PolyUtil.<C> recursive(rfac, P);
422 GenPolynomial<GenPolynomial<C>> Sr = PolyUtil.<C> recursive(rfac, S);
423
424 GenPolynomial<GenPolynomial<C>> Dr = ufd_sr.recursiveResultant(Pr, Sr);
425 GenPolynomial<C> D = PolyUtil.<C> distribute(pfac, Dr);
426 return D;
427 }
428
429
430 /**
431 * GenPolynomial co-prime list.
432 * @param A list of GenPolynomials.
433 * @return B with gcd(b,c) = 1 for all b != c in B and for all non-constant
434 * a in A there exists b in B with b|a. B does not contain zero or
435 * constant polynomials.
436 */
437 public List<GenPolynomial<C>> coPrime(List<GenPolynomial<C>> A) {
438 if (A == null || A.isEmpty()) {
439 return A;
440 }
441 List<GenPolynomial<C>> B = new ArrayList<GenPolynomial<C>>(A.size());
442 // make a coprime to rest of list
443 GenPolynomial<C> a = A.get(0);
444 //System.out.println("a = " + a);
445 if (!a.isZERO() && !a.isConstant()) {
446 for (int i = 1; i < A.size(); i++) {
447 GenPolynomial<C> b = A.get(i);
448 GenPolynomial<C> g = gcd(a, b).abs();
449 if (!g.isONE()) {
450 a = PolyUtil.<C> basePseudoDivide(a, g);
451 b = PolyUtil.<C> basePseudoDivide(b, g);
452 GenPolynomial<C> gp = gcd(a, g).abs();
453 while (!gp.isONE()) {
454 a = PolyUtil.<C> basePseudoDivide(a, gp);
455 g = PolyUtil.<C> basePseudoDivide(g, gp);
456 B.add(g); // gcd(a,g) == 1
457 g = gp;
458 gp = gcd(a, gp).abs();
459 }
460 if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
461 B.add(g); // gcd(a,g) == 1
462 }
463 }
464 if (!b.isZERO() && !b.isConstant()) {
465 B.add(b); // gcd(a,b) == 1
466 }
467 }
468 } else {
469 B.addAll(A.subList(1, A.size()));
470 }
471 a = a.abs();
472 // make rest coprime
473 B = coPrime(B);
474 //System.out.println("B = " + B);
475 //System.out.println("red(a) = " + a);
476 if (!a.isZERO() && !a.isConstant() /*&& !B.contains(a)*/) {
477 B.add(a);
478 }
479 return B;
480 }
481
482
483 /**
484 * GenPolynomial co-prime list.
485 * @param A list of GenPolynomials.
486 * @return B with gcd(b,c) = 1 for all b != c in B and for all non-constant
487 * a in A there exists b in B with b|a. B does not contain zero or
488 * constant polynomials.
489 */
490 public List<GenPolynomial<C>> coPrimeRec(List<GenPolynomial<C>> A) {
491 if (A == null || A.isEmpty()) {
492 return A;
493 }
494 List<GenPolynomial<C>> B = new ArrayList<GenPolynomial<C>>();
495 // make a co-prime to rest of list
496 for (GenPolynomial<C> a : A) {
497 //System.out.println("a = " + a);
498 B = coPrime(a, B);
499 //System.out.println("B = " + B);
500 }
501 return B;
502 }
503
504
505 /**
506 * GenPolynomial co-prime list.
507 * @param a GenPolynomial.
508 * @param P co-prime list of GenPolynomials.
509 * @return B with gcd(b,c) = 1 for all b != c in B and for non-constant a
510 * there exists b in P with b|a. B does not contain zero or constant
511 * polynomials.
512 */
513 public List<GenPolynomial<C>> coPrime(GenPolynomial<C> a, List<GenPolynomial<C>> P) {
514 if (a == null || a.isZERO() || a.isConstant()) {
515 return P;
516 }
517 List<GenPolynomial<C>> B = new ArrayList<GenPolynomial<C>>(P.size() + 1);
518 // make a coprime to elements of the list P
519 for (int i = 0; i < P.size(); i++) {
520 GenPolynomial<C> b = P.get(i);
521 GenPolynomial<C> g = gcd(a, b).abs();
522 if (!g.isONE()) {
523 a = PolyUtil.<C> basePseudoDivide(a, g);
524 b = PolyUtil.<C> basePseudoDivide(b, g);
525 // make g co-prime to new a, g is co-prime to c != b in P, B
526 GenPolynomial<C> gp = gcd(a, g).abs();
527 while (!gp.isONE()) {
528 a = PolyUtil.<C> basePseudoDivide(a, gp);
529 g = PolyUtil.<C> basePseudoDivide(g, gp);
530 if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
531 B.add(g); // gcd(a,g) == 1 and gcd(g,c) == 1 for c != b in P, B
532 }
533 g = gp;
534 gp = gcd(a, gp).abs();
535 }
536 // make new g co-prime to new b
537 gp = gcd(b, g).abs();
538 while (!gp.isONE()) {
539 b = PolyUtil.<C> basePseudoDivide(b, gp);
540 g = PolyUtil.<C> basePseudoDivide(g, gp);
541 if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
542 B.add(g); // gcd(a,g) == 1 and gcd(g,c) == 1 for c != b in P, B
543 }
544 g = gp;
545 gp = gcd(b, gp).abs();
546 }
547 if (!g.isZERO() && !g.isConstant() /*&& !B.contains(g)*/) {
548 B.add(g); // gcd(a,g) == 1 and gcd(g,c) == 1 for c != b in P, B
549 }
550 }
551 if (!b.isZERO() && !b.isConstant() /*&& !B.contains(b)*/) {
552 B.add(b); // gcd(a,b) == 1 and gcd(b,c) == 1 for c != b in P, B
553 }
554 }
555 if (!a.isZERO() && !a.isConstant() /*&& !B.contains(a)*/) {
556 B.add(a);
557 }
558 return B;
559 }
560
561
562 /**
563 * GenPolynomial test for co-prime list.
564 * @param A list of GenPolynomials.
565 * @return true if gcd(b,c) = 1 for all b != c in B, else false.
566 */
567 public boolean isCoPrime(List<GenPolynomial<C>> A) {
568 if (A == null || A.isEmpty()) {
569 return true;
570 }
571 if (A.size() == 1) {
572 return true;
573 }
574 for (int i = 0; i < A.size(); i++) {
575 GenPolynomial<C> a = A.get(i);
576 for (int j = i + 1; j < A.size(); j++) {
577 GenPolynomial<C> b = A.get(j);
578 GenPolynomial<C> g = gcd(a, b);
579 if (!g.isONE()) {
580 System.out.println("not co-prime, a: " + a);
581 System.out.println("not co-prime, b: " + b);
582 System.out.println("not co-prime, g: " + g);
583 return false;
584 }
585 }
586 }
587 return true;
588 }
589
590
591 /**
592 * GenPolynomial test for co-prime list of given list.
593 * @param A list of GenPolynomials.
594 * @param P list of co-prime GenPolynomials.
595 * @return true if isCoPrime(P) and for all a in A exists p in P with p | a,
596 * else false.
597 */
598 public boolean isCoPrime(List<GenPolynomial<C>> P, List<GenPolynomial<C>> A) {
599 if (!isCoPrime(P)) {
600 return false;
601 }
602 if (A == null || A.isEmpty()) {
603 return true;
604 }
605 for (GenPolynomial<C> q : A) {
606 if (q.isZERO() || q.isConstant()) {
607 continue;
608 }
609 boolean divides = false;
610 for (GenPolynomial<C> p : P) {
611 GenPolynomial<C> a = PolyUtil.<C> basePseudoRemainder(q, p);
612 if (a.isZERO()) { // p divides q
613 divides = true;
614 break;
615 }
616 }
617 if (!divides) {
618 System.out.println("no divisor for: " + q);
619 return false;
620 }
621 }
622 return true;
623 }
624
625
626 /**
627 * Univariate GenPolynomial extended greatest common divisor. Uses sparse
628 * pseudoRemainder for remainder.
629 * @param P univariate GenPolynomial.
630 * @param S univariate GenPolynomial.
631 * @return [ gcd(P,S), a, b ] with a*P + b*S = gcd(P,S).
632 */
633 public GenPolynomial<C>[] baseExtendedGcd(GenPolynomial<C> P, GenPolynomial<C> S) {
634 //return P.egcd(S);
635 GenPolynomial<C>[] hegcd = baseHalfExtendedGcd(P,S);
636 GenPolynomial<C>[] ret = (GenPolynomial<C>[]) new GenPolynomial[3];
637 ret[0] = hegcd[0];
638 ret[1] = hegcd[1];
639 GenPolynomial<C> x = hegcd[0].subtract( hegcd[1].multiply(P) );
640 GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(x, S);
641 // assert qr[1].isZERO()
642 ret[2] = qr[0];
643 return ret;
644 }
645
646
647 /**
648 * Univariate GenPolynomial half extended greatest comon divisor.
649 * Uses sparse pseudoRemainder for remainder.
650 * @param S GenPolynomial.
651 * @return [ gcd(P,S), a ] with a*P + b*S = gcd(P,S).
652 */
653 public GenPolynomial<C>[] baseHalfExtendedGcd(GenPolynomial<C> P, GenPolynomial<C> S) {
654 //if ( P == null ) {
655 // throw new IllegalArgumentException("null P not allowed");
656 //}
657 GenPolynomial<C>[] ret = (GenPolynomial<C>[]) new GenPolynomial[2];
658 ret[0] = null;
659 ret[1] = null;
660 if ( S == null || S.isZERO() ) {
661 ret[0] = P;
662 ret[1] = P.ring.getONE();
663 return ret;
664 }
665 if ( P == null || P.isZERO() ) {
666 ret[0] = S;
667 ret[1] = S.ring.getZERO();
668 return ret;
669 }
670 if ( P.ring.nvar != 1 ) {
671 throw new IllegalArgumentException(this.getClass().getName()
672 + " not univariate polynomials " + P.ring);
673 }
674 GenPolynomial<C> q = P;
675 GenPolynomial<C> r = S;
676 GenPolynomial<C> c1 = P.ring.getONE().clone();
677 GenPolynomial<C> d1 = P.ring.getZERO().clone();
678 while ( !r.isZERO() ) {
679 GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(q, r);
680 //q.divideAndRemainder(r);
681 q = qr[0];
682 GenPolynomial<C> x = c1.subtract( q.multiply(d1) );
683 c1 = d1;
684 d1 = x;
685 q = r;
686 r = qr[1];
687 }
688 // normalize ldcf(q) to 1, i.e. make monic
689 C g = q.leadingBaseCoefficient();
690 if ( g.isUnit() ) {
691 C h = g.inverse();
692 q = q.multiply( h );
693 c1 = c1.multiply( h );
694 }
695 //assert ( ((c1.multiply(P)).remainder(S).equals(q) ));
696 ret[0] = q;
697 ret[1] = c1;
698 return ret;
699 }
700
701
702 /**
703 * Univariate GenPolynomial greatest common divisor diophantine version.
704 * @param P univariate GenPolynomial.
705 * @param S univariate GenPolynomial.
706 * @param c univariate GenPolynomial.
707 * @return [ a, b ] with a*P + b*S = c and deg(a) < deg(S).
708 */
709 public GenPolynomial<C>[] baseGcdDiophant(GenPolynomial<C> P, GenPolynomial<C> S, GenPolynomial<C> c) {
710 GenPolynomial<C>[] egcd = baseExtendedGcd(P,S);
711 GenPolynomial<C> g = egcd[0];
712 GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(c, g);
713 if ( !qr[1].isZERO() ) {
714 throw new ArithmeticException("not solvable, r = " + qr[1] + ", c = " + c + ", g = " + g);
715 }
716 GenPolynomial<C> q = qr[0];
717 GenPolynomial<C> a = egcd[1].multiply(q);
718 GenPolynomial<C> b = egcd[2].multiply(q);
719 if ( !a.isZERO() && a.degree(0) >= S.degree(0) ) {
720 qr = PolyUtil.<C> basePseudoQuotientRemainder(a, S);
721 a = qr[1];
722 b = b.sum( P.multiply( qr[0] ) );
723 }
724 GenPolynomial<C>[] ret = (GenPolynomial<C>[]) new GenPolynomial[2];
725 ret[0] = a;
726 ret[1] = b;
727
728 if ( true ) {
729 return ret;
730 }
731
732 GenPolynomial<C> y = ret[0].multiply(P).sum( ret[1].multiply(S) );
733 if ( !y.equals(c) ) {
734 System.out.println("P = " + P);
735 System.out.println("S = " + S);
736 System.out.println("c = " + c);
737 System.out.println("a = " + a);
738 System.out.println("b = " + b);
739 System.out.println("y = " + y);
740 throw new ArithmeticException("not diophant, x = " + y.subtract(c));
741 }
742
743 return ret;
744 }
745
746
747 /**
748 * Univariate GenPolynomial partial fraction decomposition.
749 * @param A univariate GenPolynomial.
750 * @param P univariate GenPolynomial.
751 * @param S univariate GenPolynomial.
752 * @return [ A0, Ap, As ] with A/(P*S) = A0 + Ap/P + As/S with deg(Ap) < deg(P) and deg(As) < deg(S).
753 */
754 public GenPolynomial<C>[] basePartialFraction(GenPolynomial<C> A, GenPolynomial<C> P, GenPolynomial<C> S) {
755 GenPolynomial<C>[] ret = (GenPolynomial<C>[]) new GenPolynomial[3];
756 ret[0] = null;
757 ret[1] = null;
758 ret[2] = null;
759 GenPolynomial<C> ps = P.multiply(S);
760 GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(A, ps);
761 ret[0] = qr[0];
762 GenPolynomial<C> r = qr[1];
763 GenPolynomial<C>[] diop = baseGcdDiophant(S,P,r); // switch arguments
764
765 // GenPolynomial<C> x = diop[0].multiply(S).sum( diop[1].multiply(P) );
766 // if ( !x.equals(r) ) {
767 // System.out.println("r = " + r);
768 // System.out.println("x = " + x);
769 // throw new RuntimeException("not partial fraction, x = " + x);
770 // }
771
772 ret[1] = diop[0];
773 ret[2] = diop[1];
774 if ( ret[1].degree(0) >= P.degree(0) ) {
775 qr = PolyUtil.<C> basePseudoQuotientRemainder(ret[1], P);
776 ret[0] = ret[0].sum( qr[0] );
777 ret[1] = qr[1];
778 }
779 if ( ret[2].degree(0) >= S.degree(0) ) {
780 qr = PolyUtil.<C> basePseudoQuotientRemainder(ret[2], S);
781 ret[0] = ret[0].sum( qr[0] );
782 ret[2] = qr[1];
783 }
784 return ret;
785 }
786
787
788 /**
789 * Univariate GenPolynomial partial fraction decomposition.
790 * @param A univariate GenPolynomial.
791 * @param P univariate GenPolynomial.
792 * @param e exponent for P.
793 * @return [ F0, F1, ..., Fe ] with A/(P^e) = sum( Fi / P^i ) with deg(Fi) < deg(P).
794 */
795 public List<GenPolynomial<C>> basePartialFraction(GenPolynomial<C> A, GenPolynomial<C> P, int e) {
796 if ( A == null || P == null || e == 0 ) {
797 throw new IllegalArgumentException("null A, P or e = 0 not allowed");
798 }
799 List<GenPolynomial<C>> pf = new ArrayList<GenPolynomial<C>>( e );
800 if ( A.isZERO() ) {
801 for ( int i = 0; i < e; i++ ) {
802 pf.add(A);
803 }
804 return pf;
805 }
806 if ( e == 1 ) {
807 GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(A, P);
808 pf.add(qr[0]);
809 pf.add(qr[1]);
810 return pf;
811 }
812 GenPolynomial<C> a = A;
813 for ( int j = e; j > 0; j-- ) {
814 GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(a, P);
815 a = qr[0];
816 pf.add(0,qr[1]);
817 }
818 pf.add(0,a);
819 return pf;
820 }
821
822
823 /**
824 * Univariate GenPolynomial partial fraction decomposition.
825 * @param A univariate GenPolynomial.
826 * @param D list of co-prime univariate GenPolynomials.
827 * @return [ A0, A1,..., An ] with A/prod(D) = A0 + sum( Ai/Di ) with deg(Ai) < deg(Di).
828 */
829 public List<GenPolynomial<C>> basePartialFraction(GenPolynomial<C> A, List<GenPolynomial<C>> D) {
830 if ( D == null || A == null ) {
831 throw new IllegalArgumentException("null A or D not allowed");
832 }
833 List<GenPolynomial<C>> pf = new ArrayList<GenPolynomial<C>>( D.size()+1 );
834 if ( A.isZERO() || D.size() == 0 ) {
835 pf.add(A);
836 for ( int i = 0; i < D.size(); i++ ) {
837 pf.add(A);
838 }
839 return pf;
840 }
841 List<GenPolynomial<C>> Dp = new ArrayList<GenPolynomial<C>>( D.size()-1 );
842 GenPolynomial<C> P = A.ring.getONE();
843 GenPolynomial<C> d1 = null;
844 for ( GenPolynomial<C> d : D ) {
845 if ( d1 == null ) {
846 d1 = d;
847 } else {
848 P = P.multiply(d);
849 Dp.add(d);
850 }
851 }
852 GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(A, P.multiply(d1));
853 GenPolynomial<C> A0 = qr[0];
854 GenPolynomial<C> r = qr[1];
855 if ( D.size() == 1 ) {
856 pf.add(A0);
857 pf.add(r);
858 return pf;
859 }
860 GenPolynomial<C>[] diop = baseGcdDiophant(P,d1,r); // switch arguments
861 GenPolynomial<C> A1 = diop[0];
862 GenPolynomial<C> S = diop[1];
863 List<GenPolynomial<C>> Fr = basePartialFraction(S,Dp);
864 A0 = A0.sum( Fr.remove(0) );
865 pf.add(A0);
866 pf.add(A1);
867 pf.addAll(Fr);
868 return pf;
869 }
870
871
872 /**
873 * Test for Univariate GenPolynomial partial fraction decomposition.
874 * @param A univariate GenPolynomial.
875 * @param D list of (co-prime) univariate GenPolynomials.
876 * @param F list of univariate GenPolynomials from a partial fraction computation.
877 * @return true if A/prod(D) = F0 + sum( Fi/Di ) with deg(Fi) < deg(Di), Fi in F,
878 else false.
879 */
880 public boolean isBasePartialFraction(GenPolynomial<C> A, List<GenPolynomial<C>> D, List<GenPolynomial<C>> F) {
881 if ( D == null || A == null || F == null ) {
882 throw new IllegalArgumentException("null A, F or D not allowed");
883 }
884 if ( D.size() != F.size()-1 ) {
885 return false;
886 }
887 // A0*prod(D) + sum( Ai * Dip ), Dip = prod(D,j!=i)
888 GenPolynomial<C> P = A.ring.getONE();
889 for ( GenPolynomial<C> d : D ) {
890 P = P.multiply(d);
891 }
892 List<GenPolynomial<C>> Fp = new ArrayList<GenPolynomial<C>>( F );
893 GenPolynomial<C> A0 = Fp.remove(0).multiply(P);
894 //System.out.println("A0 = " + A0);
895 int j = 0;
896 for ( GenPolynomial<C> Fi : Fp ) {
897 P = A.ring.getONE();
898 int i = 0;
899 for ( GenPolynomial<C> d : D ) {
900 if ( i != j ) {
901 P = P.multiply(d);
902 }
903 i++;
904 }
905 //System.out.println("Fi = " + Fi);
906 //System.out.println("P = " + P);
907 A0 = A0.sum( Fi.multiply(P) );
908 //System.out.println("A0 = " + A0);
909 j++;
910 }
911 boolean t = A.equals(A0);
912 if ( ! t ) {
913 System.out.println("not isPartFrac = " + A0);
914 }
915 return t;
916 }
917
918
919 /**
920 * Test for Univariate GenPolynomial partial fraction decomposition.
921 * @param A univariate GenPolynomial.
922 * @param P univariate GenPolynomial.
923 * @param e exponent for P.
924 * @param F list of univariate GenPolynomials from a partial fraction computation.
925 * @return true if A/(P^e) = F0 + sum( Fi / P^i ) with deg(Fi) < deg(P), Fi in F,
926 else false.
927 */
928 public boolean isBasePartialFraction(GenPolynomial<C> A, GenPolynomial<C> P, int e, List<GenPolynomial<C>> F) {
929 if ( A == null || P == null || F == null || e == 0 ) {
930 throw new IllegalArgumentException("null A, P, F or e = 0 not allowed");
931 }
932 GenPolynomial<C> A0 = basePartialFractionValue(P,e,F);
933 // A.ring.getZERO();
934 // for ( GenPolynomial<C> Fi : F ) {
935 // A0 = A0.multiply(P);
936 // A0 = A0.sum(Fi);
937 // //System.out.println("A0 = " + A0);
938 // }
939 boolean t = A.equals(A0);
940 if ( ! t ) {
941 System.out.println("not isPartFrac = " + A0);
942 }
943 return t;
944 }
945
946
947 /**
948 * Test for Univariate GenPolynomial partial fraction decomposition.
949 * @param P univariate GenPolynomial.
950 * @param e exponent for P.
951 * @param F list of univariate GenPolynomials from a partial fraction computation.
952 * @return (F0 + sum( Fi / P^i )) * P^e.
953 */
954 public GenPolynomial<C> basePartialFractionValue(GenPolynomial<C> P, int e, List<GenPolynomial<C>> F) {
955 if ( P == null || F == null || e == 0 ) {
956 throw new IllegalArgumentException("null P, F or e = 0 not allowed");
957 }
958 GenPolynomial<C> A0 = P.ring.getZERO();
959 for ( GenPolynomial<C> Fi : F ) {
960 A0 = A0.multiply(P);
961 A0 = A0.sum(Fi);
962 //System.out.println("A0 = " + A0);
963 }
964 return A0;
965 }
966
967 }