001/* 002 * $Id: FactorInteger.java 4115 2012-08-19 13:18:59Z kredel $ 003 */ 004 005package edu.jas.ufd; 006 007 008import java.util.ArrayList; 009import java.util.BitSet; 010import java.util.Iterator; 011import java.util.List; 012import java.util.SortedMap; 013 014import org.apache.log4j.Logger; 015 016import edu.jas.arith.BigInteger; 017import edu.jas.arith.ModIntegerRing; 018import edu.jas.arith.ModLongRing; 019import edu.jas.arith.Modular; 020import edu.jas.arith.ModularRingFactory; 021import edu.jas.arith.PrimeList; 022import edu.jas.poly.ExpVector; 023import edu.jas.poly.GenPolynomial; 024import edu.jas.poly.GenPolynomialRing; 025import edu.jas.poly.OptimizedPolynomialList; 026import edu.jas.poly.PolyUtil; 027import edu.jas.poly.TermOrderOptimization; 028import edu.jas.structure.GcdRingElem; 029import edu.jas.structure.Power; 030import edu.jas.structure.RingElem; 031import edu.jas.structure.RingFactory; 032import edu.jas.util.KsubSet; 033 034 035/** 036 * Integer coefficients factorization algorithms. This class implements 037 * factorization methods for polynomials over integers. 038 * @author Heinz Kredel 039 */ 040 041/** 042 * @author kredel 043 * 044 * @param <MOD> 045 */ 046public class FactorInteger<MOD extends GcdRingElem<MOD> & Modular> extends FactorAbstract<BigInteger> { 047 048 049 private static final Logger logger = Logger.getLogger(FactorInteger.class); 050 051 052 private final boolean debug = logger.isDebugEnabled(); 053 054 055 /** 056 * Factorization engine for modular base coefficients. 057 */ 058 protected final FactorAbstract<MOD> mfactor; 059 060 061 /** 062 * Gcd engine for modular base coefficients. 063 */ 064 protected final GreatestCommonDivisorAbstract<MOD> mengine; 065 066 067 /** 068 * No argument constructor. 069 */ 070 public FactorInteger() { 071 this(BigInteger.ONE); 072 } 073 074 075 /** 076 * Constructor. 077 * @param cfac coefficient ring factory. 078 */ 079 public FactorInteger(RingFactory<BigInteger> cfac) { 080 super(cfac); 081 ModularRingFactory<MOD> mcofac = (ModularRingFactory<MOD>) (Object) new ModLongRing(13, true); // hack 082 mfactor = FactorFactory.getImplementation(mcofac); //new FactorModular(mcofac); 083 mengine = GCDFactory.getImplementation(mcofac); 084 //mengine = GCDFactory.getProxy(mcofac); 085 } 086 087 088 /** 089 * GenPolynomial base factorization of a squarefree polynomial. 090 * @param P squarefree and primitive! GenPolynomial. 091 * @return [p_1,...,p_k] with P = prod_{i=1, ..., k} p_i. 092 */ 093 @SuppressWarnings("unchecked") 094 @Override 095 public List<GenPolynomial<BigInteger>> baseFactorsSquarefree(GenPolynomial<BigInteger> P) { 096 if (P == null) { 097 throw new IllegalArgumentException(this.getClass().getName() + " P == null"); 098 } 099 List<GenPolynomial<BigInteger>> factors = new ArrayList<GenPolynomial<BigInteger>>(); 100 if (P.isZERO()) { 101 return factors; 102 } 103 if (P.isONE()) { 104 factors.add(P); 105 return factors; 106 } 107 GenPolynomialRing<BigInteger> pfac = P.ring; 108 if (pfac.nvar > 1) { 109 throw new IllegalArgumentException(this.getClass().getName() + " only for univariate polynomials"); 110 } 111 if (!engine.baseContent(P).isONE()) { 112 throw new IllegalArgumentException(this.getClass().getName() + " P not primitive"); 113 } 114 if (P.degree(0) <= 1L) { // linear is irreducible 115 factors.add(P); 116 return factors; 117 } 118 // compute norm 119 BigInteger an = P.maxNorm(); 120 BigInteger ac = P.leadingBaseCoefficient(); 121 //compute factor coefficient bounds 122 ExpVector degv = P.degreeVector(); 123 int degi = (int) P.degree(0); 124 BigInteger M = an.multiply(PolyUtil.factorBound(degv)); 125 M = M.multiply(ac.abs().multiply(ac.fromInteger(8))); 126 //System.out.println("M = " + M); 127 //M = M.multiply(M); // test 128 129 //initialize prime list and degree vector 130 PrimeList primes = new PrimeList(PrimeList.Range.small); 131 int pn = 30; //primes.size(); 132 ModularRingFactory<MOD> cofac = null; 133 GenPolynomial<MOD> am = null; 134 GenPolynomialRing<MOD> mfac = null; 135 final int TT = 5; // 7 136 List<GenPolynomial<MOD>>[] modfac = new List[TT]; 137 List<GenPolynomial<BigInteger>>[] intfac = new List[TT]; 138 BigInteger[] plist = new BigInteger[TT]; 139 List<GenPolynomial<MOD>> mlist = null; 140 List<GenPolynomial<BigInteger>> ilist = null; 141 int i = 0; 142 if (debug) { 143 logger.debug("an = " + an); 144 logger.debug("ac = " + ac); 145 logger.debug("M = " + M); 146 logger.info("degv = " + degv); 147 } 148 Iterator<java.math.BigInteger> pit = primes.iterator(); 149 pit.next(); // skip p = 2 150 pit.next(); // skip p = 3 151 MOD nf = null; 152 for (int k = 0; k < TT; k++) { 153 if (k == TT - 1) { // -2 154 primes = new PrimeList(PrimeList.Range.medium); 155 pit = primes.iterator(); 156 } 157 if (k == TT + 1) { // -1 158 primes = new PrimeList(PrimeList.Range.large); 159 pit = primes.iterator(); 160 } 161 while (pit.hasNext()) { 162 java.math.BigInteger p = pit.next(); 163 //System.out.println("next run ++++++++++++++++++++++++++++++++++"); 164 if (++i >= pn) { 165 logger.error("prime list exhausted, pn = " + pn); 166 throw new ArithmeticException("prime list exhausted"); 167 } 168 if (ModLongRing.MAX_LONG.compareTo(p) > 0) { 169 cofac = (ModularRingFactory) new ModLongRing(p, true); 170 } else { 171 cofac = (ModularRingFactory) new ModIntegerRing(p, true); 172 } 173 logger.info("prime = " + cofac); 174 nf = cofac.fromInteger(ac.getVal()); 175 if (nf.isZERO()) { 176 logger.info("unlucky prime (nf) = " + p); 177 //System.out.println("unlucky prime (nf) = " + p); 178 continue; 179 } 180 // initialize polynomial factory and map polynomial 181 mfac = new GenPolynomialRing<MOD>(cofac, pfac); 182 am = PolyUtil.<MOD> fromIntegerCoefficients(mfac, P); 183 if (!am.degreeVector().equals(degv)) { // allways true 184 logger.info("unlucky prime (deg) = " + p); 185 //System.out.println("unlucky prime (deg) = " + p); 186 continue; 187 } 188 GenPolynomial<MOD> ap = PolyUtil.<MOD> baseDeriviative(am); 189 if (ap.isZERO()) { 190 logger.info("unlucky prime (a')= " + p); 191 //System.out.println("unlucky prime (a')= " + p); 192 continue; 193 } 194 GenPolynomial<MOD> g = mengine.baseGcd(am, ap); 195 if (g.isONE()) { 196 logger.info("**lucky prime = " + p); 197 //System.out.println("**lucky prime = " + p); 198 break; 199 } 200 } 201 // now am is squarefree mod p, make monic and factor mod p 202 if (!nf.isONE()) { 203 //System.out.println("nf = " + nf); 204 am = am.divide(nf); // make monic 205 } 206 mlist = mfactor.baseFactorsSquarefree(am); 207 if (logger.isInfoEnabled()) { 208 logger.info("modlist = " + mlist); 209 } 210 if (mlist.size() <= 1) { 211 factors.add(P); 212 return factors; 213 } 214 if (!nf.isONE()) { 215 GenPolynomial<MOD> mp = mfac.getONE(); //mlist.get(0); 216 //System.out.println("mp = " + mp); 217 mp = mp.multiply(nf); 218 //System.out.println("mp = " + mp); 219 mlist.add(0, mp); // set(0,mp); 220 } 221 modfac[k] = mlist; 222 plist[k] = cofac.getIntegerModul(); // p 223 } 224 225 // search shortest factor list 226 int min = Integer.MAX_VALUE; 227 BitSet AD = null; 228 for (int k = 0; k < TT; k++) { 229 List<ExpVector> ev = PolyUtil.<MOD> leadingExpVector(modfac[k]); 230 BitSet D = factorDegrees(ev, degi); 231 if (AD == null) { 232 AD = D; 233 } else { 234 AD.and(D); 235 } 236 int s = modfac[k].size(); 237 logger.info("mod(" + plist[k] + ") #s = " + s + ", D = " + D /*+ ", lt = " + ev*/); 238 //System.out.println("mod s = " + s); 239 if (s < min) { 240 min = s; 241 mlist = modfac[k]; 242 } 243 } 244 logger.info("min = " + min + ", AD = " + AD); 245 if (mlist.size() <= 1) { 246 logger.info("mlist.size() = 1"); 247 factors.add(P); 248 return factors; 249 } 250 if (AD.cardinality() <= 2) { // only one possible factor 251 logger.info("degree set cardinality = " + AD.cardinality()); 252 factors.add(P); 253 return factors; 254 } 255 256 boolean allLists = false; //true; //false; 257 if (allLists) { 258 // try each factor list 259 for (int k = 0; k < TT; k++) { 260 mlist = modfac[k]; 261 if (debug) { 262 logger.info("lifting from " + mlist); 263 } 264 if (P.leadingBaseCoefficient().isONE()) { // monic case 265 factors = searchFactorsMonic(P, M, mlist, AD); // does now work in all cases 266 if (factors.size() == 1) { 267 factors = searchFactorsNonMonic(P, M, mlist, AD); 268 } 269 } else { 270 factors = searchFactorsNonMonic(P, M, mlist, AD); 271 } 272 intfac[k] = factors; 273 } 274 } else { 275 // try only shortest factor list 276 if (debug) { 277 logger.info("lifting shortest from " + mlist); 278 } 279 if (true && P.leadingBaseCoefficient().isONE()) { 280 long t = System.currentTimeMillis(); 281 try { 282 mlist = PolyUtil.<MOD> monic(mlist); 283 factors = searchFactorsMonic(P, M, mlist, AD); // does now work in all cases 284 t = System.currentTimeMillis() - t; 285 //System.out.println("monic time = " + t); 286 if (debug) { 287 t = System.currentTimeMillis(); 288 List<GenPolynomial<BigInteger>> fnm = searchFactorsNonMonic(P, M, mlist, AD); 289 t = System.currentTimeMillis() - t; 290 System.out.println("non monic time = " + t); 291 if (!factors.equals(fnm)) { 292 System.out.println("monic factors = " + factors); 293 System.out.println("non monic factors = " + fnm); 294 } 295 } 296 } catch (RuntimeException e) { 297 t = System.currentTimeMillis(); 298 factors = searchFactorsNonMonic(P, M, mlist, AD); 299 t = System.currentTimeMillis() - t; 300 //System.out.println("only non monic time = " + t); 301 } 302 } else { 303 long t = System.currentTimeMillis(); 304 factors = searchFactorsNonMonic(P, M, mlist, AD); 305 t = System.currentTimeMillis() - t; 306 //System.out.println("non monic time = " + t); 307 } 308 return normalizeFactorization(factors); 309 } 310 311 // search longest factor list 312 int max = 0; 313 for (int k = 0; k < TT; k++) { 314 int s = intfac[k].size(); 315 logger.info("int s = " + s); 316 //System.out.println("int s = " + s); 317 if (s > max) { 318 max = s; 319 ilist = intfac[k]; 320 } 321 } 322 factors = normalizeFactorization(ilist); 323 return factors; 324 } 325 326 327 /** 328 * BitSet for factor degree list. 329 * @param E exponent vector list. 330 * @return b_0,...,b_k} a BitSet of possible factor degrees. 331 */ 332 public BitSet factorDegrees(List<ExpVector> E, int deg) { 333 BitSet D = new BitSet(deg + 1); 334 D.set(0); // constant factor 335 for (ExpVector e : E) { 336 int i = (int) e.getVal(0); 337 BitSet s = new BitSet(deg + 1); 338 for (int k = 0; k < deg + 1 - i; k++) { // shift by i places 339 s.set(i + k, D.get(k)); 340 } 341 //System.out.println("s = " + s); 342 D.or(s); 343 //System.out.println("D = " + D); 344 } 345 return D; 346 } 347 348 349 /** 350 * Sum of all degrees. 351 * @param L univariate polynomial list. 352 * @return sum deg(p) for p in L. 353 */ 354 public static <C extends RingElem<C>> long degreeSum(List<GenPolynomial<C>> L) { 355 long s = 0L; 356 for (GenPolynomial<C> p : L) { 357 ExpVector e = p.leadingExpVector(); 358 long d = e.getVal(0); 359 s += d; 360 } 361 return s; 362 } 363 364 365 /** 366 * Factor search with modular Hensel lifting algorithm. Let p = 367 * f_i.ring.coFac.modul() i = 0, ..., n-1 and assume C == prod_{0,...,n-1} 368 * f_i mod p with ggt(f_i,f_j) == 1 mod p for i != j 369 * @param C GenPolynomial. 370 * @param M bound on the coefficients of g_i as factors of C. 371 * @param F = [f_0,...,f_{n-1}] List<GenPolynomial>. 372 * @param D bit set of possible factor degrees. 373 * @return [g_0,...,g_{n-1}] = lift(C,F), with C = prod_{0,...,n-1} g_i mod 374 * p**e. <b>Note:</b> does not work in all cases. 375 */ 376 List<GenPolynomial<BigInteger>> searchFactorsMonic(GenPolynomial<BigInteger> C, BigInteger M, 377 List<GenPolynomial<MOD>> F, BitSet D) { 378 //System.out.println("*** monic factor combination ***"); 379 if (C == null || C.isZERO() || F == null || F.size() == 0) { 380 throw new IllegalArgumentException("C must be nonzero and F must be nonempty"); 381 } 382 GenPolynomialRing<BigInteger> pfac = C.ring; 383 if (pfac.nvar != 1) { // todo assert 384 throw new IllegalArgumentException("polynomial ring not univariate"); 385 } 386 List<GenPolynomial<BigInteger>> factors = new ArrayList<GenPolynomial<BigInteger>>(F.size()); 387 List<GenPolynomial<MOD>> mlist = F; 388 List<GenPolynomial<MOD>> lift; 389 390 //MOD nf = null; 391 GenPolynomial<MOD> ct = mlist.get(0); 392 if (ct.isConstant()) { 393 //nf = ct.leadingBaseCoefficient(); 394 mlist.remove(ct); 395 //System.out.println("=== nf = " + nf); 396 if (mlist.size() <= 1) { 397 factors.add(C); 398 return factors; 399 } 400 } else { 401 //nf = ct.ring.coFac.getONE(); 402 } 403 //System.out.println("modlist = " + mlist); // includes not ldcf 404 ModularRingFactory<MOD> mcfac = (ModularRingFactory<MOD>) ct.ring.coFac; 405 BigInteger m = mcfac.getIntegerModul(); 406 long k = 1; 407 BigInteger pi = m; 408 while (pi.compareTo(M) < 0) { 409 k++; 410 pi = pi.multiply(m); 411 } 412 logger.info("p^k = " + m + "^" + k); 413 GenPolynomial<BigInteger> PP = C, P = C; 414 // lift via Hensel 415 try { 416 lift = HenselUtil.<MOD> liftHenselMonic(PP, mlist, k); 417 //System.out.println("lift = " + lift); 418 } catch (NoLiftingException e) { 419 throw new RuntimeException(e); 420 } 421 if (logger.isInfoEnabled()) { 422 logger.info("lifted modlist = " + lift); 423 } 424 GenPolynomialRing<MOD> mpfac = lift.get(0).ring; 425 426 // combine trial factors 427 int dl = (lift.size() + 1) / 2; 428 //System.out.println("dl = " + dl); 429 GenPolynomial<BigInteger> u = PP; 430 long deg = (u.degree(0) + 1L) / 2L; 431 //System.out.println("deg = " + deg); 432 //BigInteger ldcf = u.leadingBaseCoefficient(); 433 //System.out.println("ldcf = " + ldcf); 434 for (int j = 1; j <= dl; j++) { 435 //System.out.println("j = " + j + ", dl = " + dl + ", lift = " + lift); 436 KsubSet<GenPolynomial<MOD>> ps = new KsubSet<GenPolynomial<MOD>>(lift, j); 437 for (List<GenPolynomial<MOD>> flist : ps) { 438 //System.out.println("degreeSum = " + degreeSum(flist)); 439 if (!D.get((int) FactorInteger.<MOD> degreeSum(flist))) { 440 logger.info("skipped by degree set " + D + ", deg = " + degreeSum(flist)); 441 continue; 442 } 443 GenPolynomial<MOD> mtrial = Power.<GenPolynomial<MOD>> multiply(mpfac, flist); 444 //GenPolynomial<MOD> mtrial = mpfac.getONE(); 445 //for (int kk = 0; kk < flist.size(); kk++) { 446 // GenPolynomial<MOD> fk = flist.get(kk); 447 // mtrial = mtrial.multiply(fk); 448 //} 449 //System.out.println("+flist = " + flist + ", mtrial = " + mtrial); 450 if (mtrial.degree(0) > deg) { // this test is sometimes wrong 451 logger.info("degree " + mtrial.degree(0) + " > deg " + deg); 452 //continue; 453 } 454 //System.out.println("+flist = " + flist); 455 GenPolynomial<BigInteger> trial = PolyUtil.integerFromModularCoefficients(pfac, mtrial); 456 //System.out.println("+trial = " + trial); 457 //trial = engine.basePrimitivePart( trial.multiply(ldcf) ); 458 trial = engine.basePrimitivePart(trial); 459 //System.out.println("pp(trial)= " + trial); 460 if (PolyUtil.<BigInteger> baseSparsePseudoRemainder(u, trial).isZERO()) { 461 logger.info("successful trial = " + trial); 462 //System.out.println("trial = " + trial); 463 //System.out.println("flist = " + flist); 464 //trial = engine.basePrimitivePart(trial); 465 //System.out.println("pp(trial)= " + trial); 466 factors.add(trial); 467 u = PolyUtil.<BigInteger> basePseudoDivide(u, trial); //u.divide( trial ); 468 //System.out.println("u = " + u); 469 //if (lift.removeAll(flist)) { 470 lift = removeOnce(lift, flist); 471 logger.info("new lift= " + lift); 472 dl = (lift.size() + 1) / 2; 473 //System.out.println("dl = " + dl); 474 j = 0; // since j++ 475 break; 476 //} logger.error("error removing flist from lift = " + lift); 477 } 478 } 479 } 480 if (!u.isONE() && !u.equals(P)) { 481 logger.info("rest u = " + u); 482 //System.out.println("rest u = " + u); 483 factors.add(u); 484 } 485 if (factors.size() == 0) { 486 logger.info("irred u = " + u); 487 //System.out.println("irred u = " + u); 488 factors.add(PP); 489 } 490 return normalizeFactorization(factors); 491 } 492 493 494 /** 495 * Factor search with modular Hensel lifting algorithm. Let p = 496 * f_i.ring.coFac.modul() i = 0, ..., n-1 and assume C == prod_{0,...,n-1} 497 * f_i mod p with ggt(f_i,f_j) == 1 mod p for i != j 498 * @param C GenPolynomial. 499 * @param M bound on the coefficients of g_i as factors of C. 500 * @param F = [f_0,...,f_{n-1}] List<GenPolynomial>. 501 * @param D bit set of possible factor degrees. 502 * @return [g_0,...,g_{n-1}] = lift(C,F), with C = prod_{0,...,n-1} g_i mod 503 * p**e. 504 */ 505 List<GenPolynomial<BigInteger>> searchFactorsNonMonic(GenPolynomial<BigInteger> C, BigInteger M, 506 List<GenPolynomial<MOD>> F, BitSet D) { 507 //System.out.println("*** non monic factor combination ***"); 508 if (C == null || C.isZERO() || F == null || F.size() == 0) { 509 throw new IllegalArgumentException("C must be nonzero and F must be nonempty"); 510 } 511 GenPolynomialRing<BigInteger> pfac = C.ring; 512 if (pfac.nvar != 1) { // todo assert 513 throw new IllegalArgumentException("polynomial ring not univariate"); 514 } 515 List<GenPolynomial<BigInteger>> factors = new ArrayList<GenPolynomial<BigInteger>>(F.size()); 516 List<GenPolynomial<MOD>> mlist = F; 517 518 MOD nf = null; 519 GenPolynomial<MOD> ct = mlist.get(0); 520 if (ct.isConstant()) { 521 nf = ct.leadingBaseCoefficient(); 522 mlist.remove(ct); 523 //System.out.println("=== nf = " + nf); 524 //System.out.println("=== ldcf = " + C.leadingBaseCoefficient()); 525 if (mlist.size() <= 1) { 526 factors.add(C); 527 return factors; 528 } 529 } else { 530 nf = ct.ring.coFac.getONE(); 531 } 532 //System.out.println("modlist = " + mlist); // includes not ldcf 533 GenPolynomialRing<MOD> mfac = ct.ring; 534 GenPolynomial<MOD> Pm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, C); 535 GenPolynomial<BigInteger> PP = C, P = C; 536 537 // combine trial factors 538 int dl = (mlist.size() + 1) / 2; 539 GenPolynomial<BigInteger> u = PP; 540 long deg = (u.degree(0) + 1L) / 2L; 541 GenPolynomial<MOD> um = Pm; 542 //BigInteger ldcf = u.leadingBaseCoefficient(); 543 //System.out.println("ldcf = " + ldcf); 544 HenselApprox<MOD> ilist = null; 545 for (int j = 1; j <= dl; j++) { 546 //System.out.println("j = " + j + ", dl = " + dl + ", ilist = " + ilist); 547 KsubSet<GenPolynomial<MOD>> ps = new KsubSet<GenPolynomial<MOD>>(mlist, j); 548 for (List<GenPolynomial<MOD>> flist : ps) { 549 //System.out.println("degreeSum = " + degreeSum(flist)); 550 if (!D.get((int) FactorInteger.<MOD> degreeSum(flist))) { 551 logger.info("skipped by degree set " + D + ", deg = " + degreeSum(flist)); 552 continue; 553 } 554 GenPolynomial<MOD> trial = mfac.getONE().multiply(nf); 555 for (int kk = 0; kk < flist.size(); kk++) { 556 GenPolynomial<MOD> fk = flist.get(kk); 557 trial = trial.multiply(fk); 558 } 559 if (trial.degree(0) > deg) { // this test is sometimes wrong 560 logger.info("degree > deg " + deg + ", degree = " + trial.degree(0)); 561 //continue; 562 } 563 GenPolynomial<MOD> cofactor = um.divide(trial); 564 //System.out.println("trial = " + trial); 565 //System.out.println("cofactor = " + cofactor); 566 567 // lift via Hensel 568 try { 569 // ilist = HenselUtil.liftHenselQuadraticFac(PP, M, trial, cofactor); 570 ilist = HenselUtil.<MOD> liftHenselQuadratic(PP, M, trial, cofactor); 571 //ilist = HenselUtil.<MOD> liftHensel(PP, M, trial, cofactor); 572 } catch (NoLiftingException e) { 573 // no liftable factors 574 if ( /*debug*/logger.isDebugEnabled()) { 575 logger.info("no liftable factors " + e); 576 //e.printStackTrace(); 577 } 578 continue; 579 } 580 GenPolynomial<BigInteger> itrial = ilist.A; 581 GenPolynomial<BigInteger> icofactor = ilist.B; 582 if (logger.isDebugEnabled()) { 583 logger.info(" modlist = " + trial + ", cofactor " + cofactor); 584 logger.info("lifted intlist = " + itrial + ", cofactor " + icofactor); 585 } 586 //System.out.println("lifted intlist = " + itrial + ", cofactor " + icofactor); 587 588 itrial = engine.basePrimitivePart(itrial); 589 //System.out.println("pp(trial)= " + itrial); 590 if (PolyUtil.<BigInteger> baseSparsePseudoRemainder(u, itrial).isZERO()) { 591 logger.info("successful trial = " + itrial); 592 //System.out.println("trial = " + itrial); 593 //System.out.println("cofactor = " + icofactor); 594 //System.out.println("flist = " + flist); 595 //itrial = engine.basePrimitivePart(itrial); 596 //System.out.println("pp(itrial)= " + itrial); 597 factors.add(itrial); 598 //u = PolyUtil.<BigInteger> basePseudoDivide(u, itrial); //u.divide( trial ); 599 u = icofactor; 600 PP = u; // fixed finally on 2009-05-03 601 um = cofactor; 602 //System.out.println("u = " + u); 603 //System.out.println("um = " + um); 604 //if (mlist.removeAll(flist)) { 605 mlist = removeOnce(mlist, flist); 606 logger.info("new mlist= " + mlist); 607 dl = (mlist.size() + 1) / 2; 608 j = 0; // since j++ 609 break; 610 //} logger.error("error removing flist from ilist = " + mlist); 611 } 612 } 613 } 614 if (!u.isONE() && !u.equals(P)) { 615 logger.info("rest u = " + u); 616 //System.out.println("rest u = " + u); 617 factors.add(u); 618 } 619 if (factors.size() == 0) { 620 logger.info("irred u = " + u); 621 //System.out.println("irred u = " + u); 622 factors.add(PP); 623 } 624 return normalizeFactorization(factors); 625 } 626 627 628 /** 629 * GenPolynomial factorization of a multivariate squarefree polynomial, 630 * using Hensel lifting if possible. 631 * @param P squarefree and primitive! (respectively monic) multivariate 632 * GenPolynomial over the integers. 633 * @return [p_1,...,p_k] with P = prod_{i=1,...,r} p_i. 634 */ 635 @Override 636 public List<GenPolynomial<BigInteger>> factorsSquarefree(GenPolynomial<BigInteger> P) { 637 GenPolynomialRing<BigInteger> pfac = P.ring; 638 if (pfac.nvar <= 1) { 639 return baseFactorsSquarefree(P); 640 } 641 List<GenPolynomial<BigInteger>> topt = new ArrayList<GenPolynomial<BigInteger>>(1); 642 topt.add(P); 643 OptimizedPolynomialList<BigInteger> opt = TermOrderOptimization.<BigInteger> optimizeTermOrder(pfac, 644 topt); 645 P = opt.list.get(0); 646 logger.info("optimized polynomial: " + P); 647 List<Integer> iperm = TermOrderOptimization.inversePermutation(opt.perm); 648 logger.info("optimize perm: " + opt.perm + ", de-optimize perm: " + iperm); 649 650 ExpVector degv = P.degreeVector(); 651 int[] donv = degv.dependencyOnVariables(); 652 List<GenPolynomial<BigInteger>> facs = null; 653 if (degv.length() == donv.length) { // all variables appear, hack for Hensel, TODO check 654 try { 655 logger.info("try factorsSquarefreeHensel: " + P); 656 facs = factorsSquarefreeHensel(P); 657 } catch (Exception e) { 658 logger.warn("exception " + e); 659 //e.printStackTrace(); 660 } 661 } else { // not all variables appear, remove unused variables, hack for Hensel, TODO check 662 GenPolynomial<BigInteger> pu = PolyUtil.<BigInteger> removeUnusedUpperVariables(P); 663 GenPolynomial<BigInteger> pl = PolyUtil.<BigInteger> removeUnusedLowerVariables(pu); // not useful 664 try { 665 logger.info("try factorsSquarefreeHensel: " + pl); 666 facs = factorsSquarefreeHensel(pu); 667 List<GenPolynomial<BigInteger>> fs = new ArrayList<GenPolynomial<BigInteger>>(facs.size()); 668 GenPolynomialRing<BigInteger> pf = P.ring; 669 GenPolynomialRing<BigInteger> pfu = pu.ring; 670 for (GenPolynomial<BigInteger> p : facs) { 671 GenPolynomial<BigInteger> pel = p.extendLower(pfu, 0, 0L); 672 GenPolynomial<BigInteger> pe = pel.extend(pf, 0, 0L); 673 fs.add(pe); 674 } 675 //System.out.println("fs = " + fs); 676 facs = fs; 677 } catch (Exception e) { 678 logger.warn("exception " + e); 679 //e.printStackTrace(); 680 } 681 } 682 if (facs == null) { 683 logger.info("factorsSquarefreeHensel not applicable or failed, reverting to Kronecker for: " + P); 684 facs = super.factorsSquarefree(P); 685 } 686 List<GenPolynomial<BigInteger>> iopt = TermOrderOptimization.<BigInteger> permutation(iperm, pfac, 687 facs); 688 logger.info("de-optimized polynomials: " + iopt); 689 facs = normalizeFactorization(iopt); 690 return facs; 691 } 692 693 694 /** 695 * GenPolynomial factorization of a multivariate squarefree polynomial, 696 * using Hensel lifting. 697 * @param P squarefree and primitive! (respectively monic) multivariate 698 * GenPolynomial over the integers. 699 * @return [p_1,...,p_k] with P = prod_{i=1,...,r} p_i. 700 */ 701 public List<GenPolynomial<BigInteger>> factorsSquarefreeHensel(GenPolynomial<BigInteger> P) { 702 if (P == null) { 703 throw new IllegalArgumentException(this.getClass().getName() + " P != null"); 704 } 705 GenPolynomialRing<BigInteger> pfac = P.ring; 706 if (pfac.nvar == 1) { 707 return baseFactorsSquarefree(P); 708 } 709 List<GenPolynomial<BigInteger>> factors = new ArrayList<GenPolynomial<BigInteger>>(); 710 if (P.isZERO()) { 711 return factors; 712 } 713 if (P.degreeVector().totalDeg() <= 1L) { 714 factors.add(P); 715 return factors; 716 } 717 GenPolynomial<BigInteger> pd = P; 718 //System.out.println("pd = " + pd); 719 // ldcf(pd) 720 BigInteger ac = pd.leadingBaseCoefficient(); 721 722 // factor leading coefficient as polynomial in the lowest! variable 723 GenPolynomialRing<GenPolynomial<BigInteger>> rnfac = pfac.recursive(pfac.nvar - 1); 724 GenPolynomial<GenPolynomial<BigInteger>> pr = PolyUtil.<BigInteger> recursive(rnfac, pd); 725 GenPolynomial<GenPolynomial<BigInteger>> prr = PolyUtil.<BigInteger> switchVariables(pr); 726 727 GenPolynomial<BigInteger> prrc = engine.recursiveContent(prr); // can have content wrt this variable 728 List<GenPolynomial<BigInteger>> cfactors = null; 729 if (!prrc.isONE()) { 730 prr = PolyUtil.<BigInteger> recursiveDivide(prr, prrc); 731 GenPolynomial<BigInteger> prrcu = prrc.extendLower(pfac, 0, 0L); // since switched vars 732 pd = PolyUtil.<BigInteger> basePseudoDivide(pd, prrcu); 733 logger.info("recursive content = " + prrc + ", new P = " + pd); 734 cfactors = factorsSquarefree(prrc); 735 List<GenPolynomial<BigInteger>> cff = new ArrayList<GenPolynomial<BigInteger>>(cfactors.size()); 736 for (GenPolynomial<BigInteger> fs : cfactors) { 737 GenPolynomial<BigInteger> fsp = fs.extendLower(pfac, 0, 0L); // since switched vars 738 cff.add(fsp); 739 } 740 cfactors = cff; 741 logger.info("cfactors = " + cfactors); 742 } 743 GenPolynomial<BigInteger> lprr = prr.leadingBaseCoefficient(); 744 //System.out.println("prr = " + prr); 745 logger.info("leading coeffcient = " + lprr); 746 boolean isMonic = false; // multivariate monic 747 if (lprr.isConstant()) { // isONE ? 748 isMonic = true; 749 } 750 SortedMap<GenPolynomial<BigInteger>, Long> lfactors = factors(lprr); 751 //System.out.println("lfactors = " + lfactors); 752 List<GenPolynomial<BigInteger>> lfacs = new ArrayList<GenPolynomial<BigInteger>>(lfactors.keySet()); 753 logger.info("leading coefficient factors = " + lfacs); 754 755 // search evaluation point and evaluate 756 GenPolynomialRing<BigInteger> cpfac = pfac; 757 GenPolynomial<BigInteger> pe = pd; 758 GenPolynomial<BigInteger> pep; 759 GenPolynomialRing<BigInteger> ccpfac = lprr.ring; 760 List<GenPolynomial<BigInteger>> ce = lfacs; 761 List<GenPolynomial<BigInteger>> cep = null; 762 List<BigInteger> cei = null; 763 List<BigInteger> dei = new ArrayList<BigInteger>(); 764 BigInteger pec = null; 765 BigInteger pecw = null; 766 BigInteger ped = null; 767 768 List<GenPolynomial<BigInteger>> ufactors = null; 769 List<TrialParts> tParts = new ArrayList<TrialParts>(); 770 List<GenPolynomial<BigInteger>> lf = null; 771 GenPolynomial<BigInteger> lpx = null; 772 List<GenPolynomial<BigInteger>> ln = null; 773 List<GenPolynomial<BigInteger>> un = null; 774 GenPolynomial<BigInteger> pes = null; 775 776 List<BigInteger> V = null; 777 long evStart = 0L; //3L * 5L; 778 List<Long> Evs = new ArrayList<Long>(pfac.nvar + 1); // Evs(0), Evs(1) unused 779 for (int j = 0; j <= pfac.nvar; j++) { 780 Evs.add(evStart); 781 } 782 final int trials = 4; 783 int countSeparate = 0; 784 final int COUNT_MAX = 50; 785 double ran = 1.001; // higher values not good 786 boolean isPrimitive = true; 787 boolean notLucky = true; 788 while (notLucky) { // for Wang's test 789 if (Math.abs(evStart) > 371L) { 790 logger.warn("no lucky evaluation point for: P = " + P + ", lprr = " + lprr + ", lfacs = " 791 + lfacs); 792 throw new RuntimeException("no lucky evaluation point found after " + Math.abs(evStart) 793 + " iterations"); 794 } 795 if (Math.abs(evStart) % 100L <= 3L) { 796 ran = ran * (Math.PI - 2.14); 797 } 798 //System.out.println("-------------------------------------------- Evs = " + Evs); 799 notLucky = false; 800 V = new ArrayList<BigInteger>(); 801 cpfac = pfac; 802 pe = pd; 803 ccpfac = lprr.ring; 804 ce = lfacs; 805 cep = null; 806 cei = null; 807 pec = null; 808 ped = null; 809 long vi = 0L; 810 for (int j = pfac.nvar; j > 1; j--) { 811 // evaluation up to univariate case 812 long degp = pe.degree(cpfac.nvar - 2); 813 cpfac = cpfac.contract(1); 814 ccpfac = ccpfac.contract(1); 815 //vi = evStart; // + j;//0L; //(long)(pfac.nvar-j); // 1L; 0 not so good for small p 816 vi = Evs.get(j); //evStart + j;//0L; //(long)(pfac.nvar-j); // 1L; 0 not so good for small p 817 BigInteger Vi; 818 819 // search evaluation point 820 boolean doIt = true; 821 Vi = null; 822 pep = null; 823 while (doIt) { 824 logger.info("vi(" + j + ") = " + vi); 825 Vi = new BigInteger(vi); 826 pep = PolyUtil.<BigInteger> evaluateMain(cpfac, pe, Vi); 827 //System.out.println("pep = " + pep); 828 // check lucky evaluation point 829 if (degp == pep.degree(cpfac.nvar - 1)) { 830 logger.info("pep = " + pep); 831 //System.out.println("deg(pe) = " + degp + ", deg(pep) = " + pep.degree(cpfac.nvar-1)); 832 // check squarefree 833 if (sengine.isSquarefree(pep)) { // cpfac.nvar == 1 && ?? no, must test on each variable 834 //if ( isNearlySquarefree(pep) ) { 835 //System.out.println("squarefeee(pep)"); // + pep); 836 doIt = false; //break; 837 } 838 } 839 if (vi > 0L) { 840 vi = -vi; 841 } else { 842 vi = 1L - vi; 843 } 844 } 845 //if ( !isMonic ) { 846 if (ccpfac.nvar >= 1) { 847 cep = PolyUtil.<BigInteger> evaluateMain(ccpfac, ce, Vi); 848 } else { 849 cei = PolyUtil.<BigInteger> evaluateMain(ccpfac.coFac, ce, Vi); 850 } 851 //} 852 int jj = (int) Math.round(ran + 0.52 * Math.random()); // j, random increment 853 //jj = 1; // ...4 test 854 //System.out.println("minimal jj = " + jj + ", vi " + vi); 855 if (vi > 0L) { 856 Evs.set(j, vi + jj); // record last tested value plus increment 857 evStart = vi + jj; 858 } else { 859 Evs.set(j, vi - jj); // record last tested value minus increment 860 evStart = vi - jj; 861 } 862 //evStart = vi+1L; 863 V.add(Vi); 864 pe = pep; 865 ce = cep; 866 } 867 //System.out.println("ce = " + ce + ", pe = " + pe); 868 pecw = engine.baseContent(pe); // original Wang 869 isPrimitive = pecw.isONE(); 870 ped = ccpfac.coFac.getONE(); 871 pec = pe.ring.coFac.getONE(); 872 //System.out.println("cei = " + cei + ", pecw = " + pecw); 873 if (!isMonic) { 874 if (countSeparate > COUNT_MAX) { 875 pec = pe.ring.coFac.getONE(); // hack is sometimes better 876 } else { 877 pec = pecw; 878 } 879 //pec = pecw; 880 //System.out.println("cei = " + cei + ", pec = " + pec + ", pe = " + pe); 881 if (lfacs.get(0).isConstant()) { 882 ped = cei.remove(0); 883 //lfacs.remove(0); // later 884 } 885 //System.out.println("lfacs = " + lfacs + ", cei = " + cei + ", ped = " + ped + ", pecw = " + pecw); 886 // test Wang's condition 887 dei = new ArrayList<BigInteger>(); 888 dei.add(pec.multiply(ped).abs()); // .abs() 889 int i = 1; 890 for (BigInteger ci : cei) { 891 if (ci.isZERO()) { 892 logger.info("condition (0) not met for cei = " + cei); // + ", dei = " + dei); 893 notLucky = true; 894 break; 895 } 896 BigInteger q = ci.abs(); 897 //System.out.println("q = " + q); 898 for (int ii = i - 1; ii >= 0; ii--) { 899 BigInteger r = dei.get(ii); 900 //System.out.println("r = " + r); 901 while (!r.isONE()) { 902 r = r.gcd(q); 903 q = q.divide(r); 904 //System.out.println("r = " + r + ", q = " + q); 905 } 906 } 907 dei.add(q); 908 if (q.isONE()) { 909 logger.info("condition (1) not met for dei = " + dei + ", cei = " + cei); 910 if (!testSeparate(cei, pecw)) { 911 countSeparate++; 912 if (countSeparate > COUNT_MAX) { 913 logger.info("too many inseparable evaluation points: " + countSeparate 914 + ", removing " + pecw); 915 } 916 } 917 notLucky = true; 918 break; 919 } 920 i++; 921 } 922 //System.out.println("dei = " + dei); 923 } 924 if (notLucky) { 925 continue; 926 } 927 logger.info("evaluation points = " + V + ", dei = " + dei); 928 //System.out.println("Evs = " + Evs); 929 logger.info("univariate polynomial = " + pe + ", pecw = " + pecw); 930 //pe = pe.abs(); 931 //ufactors = baseFactorsRadical(pe); //baseFactorsSquarefree(pe); wrong since not primitive 932 ufactors = baseFactorsSquarefree(pe.divide(pecw)); //wrong if not primitive 933 if (!pecw.isONE()) { 934 ufactors.add(0, cpfac.getONE().multiply(pecw)); 935 } 936 if (ufactors.size() <= 1) { 937 logger.info("irreducible univariate polynomial"); 938 factors.add(pd); // P 939 if (cfactors != null) { 940 cfactors.addAll(factors); 941 factors = cfactors; 942 } 943 return factors; 944 } 945 logger.info("univariate factors = " + ufactors); // + ", of " + pe); 946 //System.out.println("lfacs = " + lfacs); 947 //System.out.println("cei = " + cei); 948 //System.out.println("pecw = " + pecw); 949 950 // determine leading coefficient polynomials for factors 951 lf = new ArrayList<GenPolynomial<BigInteger>>(); 952 lpx = lprr.ring.getONE(); 953 for (GenPolynomial<BigInteger> unused : ufactors) { 954 lf.add(lprr.ring.getONE()); 955 } 956 //System.out.println("lf = " + lf); 957 if (!isMonic || !pecw.isONE()) { 958 if (lfacs.size() > 0 && lfacs.get(0).isConstant()) { 959 GenPolynomial<BigInteger> unused = lfacs.remove(0); 960 //BigInteger xxi = xx.leadingBaseCoefficient(); 961 //System.out.println("xx = " + xx + " == ped = " +ped); 962 } 963 for (int i = ufactors.size() - 1; i >= 0; i--) { 964 GenPolynomial<BigInteger> pp = ufactors.get(i); 965 BigInteger ppl = pp.leadingBaseCoefficient(); 966 //System.out.println("ppl = " + ppl + ", pp = " + pp); 967 ppl = ppl.multiply(pec); // content 968 GenPolynomial<BigInteger> lfp = lf.get(i); 969 int ii = 0; 970 for (BigInteger ci : cei) { 971 //System.out.println("ci = " + ci + ", lfp = " + lfp + ", lfacs.get(ii) = " + lfacs.get(ii)); 972 if (ci.abs().isONE()) { 973 System.out.println("ppl = " + ppl + ", ci = " + ci + ", lfp = " + lfp 974 + ", lfacs.get(ii) = " + lfacs.get(ii)); 975 throw new RuntimeException("something is wrong, ci is a unit"); 976 //notLucky = true; 977 } 978 while (ppl.remainder(ci).isZERO() && lfacs.size() > ii) { 979 ppl = ppl.divide(ci); 980 lfp = lfp.multiply(lfacs.get(ii)); 981 } 982 ii++; 983 } 984 //System.out.println("ppl = " + ppl + ", lfp = " + lfp); 985 lfp = lfp.multiply(ppl); 986 lf.set(i, lfp); 987 } 988 // adjust if pec != 1 989 pec = pecw; 990 lpx = Power.<GenPolynomial<BigInteger>> multiply(lprr.ring, lf); // test only, not used 991 //System.out.println("lpx = " + lpx); 992 if (!lprr.degreeVector().equals(lpx.degreeVector())) { 993 logger.info("deg(lprr) != deg(lpx): lprr = " + lprr + ", lpx = " + lpx); 994 notLucky = true; 995 continue; 996 } 997 if (!pec.isONE()) { // content, was always false by hack 998 // evaluate factors of ldcf 999 List<GenPolynomial<BigInteger>> lfe = lf; 1000 List<BigInteger> lfei = null; 1001 ccpfac = lprr.ring; 1002 for (int j = lprr.ring.nvar; j > 0; j--) { 1003 ccpfac = ccpfac.contract(1); 1004 BigInteger Vi = V.get(lprr.ring.nvar - j); 1005 if (ccpfac.nvar >= 1) { 1006 lfe = PolyUtil.<BigInteger> evaluateMain(ccpfac, lfe, Vi); 1007 } else { 1008 lfei = PolyUtil.<BigInteger> evaluateMain(ccpfac.coFac, lfe, Vi); 1009 } 1010 } 1011 //System.out.println("lfe = " + lfe + ", lfei = " + lfei + ", V = " + V); 1012 1013 ln = new ArrayList<GenPolynomial<BigInteger>>(lf.size()); 1014 un = new ArrayList<GenPolynomial<BigInteger>>(lf.size()); 1015 for (int jj = 0; jj < lf.size(); jj++) { 1016 GenPolynomial<BigInteger> up = ufactors.get(jj); 1017 BigInteger ui = up.leadingBaseCoefficient(); 1018 BigInteger li = lfei.get(jj); 1019 BigInteger di = ui.gcd(li).abs(); 1020 BigInteger udi = ui.divide(di); 1021 BigInteger ldi = li.divide(di); 1022 GenPolynomial<BigInteger> lp = lf.get(jj); 1023 GenPolynomial<BigInteger> lpd = lp.multiply(udi); 1024 GenPolynomial<BigInteger> upd = up.multiply(ldi); 1025 if (pec.isONE()) { 1026 ln.add(lp); 1027 un.add(up); 1028 } else { 1029 ln.add(lpd); 1030 un.add(upd); 1031 BigInteger pec1 = pec.divide(ldi); 1032 //System.out.println("pec = " + pec + ", pec1 = " + pec1); 1033 pec = pec1; 1034 } 1035 } 1036 if (!lf.equals(ln) || !un.equals(ufactors)) { 1037 logger.debug("!lf.equals(ln) || !un.equals(ufactors)"); 1038 //System.out.println("pe = " + pe); 1039 //System.out.println("#ln = " + ln + ", #lf = " + lf); 1040 //System.out.println("#un = " + un + ", #ufactors = " + ufactors); 1041 //lf = ln; 1042 //ufactors = un; 1043 // adjust pe 1044 } 1045 if (!pec.isONE()) { // still not 1 1046 ln = new ArrayList<GenPolynomial<BigInteger>>(lf.size()); 1047 un = new ArrayList<GenPolynomial<BigInteger>>(lf.size()); 1048 pes = pe; 1049 for (int jj = 0; jj < lf.size(); jj++) { 1050 GenPolynomial<BigInteger> up = ufactors.get(jj); 1051 GenPolynomial<BigInteger> lp = lf.get(jj); 1052 //System.out.println("up = " + up + ", lp = " + lp); 1053 if (!up.isConstant()) { 1054 up = up.multiply(pec); 1055 } 1056 lp = lp.multiply(pec); 1057 if (jj != 0) { 1058 pes = pes.multiply(pec); 1059 } 1060 un.add(up); 1061 ln.add(lp); 1062 } 1063 if (pes.equals(Power.<GenPolynomial<BigInteger>> multiply(pe.ring, un))) { 1064 //System.out.println("*pe = " + pes + ", pec = " + pec); 1065 //ystem.out.println("*ln = " + ln + ", *lf = " + lf); 1066 //System.out.println("*un = " + un + ", *ufactors = " + ufactors); 1067 //System.out.println("*pe == prod(un) "); 1068 isPrimitive = false; 1069 //pe = pes; 1070 //lf = ln; 1071 //ufactors = un; 1072 } else { 1073 //System.out.println("*pe != prod(un): " + Power.<GenPolynomial<BigInteger>> multiply(pe.ring,un)); 1074 } 1075 } 1076 } 1077 if (notLucky) { 1078 continue; 1079 } 1080 logger.info("distributed factors of leading coefficient = " + lf); 1081 lpx = Power.<GenPolynomial<BigInteger>> multiply(lprr.ring, lf); 1082 if (!lprr.abs().equals(lpx.abs())) { // not correctly distributed 1083 if (!lprr.degreeVector().equals(lpx.degreeVector())) { 1084 logger.info("lprr != lpx: lprr = " + lprr + ", lpx = " + lpx); 1085 notLucky = true; 1086 } 1087 } 1088 } // end determine leading coefficients for factors 1089 1090 if (!notLucky) { 1091 TrialParts tp = null; 1092 if (isPrimitive) { 1093 tp = new TrialParts(V, pe, ufactors, cei, lf); 1094 } else { 1095 tp = new TrialParts(V, pes, un, cei, ln); 1096 } 1097 //System.out.println("trialParts = " + tp); 1098 if (tp.univPoly != null) { 1099 if (tp.ldcfEval.size() != 0) { 1100 tParts.add(tp); 1101 } 1102 } 1103 if (tParts.size() < trials) { 1104 notLucky = true; 1105 } 1106 } 1107 } // end notLucky loop 1108 1109 // search TrialParts with shortest factorization of univariate polynomial 1110 int min = Integer.MAX_VALUE; 1111 TrialParts tpmin = null; 1112 for (TrialParts tp : tParts) { 1113 logger.info("tp.univFactors.size() = " + tp.univFactors.size()); 1114 if (tp.univFactors.size() < min) { 1115 min = tp.univFactors.size(); 1116 tpmin = tp; 1117 } 1118 } 1119 for (TrialParts tp : tParts) { 1120 //logger.info("tp.univFactors.get(0) = " + tp.univFactors.get(0)); 1121 if (tp.univFactors.size() == min) { 1122 if (!tp.univFactors.get(0).isConstant()) { 1123 tpmin = tp; 1124 break; 1125 } 1126 } 1127 } 1128 // set to (first) shortest 1129 V = tpmin.evalPoints; 1130 pe = tpmin.univPoly; 1131 ufactors = tpmin.univFactors; 1132 cei = tpmin.ldcfEval; // unused 1133 lf = tpmin.ldcfFactors; 1134 logger.info("iterations = " + Math.abs(evStart)); 1135 logger.info("minimal trial = " + tpmin); 1136 1137 GenPolynomialRing<BigInteger> ufac = pe.ring; 1138 1139 //initialize prime list 1140 PrimeList primes = new PrimeList(PrimeList.Range.medium); // PrimeList.Range.medium); 1141 Iterator<java.math.BigInteger> primeIter = primes.iterator(); 1142 int pn = 50; //primes.size(); 1143 BigInteger ae = pe.leadingBaseCoefficient(); 1144 GenPolynomial<MOD> Pm = null; 1145 ModularRingFactory<MOD> cofac = null; 1146 GenPolynomialRing<MOD> mufac = null; 1147 1148 // search lucky prime 1149 for (int i = 0; i < 11; i++) { // prime meta loop 1150 //for ( int i = 0; i < 1; i++ ) { // meta loop 1151 java.math.BigInteger p = null; //new java.math.BigInteger("19"); //primes.next(); 1152 // 2 small, 5 medium and 4 large size primes 1153 if (i == 0) { // medium size 1154 primes = new PrimeList(PrimeList.Range.medium); 1155 primeIter = primes.iterator(); 1156 } 1157 if (i == 5) { // small size 1158 primes = new PrimeList(PrimeList.Range.small); 1159 primeIter = primes.iterator(); 1160 p = primeIter.next(); // 2 1161 p = primeIter.next(); // 3 1162 p = primeIter.next(); // 5 1163 p = primeIter.next(); // 7 1164 } 1165 if (i == 7) { // large size 1166 primes = new PrimeList(PrimeList.Range.large); 1167 primeIter = primes.iterator(); 1168 } 1169 int pi = 0; 1170 while (pi < pn && primeIter.hasNext()) { 1171 p = primeIter.next(); 1172 logger.info("prime = " + p); 1173 // initialize coefficient factory and map normalization factor and polynomials 1174 ModularRingFactory<MOD> cf = null; 1175 if (ModLongRing.MAX_LONG.compareTo(p) > 0) { 1176 cf = (ModularRingFactory) new ModLongRing(p, true); 1177 } else { 1178 cf = (ModularRingFactory) new ModIntegerRing(p, true); 1179 } 1180 MOD nf = cf.fromInteger(ae.getVal()); 1181 if (nf.isZERO()) { 1182 continue; 1183 } 1184 mufac = new GenPolynomialRing<MOD>(cf, ufac); 1185 //System.out.println("mufac = " + mufac.toScript()); 1186 Pm = PolyUtil.<MOD> fromIntegerCoefficients(mufac, pe); 1187 //System.out.println("Pm = " + Pm); 1188 if (!mfactor.isSquarefree(Pm)) { 1189 continue; 1190 } 1191 cofac = cf; 1192 break; 1193 } 1194 if (cofac != null) { 1195 break; 1196 } 1197 } // end prime meta loop 1198 if (cofac == null) { // no lucky prime found 1199 throw new RuntimeException("giving up on Hensel preparation, no lucky prime found"); 1200 } 1201 logger.info("lucky prime = " + cofac.getIntegerModul()); 1202 if (logger.isDebugEnabled()) { 1203 logger.debug("univariate modulo p: = " + Pm); 1204 } 1205 1206 // coefficient bound 1207 BigInteger an = pd.maxNorm(); 1208 BigInteger mn = an.multiply(ac.abs()).multiply(new BigInteger(2L)); 1209 long k = Power.logarithm(cofac.getIntegerModul(), mn) + 1L; 1210 //System.out.println("mn = " + mn + ", k = " +k); 1211 1212 BigInteger q = Power.positivePower(cofac.getIntegerModul(), k); 1213 ModularRingFactory<MOD> muqfac; 1214 if (ModLongRing.MAX_LONG.compareTo(q.getVal()) > 0) { 1215 muqfac = (ModularRingFactory) new ModLongRing(q.getVal()); 1216 } else { 1217 muqfac = (ModularRingFactory) new ModIntegerRing(q.getVal()); 1218 } 1219 //System.out.println("muqfac = " + muqfac); 1220 GenPolynomialRing<MOD> mucpfac = new GenPolynomialRing<MOD>(muqfac, ufac); 1221 1222 List<GenPolynomial<MOD>> muqfactors = PolyUtil.<MOD> fromIntegerCoefficients(mucpfac, ufactors); 1223 GenPolynomial<MOD> peqq = PolyUtil.<MOD> fromIntegerCoefficients(mucpfac, pe); 1224 if (debug) { 1225 if (!mfactor.isFactorization(peqq, muqfactors)) { // should not happen 1226 System.out.println("muqfactors = " + muqfactors); 1227 System.out.println("peqq = " + peqq); 1228 throw new RuntimeException("something is wrong, no modular p^k factorization"); 1229 } 1230 } 1231 logger.info("univariate modulo p^k: " + peqq + " = " + muqfactors); 1232 1233 // convert C from Z[...] to Z_q[...] 1234 GenPolynomialRing<MOD> qcfac = new GenPolynomialRing<MOD>(muqfac, pd.ring); 1235 GenPolynomial<MOD> pq = PolyUtil.<MOD> fromIntegerCoefficients(qcfac, pd); 1236 //System.out.println("pd = " + pd); 1237 logger.info("multivariate modulo p^k: " + pq); 1238 1239 List<MOD> Vm = new ArrayList<MOD>(V.size()); 1240 for (BigInteger v : V) { 1241 MOD vm = muqfac.fromInteger(v.getVal()); 1242 Vm.add(vm); 1243 } 1244 //System.out.println("Vm = " + Vm); 1245 1246 // Hensel lifting of factors 1247 List<GenPolynomial<MOD>> mlift; 1248 try { 1249 mlift = HenselMultUtil.<MOD> liftHensel(pd, pq, muqfactors, Vm, k, lf); 1250 logger.info("mlift = " + mlift); 1251 } catch (NoLiftingException nle) { 1252 //System.out.println("exception : " + nle); 1253 //nle.printStackTrace(); 1254 mlift = new ArrayList<GenPolynomial<MOD>>(); 1255 throw new RuntimeException(nle); 1256 } catch (ArithmeticException aex) { 1257 //System.out.println("exception : " + aex); 1258 //aex.printStackTrace(); 1259 mlift = new ArrayList<GenPolynomial<MOD>>(); 1260 throw aex; 1261 } 1262 if (mlift.size() <= 1) { // irreducible mod I, p^k, can this happen? 1263 logger.info("modular lift size == 1: " + mlift); 1264 factors.add(pd); // P 1265 if (cfactors != null) { 1266 cfactors.addAll(factors); 1267 factors = cfactors; 1268 } 1269 return factors; 1270 } 1271 1272 // combine trial factors 1273 GenPolynomialRing<MOD> mfac = mlift.get(0).ring; 1274 int dl = (mlift.size() + 1) / 2; 1275 GenPolynomial<BigInteger> u = P; 1276 long deg = (u.degree() + 1L) / 2L; 1277 1278 GenPolynomial<BigInteger> ui = pd; 1279 for (int j = 1; j <= dl; j++) { 1280 //System.out.println("j = " + j + ", dl = " + dl + ", mlift = " + mlift); 1281 KsubSet<GenPolynomial<MOD>> subs = new KsubSet<GenPolynomial<MOD>>(mlift, j); 1282 for (List<GenPolynomial<MOD>> flist : subs) { 1283 //System.out.println("degreeSum = " + degreeSum(flist)); 1284 GenPolynomial<MOD> mtrial = Power.<GenPolynomial<MOD>> multiply(mfac, flist); 1285 if (mtrial.degree() > deg) { // this test is sometimes wrong 1286 logger.info("degree > deg " + deg + ", degree = " + mtrial.degree()); 1287 //continue; 1288 } 1289 GenPolynomial<BigInteger> trial = PolyUtil.integerFromModularCoefficients(pfac, mtrial); 1290 trial = engine.basePrimitivePart(trial); 1291 //if ( ! isPrimitive ) { 1292 //} 1293 if (debug) { 1294 logger.info("trial = " + trial); // + ", mtrial = " + mtrial); 1295 } 1296 if (PolyUtil.<BigInteger> baseSparsePseudoRemainder(ui, trial).isZERO()) { 1297 logger.info("successful trial = " + trial); 1298 factors.add(trial); 1299 ui = PolyUtil.<BigInteger> basePseudoDivide(ui, trial); 1300 //System.out.println("ui = " + ui); 1301 mlift = removeOnce(mlift, flist); 1302 logger.info("new mlift= " + mlift); 1303 //System.out.println("dl = " + dl); 1304 if (mlift.size() > 1) { 1305 dl = (mlift.size() + 1) / 2; 1306 j = 0; // since j++ 1307 break; 1308 } 1309 logger.info("last factor = " + ui); 1310 factors.add(ui); 1311 if (cfactors != null) { 1312 cfactors.addAll(factors); 1313 factors = cfactors; 1314 } 1315 return normalizeFactorization(factors); 1316 } 1317 } 1318 } 1319 if (!ui.isONE() && !ui.equals(pd)) { 1320 logger.info("rest factor = " + ui); 1321 // pp(ui) ?? no ?? 1322 factors.add(ui); 1323 } 1324 if (factors.size() == 0) { 1325 logger.info("irreducible P = " + P); 1326 factors.add(pd); // P 1327 } 1328 if (cfactors != null) { 1329 cfactors.addAll(factors); 1330 factors = cfactors; 1331 } 1332 return normalizeFactorization(factors); 1333 } 1334 1335 1336 /** 1337 * Test if b has a prime factor different to the elements of A. 1338 * @param A list of integer with at least one different prime factor. 1339 * @param b integer to test with A. 1340 * @return true, if b hase a prime factor different to elements of A 1341 */ 1342 boolean testSeparate(List<BigInteger> A, BigInteger b) { 1343 int i = 0; 1344 List<BigInteger> gei = new ArrayList<BigInteger>(A.size()); 1345 for (BigInteger c : A) { 1346 BigInteger g = c.gcd(b).abs(); 1347 gei.add(g); 1348 if (!g.isONE()) { 1349 i++; 1350 } 1351 } 1352 //if ( i >= 1 ) { 1353 //System.out.println("gei = " + gei + ", cei = " + cei + ", pec(w) = " + pec); 1354 //} 1355 return (i <= 1); 1356 } 1357 1358 1359 // not useable 1360 boolean isNearlySquarefree(GenPolynomial<BigInteger> P) { // unused 1361 // in main variable 1362 GenPolynomialRing<BigInteger> pfac = P.ring; 1363 if (pfac.nvar >= 0) { // allways true 1364 return sengine.isSquarefree(P); 1365 } 1366 GenPolynomialRing<GenPolynomial<BigInteger>> rfac = pfac.recursive(1); 1367 GenPolynomial<GenPolynomial<BigInteger>> Pr = PolyUtil.<BigInteger> recursive(rfac, P); 1368 GenPolynomial<GenPolynomial<BigInteger>> Ps = PolyUtil.<BigInteger> recursiveDeriviative(Pr); 1369 System.out.println("Pr = " + Pr); 1370 System.out.println("Ps = " + Ps); 1371 GenPolynomial<GenPolynomial<BigInteger>> g = engine.recursiveUnivariateGcd(Pr, Ps); 1372 System.out.println("g_m = " + g); 1373 if (!g.isONE()) { 1374 return false; 1375 } 1376 // in lowest variable 1377 rfac = pfac.recursive(pfac.nvar - 1); 1378 Pr = PolyUtil.<BigInteger> recursive(rfac, P); 1379 Pr = PolyUtil.<BigInteger> switchVariables(Pr); 1380 Ps = PolyUtil.<BigInteger> recursiveDeriviative(Pr); 1381 System.out.println("Pr = " + Pr); 1382 System.out.println("Ps = " + Ps); 1383 g = engine.recursiveUnivariateGcd(Pr, Ps); 1384 System.out.println("g_1 = " + g); 1385 if (!g.isONE()) { 1386 return false; 1387 } 1388 return true; 1389 } 1390 1391} 1392 1393 1394/** 1395 * Container for factorization trial lifting parameters. 1396 */ 1397class TrialParts { 1398 1399 1400 /** 1401 * evaluation points 1402 */ 1403 public final List<BigInteger> evalPoints; 1404 1405 1406 /** 1407 * univariate polynomial 1408 */ 1409 public final GenPolynomial<BigInteger> univPoly; 1410 1411 1412 /** 1413 * irreducible factors of univariate polynomial 1414 */ 1415 public final List<GenPolynomial<BigInteger>> univFactors; 1416 1417 1418 /** 1419 * irreducible factors of leading coefficient 1420 */ 1421 public final List<GenPolynomial<BigInteger>> ldcfFactors; 1422 1423 1424 /** 1425 * evaluated factors of leading coefficient factors by evaluation points 1426 */ 1427 public final List<BigInteger> ldcfEval; 1428 1429 1430 /** 1431 * Constructor. 1432 * @param ev evaluation points. 1433 * @param up univariate polynomial. 1434 * @param uf irreducible factors of up. 1435 * @param le irreducible factors of leading coefficient. 1436 * @param lf evaluated le by evaluation points. 1437 */ 1438 public TrialParts(List<BigInteger> ev, GenPolynomial<BigInteger> up, List<GenPolynomial<BigInteger>> uf, 1439 List<BigInteger> le, List<GenPolynomial<BigInteger>> lf) { 1440 evalPoints = ev; 1441 univPoly = up; 1442 univFactors = uf; 1443 //ldcfPoly = lp; 1444 ldcfFactors = lf; 1445 ldcfEval = le; 1446 } 1447 1448 1449 /** 1450 * @see java.lang.Object#toString() 1451 */ 1452 @Override 1453 public String toString() { 1454 StringBuffer sb = new StringBuffer(); 1455 sb.append("TrialParts["); 1456 sb.append("evalPoints = " + evalPoints); 1457 sb.append(", univPoly = " + univPoly); 1458 sb.append(", univFactors = " + univFactors); 1459 sb.append(", ldcfEval = " + ldcfEval); 1460 sb.append(", ldcfFactors = " + ldcfFactors); 1461 sb.append("]"); 1462 return sb.toString(); 1463 } 1464 1465}