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