001/* 002 * $Id: GreatestCommonDivisorModular.java 4965 2014-10-17 20:07:51Z kredel $ 003 */ 004 005package edu.jas.ufd; 006 007 008import org.apache.log4j.Logger; 009 010import edu.jas.arith.BigInteger; 011import edu.jas.arith.Combinatoric; 012import edu.jas.arith.ModIntegerRing; 013import edu.jas.arith.ModLongRing; 014import edu.jas.arith.Modular; 015import edu.jas.arith.ModularRingFactory; 016import edu.jas.arith.PrimeList; 017import edu.jas.poly.ExpVector; 018import edu.jas.poly.GenPolynomial; 019import edu.jas.poly.GenPolynomialRing; 020import edu.jas.poly.PolyUtil; 021import edu.jas.structure.GcdRingElem; 022import edu.jas.structure.Power; 023 024 025/** 026 * Greatest common divisor algorithms with modular computation and chinese 027 * remainder algorithm. 028 * @author Heinz Kredel 029 */ 030 031public class GreatestCommonDivisorModular<MOD extends GcdRingElem<MOD> & Modular> extends 032 GreatestCommonDivisorAbstract<BigInteger> { 033 034 035 private static final Logger logger = Logger.getLogger(GreatestCommonDivisorModular.class); 036 037 038 private final boolean debug = logger.isDebugEnabled(); //logger.isInfoEnabled(); 039 040 041 /* 042 * Modular gcd algorithm to use. 043 */ 044 protected final GreatestCommonDivisorAbstract<MOD> mufd; 045 046 047 /* 048 * Integer gcd algorithm for fall back. 049 */ 050 protected final GreatestCommonDivisorAbstract<BigInteger> iufd = new GreatestCommonDivisorSubres<BigInteger>(); 051 052 053 /** 054 * Constructor to set recursive algorithm. Use modular evaluation GCD 055 * algorithm. 056 */ 057 public GreatestCommonDivisorModular() { 058 this(false); 059 } 060 061 062 /** 063 * Constructor to set recursive algorithm. 064 * @param simple , true if the simple PRS should be used. 065 */ 066 public GreatestCommonDivisorModular(boolean simple) { 067 if (simple) { 068 mufd = new GreatestCommonDivisorSimple<MOD>(); 069 } else { 070 mufd = new GreatestCommonDivisorModEval<MOD>(); 071 } 072 } 073 074 075 /** 076 * Univariate GenPolynomial greatest comon divisor. Delegate to subresultant 077 * baseGcd, should not be needed. 078 * @param P univariate GenPolynomial. 079 * @param S univariate GenPolynomial. 080 * @return gcd(P,S). 081 */ 082 @Override 083 public GenPolynomial<BigInteger> baseGcd(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) { 084 return iufd.baseGcd(P, S); 085 } 086 087 088 /** 089 * Univariate GenPolynomial recursive greatest comon divisor. Delegate to 090 * subresultant recursiveGcd, should not be needed. 091 * @param P univariate recursive GenPolynomial. 092 * @param S univariate recursive GenPolynomial. 093 * @return gcd(P,S). 094 */ 095 @Override 096 public GenPolynomial<GenPolynomial<BigInteger>> recursiveUnivariateGcd( 097 GenPolynomial<GenPolynomial<BigInteger>> P, GenPolynomial<GenPolynomial<BigInteger>> S) { 098 return iufd.recursiveUnivariateGcd(P, S); 099 } 100 101 102 /** 103 * GenPolynomial greatest comon divisor, modular algorithm. 104 * @param P GenPolynomial. 105 * @param S GenPolynomial. 106 * @return gcd(P,S). 107 */ 108 @Override 109 @SuppressWarnings("unchecked") 110 public GenPolynomial<BigInteger> gcd(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) { 111 if (S == null || S.isZERO()) { 112 return P; 113 } 114 if (P == null || P.isZERO()) { 115 return S; 116 } 117 GenPolynomialRing<BigInteger> fac = P.ring; 118 // special case for univariate polynomials 119 if (fac.nvar <= 1) { 120 GenPolynomial<BigInteger> T = baseGcd(P, S); 121 return T; 122 } 123 long e = P.degree(0); 124 long f = S.degree(0); 125 GenPolynomial<BigInteger> q; 126 GenPolynomial<BigInteger> r; 127 if (f > e) { 128 r = P; 129 q = S; 130 long g = f; 131 f = e; 132 e = g; 133 } else { 134 q = P; 135 r = S; 136 } 137 if (debug) { 138 logger.debug("degrees: e = " + e + ", f = " + f); 139 } 140 r = r.abs(); 141 q = q.abs(); 142 // compute contents and primitive parts 143 BigInteger a = baseContent(r); 144 BigInteger b = baseContent(q); 145 // gcd of coefficient contents 146 BigInteger c = gcd(a, b); // indirection 147 r = divide(r, a); // indirection 148 q = divide(q, b); // indirection 149 if (r.isONE()) { 150 return r.multiply(c); 151 } 152 if (q.isONE()) { 153 return q.multiply(c); 154 } 155 // compute normalization factor 156 BigInteger ac = r.leadingBaseCoefficient(); 157 BigInteger bc = q.leadingBaseCoefficient(); 158 BigInteger cc = gcd(ac, bc); // indirection 159 // compute norms 160 BigInteger an = r.maxNorm(); 161 BigInteger bn = q.maxNorm(); 162 BigInteger n = (an.compareTo(bn) < 0 ? bn : an); 163 n = n.multiply(cc).multiply(n.fromInteger(2)); 164 // compute degree vectors 165 ExpVector rdegv = r.degreeVector(); 166 ExpVector qdegv = q.degreeVector(); 167 //compute factor coefficient bounds 168 BigInteger af = an.multiply(PolyUtil.factorBound(rdegv)); 169 BigInteger bf = bn.multiply(PolyUtil.factorBound(qdegv)); 170 BigInteger cf = (af.compareTo(bf) < 0 ? bf : af); 171 cf = cf.multiply(cc.multiply(cc.fromInteger(8))); 172 //initialize prime list and degree vector 173 PrimeList primes = new PrimeList(); 174 int pn = 10; //primes.size(); 175 ExpVector wdegv = rdegv.subst(0, rdegv.getVal(0) + 1); 176 // +1 seems to be a hack for the unlucky prime test 177 ModularRingFactory<MOD> cofac; 178 ModularRingFactory<MOD> cofacM = null; 179 GenPolynomial<MOD> qm; 180 GenPolynomial<MOD> rm; 181 GenPolynomialRing<MOD> mfac; 182 GenPolynomialRing<MOD> rfac = null; 183 int i = 0; 184 BigInteger M = null; 185 BigInteger cfe = null; 186 GenPolynomial<MOD> cp = null; 187 GenPolynomial<MOD> cm = null; 188 GenPolynomial<BigInteger> cpi = null; 189 if (debug) { 190 logger.debug("c = " + c); 191 logger.debug("cc = " + cc); 192 logger.debug("n = " + n); 193 logger.debug("cf = " + cf); 194 logger.info("wdegv = " + wdegv); 195 } 196 for (java.math.BigInteger p : primes) { 197 //System.out.println("next run ++++++++++++++++++++++++++++++++++"); 198 if (p.longValue() == 2L) { // skip 2 199 continue; 200 } 201 if (++i >= pn) { 202 logger.warn("prime list exhausted, pn = " + pn); 203 return iufd.gcd(P, S); 204 //throw new ArithmeticException("prime list exhausted"); 205 } 206 // initialize coefficient factory and map normalization factor 207 if (ModLongRing.MAX_LONG.compareTo(p) > 0) { 208 cofac = (ModularRingFactory) new ModLongRing(p, true); 209 } else { 210 cofac = (ModularRingFactory) new ModIntegerRing(p, true); 211 } 212 MOD nf = cofac.fromInteger(cc.getVal()); 213 if (nf.isZERO()) { 214 continue; 215 } 216 // initialize polynomial factory and map polynomials 217 mfac = new GenPolynomialRing<MOD>(cofac, fac.nvar, fac.tord, fac.getVars()); 218 qm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, q); 219 if (qm.isZERO() || !qm.degreeVector().equals(qdegv)) { 220 continue; 221 } 222 rm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, r); 223 if (rm.isZERO() || !rm.degreeVector().equals(rdegv)) { 224 continue; 225 } 226 if (debug) { 227 logger.info("cofac = " + cofac.getIntegerModul()); 228 } 229 // compute modular gcd 230 cm = mufd.gcd(rm, qm); 231 // test for constant g.c.d 232 if (cm.isConstant()) { 233 logger.debug("cm, constant = " + cm); 234 return fac.getONE().multiply(c); 235 //return cm.abs().multiply( c ); 236 } 237 // test for unlucky prime 238 ExpVector mdegv = cm.degreeVector(); 239 if (wdegv.equals(mdegv)) { // TL = 0 240 // prime ok, next round 241 if (M != null) { 242 if (M.compareTo(cfe) > 0) { 243 System.out.println("M > cfe: " + M + " > " + cfe); 244 // continue; // why should this be required? 245 } 246 } 247 } else { // TL = 3 248 boolean ok = false; 249 if (wdegv.multipleOf(mdegv)) { // TL = 2 // EVMT(wdegv,mdegv) 250 M = null; // init chinese remainder 251 ok = true; // prime ok 252 } 253 if (mdegv.multipleOf(wdegv)) { // TL = 1 // EVMT(mdegv,wdegv) 254 continue; // skip this prime 255 } 256 if (!ok) { 257 M = null; // discard chinese remainder and previous work 258 continue; // prime not ok 259 } 260 } 261 //--wdegv = mdegv; 262 // prepare chinese remainder algorithm 263 cm = cm.multiply(nf); 264 if (M == null) { 265 // initialize chinese remainder algorithm 266 M = new BigInteger(p); 267 cofacM = cofac; 268 rfac = mfac; 269 cp = cm; 270 wdegv = wdegv.gcd(mdegv); //EVGCD(wdegv,mdegv); 271 cfe = cf; 272 for (int k = 0; k < wdegv.length(); k++) { 273 cfe = cfe.multiply(new BigInteger(wdegv.getVal(k) + 1)); 274 } 275 } else { 276 // apply chinese remainder algorithm 277 BigInteger Mp = M; 278 MOD mi = cofac.fromInteger(Mp.getVal()); 279 mi = mi.inverse(); // mod p 280 M = M.multiply(new BigInteger(p)); 281 if (ModLongRing.MAX_LONG.compareTo(M.getVal()) > 0) { 282 cofacM = (ModularRingFactory) new ModLongRing(M.getVal()); 283 } else { 284 cofacM = (ModularRingFactory) new ModIntegerRing(M.getVal()); 285 } 286 rfac = new GenPolynomialRing<MOD>(cofacM, fac); 287 if (!cofac.getClass().equals(cofacM.getClass())) { 288 logger.info("adjusting coefficents: cofacM = " + cofacM.getClass() + ", cofacP = " 289 + cofac.getClass()); 290 cofac = (ModularRingFactory) new ModIntegerRing(p); 291 mfac = new GenPolynomialRing<MOD>(cofac, fac); 292 GenPolynomial<BigInteger> mm = PolyUtil.<MOD> integerFromModularCoefficients(fac, cm); 293 cm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, mm); 294 mi = cofac.fromInteger(Mp.getVal()); 295 mi = mi.inverse(); // mod p 296 } 297 if (!cp.ring.coFac.getClass().equals(cofacM.getClass())) { 298 logger.info("adjusting coefficents: cofacM = " + cofacM.getClass() + ", cofacM' = " 299 + cp.ring.coFac.getClass()); 300 ModularRingFactory cop = (ModularRingFactory) cp.ring.coFac; 301 cofac = (ModularRingFactory) new ModIntegerRing(cop.getIntegerModul().getVal()); 302 mfac = new GenPolynomialRing<MOD>(cofac, fac); 303 GenPolynomial<BigInteger> mm = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp); 304 cp = PolyUtil.<MOD> fromIntegerCoefficients(mfac, mm); 305 } 306 cp = PolyUtil.<MOD> chineseRemainder(rfac, cp, mi, cm); 307 } 308 // test for completion 309 if (n.compareTo(M) <= 0) { 310 break; 311 } 312 // must use integer.sumNorm 313 cpi = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp); 314 BigInteger cmn = cpi.sumNorm(); 315 cmn = cmn.multiply(cmn.fromInteger(4)); 316 //if ( cmn.compareTo( M ) <= 0 ) { 317 // does not work: only if cofactors are also considered? 318 // break; 319 //} 320 if (i % 2 != 0 && !cp.isZERO()) { 321 // check if done on every second prime 322 GenPolynomial<BigInteger> x; 323 x = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp); 324 x = basePrimitivePart(x); 325 if (!PolyUtil.<BigInteger> baseSparsePseudoRemainder(q, x).isZERO()) { 326 continue; 327 } 328 if (!PolyUtil.<BigInteger> baseSparsePseudoRemainder(r, x).isZERO()) { 329 continue; 330 } 331 logger.info("done on exact division, #primes = " + i); 332 break; 333 } 334 } 335 if (debug) { 336 logger.info("done on M = " + M + ", #primes = " + i); 337 } 338 // remove normalization 339 q = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp); 340 q = basePrimitivePart(q); 341 return q.abs().multiply(c); 342 } 343 344 345 /** 346 * Univariate GenPolynomial resultant. 347 * @param P univariate GenPolynomial. 348 * @param S univariate GenPolynomial. 349 * @return res(P,S). 350 */ 351 @Override 352 public GenPolynomial<BigInteger> baseResultant(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) { 353 // not a special case here 354 return resultant(P, S); 355 } 356 357 358 /** 359 * Univariate GenPolynomial recursive resultant. 360 * @param P univariate recursive GenPolynomial. 361 * @param S univariate recursive GenPolynomial. 362 * @return res(P,S). 363 */ 364 @Override 365 public GenPolynomial<GenPolynomial<BigInteger>> recursiveUnivariateResultant( 366 GenPolynomial<GenPolynomial<BigInteger>> P, GenPolynomial<GenPolynomial<BigInteger>> S) { 367 // only in this class 368 return recursiveResultant(P, S); 369 } 370 371 372 /** 373 * GenPolynomial resultant, modular algorithm. 374 * @param P GenPolynomial. 375 * @param S GenPolynomial. 376 * @return res(P,S). 377 */ 378 @Override 379 @SuppressWarnings("unchecked") 380 public GenPolynomial<BigInteger> resultant(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) { 381 if (S == null || S.isZERO()) { 382 return S; 383 } 384 if (P == null || P.isZERO()) { 385 return P; 386 } 387 GenPolynomialRing<BigInteger> fac = P.ring; 388 // no special case for univariate polynomials in this class ! 389 //if (fac.nvar <= 1) { 390 // GenPolynomial<BigInteger> T = iufd.baseResultant(P, S); 391 // return T; 392 //} 393 long e = P.degree(0); 394 long f = S.degree(0); 395 GenPolynomial<BigInteger> q; 396 GenPolynomial<BigInteger> r; 397 if (f > e) { 398 r = P; 399 q = S; 400 long g = f; 401 f = e; 402 e = g; 403 } else { 404 q = P; 405 r = S; 406 } 407 // compute norms 408 BigInteger an = r.maxNorm(); 409 BigInteger bn = q.maxNorm(); 410 an = Power.<BigInteger> power(fac.coFac, an, f); 411 bn = Power.<BigInteger> power(fac.coFac, bn, e); 412 BigInteger cn = Combinatoric.factorial(e + f); 413 BigInteger n = cn.multiply(an).multiply(bn); 414 415 // compute degree vectors 416 ExpVector rdegv = r.leadingExpVector(); //degreeVector(); 417 ExpVector qdegv = q.leadingExpVector(); //degreeVector(); 418 419 //initialize prime list and degree vector 420 PrimeList primes = new PrimeList(); 421 int pn = 30; //primes.size(); 422 ModularRingFactory<MOD> cofac; 423 ModularRingFactory<MOD> cofacM = null; 424 GenPolynomial<MOD> qm; 425 GenPolynomial<MOD> rm; 426 GenPolynomialRing<MOD> mfac; 427 GenPolynomialRing<MOD> rfac = null; 428 int i = 0; 429 BigInteger M = null; 430 GenPolynomial<MOD> cp = null; 431 GenPolynomial<MOD> cm = null; 432 //GenPolynomial<BigInteger> cpi = null; 433 if (debug) { 434 logger.debug("an = " + an); 435 logger.debug("bn = " + bn); 436 logger.debug("e+f = " + (e + f)); 437 logger.debug("cn = " + cn); 438 logger.info("n = " + n); 439 //logger.info("q = " + q); 440 //logger.info("r = " + r); 441 //logger.info("rdegv = " + rdegv.toString(fac.getVars())); 442 //logger.info("qdegv = " + qdegv.toString(fac.getVars())); 443 } 444 for (java.math.BigInteger p : primes) { 445 //System.out.println("next run ++++++++++++++++++++++++++++++++++"); 446 if (p.longValue() == 2L) { // skip 2 447 continue; 448 } 449 if (++i >= pn) { 450 logger.warn("prime list exhausted, pn = " + pn); 451 return iufd.resultant(P, S); 452 //throw new ArithmeticException("prime list exhausted"); 453 } 454 // initialize coefficient factory and map normalization factor 455 if (ModLongRing.MAX_LONG.compareTo(p) > 0) { 456 cofac = (ModularRingFactory) new ModLongRing(p, true); 457 } else { 458 cofac = (ModularRingFactory) new ModIntegerRing(p, true); 459 } 460 // initialize polynomial factory and map polynomials 461 mfac = new GenPolynomialRing<MOD>(cofac, fac); 462 qm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, q); 463 if (qm.isZERO() || !qm.leadingExpVector().equals(qdegv)) { //degreeVector() 464 //logger.info("qm = " + qm); 465 if (debug) { 466 logger.info("unlucky prime = " + cofac.getIntegerModul() + ", degv = " 467 + qm.leadingExpVector()); 468 } 469 continue; 470 } 471 rm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, r); 472 if (rm.isZERO() || !rm.leadingExpVector().equals(rdegv)) { //degreeVector() 473 //logger.info("rm = " + rm); 474 if (debug) { 475 logger.info("unlucky prime = " + cofac.getIntegerModul() + ", degv = " 476 + rm.leadingExpVector()); 477 } 478 continue; 479 } 480 logger.info("lucky prime = " + cofac.getIntegerModul()); 481 482 // compute modular resultant 483 cm = mufd.resultant(qm, rm); 484 if (debug) { 485 logger.info("res_p = " + cm); 486 } 487 488 // prepare chinese remainder algorithm 489 if (M == null) { 490 // initialize chinese remainder algorithm 491 M = new BigInteger(p); 492 cofacM = cofac; 493 //rfac = mfac; 494 cp = cm; 495 } else { 496 // apply chinese remainder algorithm 497 BigInteger Mp = M; 498 MOD mi = cofac.fromInteger(Mp.getVal()); 499 mi = mi.inverse(); // mod p 500 M = M.multiply(new BigInteger(p)); 501 if (ModLongRing.MAX_LONG.compareTo(M.getVal()) > 0) { 502 cofacM = (ModularRingFactory) new ModLongRing(M.getVal()); 503 } else { 504 cofacM = (ModularRingFactory) new ModIntegerRing(M.getVal()); 505 } 506 rfac = new GenPolynomialRing<MOD>(cofacM, fac); 507 if (!cofac.getClass().equals(cofacM.getClass())) { 508 logger.info("adjusting coefficents: cofacM = " + cofacM.getClass() + ", cofacP = " 509 + cofac.getClass()); 510 cofac = (ModularRingFactory) new ModIntegerRing(p); 511 mfac = new GenPolynomialRing<MOD>(cofac, fac); 512 GenPolynomial<BigInteger> mm = PolyUtil.<MOD> integerFromModularCoefficients(fac, cm); 513 cm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, mm); 514 mi = cofac.fromInteger(Mp.getVal()); 515 mi = mi.inverse(); // mod p 516 } 517 if (!cp.ring.coFac.getClass().equals(cofacM.getClass())) { 518 logger.info("adjusting coefficents: cofacM = " + cofacM.getClass() + ", cofacM' = " 519 + cp.ring.coFac.getClass()); 520 ModularRingFactory cop = (ModularRingFactory) cp.ring.coFac; 521 cofac = (ModularRingFactory) new ModIntegerRing(cop.getIntegerModul().getVal()); 522 mfac = new GenPolynomialRing<MOD>(cofac, fac); 523 GenPolynomial<BigInteger> mm = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp); 524 cp = PolyUtil.<MOD> fromIntegerCoefficients(mfac, mm); 525 } 526 cp = PolyUtil.<MOD> chineseRemainder(rfac, cp, mi, cm); 527 } 528 // test for completion 529 if (n.compareTo(M) <= 0) { 530 break; 531 } 532 } 533 if (debug) { 534 logger.info("done on M = " + M + ", #primes = " + i); 535 } 536 // convert to integer polynomial 537 q = PolyUtil.<MOD> integerFromModularCoefficients(fac, cp); 538 return q; 539 } 540 541}