Source code

001/*
002 * $Id: GroebnerBaseWalk.java 5870 2018-07-20 15:56:23Z kredel $
003 */
004
005package edu.jas.gbufd;
006
007
008import java.util.ArrayList;
009import java.util.List;
010import java.util.Set;
011import java.util.SortedSet;
012
013import org.apache.logging.log4j.Logger;
014import org.apache.logging.log4j.LogManager; 
015
016import edu.jas.gb.GroebnerBaseAbstract;
017import edu.jas.gb.ReductionAbstract;
018import edu.jas.poly.ExpVector;
019import edu.jas.poly.GenPolynomial;
020import edu.jas.poly.GenPolynomialRing;
021import edu.jas.poly.Monomial;
022import edu.jas.poly.PolyUtil;
023import edu.jas.poly.PolynomialList;
024import edu.jas.poly.TermOrder;
025import edu.jas.poly.TermOrderByName;
026import edu.jas.structure.GcdRingElem;
027import edu.jas.structure.RingFactory;
028
029
030/**
031 * Groebner Base sequential Groebner Walk algorithm. Implements Groebner base
032 * computation via Groebner Walk algorithm. See "The generic Groebner walk" by
033 * Fukuda, Jensen, Lauritzen, Thomas, 2005.
034 *
035 * The start term order t1 can be set by a constructor. The target term order t2
036 * is taken from the input polynomials term order.
037 *
038 * @param <C> coefficient type
039 * @author Heinz Kredel
040 * 
041 * @see edu.jas.application.GBAlgorithmBuilder
042 * @see edu.jas.gbufd.GBFactory
043 */
044public class GroebnerBaseWalk<C extends GcdRingElem<C>> extends GroebnerBaseAbstract<C> {
045
046
047    private static final Logger logger = LogManager.getLogger(GroebnerBaseWalk.class);
048
049
050    private static final boolean debug = logger.isDebugEnabled();
051
052
053    /**
054     * The backing GB algorithm implementation.
055     */
056    protected GroebnerBaseAbstract<C> sgb;
057
058
059    /**
060     * The start term order t1.
061     */
062    //protected TermOrder startTO = TermOrderByName.IGRLEX.blockOrder(2); 
063    protected TermOrder startTO = TermOrderByName.IGRLEX;
064
065
066    /**
067     * Print intermediate GB after this number of iterations.
068     */
069    int iterPrint = 100;
070
071
072    /**
073     * Constructor.
074     */
075    public GroebnerBaseWalk() {
076        super();
077        sgb = null;
078    }
079
080
081    /**
082     * Constructor.
083     * @param coFac coefficient ring of polynomial ring.
084     */
085    public GroebnerBaseWalk(RingFactory<C> coFac) {
086        this(GBFactory.<C> getImplementation(coFac));
087    }
088
089
090    /**
091     * Constructor.
092     * @param coFac coefficient ring of polynomial ring.
093     * @param t1 start term order.
094     */
095    public GroebnerBaseWalk(RingFactory<C> coFac, TermOrder t1) {
096        this(GBFactory.<C> getImplementation(coFac), t1);
097    }
098
099
100    /**
101     * Constructor.
102     * @param gb backing GB algorithm.
103     */
104    public GroebnerBaseWalk(GroebnerBaseAbstract<C> gb) {
105        super(gb.red, gb.strategy);
106        sgb = gb;
107    }
108
109
110    /**
111     * Constructor.
112     * @param gb backing GB algorithm.
113     * @param t1 start term order.
114     */
115    public GroebnerBaseWalk(GroebnerBaseAbstract<C> gb, TermOrder t1) {
116        this(gb);
117        startTO = t1;
118    }
119
120
121    /**
122     * Get the String representation with GB engine.
123     * @see java.lang.Object#toString()
124     */
125    @Override
126    public String toString() {
127        if (sgb == null) {
128            return "GroebnerBaseWalk(" + startTO.toScript() + ")";
129        }
130        return "GroebnerBaseWalk( " + sgb.toString() + ", " + startTO.toScript() + " )";
131    }
132
133
134    /**
135     * Groebner base using Groebner Walk algorithm.
136     * @param modv module variable number.
137     * @param F polynomial list in target term order.
138     * @return GB(F) a INVLEX / target term order Groebner base of F.
139     */
140    public List<GenPolynomial<C>> GB(int modv, List<GenPolynomial<C>> F) {
141        List<GenPolynomial<C>> G = normalizeZerosOnes(F);
142        G = PolyUtil.<C> monic(G);
143        if (G == null || G.size() == 0) {
144            return G;
145        }
146        GenPolynomialRing<C> pfac = G.get(0).ring;
147        if (!pfac.coFac.isField()) {
148            throw new IllegalArgumentException("coefficients not from a field: " + pfac.coFac);
149        }
150        if (G.size() <= 1) {
151            GenPolynomial<C> p = pfac.copy(G.get(0)); // change term order
152            G.clear();
153            G.add(p);
154            return G;
155        }
156        // compute in graded / start term order
157        TermOrder grord = startTO; //new TermOrder(TermOrder.IGRLEX);
158        GenPolynomialRing<C> gfac = new GenPolynomialRing<C>(pfac, grord);
159        if (debug) {
160            logger.info("gfac = " + gfac.toScript());
161        }
162        List<GenPolynomial<C>> Fp = gfac.copy(F);
163        logger.info("Term order: graded = " + grord + ", target = " + pfac.tord);
164
165        // compute graded / start term order Groebner base
166        if (sgb == null) {
167            sgb = GBFactory.<C> getImplementation(pfac.coFac, strategy);
168        }
169        List<GenPolynomial<C>> Gp = sgb.GB(modv, Fp);
170        logger.info("graded / start GB = " + Gp);
171        if (grord.equals(pfac.tord)) {
172            return Gp;
173        }
174        if (Gp.size() == 1) { // also dimension -1
175            GenPolynomial<C> p = pfac.copy(Gp.get(0)); // change term order
176            G.clear();
177            G.add(p);
178            return G;
179        }
180        // info on FGLM
181        int z = commonZeroTest(Gp);
182        if (z == 0) {
183            logger.info("ideal zero dimensional, can use also FGLM algorithm");
184        }
185        // compute target term order Groebner base via Groebner Walk
186        G = walkGroebnerToTarget(modv, Gp, pfac);
187        return G;
188    }
189
190
191    /**
192     * Converts Groebner bases w.r.t. total degree / start term order to
193     * Groebner base w.r.t to inverse lexicographical / target term order.
194     * @param modv module variable number.
195     * @param Gl Groebner base with respect to graded / start term order.
196     * @param ufac target polynomial ring and term order.
197     * @return Groebner base w.r.t inverse lexicographical / target term order
198     */
199    public List<GenPolynomial<C>> walkGroebnerToTarget(int modv, List<GenPolynomial<C>> Gl,
200                    GenPolynomialRing<C> ufac) {
201        if (Gl == null || Gl.size() == 0) {
202            throw new IllegalArgumentException("G may not be null or empty");
203        }
204        //Polynomial ring of input Groebner basis Gl
205        GenPolynomialRing<C> ring = Gl.get(0).ring;
206        if (debug) {
207            logger.info("ring = " + ring.toScript());
208        }
209        //Polynomial ring of newGB with INVLEX / target term order
210        //TermOrder lexi = ufac.tord; //new TermOrder(TermOrder.INVLEX);
211        logger.info("G walk from ev1 = " + ring.tord.toScript() + " to ev2 = " + ufac.tord.toScript());
212        List<GenPolynomial<C>> Giter = Gl;
213        // determine initial marks
214        List<ExpVector> marks = new ArrayList<ExpVector>();
215        List<Monomial<C>> M = new ArrayList<Monomial<C>>();
216        for (GenPolynomial<C> p : Giter) {
217            marks.add(p.leadingExpVector());
218            M.add(new Monomial<C>(p.leadingMonomial()));
219        }
220        logger.info("marks = " + marks);
221        List<Monomial<C>> Mp = new ArrayList<Monomial<C>>(M);
222        // weight matrix for target term order
223        long[][] ufweight = TermOrderByName.weightForOrder(ufac.tord, ring.nvar);
224        //TermOrder word = TermOrder.reverseWeight(ufweight);
225        TermOrder word = new TermOrder(ufweight);
226        logger.info("weight order: " + word);
227        ufweight = word.getWeight(); // because of weightDeg usage
228
229        // loop throught term orders
230        ExpVector w = null;
231        int iter = 0; // count #loops
232        boolean done = false;
233        while (!done) {
234            iter++;
235            // determine V and w
236            PolynomialList<C> Pl = new PolynomialList<C>(ring, Giter);
237            SortedSet<ExpVector> delta = Pl.deltaExpVectors(marks);
238            //logger.info("delta(marks) = " + delta);
239            logger.info("w_old = " + w);
240            ExpVector v = facetNormal(ring.tord, ufac.tord, delta, ring.evzero, ufweight);
241            logger.info("minimal v = " + v);
242            if (v == null) {
243                done = true;
244                break; // finished
245            }
246            w = v;
247            // determine facet polynomials for w
248            List<GenPolynomial<C>> iG = new ArrayList<GenPolynomial<C>>();
249            int i = 0;
250            for (GenPolynomial<C> f : Giter) {
251                ExpVector h = marks.get(i++);
252                GenPolynomial<C> ing = f.leadingFacetPolynomial(h, w);
253                if (debug) {
254                    logger.info("ing_g = [" + ing + "], lt(ing) = "
255                                    + ing.ring.toScript(ing.leadingExpVector()) + ", f = "
256                                    + f.ring.toScript(f.leadingExpVector()));
257                }
258                iG.add(ing);
259            }
260            List<GenPolynomial<C>> inOmega = ufac.copy(iG);
261            if (debug) {
262                logger.info("inOmega = " + inOmega);
263                logger.info("inOmega.ring: " + inOmega.get(0).ring.toScript());
264            }
265
266            // INVLEX / target term order GB of inOmega
267            List<GenPolynomial<C>> inOG = sgb.GB(modv, inOmega);
268            if (debug) {
269                logger.info("GB(inOmega) = " + inOG);
270            }
271            // remark polynomials 
272            marks.clear();
273            M.clear();
274            for (GenPolynomial<C> p : inOG) {
275                marks.add(p.leadingExpVector());
276                M.add(new Monomial<C>(p.leadingMonomial()));
277            }
278            logger.info("new marks/M = " + marks);
279            // lift and reduce
280            List<GenPolynomial<C>> G = liftReductas(M, Mp, Giter, inOG);
281            if (debug || (iter % iterPrint == 0)) {
282                logger.info("lift(" + iter + ") inOG, new GB: " + G);
283            }
284            if (G.size() == 1) { // will not happen
285                GenPolynomial<C> p = ufac.copy(G.get(0)); // change term order
286                G.clear();
287                G.add(p);
288                return G;
289            }
290            // iterate
291            Giter = G;
292            Mp.clear();
293            Mp.addAll(M);
294        }
295        return Giter;
296    }
297
298
299    /**
300     * Determine new facet normal.
301     * @param t1 old term order.
302     * @param t2 new term order.
303     * @param delta exponent vectors deltas.
304     * @param zero exponent vector.
305     * @param t2weight weight representation of t2.
306     * @return new facet normal v or null if no new facet normal exists.
307     */
308    public ExpVector facetNormal(TermOrder t1, TermOrder t2, Set<ExpVector> delta, ExpVector zero,
309                    long[][] t2weight) {
310        TermOrder.EVComparator ev1 = t1.getAscendComparator();
311        TermOrder.EVComparator ev2 = t2.getAscendComparator();
312        ExpVector v = null;
313        long d = 0; // = Long.MIN_VALUE;
314        for (ExpVector e : delta) {
315            if (ev1.compare(zero, e) >= 0 || ev2.compare(zero, e) <= 0) { //ring.evzero
316                //logger.info("skip e = " + e);
317                continue;
318            }
319            int s = 0;
320            long d2 = 0;
321            ExpVector vt = null;
322            ExpVector et = null;
323            if (v == null) {
324                v = e;
325                logger.info("init v = " + v);
326                continue;
327            }
328            for (long[] tau : t2weight) {
329                //logger.info("current tau = " + Arrays.toString(tau));
330                //compare
331                d = v.weightDeg(tau);
332                d2 = e.weightDeg(tau);
333                vt = v.scalarMultiply(d2);
334                et = e.scalarMultiply(d);
335                s = ev1.compare(et, vt);
336                if (s == 0) {
337                    continue; // next tau
338                } else if (s > 0) { // <, (> by example)
339                    v = e;
340                    break;
341                } else {
342                    break;
343                }
344            }
345            //logger.info("step s  = " + s + ", d = " + d + ", d2 = " + d2 + ", e = " + e); 
346            //   + ", v = " + v + ", et = " + et + ", vt = " + vt);
347        }
348        return v;
349    }
350
351
352    /**
353     * Cleanup and terminate ThreadPool.
354     */
355    @Override
356    public void terminate() {
357        if (sgb == null) {
358            return;
359        }
360        sgb.terminate();
361    }
362
363
364    /**
365     * Cancel ThreadPool.
366     */
367    @Override
368    public int cancel() {
369        if (sgb == null) {
370            return 0;
371        }
372        return sgb.cancel();
373    }
374
375
376    /**
377     * Lift leading polynomials to full Groebner base with respect to term
378     * order.
379     * @param M new leading monomial list of polynomials as marks.
380     * @param Mp old leading monomial list of polynomials as marks.
381     * @param G Groebner base polynomials for lift.
382     * @param A polynomial list of leading omega polynomials to lift.
383     * @return lift(A) a Groebner base wrt M of ideal(G).
384     */
385    public List<GenPolynomial<C>> liftReductas(List<Monomial<C>> M, List<Monomial<C>> Mp,
386                    List<GenPolynomial<C>> G, List<GenPolynomial<C>> A) {
387        if (G == null || M == null || Mp == null || G.isEmpty()) {
388            throw new IllegalArgumentException("null or empty lists not allowed");
389        }
390        if (A == null || A.isEmpty()) {
391            return A;
392        }
393        if (G.size() != Mp.size() || A.size() != M.size()) {
394            throw new IllegalArgumentException("equal sized lists required");
395        }
396        GenPolynomial<C> pol = G.get(0);
397        if (pol == null) {
398            throw new IllegalArgumentException("null polynomial not allowed");
399        }
400        // remove mark monomials
401        List<GenPolynomial<C>> Gp = new ArrayList<GenPolynomial<C>>(G.size());
402        int i = 0;
403        int len = G.size();
404        for (i = 0; i < len; i++) {
405            Monomial<C> mon = Mp.get(i);
406            GenPolynomial<C> s = G.get(i).subtract(mon);
407            Gp.add(s);
408        }
409        if (debug) {
410            logger.info("lifter GB: Gp  = " + Gp + ", Mp = " + Mp);
411        }
412        // compute f^Gp for f in A
413        //GenPolynomialRing<C> oring = pol.ring;
414        logger.info("liftReductas: G = " + G.size() + ", Mp = " + Mp.size()); // + ", old = " + oring.toScript());
415        List<GenPolynomial<C>> Ap = A; //oring.copy(A);
416        //logger.info("to lift Ap = " + Ap);
417        ReductionAbstract<C> sred = (ReductionAbstract<C>) sgb.red; //new ReductionSeq<C>();
418        List<GenPolynomial<C>> red = new ArrayList<GenPolynomial<C>>();
419        GenPolynomialRing<C> tring = A.get(0).ring;
420        len = Ap.size();
421        for (i = 0; i < len; i++) {
422            GenPolynomial<C> a = Ap.get(i);
423            GenPolynomial<C> r = sred.normalformMarked(Mp, Gp, a);
424            red.add(r);
425        }
426        logger.info("liftReductas: red(A) = " + red.size());
427        // combine f - f^Gp in tring
428        if (debug) {
429            logger.info("tring = " + tring.toScript()); // + ", M = " + M);
430        }
431        List<GenPolynomial<C>> nb = new ArrayList<GenPolynomial<C>>(red.size());
432        for (i = 0; i < A.size(); i++) {
433            GenPolynomial<C> x = tring.copy(A.get(i)); // Ap? A!
434            GenPolynomial<C> r = tring.copy(red.get(i));
435            GenPolynomial<C> s = x.subtract(r);
436            Monomial<C> m = M.get(i);
437            s.doAddTo(m.coefficient().negate(), m.exponent()); // remove marked term
438            if (!s.coefficient(m.exponent()).isZERO()) {
439                System.out.println("L-M: x = " + x + ", r = " + r);
440                throw new IllegalArgumentException("mark not removed: " + s + ", m = " + m);
441            }
442            nb.add(s);
443        }
444        if (debug) {
445            logger.info("lifted-M, nb = " + nb.size());
446        }
447        // minimal GB with preserved marks
448        //Collections.reverse(nb); // important for lex GB
449        len = nb.size();
450        i = 0;
451        while (i < len) {
452            GenPolynomial<C> a = nb.remove(0);
453            Monomial<C> m = M.remove(0); // in step with element from nb
454            if (debug) {
455                logger.info("doing " + a + ", lt = " + tring.toScript(m.exponent()));
456            }
457            //a = sgb.red.normalform(nb, a);
458            a = sred.normalformMarked(M, nb, a);
459            if (debug) {
460                logger.info("done, a = " + a + ", lt = " + tring.toScript(a.leadingExpVector()));
461            }
462            nb.add(a); // adds as last
463            M.add(m);
464            i++;
465        }
466        // re-add mark after minimal
467        for (i = 0; i < len; i++) {
468            GenPolynomial<C> a = nb.get(i);
469            Monomial<C> m = M.get(i);
470            a.doAddTo(m.coefficient(), m.exponent()); // re-add marked term
471            nb.set(i, a);
472        }
473        logger.info("liftReductas: nb = " + nb.size() + ", M = " + M.size());
474        //Collections.reverse(nb); // undo reverse
475        return nb;
476    }
477
478}