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}