001/* 002 * $Id: GreatestCommonDivisorSimple.java 5105 2015-02-07 14:35:17Z kredel $ 003 */ 004 005package edu.jas.fd; 006 007 008import org.apache.log4j.Logger; 009 010import edu.jas.poly.GenPolynomial; 011import edu.jas.poly.GenSolvablePolynomial; 012import edu.jas.poly.PolyUtil; 013import edu.jas.structure.GcdRingElem; 014import edu.jas.structure.RingFactory; 015 016 017/** 018 * (Non-unique) factorization domain greatest common divisor common algorithms 019 * with monic polynomial remainder sequence. If C is a field, then the monic PRS 020 * (on coefficients) is computed otherwise no simplifications in the reduction 021 * are made. 022 * @param <C> coefficient type 023 * @author Heinz Kredel 024 */ 025 026public class GreatestCommonDivisorSimple<C extends GcdRingElem<C>> extends GreatestCommonDivisorAbstract<C> { 027 028 029 private static final Logger logger = Logger.getLogger(GreatestCommonDivisorSimple.class); 030 031 032 private final boolean debug = true; //logger.isDebugEnabled(); 033 034 035 /** 036 * Constructor. 037 * @param cf coefficient ring. 038 */ 039 public GreatestCommonDivisorSimple(RingFactory<C> cf) { 040 super(cf); 041 } 042 043 044 /** 045 * Univariate GenSolvablePolynomial greatest common divisor. Uses 046 * pseudoRemainder for remainder. 047 * @param P univariate GenSolvablePolynomial. 048 * @param S univariate GenSolvablePolynomial. 049 * @return gcd(P,S) with P = P'*gcd(P,S) and S = S'*gcd(P,S). 050 */ 051 @Override 052 public GenSolvablePolynomial<C> leftBaseGcd(GenSolvablePolynomial<C> P, GenSolvablePolynomial<C> S) { 053 if (S == null || S.isZERO()) { 054 return P; 055 } 056 if (P == null || P.isZERO()) { 057 return S; 058 } 059 if (P.ring.nvar > 1) { 060 throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial"); 061 } 062 boolean field = P.ring.coFac.isField(); 063 long e = P.degree(0); 064 long f = S.degree(0); 065 GenSolvablePolynomial<C> q; 066 GenSolvablePolynomial<C> r; 067 if (f > e) { 068 r = P; 069 q = S; 070 long g = f; 071 f = e; 072 e = g; 073 } else { 074 q = P; 075 r = S; 076 } 077 if (debug) { 078 logger.debug("degrees: e = " + e + ", f = " + f); 079 } 080 C c; 081 if (field) { 082 r = r.monic(); 083 q = q.monic(); 084 c = P.ring.getONECoefficient(); 085 } else { 086 r = (GenSolvablePolynomial<C>) r.abs(); 087 q = (GenSolvablePolynomial<C>) q.abs(); 088 C a = rightBaseContent(r); 089 C b = rightBaseContent(q); 090 r = divide(r, a); // indirection 091 q = divide(q, b); // indirection 092 c = gcd(a, b); // indirection 093 } 094 //System.out.println("baseCont: gcd(cont) = " + b); 095 if (r.isONE()) { 096 return r.multiply(c); 097 } 098 if (q.isONE()) { 099 return q.multiply(c); 100 } 101 GenSolvablePolynomial<C> x; 102 //System.out.println("baseGCD: q = " + q); 103 //System.out.println("baseGCD: r = " + r); 104 while (!r.isZERO()) { 105 x = FDUtil.<C> leftBaseSparsePseudoRemainder(q, r); 106 q = r; 107 if (field) { 108 r = x.monic(); 109 } else { 110 r = x; 111 } 112 //System.out.println("baseGCD: q = " + q); 113 //System.out.println("baseGCD: r = " + r); 114 } 115 q = leftBasePrimitivePart(q); 116 return (GenSolvablePolynomial<C>) (q.multiply(c)).abs(); 117 } 118 119 120 /** 121 * Univariate GenSolvablePolynomial right greatest common divisor. Uses 122 * pseudoRemainder for remainder. 123 * @param P univariate GenSolvablePolynomial. 124 * @param S univariate GenSolvablePolynomial. 125 * @return gcd(P,S) with P = gcd(P,S)*P' and S = gcd(P,S)*S'. 126 */ 127 @Override 128 public GenSolvablePolynomial<C> rightBaseGcd(GenSolvablePolynomial<C> P, GenSolvablePolynomial<C> S) { 129 if (S == null || S.isZERO()) { 130 return P; 131 } 132 if (P == null || P.isZERO()) { 133 return S; 134 } 135 if (P.ring.nvar > 1) { 136 throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial"); 137 } 138 boolean field = P.ring.coFac.isField(); 139 long e = P.degree(0); 140 long f = S.degree(0); 141 GenSolvablePolynomial<C> q; 142 GenSolvablePolynomial<C> r; 143 if (f > e) { 144 r = P; 145 q = S; 146 long g = f; 147 f = e; 148 e = g; 149 } else { 150 q = P; 151 r = S; 152 } 153 if (debug) { 154 logger.debug("degrees: e = " + e + ", f = " + f); 155 } 156 C c; 157 if (field) { 158 r = r.monic(); 159 q = q.monic(); 160 c = P.ring.getONECoefficient(); 161 } else { 162 r = (GenSolvablePolynomial<C>) r.abs(); 163 q = (GenSolvablePolynomial<C>) q.abs(); 164 C a = leftBaseContent(r); 165 C b = leftBaseContent(q); 166 r = divide(r, a); // indirection 167 q = divide(q, b); // indirection 168 c = gcd(a, b); // indirection 169 } 170 //System.out.println("baseCont: gcd(cont) = " + b); 171 if (r.isONE()) { 172 return r.multiply(c); 173 } 174 if (q.isONE()) { 175 return q.multiply(c); 176 } 177 GenSolvablePolynomial<C> x; 178 //System.out.println("baseGCD: q = " + q); 179 //System.out.println("baseGCD: r = " + r); 180 while (!r.isZERO()) { 181 x = FDUtil.<C> rightBaseSparsePseudoRemainder(q, r); 182 q = r; 183 if (field) { 184 r = x.monic(); 185 } else { 186 r = x; 187 } 188 //System.out.println("baseGCD: q = " + q); 189 //System.out.println("baseGCD: r = " + r); 190 } 191 q = rightBasePrimitivePart(q); 192 return (GenSolvablePolynomial<C>) (q.multiplyLeft(c)).abs(); 193 } 194 195 196 /** 197 * Univariate GenSolvablePolynomial left recursive greatest common divisor. 198 * Uses pseudoRemainder for remainder. 199 * @param P univariate recursive GenSolvablePolynomial. 200 * @param S univariate recursive GenSolvablePolynomial. 201 * @return gcd(P,S) with P = P'*gcd(P,S)*p and S = S'*gcd(P,S)*s, where 202 * deg_main(p) = deg_main(s) == 0. 203 */ 204 @Override 205 public GenSolvablePolynomial<GenPolynomial<C>> leftRecursiveUnivariateGcd( 206 GenSolvablePolynomial<GenPolynomial<C>> P, GenSolvablePolynomial<GenPolynomial<C>> S) { 207 if (S == null || S.isZERO()) { 208 return P; 209 } 210 if (P == null || P.isZERO()) { 211 return S; 212 } 213 if (P.ring.nvar > 1) { 214 throw new IllegalArgumentException("no univariate polynomial"); 215 } 216 boolean field = P.leadingBaseCoefficient().ring.coFac.isField(); 217 long e = P.degree(0); 218 long f = S.degree(0); 219 GenSolvablePolynomial<GenPolynomial<C>> q, r, x, qs, rs, qp, rp; 220 if (f > e) { 221 r = P; 222 q = S; 223 long g = f; 224 f = e; 225 e = g; 226 } else if (f < e) { 227 q = P; 228 r = S; 229 } else { // f == e 230 if (P.leadingBaseCoefficient().degree() > S.leadingBaseCoefficient().degree()) { 231 q = P; 232 r = S; 233 } else { 234 r = P; 235 q = S; 236 } 237 } 238 if (debug) { 239 logger.debug("degrees: e = " + e + ", f = " + f); 240 } 241 if (field) { 242 r = PolyUtil.<C> monic(r); 243 q = PolyUtil.<C> monic(q); 244 } else { 245 r = (GenSolvablePolynomial<GenPolynomial<C>>) r.abs(); 246 q = (GenSolvablePolynomial<GenPolynomial<C>>) q.abs(); 247 } 248 GenSolvablePolynomial<C> a = rightRecursiveContent(r); 249 rs = FDUtil.<C> recursiveDivideRightEval(r, a); 250 if (debug) { 251 logger.info("recCont a = " + a + ", r = " + r); 252 logger.info("recCont r/a = " + rs + ", r%a = " + r.subtract(rs.multiply(a))); 253 if (!r.equals(rs.multiply(a))) { 254 if (!rs.multiplyLeft(a).equals(r)) { 255 System.out.println("recGcd, r = " + r); 256 System.out.println("recGcd, cont(r) = " + a); 257 System.out.println("recGcd, pp(r) = " + rs); 258 System.out.println("recGcd, pp(r)c(r) = " + rs.multiply(a)); 259 System.out.println("recGcd, c(r)pp(r) = " + rs.multiplyLeft(a)); 260 throw new RuntimeException("recGcd, pp: not divisible"); 261 } 262 } 263 } 264 r = rs; 265 GenSolvablePolynomial<C> b = rightRecursiveContent(q); 266 qs = FDUtil.<C> recursiveDivideRightEval(q, b); 267 if (debug) { 268 logger.info("recCont b = " + b + ", q = " + q); 269 logger.info("recCont q/b = " + qs + ", q%b = " + q.subtract(qs.multiply(b))); 270 if (!q.equals(qs.multiply(b))) { 271 if (!qs.multiplyLeft(b).equals(q)) { 272 System.out.println("recGcd, q = " + q); 273 System.out.println("recGcd, cont(q) = " + b); 274 System.out.println("recGcd, pp(q) = " + qs); 275 System.out.println("recGcd, pp(q)c(q) = " + qs.multiply(b)); 276 System.out.println("recGcd, c(q)pp(q) = " + qs.multiplyLeft(b)); 277 throw new RuntimeException("recGcd, pp: not divisible"); 278 } 279 } 280 } 281 q = qs; 282 logger.info("Gcd(content).ring = " + a.ring.toScript() + ", a = " + a + ", b = " + b); 283 //no: GenSolvablePolynomial<C> c = leftGcd(a, b); // go to recursion 284 GenSolvablePolynomial<C> c = rightGcd(a, b); // go to recursion 285 logger.info("Gcd(contents) c = " + c + ", a = " + a + ", b = " + b); 286 if (r.isONE()) { 287 return r.multiply(c); 288 } 289 if (q.isONE()) { 290 return q.multiply(c); 291 } 292 rs = r; 293 qs = q; 294 if (debug) { 295 logger.info("r.ring = " + r.ring.toScript()); 296 logger.info("gcd-loop, start: q = " + q + ", r = " + r); 297 } 298 while (!r.isZERO()) { 299 x = FDUtil.<C> recursiveSparsePseudoRemainder(q, r); 300 q = r; 301 if (field) { 302 r = PolyUtil.<C> monic(x); 303 } else { 304 r = x; 305 } 306 if (debug) { 307 logger.info("gcd-loop, rem: q = " + q + ", r = " + r); 308 } 309 } 310 if (debug) { 311 logger.info("gcd(div) = " + q + ", rs = " + rs + ", qs = " + qs); 312 rp = FDUtil.<C> recursiveSparsePseudoRemainder(rs, q); 313 qp = FDUtil.<C> recursiveSparsePseudoRemainder(qs, q); 314 if (!qp.isZERO() || !rp.isZERO()) { 315 logger.info("gcd(div): rem(r,g) = " + rp + ", rem(q,g) = " + qp); 316 rp = FDUtil.<C> recursivePseudoQuotient(rs, q); 317 qp = FDUtil.<C> recursivePseudoQuotient(qs, q); 318 logger.info("gcd(div): r/g = " + rp + ", q/g = " + qp); 319 throw new RuntimeException("recGcd, div: not divisible"); 320 } 321 } 322 qp = q; 323 if (false) { // not checkable 324 System.out.println("recGcd, qs = " + qs + ", P = " + P); 325 System.out.println("recGcd, rs = " + rs + ", S = " + S); 326 qs = FDUtil.<C> recursiveSparsePseudoRemainder(P, qp); // P 327 rs = FDUtil.<C> recursiveSparsePseudoRemainder(S, qp); // S 328 if (!qs.isZERO() || !rs.isZERO()) { 329 System.out.println("recGcd, P = " + P); 330 System.out.println("recGcd, S = " + S); 331 System.out.println("recGcd, qp = " + qp); 332 System.out.println("recGcd, qs = " + qs); 333 System.out.println("recGcd, rs = " + rs); 334 rs = FDUtil.<C> recursiveRightSparsePseudoRemainder(S, qp); 335 System.out.println("recGcd, rs = " + rs); 336 if (!rs.isZERO()) { 337 throw new RuntimeException("recGcd, pre pp: not divisible"); 338 } 339 } 340 logger.info("recGcd(P,S) okay: qp = " + qp); 341 } 342 q = leftRecursivePrimitivePart(q); 343 if (!qp.equals(q)) { 344 logger.info("gcd(pp) = " + q + ", qp = " + qp); 345 } 346 q = (GenSolvablePolynomial<GenPolynomial<C>>) q.multiply(c).abs(); 347 if (false) { // not checkable 348 logger.info("gcd(pp) = " + q + ", c = " + c); 349 qs = FDUtil.<C> recursiveSparsePseudoRemainder(P, q); 350 rs = FDUtil.<C> recursiveSparsePseudoRemainder(S, q); 351 if (!qs.isZERO() || !rs.isZERO()) { 352 System.out.println("recGcd-pp, P = " + P); 353 System.out.println("recGcd-pp, S = " + S); 354 System.out.println("recGcd-pp, q = " + q); 355 System.out.println("recGcd-pp, qs = " + qs); 356 System.out.println("recGcd-pp, rs = " + rs); 357 rs = FDUtil.<C> recursiveRightSparsePseudoRemainder(S, q); 358 System.out.println("recGcd-pp, rs = " + rs); 359 if (!rs.isZERO()) { 360 throw new RuntimeException("recGcd: not divisible"); 361 } 362 } 363 logger.info("recGcd(P,S) okay: q = " + q); 364 } 365 return q; 366 } 367 368 369 /** 370 * Univariate GenSolvablePolynomial right recursive greatest common divisor. 371 * Uses pseudoRemainder for remainder. 372 * @param P univariate recursive GenSolvablePolynomial. 373 * @param S univariate recursive GenSolvablePolynomial. 374 * @return gcd(P,S) with P = p*gcd(P,S)*P' and S = s*gcd(P,S)*S', where 375 * deg_main(p) = deg_main(s) == 0. 376 */ 377 @Override 378 public GenSolvablePolynomial<GenPolynomial<C>> rightRecursiveUnivariateGcd( 379 GenSolvablePolynomial<GenPolynomial<C>> P, GenSolvablePolynomial<GenPolynomial<C>> S) { 380 if (S == null || S.isZERO()) { 381 return P; 382 } 383 if (P == null || P.isZERO()) { 384 return S; 385 } 386 if (P.ring.nvar > 1) { 387 throw new IllegalArgumentException("no univariate polynomial"); 388 } 389 boolean field = P.leadingBaseCoefficient().ring.coFac.isField(); 390 long e = P.degree(0); 391 long f = S.degree(0); 392 GenSolvablePolynomial<GenPolynomial<C>> q, r, x, qs, rs, qp, rp; 393 if (f > e) { 394 r = P; 395 q = S; 396 long g = f; 397 f = e; 398 e = g; 399 } else if (f < e) { 400 q = P; 401 r = S; 402 } else { // f == e 403 if (P.leadingBaseCoefficient().degree() > S.leadingBaseCoefficient().degree()) { 404 q = P; 405 r = S; 406 } else { 407 r = P; 408 q = S; 409 } 410 } 411 if (debug) { 412 logger.debug("RI-degrees: e = " + e + ", f = " + f); 413 } 414 if (field) { 415 r = PolyUtil.<C> monic(r); 416 q = PolyUtil.<C> monic(q); 417 } else { 418 r = (GenSolvablePolynomial<GenPolynomial<C>>) r.abs(); 419 q = (GenSolvablePolynomial<GenPolynomial<C>>) q.abs(); 420 } 421 GenSolvablePolynomial<C> a = leftRecursiveContent(r); 422 rs = FDUtil.<C> recursiveDivide(r, a); 423 if (debug) { 424 logger.info("RI-recCont a = " + a + ", r = " + r); 425 logger.info("RI-recCont r/a = " + r + ", r%a = " + r.subtract(rs.multiplyLeft(a))); 426 if (!r.equals(rs.multiplyLeft(a))) { 427 System.out.println("RI-recGcd, r = " + r); 428 System.out.println("RI-recGcd, cont(r) = " + a); 429 System.out.println("RI-recGcd, pp(r) = " + rs); 430 System.out.println("RI-recGcd, pp(r)c(r) = " + rs.multiplyLeft(a)); 431 throw new RuntimeException("RI-recGcd, pp: not divisible"); 432 } 433 } 434 r = rs; 435 GenSolvablePolynomial<C> b = leftRecursiveContent(q); 436 qs = FDUtil.<C> recursiveDivide(q, b); 437 if (debug) { 438 logger.info("RI-recCont b = " + b + ", q = " + q); 439 logger.info("RI-recCont q/b = " + qs + ", q%b = " + q.subtract(qs.multiplyLeft(b))); 440 if (!q.equals(qs.multiplyLeft(b))) { 441 System.out.println("RI-recGcd, q = " + q); 442 System.out.println("RI-recGcd, cont(q) = " + b); 443 System.out.println("RI-recGcd, pp(q) = " + qs); 444 System.out.println("RI-recGcd, pp(q)c(q) = " + qs.multiplyLeft(b)); 445 throw new RuntimeException("RI-recGcd, pp: not divisible"); 446 } 447 } 448 q = qs; 449 //no: GenSolvablePolynomial<C> c = rightGcd(a, b); // go to recursion 450 GenSolvablePolynomial<C> c = leftGcd(a, b); // go to recursion 451 logger.info("RI-Gcd(contents) c = " + c + ", a = " + a + ", b = " + b); 452 if (r.isONE()) { 453 return r.multiplyLeft(c); 454 } 455 if (q.isONE()) { 456 return q.multiplyLeft(c); 457 } 458 rs = r; 459 qs = q; 460 if (debug) { 461 logger.info("RI-r.ring = " + r.ring.toScript()); 462 logger.info("RI-gcd-loop, start: q = " + q + ", r = " + r); 463 } 464 while (!r.isZERO()) { 465 x = FDUtil.<C> recursiveRightSparsePseudoRemainder(q, r); 466 q = r; 467 if (field) { 468 r = PolyUtil.<C> monic(x); 469 } else { 470 r = x; 471 } 472 if (debug) { 473 logger.info("RI-gcd-loop, rem: q = " + q + ", r = " + r); 474 } 475 } 476 if (debug) { 477 logger.info("RI-gcd(div) = " + q + ", rs = " + rs + ", qs = " + qs); 478 rp = FDUtil.<C> recursiveRightSparsePseudoRemainder(rs, q); 479 qp = FDUtil.<C> recursiveRightSparsePseudoRemainder(qs, q); 480 if (!qp.isZERO() || !rp.isZERO()) { 481 logger.info("RI-gcd(div): rem(r,g) = " + rp + ", rem(q,g) = " + qp); 482 rp = FDUtil.<C> recursivePseudoQuotient(rs, q); 483 qp = FDUtil.<C> recursivePseudoQuotient(qs, q); 484 logger.info("RI-gcd(div): r/g = " + rp + ", q/g = " + qp); 485 //logger.info("gcd(div): rp*g = " + rp.multiply(q) + ", qp*g = " + qp.multiply(q)); 486 throw new RuntimeException("recGcd, div: not divisible"); 487 } 488 } 489 qp = q; 490 if (false) { // not checkable 491 System.out.println("RI-recGcd, qs = " + qs + ", P = " + P); 492 System.out.println("RI-recGcd, rs = " + rs + ", S = " + S); 493 qs = FDUtil.<C> recursiveRightSparsePseudoRemainder(P, qp); // P 494 rs = FDUtil.<C> recursiveRightSparsePseudoRemainder(S, qp); // S 495 if (!qs.isZERO() || !rs.isZERO()) { 496 System.out.println("RI-recGcd, P = " + P); 497 System.out.println("RI-recGcd, S = " + S); 498 System.out.println("RI-recGcd, qp = " + qp); 499 System.out.println("RI-recGcd, qs = " + qs); 500 System.out.println("RI-recGcd, rs = " + rs); 501 throw new RuntimeException("RI-recGcd, pre pp: not divisible"); 502 } 503 logger.info("RI-recGcd(P,S) pre pp okay: qp = " + qp); 504 } 505 q = rightRecursivePrimitivePart(q); 506 if (!qp.equals(q)) { 507 logger.info("RI-gcd(pp) = " + q + ", qp = " + qp); // + ", ring = " + P.ring.toScript()); 508 } 509 q = (GenSolvablePolynomial<GenPolynomial<C>>) q.multiplyLeft(c).abs(); 510 if (false) { // not checkable 511 qs = FDUtil.<C> recursiveRightSparsePseudoRemainder(P, q); 512 rs = FDUtil.<C> recursiveRightSparsePseudoRemainder(S, q); 513 if (!qs.isZERO() || !rs.isZERO()) { 514 System.out.println("RI-recGcd-pp, P = " + P); 515 System.out.println("RI-recGcd-pp, S = " + S); 516 System.out.println("RI-recGcd-pp, q = " + q); 517 System.out.println("RI-recGcd-pp, qs = " + qs); 518 System.out.println("RI-recGcd-pp, rs = " + rs); 519 throw new RuntimeException("RI-recGcd: not divisible"); 520 } 521 logger.info("RI-recGcd(P,S) post pp okay: q = " + q); 522 } 523 return q; 524 } 525 526}