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