001 /*
002 * $Id: FactorInteger.java 3753 2011-08-27 20:34:30Z kredel $
003 */
004
005 package edu.jas.ufd;
006
007
008 import java.util.ArrayList;
009 import java.util.BitSet;
010 import java.util.Iterator;
011 import java.util.List;
012 import java.util.SortedMap;
013
014 import org.apache.log4j.Logger;
015
016 import edu.jas.arith.BigInteger;
017 import edu.jas.arith.ModIntegerRing;
018 import edu.jas.arith.ModLongRing;
019 import edu.jas.arith.Modular;
020 import edu.jas.arith.ModularRingFactory;
021 import edu.jas.arith.PrimeList;
022 import edu.jas.poly.ExpVector;
023 import edu.jas.poly.GenPolynomial;
024 import edu.jas.poly.GenPolynomialRing;
025 import edu.jas.poly.PolyUtil;
026 import edu.jas.structure.GcdRingElem;
027 import edu.jas.structure.Power;
028 import edu.jas.structure.RingElem;
029 import edu.jas.structure.RingFactory;
030 import edu.jas.util.KsubSet;
031
032
033 /**
034 * Integer coefficients factorization algorithms. This class implements
035 * factorization methods for polynomials over integers.
036 * @author Heinz Kredel
037 */
038
039 /**
040 * @author kredel
041 *
042 * @param <MOD>
043 */
044 public class FactorInteger<MOD extends GcdRingElem<MOD> & Modular> extends FactorAbstract<BigInteger> {
045
046
047 private static final Logger logger = Logger.getLogger(FactorInteger.class);
048
049
050 private final boolean debug = true || logger.isDebugEnabled();
051
052
053 /**
054 * Factorization engine for modular base coefficients.
055 */
056 protected final FactorAbstract<MOD> mfactor;
057
058
059 /**
060 * Gcd engine for modular base coefficients.
061 */
062 protected final GreatestCommonDivisorAbstract<MOD> mengine;
063
064
065 /**
066 * No argument constructor.
067 */
068 public FactorInteger() {
069 this(BigInteger.ONE);
070 }
071
072
073 /**
074 * Constructor.
075 * @param cfac coefficient ring factory.
076 */
077 public FactorInteger(RingFactory<BigInteger> cfac) {
078 super(cfac);
079 ModularRingFactory<MOD> mcofac = (ModularRingFactory<MOD>) (Object) new ModLongRing(13, true); // hack
080 mfactor = FactorFactory.getImplementation(mcofac); //new FactorModular(mcofac);
081 mengine = GCDFactory.getImplementation(mcofac);
082 //mengine = GCDFactory.getProxy(mcofac);
083 }
084
085
086 /**
087 * GenPolynomial base factorization of a squarefree polynomial.
088 * @param P squarefree and primitive! GenPolynomial.
089 * @return [p_1,...,p_k] with P = prod_{i=1, ..., k} p_i.
090 */
091 @SuppressWarnings("unchecked")
092 @Override
093 public List<GenPolynomial<BigInteger>> baseFactorsSquarefree(GenPolynomial<BigInteger> P) {
094 if (P == null) {
095 throw new IllegalArgumentException(this.getClass().getName() + " P == null");
096 }
097 List<GenPolynomial<BigInteger>> factors = new ArrayList<GenPolynomial<BigInteger>>();
098 if (P.isZERO()) {
099 return factors;
100 }
101 if (P.isONE()) {
102 factors.add(P);
103 return factors;
104 }
105 GenPolynomialRing<BigInteger> pfac = P.ring;
106 if (pfac.nvar > 1) {
107 throw new IllegalArgumentException(this.getClass().getName() + " only for univariate polynomials");
108 }
109 if (!engine.baseContent(P).isONE()) {
110 throw new IllegalArgumentException(this.getClass().getName() + " P not primitive");
111 }
112 if (P.degree(0) <= 1L) { // linear is irreducible
113 factors.add(P);
114 return factors;
115 }
116 // compute norm
117 BigInteger an = P.maxNorm();
118 BigInteger ac = P.leadingBaseCoefficient();
119 //compute factor coefficient bounds
120 ExpVector degv = P.degreeVector();
121 int degi = (int) P.degree(0);
122 BigInteger M = an.multiply(PolyUtil.factorBound(degv));
123 M = M.multiply(ac.abs().multiply(ac.fromInteger(8)));
124 //System.out.println("M = " + M);
125 //M = M.multiply(M); // test
126
127 //initialize prime list and degree vector
128 PrimeList primes = new PrimeList(PrimeList.Range.small);
129 int pn = 30; //primes.size();
130 ModularRingFactory<MOD> cofac = null;
131 GenPolynomial<MOD> am = null;
132 GenPolynomialRing<MOD> mfac = null;
133 final int TT = 5; // 7
134 List<GenPolynomial<MOD>>[] modfac = new List[TT];
135 List<GenPolynomial<BigInteger>>[] intfac = new List[TT];
136 BigInteger[] plist = new BigInteger[TT];
137 List<GenPolynomial<MOD>> mlist = null;
138 List<GenPolynomial<BigInteger>> ilist = null;
139 int i = 0;
140 if (debug) {
141 logger.debug("an = " + an);
142 logger.debug("ac = " + ac);
143 logger.debug("M = " + M);
144 logger.info("degv = " + degv);
145 }
146 Iterator<java.math.BigInteger> pit = primes.iterator();
147 pit.next(); // skip p = 2
148 pit.next(); // skip p = 3
149 MOD nf = null;
150 for (int k = 0; k < TT; k++) {
151 if (k == TT - 1) { // -2
152 primes = new PrimeList(PrimeList.Range.medium);
153 pit = primes.iterator();
154 }
155 if (k == TT + 1) { // -1
156 primes = new PrimeList(PrimeList.Range.large);
157 pit = primes.iterator();
158 }
159 while (pit.hasNext()) {
160 java.math.BigInteger p = pit.next();
161 //System.out.println("next run ++++++++++++++++++++++++++++++++++");
162 if (++i >= pn) {
163 logger.error("prime list exhausted, pn = " + pn);
164 throw new ArithmeticException("prime list exhausted");
165 }
166 if (ModLongRing.MAX_LONG.compareTo(p) > 0) {
167 cofac = (ModularRingFactory) new ModLongRing(p, true);
168 } else {
169 cofac = (ModularRingFactory) new ModIntegerRing(p, true);
170 }
171 logger.info("prime = " + cofac);
172 nf = cofac.fromInteger(ac.getVal());
173 if (nf.isZERO()) {
174 logger.info("unlucky prime (nf) = " + p);
175 //System.out.println("unlucky prime (nf) = " + p);
176 continue;
177 }
178 // initialize polynomial factory and map polynomial
179 mfac = new GenPolynomialRing<MOD>(cofac, pfac);
180 am = PolyUtil.<MOD> fromIntegerCoefficients(mfac, P);
181 if (!am.degreeVector().equals(degv)) { // allways true
182 logger.info("unlucky prime (deg) = " + p);
183 //System.out.println("unlucky prime (deg) = " + p);
184 continue;
185 }
186 GenPolynomial<MOD> ap = PolyUtil.<MOD> baseDeriviative(am);
187 if (ap.isZERO()) {
188 logger.info("unlucky prime (a')= " + p);
189 //System.out.println("unlucky prime (a')= " + p);
190 continue;
191 }
192 GenPolynomial<MOD> g = mengine.baseGcd(am, ap);
193 if (g.isONE()) {
194 logger.info("**lucky prime = " + p);
195 //System.out.println("**lucky prime = " + p);
196 break;
197 }
198 }
199 // now am is squarefree mod p, make monic and factor mod p
200 if (!nf.isONE()) {
201 //System.out.println("nf = " + nf);
202 am = am.divide(nf); // make monic
203 }
204 mlist = mfactor.baseFactorsSquarefree(am);
205 if (logger.isInfoEnabled()) {
206 logger.info("modlist = " + mlist);
207 }
208 if (mlist.size() <= 1) {
209 factors.add(P);
210 return factors;
211 }
212 if (!nf.isONE()) {
213 GenPolynomial<MOD> mp = mfac.getONE(); //mlist.get(0);
214 //System.out.println("mp = " + mp);
215 mp = mp.multiply(nf);
216 //System.out.println("mp = " + mp);
217 mlist.add(0, mp); // set(0,mp);
218 }
219 modfac[k] = mlist;
220 plist[k] = cofac.getIntegerModul(); // p
221 }
222
223 // search shortest factor list
224 int min = Integer.MAX_VALUE;
225 BitSet AD = null;
226 for (int k = 0; k < TT; k++) {
227 List<ExpVector> ev = PolyUtil.<MOD> leadingExpVector(modfac[k]);
228 BitSet D = factorDegrees(ev, degi);
229 if (AD == null) {
230 AD = D;
231 } else {
232 AD.and(D);
233 }
234 int s = modfac[k].size();
235 logger.info("mod(" + plist[k] + ") #s = " + s + ", D = " + D /*+ ", lt = " + ev*/);
236 //System.out.println("mod s = " + s);
237 if (s < min) {
238 min = s;
239 mlist = modfac[k];
240 }
241 }
242 logger.info("min = " + min + ", AD = " + AD);
243 if (mlist.size() <= 1) {
244 logger.info("mlist.size() = 1");
245 factors.add(P);
246 return factors;
247 }
248 if (AD.cardinality() <= 2) { // only one possible factor
249 logger.info("degree set cardinality = " + AD.cardinality());
250 factors.add(P);
251 return factors;
252 }
253
254 boolean allLists = false; //true; //false;
255 if (allLists) {
256 // try each factor list
257 for (int k = 0; k < TT; k++) {
258 mlist = modfac[k];
259 if (debug) {
260 logger.info("lifting from " + mlist);
261 }
262 if (false && P.leadingBaseCoefficient().isONE()) {
263 factors = searchFactorsMonic(P, M, mlist, AD); // does now work in all cases
264 if (factors.size() == 1) {
265 factors = searchFactorsNonMonic(P, M, mlist, AD);
266 }
267 } else {
268 factors = searchFactorsNonMonic(P, M, mlist, AD);
269 }
270 intfac[k] = factors;
271 }
272 } else {
273 // try only shortest factor list
274 if (debug) {
275 logger.info("lifting shortest from " + mlist);
276 }
277 if (true && P.leadingBaseCoefficient().isONE()) {
278 long t = System.currentTimeMillis();
279 try {
280 mlist = PolyUtil.<MOD> monic(mlist);
281 factors = searchFactorsMonic(P, M, mlist, AD); // does now work in all cases
282 t = System.currentTimeMillis() - t;
283 //System.out.println("monic time = " + t);
284 if (false && debug) {
285 t = System.currentTimeMillis();
286 List<GenPolynomial<BigInteger>> fnm = searchFactorsNonMonic(P, M, mlist, AD);
287 t = System.currentTimeMillis() - t;
288 System.out.println("non monic time = " + t);
289 if (debug) {
290 if (!factors.equals(fnm)) {
291 System.out.println("monic factors = " + factors);
292 System.out.println("non monic factors = " + fnm);
293 }
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 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 = 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> basePseudoRemainder(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 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> basePseudoRemainder(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 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 ExpVector degv = P.degreeVector();
638 int[] donv = degv.dependencyOnVariables();
639 List<GenPolynomial<BigInteger>> facs = null;
640 if (degv.length() == donv.length) { // all variables appear, hack for Hensel, TODO check
641 try {
642 logger.info("try factorsSquarefreeHensel: " + P);
643 facs = factorsSquarefreeHensel(P);
644 } catch (Exception e) {
645 logger.warn("exception " + e);
646 //e.printStackTrace();
647 }
648 } else { // not all variables appear, remove unused variables, hack for Hensel, TODO check
649 GenPolynomial<BigInteger> pu = PolyUtil.<BigInteger> removeUnusedUpperVariables(P);
650 GenPolynomial<BigInteger> pl = PolyUtil.<BigInteger> removeUnusedLowerVariables(pu);
651 try {
652 logger.info("try factorsSquarefreeHensel: " + pl);
653 facs = factorsSquarefreeHensel(pu);
654 List<GenPolynomial<BigInteger>> fs = new ArrayList<GenPolynomial<BigInteger>>(facs.size());
655 GenPolynomialRing<BigInteger> pf = P.ring;
656 GenPolynomialRing<BigInteger> pfu = pu.ring;
657 for (GenPolynomial<BigInteger> p : facs) {
658 GenPolynomial<BigInteger> pel = p.extendLower(pfu, 0, 0L);
659 GenPolynomial<BigInteger> pe = pel.extend(pf, 0, 0L);
660 fs.add(pe);
661 }
662 //System.out.println("fs = " + fs);
663 facs = fs;
664 } catch (Exception e) {
665 logger.warn("exception " + e);
666 //e.printStackTrace();
667 }
668 }
669 if (facs == null) {
670 logger.info("factorsSquarefreeHensel not applicable or failed, reverting to Kronecker for: " + P);
671 facs = super.factorsSquarefree(P);
672 }
673 return facs;
674 }
675
676
677 /**
678 * GenPolynomial factorization of a multivariate squarefree polynomial,
679 * using Hensel lifting.
680 * @param P squarefree and primitive! (respectively monic) multivariate
681 * GenPolynomial over the integers.
682 * @return [p_1,...,p_k] with P = prod_{i=1,...,r} p_i.
683 */
684 public List<GenPolynomial<BigInteger>> factorsSquarefreeHensel(GenPolynomial<BigInteger> P) {
685 if (P == null) {
686 throw new IllegalArgumentException(this.getClass().getName() + " P != null");
687 }
688 GenPolynomialRing<BigInteger> pfac = P.ring;
689 if (pfac.nvar == 1) {
690 return baseFactorsSquarefree(P);
691 }
692 List<GenPolynomial<BigInteger>> factors = new ArrayList<GenPolynomial<BigInteger>>();
693 if (P.isZERO()) {
694 return factors;
695 }
696 if (P.degreeVector().totalDeg() <= 1L) {
697 factors.add(P);
698 return factors;
699 }
700 GenPolynomial<BigInteger> pd = P;
701 //System.out.println("pd = " + pd);
702 // ldcf(pd)
703 BigInteger ac = pd.leadingBaseCoefficient();
704
705 // factor leading coefficient as polynomial in the lowest! variable
706 GenPolynomialRing<GenPolynomial<BigInteger>> rnfac = pfac.recursive(pfac.nvar - 1);
707 GenPolynomial<GenPolynomial<BigInteger>> pr = PolyUtil.<BigInteger> recursive(rnfac, pd);
708 GenPolynomial<GenPolynomial<BigInteger>> prr = PolyUtil.<BigInteger> switchVariables(pr);
709
710 GenPolynomial<BigInteger> prrc = engine.recursiveContent(prr); // can have content wrt this variable
711 List<GenPolynomial<BigInteger>> cfactors = null;
712 if (!prrc.isONE()) {
713 prr = PolyUtil.<BigInteger> recursiveDivide(prr, prrc);
714 GenPolynomial<BigInteger> prrcu = prrc.extendLower(pfac, 0, 0L); // since switched vars
715 pd = PolyUtil.<BigInteger> basePseudoDivide(pd, prrcu);
716 logger.info("recursive content = " + prrc + ", new P = " + pd);
717 cfactors = factorsSquarefree(prrc);
718 List<GenPolynomial<BigInteger>> cff = new ArrayList<GenPolynomial<BigInteger>>(cfactors.size());
719 for (GenPolynomial<BigInteger> fs : cfactors) {
720 GenPolynomial<BigInteger> fsp = fs.extendLower(pfac, 0, 0L); // since switched vars
721 cff.add(fsp);
722 }
723 cfactors = cff;
724 logger.info("cfactors = " + cfactors);
725 }
726 GenPolynomial<BigInteger> lprr = prr.leadingBaseCoefficient();
727 //System.out.println("prr = " + prr);
728 logger.info("leading coeffcient = " + lprr);
729 boolean isMonic = false; // multivariate monic
730 if (lprr.isConstant()) { // isONE ?
731 isMonic = true;
732 }
733 SortedMap<GenPolynomial<BigInteger>, Long> lfactors = factors(lprr);
734 //System.out.println("lfactors = " + lfactors);
735 List<GenPolynomial<BigInteger>> lfacs = new ArrayList<GenPolynomial<BigInteger>>(lfactors.keySet());
736 logger.info("leading coefficient factors = " + lfacs);
737
738 // search evaluation point and evaluate
739 GenPolynomialRing<BigInteger> cpfac = pfac;
740 GenPolynomial<BigInteger> pe = pd;
741 GenPolynomial<BigInteger> pep;
742 GenPolynomialRing<BigInteger> ccpfac = lprr.ring;
743 List<GenPolynomial<BigInteger>> ce = lfacs;
744 List<GenPolynomial<BigInteger>> cep = null;
745 List<BigInteger> cei = null;
746 List<BigInteger> dei = new ArrayList<BigInteger>();
747 BigInteger pec = null;
748 BigInteger pecw = null;
749 BigInteger ped = null;
750
751 List<GenPolynomial<BigInteger>> ufactors = null;
752 List<TrialParts> tParts = new ArrayList<TrialParts>();
753 List<GenPolynomial<BigInteger>> lf = null;
754 GenPolynomial<BigInteger> lpx = null;
755 List<GenPolynomial<BigInteger>> ln = null;
756 List<GenPolynomial<BigInteger>> un = null;
757 GenPolynomial<BigInteger> pes = null;
758
759 List<BigInteger> V = null;
760 long evStart = 0L; //3L * 5L;
761 List<Long> Evs = new ArrayList<Long>(pfac.nvar + 1); // Evs(0), Evs(1) unused
762 for (int j = 0; j <= pfac.nvar; j++) {
763 Evs.add(evStart);
764 }
765 final int trials = 4;
766 int countSeparate = 0;
767 final int COUNT_MAX = 50;
768 double ran = 1.001; // higher values not good
769 boolean isPrimitive = true;
770 boolean notLucky = true;
771 while (notLucky) { // for Wang's test
772 if (Math.abs(evStart) > 400L) {
773 System.out.println("P = " + P + ", lprr = " + lprr + ", lfacs = " + lfacs);
774 throw new RuntimeException("no lucky evaluation point found after " + evStart + " iterations");
775 }
776 if (Math.abs(evStart) % 100L <= 3L) {
777 ran = ran * (Math.PI - 2.14);
778 }
779 //System.out.println("-------------------------------------------- Evs = " + Evs);
780 notLucky = false;
781 V = new ArrayList<BigInteger>();
782 cpfac = pfac;
783 pe = pd;
784 ccpfac = lprr.ring;
785 ce = lfacs;
786 cep = null;
787 cei = null;
788 pec = null;
789 ped = null;
790 long vi = 0L;
791 for (int j = pfac.nvar; j > 1; j--) {
792 // evaluation up to univariate case
793 long degp = pe.degree(cpfac.nvar - 2);
794 cpfac = cpfac.contract(1);
795 ccpfac = ccpfac.contract(1);
796 //vi = evStart; // + j;//0L; //(long)(pfac.nvar-j); // 1L; 0 not so good for small p
797 vi = Evs.get(j); //evStart + j;//0L; //(long)(pfac.nvar-j); // 1L; 0 not so good for small p
798 BigInteger Vi;
799
800 // search evaluation point
801 boolean doIt = true;
802 Vi = null;
803 pep = null;
804 while (doIt) {
805 logger.info("vi(" + j + ") = " + vi);
806 Vi = new BigInteger(vi);
807 pep = PolyUtil.<BigInteger> evaluateMain(cpfac, pe, Vi);
808 //System.out.println("pep = " + pep);
809 // check lucky evaluation point
810 if (degp == pep.degree(cpfac.nvar - 1)) {
811 logger.info("pep = " + pep);
812 //System.out.println("deg(pe) = " + degp + ", deg(pep) = " + pep.degree(cpfac.nvar-1));
813 // check squarefree
814 if (sengine.isSquarefree(pep)) { // cpfac.nvar == 1 && ?? no, must test on each variable
815 //if ( isNearlySquarefree(pep) ) {
816 //System.out.println("squarefeee(pep)"); // + pep);
817 doIt = false; //break;
818 }
819 }
820 if (vi > 0L) {
821 vi = -vi;
822 } else {
823 vi = 1L - vi;
824 }
825 }
826 //if ( !isMonic ) {
827 if (ccpfac.nvar >= 1) {
828 cep = PolyUtil.<BigInteger> evaluateMain(ccpfac, ce, Vi);
829 } else {
830 cei = PolyUtil.<BigInteger> evaluateMain(ccpfac.coFac, ce, Vi);
831 }
832 //}
833 int jj = (int) Math.round(ran + 0.52 * Math.random()); // j, random increment
834 //jj = 1; // ...4 test
835 //System.out.println("minimal jj = " + jj + ", vi " + vi);
836 if (vi > 0L) {
837 Evs.set(j, vi + jj); // record last tested value plus increment
838 evStart = vi + jj;
839 } else {
840 Evs.set(j, vi - jj); // record last tested value minus increment
841 evStart = vi - jj;
842 }
843 //evStart = vi+1L;
844 V.add(Vi);
845 pe = pep;
846 ce = cep;
847 }
848 //System.out.println("ce = " + ce + ", pe = " + pe);
849 pecw = engine.baseContent(pe); // original Wang
850 isPrimitive = pecw.isONE();
851 ped = ccpfac.coFac.getONE();
852 pec = pe.ring.coFac.getONE();
853 //System.out.println("cei = " + cei + ", pecw = " + pecw);
854 if (!isMonic) {
855 if (countSeparate > COUNT_MAX) {
856 pec = pe.ring.coFac.getONE(); // hack is sometimes better
857 } else {
858 pec = pecw;
859 }
860 //pec = pecw;
861 //System.out.println("cei = " + cei + ", pec = " + pec + ", pe = " + pe);
862 if (lfacs.get(0).isConstant()) {
863 ped = cei.remove(0);
864 //lfacs.remove(0); // later
865 }
866 //System.out.println("lfacs = " + lfacs + ", cei = " + cei + ", ped = " + ped + ", pecw = " + pecw);
867 // test Wang's condition
868 dei = new ArrayList<BigInteger>();
869 dei.add(pec.multiply(ped).abs()); // .abs()
870 int i = 1;
871 for (BigInteger ci : cei) {
872 if (ci.isZERO()) {
873 logger.info("condition (0) not met for cei = " + cei); // + ", dei = " + dei);
874 notLucky = true;
875 break;
876 }
877 BigInteger q = ci.abs();
878 //System.out.println("q = " + q);
879 for (int ii = i - 1; ii >= 0; ii--) {
880 BigInteger r = dei.get(ii);
881 //System.out.println("r = " + r);
882 while (!r.isONE()) {
883 r = r.gcd(q);
884 q = q.divide(r);
885 //System.out.println("r = " + r + ", q = " + q);
886 }
887 }
888 dei.add(q);
889 if (q.isONE()) {
890 logger.info("condition (1) not met for dei = " + dei + ", cei = " + cei);
891 if (!testSeparate(cei, pecw)) {
892 countSeparate++;
893 if (countSeparate > COUNT_MAX) {
894 logger.info("too many inseparable evaluation points: " + countSeparate
895 + ", removing " + pecw);
896 }
897 }
898 notLucky = true;
899 break;
900 }
901 i++;
902 }
903 //System.out.println("dei = " + dei);
904 }
905 if (notLucky) {
906 continue;
907 }
908 logger.info("evaluation points = " + V + ", dei = " + dei);
909 //System.out.println("Evs = " + Evs);
910 logger.info("univariate polynomial = " + pe + ", pecw = " + pecw);
911 //pe = pe.abs();
912 //ufactors = baseFactorsRadical(pe); //baseFactorsSquarefree(pe); wrong since not primitive
913 ufactors = baseFactorsSquarefree(pe.divide(pecw)); //wrong if not primitive
914 if (!pecw.isONE()) {
915 ufactors.add(0, cpfac.getONE().multiply(pecw));
916 }
917 if (ufactors.size() <= 1) {
918 logger.info("irreducible univariate polynomial");
919 factors.add(pd); // P
920 if (cfactors != null) {
921 cfactors.addAll(factors);
922 factors = cfactors;
923 }
924 return factors;
925 }
926 logger.info("univariate factors = " + ufactors); // + ", of " + pe);
927 //System.out.println("lfacs = " + lfacs);
928 //System.out.println("cei = " + cei);
929 //System.out.println("pecw = " + pecw);
930
931 // determine leading coefficient polynomials for factors
932 lf = new ArrayList<GenPolynomial<BigInteger>>();
933 lpx = lprr.ring.getONE();
934 for (GenPolynomial<BigInteger> pp : ufactors) {
935 lf.add(lprr.ring.getONE());
936 }
937 //System.out.println("lf = " + lf);
938 if (!isMonic || !pecw.isONE()) {
939 if (lfacs.size() > 0 && lfacs.get(0).isConstant()) {
940 GenPolynomial<BigInteger> xx = lfacs.remove(0);
941 //BigInteger xxi = xx.leadingBaseCoefficient();
942 //System.out.println("xx = " + xx + " == ped = " +ped);
943 }
944 for (int i = ufactors.size() - 1; i >= 0; i--) {
945 GenPolynomial<BigInteger> pp = ufactors.get(i);
946 BigInteger ppl = pp.leadingBaseCoefficient();
947 //System.out.println("ppl = " + ppl + ", pp = " + pp);
948 ppl = ppl.multiply(pec); // content
949 GenPolynomial<BigInteger> lfp = lf.get(i);
950 int ii = 0;
951 for (BigInteger ci : cei) {
952 //System.out.println("ci = " + ci + ", lfp = " + lfp + ", lfacs.get(ii) = " + lfacs.get(ii));
953 if (ci.abs().isONE()) {
954 System.out.println("ppl = " + ppl + ", ci = " + ci + ", lfp = " + lfp
955 + ", lfacs.get(ii) = " + lfacs.get(ii));
956 throw new RuntimeException("something is wrong, ci is a unit");
957 //notLucky = true;
958 }
959 while (ppl.remainder(ci).isZERO() && lfacs.size() > ii) {
960 ppl = ppl.divide(ci);
961 lfp = lfp.multiply(lfacs.get(ii));
962 }
963 ii++;
964 }
965 //System.out.println("ppl = " + ppl + ", lfp = " + lfp);
966 lfp = lfp.multiply(ppl);
967 lf.set(i, lfp);
968 }
969 // adjust if pec != 1
970 pec = pecw;
971 lpx = Power.<GenPolynomial<BigInteger>> multiply(lprr.ring, lf); // test only, not used
972 //System.out.println("lpx = " + lpx);
973 if (!lprr.degreeVector().equals(lpx.degreeVector())) {
974 logger.info("deg(lprr) != deg(lpx): lprr = " + lprr + ", lpx = " + lpx);
975 notLucky = true;
976 continue;
977 }
978 if (!pec.isONE()) { // content, was always false by hack
979 // evaluate factors of ldcf
980 List<GenPolynomial<BigInteger>> lfe = lf;
981 List<BigInteger> lfei = null;
982 ccpfac = lprr.ring;
983 for (int j = lprr.ring.nvar; j > 0; j--) {
984 ccpfac = ccpfac.contract(1);
985 BigInteger Vi = V.get(lprr.ring.nvar - j);
986 if (ccpfac.nvar >= 1) {
987 lfe = PolyUtil.<BigInteger> evaluateMain(ccpfac, lfe, Vi);
988 } else {
989 lfei = PolyUtil.<BigInteger> evaluateMain(ccpfac.coFac, lfe, Vi);
990 }
991 }
992 //System.out.println("lfe = " + lfe + ", lfei = " + lfei + ", V = " + V);
993
994 ln = new ArrayList<GenPolynomial<BigInteger>>(lf.size());
995 un = new ArrayList<GenPolynomial<BigInteger>>(lf.size());
996 for (int jj = 0; jj < lf.size(); jj++) {
997 GenPolynomial<BigInteger> up = ufactors.get(jj);
998 BigInteger ui = up.leadingBaseCoefficient();
999 BigInteger li = lfei.get(jj);
1000 BigInteger di = ui.gcd(li).abs();
1001 BigInteger udi = ui.divide(di);
1002 BigInteger ldi = li.divide(di);
1003 GenPolynomial<BigInteger> lp = lf.get(jj);
1004 GenPolynomial<BigInteger> lpd = lp.multiply(udi);
1005 GenPolynomial<BigInteger> upd = up.multiply(ldi);
1006 if (pec.isONE()) {
1007 ln.add(lp);
1008 un.add(up);
1009 } else {
1010 ln.add(lpd);
1011 un.add(upd);
1012 BigInteger pec1 = pec.divide(ldi);
1013 //System.out.println("pec = " + pec + ", pec1 = " + pec1);
1014 pec = pec1;
1015 }
1016 }
1017 if (!lf.equals(ln) || !un.equals(ufactors)) {
1018 //System.out.println("pe = " + pe);
1019 //System.out.println("#ln = " + ln + ", #lf = " + lf);
1020 //System.out.println("#un = " + un + ", #ufactors = " + ufactors);
1021 //lf = ln;
1022 //ufactors = un;
1023 // adjust pe
1024 }
1025 if (!pec.isONE()) { // still not 1
1026 ln = new ArrayList<GenPolynomial<BigInteger>>(lf.size());
1027 un = new ArrayList<GenPolynomial<BigInteger>>(lf.size());
1028 pes = pe;
1029 for (int jj = 0; jj < lf.size(); jj++) {
1030 GenPolynomial<BigInteger> up = ufactors.get(jj);
1031 GenPolynomial<BigInteger> lp = lf.get(jj);
1032 //System.out.println("up = " + up + ", lp = " + lp);
1033 if (!up.isConstant()) {
1034 up = up.multiply(pec);
1035 }
1036 lp = lp.multiply(pec);
1037 if (jj != 0) {
1038 pes = pes.multiply(pec);
1039 }
1040 un.add(up);
1041 ln.add(lp);
1042 }
1043 if (pes.equals(Power.<GenPolynomial<BigInteger>> multiply(pe.ring, un))) {
1044 //System.out.println("*pe = " + pes + ", pec = " + pec);
1045 //ystem.out.println("*ln = " + ln + ", *lf = " + lf);
1046 //System.out.println("*un = " + un + ", *ufactors = " + ufactors);
1047 //System.out.println("*pe == prod(un) ");
1048 isPrimitive = false;
1049 //pe = pes;
1050 //lf = ln;
1051 //ufactors = un;
1052 } else {
1053 //System.out.println("*pe != prod(un): " + Power.<GenPolynomial<BigInteger>> multiply(pe.ring,un));
1054 }
1055 }
1056 }
1057 if (notLucky) {
1058 continue;
1059 }
1060 logger.info("distributed factors of leading coefficient = " + lf);
1061 lpx = Power.<GenPolynomial<BigInteger>> multiply(lprr.ring, lf);
1062 if (!lprr.abs().equals(lpx.abs())) { // not correctly distributed
1063 if (!lprr.degreeVector().equals(lpx.degreeVector())) {
1064 logger.info("lprr != lpx: lprr = " + lprr + ", lpx = " + lpx);
1065 notLucky = true;
1066 }
1067 }
1068 } // end determine leading coefficients for factors
1069
1070 if (!notLucky) {
1071 TrialParts tp = null;
1072 if (isPrimitive) {
1073 tp = new TrialParts(V, pe, ufactors, cei, lf);
1074 } else {
1075 tp = new TrialParts(V, pes, un, cei, ln);
1076 }
1077 //System.out.println("trialParts = " + tp);
1078 tParts.add(tp);
1079 if (tParts.size() < trials) {
1080 notLucky = true;
1081 }
1082 }
1083 } // end notLucky loop
1084
1085 // search TrialParts with shortest factorization of univariate polynomial
1086 int min = Integer.MAX_VALUE;
1087 TrialParts tpmin = null;
1088 for (TrialParts tp : tParts) {
1089 logger.info("tp.univFactors.size() = " + tp.univFactors.size());
1090 if (tp.univFactors.size() < min) {
1091 min = tp.univFactors.size();
1092 tpmin = tp;
1093 }
1094 }
1095 for (TrialParts tp : tParts) {
1096 //logger.info("tp.univFactors.get(0) = " + tp.univFactors.get(0));
1097 if (tp.univFactors.size() == min) {
1098 if (!tp.univFactors.get(0).isConstant()) {
1099 tpmin = tp;
1100 break;
1101 }
1102 }
1103 }
1104 // set to (first) shortest
1105 V = tpmin.evalPoints;
1106 pe = tpmin.univPoly;
1107 ufactors = tpmin.univFactors;
1108 cei = tpmin.ldcfEval; // unused
1109 lf = tpmin.ldcfFactors;
1110 logger.info("minimal trial = " + tpmin);
1111
1112 GenPolynomialRing<BigInteger> ufac = pe.ring;
1113
1114 //initialize prime list
1115 PrimeList primes = new PrimeList(PrimeList.Range.medium); // PrimeList.Range.medium);
1116 Iterator<java.math.BigInteger> primeIter = primes.iterator();
1117 int pn = 50; //primes.size();
1118 BigInteger ae = pe.leadingBaseCoefficient();
1119 GenPolynomial<MOD> Pm = null;
1120 ModularRingFactory<MOD> cofac = null;
1121 GenPolynomialRing<MOD> mufac = null;
1122
1123 // search lucky prime
1124 for (int i = 0; i < 11; i++) { // prime meta loop
1125 //for ( int i = 0; i < 1; i++ ) { // meta loop
1126 java.math.BigInteger p = null; //new java.math.BigInteger("19"); //primes.next();
1127 // 2 small, 5 medium and 4 large size primes
1128 if (i == 0) { // medium size
1129 primes = new PrimeList(PrimeList.Range.medium);
1130 primeIter = primes.iterator();
1131 }
1132 if (i == 5) { // small size
1133 primes = new PrimeList(PrimeList.Range.small);
1134 primeIter = primes.iterator();
1135 p = primeIter.next(); // 2
1136 p = primeIter.next(); // 3
1137 p = primeIter.next(); // 5
1138 p = primeIter.next(); // 7
1139 }
1140 if (i == 7) { // large size
1141 primes = new PrimeList(PrimeList.Range.large);
1142 primeIter = primes.iterator();
1143 }
1144 int pi = 0;
1145 while (pi < pn && primeIter.hasNext()) {
1146 p = primeIter.next();
1147 logger.info("prime = " + p);
1148 // initialize coefficient factory and map normalization factor and polynomials
1149 ModularRingFactory<MOD> cf = null;
1150 if (ModLongRing.MAX_LONG.compareTo(p) > 0) {
1151 cf = (ModularRingFactory) new ModLongRing(p, true);
1152 } else {
1153 cf = (ModularRingFactory) new ModIntegerRing(p, true);
1154 }
1155 MOD nf = cf.fromInteger(ae.getVal());
1156 if (nf.isZERO()) {
1157 continue;
1158 }
1159 mufac = new GenPolynomialRing<MOD>(cf, ufac);
1160 //System.out.println("mufac = " + mufac.toScript());
1161 Pm = PolyUtil.<MOD> fromIntegerCoefficients(mufac, pe);
1162 //System.out.println("Pm = " + Pm);
1163 if (!mfactor.isSquarefree(Pm)) {
1164 continue;
1165 }
1166 cofac = cf;
1167 break;
1168 }
1169 if (cofac != null) {
1170 break;
1171 }
1172 } // end prime meta loop
1173 if (cofac == null) { // no lucky prime found
1174 throw new RuntimeException("giving up on Hensel preparation, no lucky prime found");
1175 }
1176 logger.info("lucky prime = " + cofac.getIntegerModul());
1177 if (logger.isDebugEnabled()) {
1178 logger.debug("univariate modulo p: = " + Pm);
1179 }
1180
1181 // coefficient bound
1182 BigInteger an = pd.maxNorm();
1183 BigInteger mn = an.multiply(ac.abs()).multiply(new BigInteger(2L));
1184 long k = Power.logarithm(cofac.getIntegerModul(), mn) + 1L;
1185 //System.out.println("mn = " + mn + ", k = " +k);
1186
1187 BigInteger q = Power.positivePower(cofac.getIntegerModul(), k);
1188 ModularRingFactory<MOD> muqfac;
1189 if (ModLongRing.MAX_LONG.compareTo(q.getVal()) > 0) {
1190 muqfac = (ModularRingFactory) new ModLongRing(q.getVal());
1191 } else {
1192 muqfac = (ModularRingFactory) new ModIntegerRing(q.getVal());
1193 }
1194 //System.out.println("muqfac = " + muqfac);
1195 GenPolynomialRing<MOD> mucpfac = new GenPolynomialRing<MOD>(muqfac, ufac);
1196
1197 List<GenPolynomial<MOD>> muqfactors = PolyUtil.<MOD> fromIntegerCoefficients(mucpfac, ufactors);
1198 GenPolynomial<MOD> peqq = PolyUtil.<MOD> fromIntegerCoefficients(mucpfac, pe);
1199 if (debug) {
1200 if (!mfactor.isFactorization(peqq, muqfactors)) { // should not happen
1201 System.out.println("muqfactors = " + muqfactors);
1202 System.out.println("peqq = " + peqq);
1203 throw new RuntimeException("something is wrong, no modular p^k factorization");
1204 }
1205 }
1206 logger.info("univariate modulo p^k: " + peqq + " = " + muqfactors);
1207
1208 // convert C from Z[...] to Z_q[...]
1209 GenPolynomialRing<MOD> qcfac = new GenPolynomialRing<MOD>(muqfac, pd.ring);
1210 GenPolynomial<MOD> pq = PolyUtil.<MOD> fromIntegerCoefficients(qcfac, pd);
1211 //System.out.println("pd = " + pd);
1212 logger.info("multivariate modulo p^k: " + pq);
1213
1214 List<MOD> Vm = new ArrayList<MOD>(V.size());
1215 for (BigInteger v : V) {
1216 MOD vm = muqfac.fromInteger(v.getVal());
1217 Vm.add(vm);
1218 }
1219 //System.out.println("Vm = " + Vm);
1220
1221 // Hensel lifting of factors
1222 List<GenPolynomial<MOD>> mlift;
1223 try {
1224 mlift = HenselMultUtil.<MOD> liftHensel(pd, pq, muqfactors, Vm, k, lf);
1225 logger.info("mlift = " + mlift);
1226 } catch (NoLiftingException nle) {
1227 //System.out.println("exception : " + nle);
1228 //nle.printStackTrace();
1229 mlift = new ArrayList<GenPolynomial<MOD>>();
1230 throw new RuntimeException(nle);
1231 } catch (ArithmeticException aex) {
1232 //System.out.println("exception : " + aex);
1233 //aex.printStackTrace();
1234 mlift = new ArrayList<GenPolynomial<MOD>>();
1235 throw aex;
1236 }
1237 if (mlift.size() <= 1) { // irreducible mod I, p^k, can this happen?
1238 logger.info("modular lift size == 1: " + mlift);
1239 factors.add(pd); // P
1240 if (cfactors != null) {
1241 cfactors.addAll(factors);
1242 factors = cfactors;
1243 }
1244 return factors;
1245 }
1246
1247 // combine trial factors
1248 GenPolynomialRing<MOD> mfac = mlift.get(0).ring;
1249 int dl = (mlift.size() + 1) / 2;
1250 GenPolynomial<BigInteger> u = P;
1251 long deg = (u.degree() + 1L) / 2L;
1252
1253 GenPolynomial<BigInteger> ui = pd;
1254 for (int j = 1; j <= dl; j++) {
1255 //System.out.println("j = " + j + ", dl = " + dl + ", mlift = " + mlift);
1256 KsubSet<GenPolynomial<MOD>> subs = new KsubSet<GenPolynomial<MOD>>(mlift, j);
1257 for (List<GenPolynomial<MOD>> flist : subs) {
1258 //System.out.println("degreeSum = " + degreeSum(flist));
1259 GenPolynomial<MOD> mtrial = Power.<GenPolynomial<MOD>> multiply(mfac, flist);
1260 if (mtrial.degree() > deg) { // this test is sometimes wrong
1261 logger.info("degree > deg " + deg + ", degree = " + mtrial.degree());
1262 //continue;
1263 }
1264 GenPolynomial<BigInteger> trial = PolyUtil.integerFromModularCoefficients(pfac, mtrial);
1265 trial = engine.basePrimitivePart(trial);
1266 //if ( ! isPrimitive ) {
1267 //}
1268 if (debug) {
1269 logger.info("trial = " + trial); // + ", mtrial = " + mtrial);
1270 }
1271 if (PolyUtil.<BigInteger> basePseudoRemainder(ui, trial).isZERO()) {
1272 logger.info("successful trial = " + trial);
1273 factors.add(trial);
1274 ui = PolyUtil.<BigInteger> basePseudoDivide(ui, trial);
1275 //System.out.println("ui = " + ui);
1276 mlift = removeOnce(mlift, flist);
1277 logger.info("new mlift= " + mlift);
1278 //System.out.println("dl = " + dl);
1279 if (mlift.size() > 1) {
1280 dl = (mlift.size() + 1) / 2;
1281 j = 0; // since j++
1282 break;
1283 }
1284 logger.info("last factor = " + ui);
1285 factors.add(ui);
1286 if (cfactors != null) {
1287 cfactors.addAll(factors);
1288 factors = cfactors;
1289 }
1290 return factors;
1291 }
1292 }
1293 }
1294 if (!ui.isONE() && !ui.equals(pd)) {
1295 logger.info("rest factor = " + ui);
1296 // pp(ui) ?? no ??
1297 factors.add(ui);
1298 }
1299 if (factors.size() == 0) {
1300 logger.info("irreducible P = " + P);
1301 factors.add(pd); // P
1302 }
1303 if (cfactors != null) {
1304 cfactors.addAll(factors);
1305 factors = cfactors;
1306 }
1307 return factors;
1308 }
1309
1310
1311 /**
1312 * Test if b has a prime factor different to the elements of A.
1313 * @param A list of integer with at least one different prime factor.
1314 * @param b integer to test with A.
1315 * @return true, if b hase a prime factor different to elements of A
1316 */
1317 boolean testSeparate(List<BigInteger> A, BigInteger b) {
1318 int i = 0;
1319 List<BigInteger> gei = new ArrayList<BigInteger>(A.size());
1320 for (BigInteger c : A) {
1321 BigInteger g = c.gcd(b).abs();
1322 gei.add(g);
1323 if (!g.isONE()) {
1324 i++;
1325 }
1326 }
1327 //if ( i >= 1 ) {
1328 //System.out.println("gei = " + gei + ", cei = " + cei + ", pec(w) = " + pec);
1329 //}
1330 return (i <= 1);
1331 }
1332
1333
1334 // not useable
1335 boolean isNearlySquarefree(GenPolynomial<BigInteger> P) { // unused
1336 // in main variable
1337 GenPolynomialRing<BigInteger> pfac = P.ring;
1338 if (true || pfac.nvar <= 4) {
1339 return sengine.isSquarefree(P);
1340 }
1341 GenPolynomialRing<GenPolynomial<BigInteger>> rfac = pfac.recursive(1);
1342 GenPolynomial<GenPolynomial<BigInteger>> Pr = PolyUtil.<BigInteger> recursive(rfac, P);
1343 GenPolynomial<GenPolynomial<BigInteger>> Ps = PolyUtil.<BigInteger> recursiveDeriviative(Pr);
1344 System.out.println("Pr = " + Pr);
1345 System.out.println("Ps = " + Ps);
1346 GenPolynomial<GenPolynomial<BigInteger>> g = engine.recursiveUnivariateGcd(Pr, Ps);
1347 System.out.println("g_m = " + g);
1348 if (!g.isONE()) {
1349 return false;
1350 }
1351 // in lowest variable
1352 rfac = pfac.recursive(pfac.nvar - 1);
1353 Pr = PolyUtil.<BigInteger> recursive(rfac, P);
1354 Pr = PolyUtil.<BigInteger> switchVariables(Pr);
1355 Ps = PolyUtil.<BigInteger> recursiveDeriviative(Pr);
1356 System.out.println("Pr = " + Pr);
1357 System.out.println("Ps = " + Ps);
1358 g = engine.recursiveUnivariateGcd(Pr, Ps);
1359 System.out.println("g_1 = " + g);
1360 if (!g.isONE()) {
1361 return false;
1362 }
1363 return true;
1364 }
1365
1366 }
1367
1368
1369 /**
1370 * Container for factorization trial lifting parameters.
1371 */
1372 class TrialParts {
1373
1374
1375 /**
1376 * evaluation points
1377 */
1378 public final List<BigInteger> evalPoints;
1379
1380
1381 /**
1382 * univariate polynomial
1383 */
1384 public final GenPolynomial<BigInteger> univPoly;
1385
1386
1387 /**
1388 * irreducible factors of univariate polynomial
1389 */
1390 public final List<GenPolynomial<BigInteger>> univFactors;
1391
1392
1393 /**
1394 * irreducible factors of leading coefficient
1395 */
1396 public final List<GenPolynomial<BigInteger>> ldcfFactors;
1397
1398
1399 /**
1400 * evaluated factors of leading coefficient factors by evaluation points
1401 */
1402 public final List<BigInteger> ldcfEval;
1403
1404
1405 /**
1406 * Constructor.
1407 * @param ev evaluation points.
1408 * @param up univariate polynomial.
1409 * @param uf irreducible factors of up.
1410 * @param le irreducible factors of leading coefficient.
1411 * @param lf evaluated le by evaluation points.
1412 */
1413 public TrialParts(List<BigInteger> ev, GenPolynomial<BigInteger> up, List<GenPolynomial<BigInteger>> uf,
1414 List<BigInteger> le, List<GenPolynomial<BigInteger>> lf) {
1415 evalPoints = ev;
1416 univPoly = up;
1417 univFactors = uf;
1418 //ldcfPoly = lp;
1419 ldcfFactors = lf;
1420 ldcfEval = le;
1421 }
1422
1423
1424 /**
1425 * @see java.lang.Object#toString()
1426 */
1427 @Override
1428 public String toString() {
1429 StringBuffer sb = new StringBuffer();
1430 sb.append("TrialParts[");
1431 sb.append("evalPoints = " + evalPoints);
1432 sb.append(", univPoly = " + univPoly);
1433 sb.append(", univFactors = " + univFactors);
1434 sb.append(", ldcfEval = " + ldcfEval);
1435 sb.append(", ldcfFactors = " + ldcfFactors);
1436 sb.append("]");
1437 return sb.toString();
1438 }
1439
1440 }