Source code

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