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