Source code

001/*
002 * $Id$
003 */
004
005package edu.jas.ufd;
006
007
008import java.util.ArrayList;
009import java.util.Iterator;
010import java.util.List;
011import java.util.SortedSet;
012import java.util.TreeSet;
013
014import org.apache.logging.log4j.LogManager;
015import org.apache.logging.log4j.Logger;
016
017import edu.jas.arith.ModLongRing;
018import edu.jas.arith.ModularRingFactory;
019import edu.jas.poly.AlgebraicNumberRing;
020import edu.jas.poly.GenPolynomial;
021import edu.jas.poly.GenPolynomialRing;
022import edu.jas.poly.PolyUtil;
023import edu.jas.structure.GcdRingElem;
024import edu.jas.structure.Power;
025import edu.jas.structure.RingFactory;
026import edu.jas.vector.GenMatrix;
027import edu.jas.vector.GenMatrixRing;
028import edu.jas.vector.GenVector;
029import edu.jas.vector.GenVectorModul;
030import edu.jas.vector.LinAlg;
031
032
033/**
034 * Modular coefficients Berlekamp factorization algorithms. This class
035 * implements Berlekamp, Cantor and Zassenhaus factorization methods for
036 * polynomials over (prime) modular integers.
037 * @author Heinz Kredel
038 */
039
040public class FactorModularBerlekamp<MOD extends GcdRingElem<MOD>> extends FactorAbsolute<MOD> {
041
042
043    private static final Logger logger = LogManager.getLogger(FactorModularBerlekamp.class);
044
045
046    //private static final boolean debug = logger.isDebugEnabled();
047
048
049    /**
050     * No argument constructor, do not use.
051     */
052    @SuppressWarnings({ "unchecked", "unused" })
053    private FactorModularBerlekamp() {
054        this((RingFactory<MOD>) (Object) new ModLongRing(13, true)); // hack, 13 unimportant
055    }
056
057
058    /**
059     * Constructor.
060     * @param cfac coefficient ring factory.
061     */
062    public FactorModularBerlekamp(RingFactory<MOD> cfac) {
063        super(cfac);
064    }
065
066
067    /**
068     * GenPolynomial base factorization of a squarefree polynomial.
069     * @param P squarefree and monic! GenPolynomial.
070     * @return [p_1,...,p_k] with P = prod_{i=1,...,r} p_i.
071     */
072    @Override
073    public List<GenPolynomial<MOD>> baseFactorsSquarefree(GenPolynomial<MOD> P) {
074        if (P == null) {
075            throw new IllegalArgumentException("P == null not allowed");
076        }
077        //ModularRingFactory cfac = (ModularRingFactory) P.ring.coFac;
078        long q = P.ring.coFac.characteristic().longValueExact();
079        if (q < 100) { // todo, 32003 too high
080            return baseFactorsSquarefreeSmallPrime(P);
081        }
082        return baseFactorsSquarefreeBigPrime(P);
083    }
084
085
086    /**
087     * GenPolynomial base factorization of a squarefree polynomial. Small prime
088     * version of Berlekamps algorithm.
089     * @param P squarefree and monic! GenPolynomial.
090     * @return [p_1,...,p_k] with P = prod_{i=1,...,r} p_i.
091     */
092    @SuppressWarnings("unchecked")
093    public List<GenPolynomial<MOD>> baseFactorsSquarefreeSmallPrime(GenPolynomial<MOD> P) {
094        if (P == null) {
095            throw new IllegalArgumentException("P == null");
096        }
097        List<GenPolynomial<MOD>> factors = new ArrayList<GenPolynomial<MOD>>();
098        if (P.isZERO()) {
099            return factors;
100        }
101        GenPolynomialRing<MOD> pfac = P.ring;
102        if (pfac.nvar > 1) {
103            throw new IllegalArgumentException("only for univariate polynomials");
104        }
105        if (!P.leadingBaseCoefficient().isONE()) {
106            throw new IllegalArgumentException("ldcf(P) != 1: " + P);
107        }
108        if (P.isONE() || P.degree(0) <= 1) {
109            factors.add(P);
110            return factors;
111        }
112        ArrayList<ArrayList<MOD>> Q = PolyUfdUtil.<MOD> constructQmatrix(P);
113        //System.out.println("Q = " + Q);
114        int n = Q.size();
115        int m = Q.get(0).size();
116        GenMatrixRing<MOD> mfac = new GenMatrixRing<MOD>(pfac.coFac, n, m);
117        GenMatrix<MOD> Qm = new GenMatrix<MOD>(mfac, Q);
118        GenMatrix<MOD> Qm1 = Qm.subtract(mfac.getONE());
119        //System.out.println("Qm1 = " + Qm1);
120        LinAlg<MOD> lu = new LinAlg<MOD>();
121        List<GenVector<MOD>> Nsb = lu.nullSpaceBasis(Qm1);
122        logger.info("Nsb = {}", Nsb);
123        int k = Nsb.size();
124        if (k == 1) {
125            factors.add(P);
126            return factors;
127        }
128        //int d = (int) P.degree(0);
129        GenMatrix<MOD> Ns = mfac.fromVectors(Nsb);
130        GenMatrix<MOD> L1 = Ns.negate(); //mfac.getONE().subtract(Ns);
131        //System.out.println("L1 = " + L1);
132        List<GenPolynomial<MOD>> trials = new ArrayList<GenPolynomial<MOD>>();
133        for (int i = 0; i < L1.ring.rows; i++) {
134            GenVector<MOD> rv = L1.getRow(i);
135            GenPolynomial<MOD> rp = pfac.fromVector(rv);
136            if (!rp.isONE()) {
137                trials.add(rp);
138            }
139        }
140        logger.info("#ofFactors k = {}", k);
141        logger.info("trials = {}", trials);
142        factors.add(P);
143        for (GenPolynomial<MOD> t : trials) {
144            if (factors.size() == k) {
145                break;
146            }
147            //System.out.println("t = " + t);
148            // depth first search, since Iterator is finite
149            GenPolynomial<MOD> a = factors.remove(0);
150            //System.out.println("a = " + a);
151            MOD s = null;
152            Iterator<MOD> eit = null;
153            if (pfac.coFac instanceof ModularRingFactory) {
154                eit = ((ModularRingFactory) pfac.coFac).iterator();
155            } else if (pfac.coFac instanceof AlgebraicNumberRing) {
156                eit = ((AlgebraicNumberRing) pfac.coFac).iterator();
157            } else {
158                throw new IllegalArgumentException("no iterator for: " + pfac.coFac);
159            }
160            //System.out.println("eit = " + eit);
161            while (eit.hasNext()) {
162                s = eit.next();
163                //System.out.println("s = " + s);
164                GenPolynomial<MOD> v = t.subtract(s);
165                GenPolynomial<MOD> g = v.gcd(a);
166                if (g.isONE() || g.equals(a)) {
167                    continue;
168                }
169                factors.add(g);
170                a = a.divide(g);
171                logger.info("s = {}, g = {}, a = {}", s, g, a);
172                if (a.isONE()) {
173                    break;
174                }
175            }
176            if (!a.isONE()) {
177                factors.add(a);
178            }
179        }
180        //System.out.println("factors  = " + factors);
181        factors = PolyUtil.<MOD> monic(factors);
182        SortedSet<GenPolynomial<MOD>> ss = new TreeSet<GenPolynomial<MOD>>(factors);
183        //System.out.println("sorted   = " + ss);
184        factors.clear();
185        factors.addAll(ss);
186        return factors;
187    }
188
189
190    /**
191     * GenPolynomial base factorization of a squarefree polynomial. Big prime
192     * version of Berlekamps algorithm.
193     * @param P squarefree and monic! GenPolynomial.
194     * @return [p_1,...,p_k] with P = prod_{i=1,...,r} p_i.
195     */
196    public List<GenPolynomial<MOD>> baseFactorsSquarefreeBigPrime(GenPolynomial<MOD> P) {
197        if (P == null) {
198            throw new IllegalArgumentException("P == null");
199        }
200        List<GenPolynomial<MOD>> factors = new ArrayList<GenPolynomial<MOD>>();
201        if (P.isZERO()) {
202            return factors;
203        }
204        GenPolynomialRing<MOD> pfac = P.ring;
205        if (pfac.nvar > 1) {
206            throw new IllegalArgumentException("only for univariate polynomials");
207        }
208        if (!P.leadingBaseCoefficient().isONE()) {
209            throw new IllegalArgumentException("ldcf(P) != 1: " + P);
210        }
211        if (P.isONE() || P.degree(0) <= 1) {
212            factors.add(P);
213            return factors;
214        }
215        ArrayList<ArrayList<MOD>> Q = PolyUfdUtil.<MOD> constructQmatrix(P);
216        //System.out.println("Q = " + Q);
217        int n = Q.size();
218        int m = Q.get(0).size();
219        GenMatrixRing<MOD> mfac = new GenMatrixRing<MOD>(pfac.coFac, n, m);
220        GenMatrix<MOD> Qm = new GenMatrix<MOD>(mfac, Q);
221        GenMatrix<MOD> Qm1 = Qm.subtract(mfac.getONE());
222        //System.out.println("Qm1 = " + Qm1);
223        LinAlg<MOD> lu = new LinAlg<MOD>();
224        List<GenVector<MOD>> Nsb = lu.nullSpaceBasis(Qm1);
225        logger.info("Nsb = {}", Nsb);
226        int k = Nsb.size();
227        if (k == 1) {
228            factors.add(P);
229            return factors;
230        }
231        //int d = (int) P.degree(0);
232        GenMatrix<MOD> Ns = mfac.fromVectors(Nsb);
233        GenMatrix<MOD> L1 = Ns.negate(); //mfac.getONE().subtract(Ns);
234        //System.out.println("L1 = " + L1);
235        List<GenPolynomial<MOD>> trials = new ArrayList<GenPolynomial<MOD>>();
236        for (int i = 0; i < L1.ring.rows; i++) {
237            GenVector<MOD> rv = L1.getRow(i);
238            GenPolynomial<MOD> rp = pfac.fromVector(rv);
239            trials.add(rp);
240        }
241        logger.info("#ofFactors k = {}", k);
242        logger.info("trials = {}", trials);
243        factors.add(P);
244        GenVectorModul<MOD> vfac = new GenVectorModul<MOD>(pfac.coFac, k);
245        //long q = pfac.coFac.characteristic().longValueExact();
246        //long lq = Power.logarithm(2, q);
247        java.math.BigInteger q = pfac.coFac.characteristic(); //.longValueExact();
248        int lq = q.bitLength(); //Power.logarithm(2, q);
249        if (pfac.coFac instanceof AlgebraicNumberRing) {
250            lq = (int) ((AlgebraicNumberRing) pfac.coFac).extensionDegree();
251            q = q.pow(lq); //Power.power(q, lq);
252        }
253        logger.info("char = {}, q = {}, lq = {}", pfac.coFac.characteristic(), q, lq);
254        do {
255            // breadth first search, since some a might be irreducible
256            GenPolynomial<MOD> a = factors.remove(0);
257            //System.out.println("a = " + a);
258            if (a.degree(0) <= 1) {
259                factors.add(a);
260                continue;
261            }
262            GenVector<MOD> rv = vfac.random(8, 0.95f);
263            //System.out.println("rv = " + rv.toScript());
264            GenPolynomial<MOD> rpol = pfac.getZERO();
265            int i = 0;
266            for (GenPolynomial<MOD> t : trials) {
267                MOD c = rv.get(i++);
268                GenPolynomial<MOD> s = t.multiply(c);
269                rpol = rpol.sum(s);
270            }
271            rpol = rpol.monic();
272            if (rpol.isONE()) {
273                factors.add(a);
274                continue;
275            }
276            //System.out.println("rpol = " + rpol.toScript());
277            if (!q.testBit(0)) { // q % 2 == 0
278                long e = lq - 1;
279                //System.out.println("q = " + q + ", e = " + e);
280                GenPolynomial<MOD> pow = rpol;
281                GenPolynomial<MOD> v = rpol;
282                for (int l = 1; l < e; l++) {
283                    pow = pow.multiply(pow).remainder(a);
284                    v = v.sum(pow);
285                }
286                rpol = v.remainder(a).monic(); // automatic monic
287                //System.out.println("sum_l rpol^l = " + rpol.toScript());
288            } else {
289                //long e = (q - 1) / 2;
290                java.math.BigInteger e = q.subtract(java.math.BigInteger.ONE).shiftRight(1);
291                //System.out.println("q = " + q + ", e = " + e);
292                GenPolynomial<MOD> pow = Power.<GenPolynomial<MOD>> modPositivePower(rpol, e, a);
293                rpol = pow.subtract(pfac.getONE()).monic();
294                //System.out.println("rpol^e-1 = " + rpol.toScript());
295            }
296            if (rpol.isZERO() || rpol.isONE()) {
297                factors.add(a);
298                continue;
299            }
300            GenPolynomial<MOD> g = rpol.gcd(a);
301            if (g.isONE()) {
302                factors.add(a);
303                continue;
304            }
305            factors.add(g);
306            //System.out.println("a = " + a);
307            a = a.divide(g);
308            logger.info("rv = {}, g = {}, a/g = {}", rv, g, a);
309            if (!a.isONE()) {
310                factors.add(a);
311            }
312            //System.out.println("factors  = " + factors);
313        } while (factors.size() < k);
314
315        //System.out.println("factors  = " + factors);
316        factors = PolyUtil.<MOD> monic(factors);
317        SortedSet<GenPolynomial<MOD>> ss = new TreeSet<GenPolynomial<MOD>>(factors);
318        //System.out.println("sorted   = " + ss);
319        factors.clear();
320        factors.addAll(ss);
321        return factors;
322    }
323
324}