Source code

001/*
002 * $Id$
003 */
004
005package edu.jas.gb;
006
007
008import java.util.ArrayList;
009import java.util.List;
010import java.util.ListIterator;
011
012import org.apache.logging.log4j.LogManager;
013import org.apache.logging.log4j.Logger;
014
015import edu.jas.poly.ExpVector;
016import edu.jas.poly.GenPolynomial;
017import edu.jas.poly.GenPolynomialRing;
018import edu.jas.structure.RingElem;
019import edu.jas.structure.NotInvertibleException;
020
021
022/**
023 * D-Groebner Base sequential algorithm. Implements D-Groebner bases and GB
024 * test. <b>Note:</b> Minimal reduced GBs are not unique. see BWK, section 10.1.
025 * @param <C> coefficient type
026 * @author Heinz Kredel
027 */
028
029public class DGroebnerBaseSeq<C extends RingElem<C>> extends GroebnerBaseAbstract<C> {
030
031
032    private static final Logger logger = LogManager.getLogger(DGroebnerBaseSeq.class);
033
034
035    private static final boolean debug = logger.isDebugEnabled();
036
037
038    /**
039     * Reduction engine.
040     */
041    protected DReduction<C> dred; // shadow super.red ??
042
043
044    /**
045     * Constructor.
046     */
047    public DGroebnerBaseSeq() {
048        this(new DReductionSeq<C>());
049    }
050
051
052    /**
053     * Constructor.
054     * @param dred D-Reduction engine
055     */
056    public DGroebnerBaseSeq(DReduction<C> dred) {
057        super(dred);
058        this.dred = dred;
059        assert this.dred == super.red;
060    }
061
062
063    /**
064     * D-Groebner base test.
065     * @param modv module variable number.
066     * @param F polynomial list.
067     * @return true, if F is a D-Groebner base, else false.
068     */
069    @Override
070    public boolean isGB(int modv, List<GenPolynomial<C>> F) {
071        GenPolynomial<C> pi, pj, s, d;
072        for (int i = 0; i < F.size(); i++) {
073            pi = F.get(i);
074            for (int j = i + 1; j < F.size(); j++) {
075                pj = F.get(j);
076                if (!dred.moduleCriterion(modv, pi, pj)) {
077                    continue;
078                }
079                d = dred.GPolynomial(pi, pj);
080                if (!d.isZERO()) {
081                    // better check for top reduction only
082                    d = dred.normalform(F, d);
083                }
084                if (!d.isZERO()) {
085                    System.out.println("d-pol(" + pi + "," + pj + ") != 0: " + d);
086                    return false;
087                }
088                // works ok
089                if (!dred.criterion4(pi, pj)) {
090                    continue;
091                }
092                s = dred.SPolynomial(pi, pj);
093                if (!s.isZERO()) {
094                    s = dred.normalform(F, s);
095                }
096                if (!s.isZERO()) {
097                    System.out.println("s-pol(" + i + "," + j + ") != 0: " + s);
098                    return false;
099                }
100            }
101        }
102        return true;
103    }
104
105
106    /**
107     * D-Groebner base using pairlist class.
108     * @param modv module variable number.
109     * @param F polynomial list.
110     * @return GB(F) a D-Groebner base of F.
111     */
112    public List<GenPolynomial<C>> GB(int modv, List<GenPolynomial<C>> F) {
113        List<GenPolynomial<C>> G = normalizeZerosOnes(F);
114        if (G.size() <= 1) {
115            return G;
116        }
117        GenPolynomialRing<C> ring = G.get(0).ring;
118        OrderedDPairlist<C> pairlist = new OrderedDPairlist<C>(modv, ring);
119        pairlist.put(G);
120
121        Pair<C> pair;
122        GenPolynomial<C> pi, pj, S, D, H;
123        while (pairlist.hasNext()) {
124            pair = pairlist.removeNext();
125            //System.out.println("pair = " + pair);
126            if (pair == null)
127                continue;
128            pi = pair.pi;
129            pj = pair.pj;
130            if (debug) {
131                logger.debug("pi    = {}", pi);
132                logger.debug("pj    = {}", pj);
133            }
134            H = null;
135
136            // D-polynomial case ----------------------
137            D = dred.GPolynomial(pi, pj);
138            //System.out.println("D_d = " + D);
139            if (!D.isZERO() && !dred.isTopReducible(G, D)) {
140                H = dred.normalform(G, D);
141                if (H.isONE()) {
142                    G.clear();
143                    G.add(H);
144                    return G; // since no threads are activated
145                }
146                if (!H.isZERO()) {
147                    logger.info("Dred = {}", H);
148                    //l++;
149                    G.add(H);
150                    pairlist.put(H);
151                }
152            }
153
154            // S-polynomial case -----------------------
155            if (pair.getUseCriterion3() && pair.getUseCriterion4()) {
156                S = dred.SPolynomial(pi, pj);
157                //System.out.println("S_d = " + S);
158                if (S.isZERO()) {
159                    pair.setZero();
160                    //continue;
161                }
162                if (logger.isDebugEnabled()) {
163                    logger.debug("ht(S) = {}", S.leadingExpVector());
164                }
165                H = dred.normalform(G, S);
166                if (H.isZERO()) {
167                    pair.setZero();
168                    //continue;
169                }
170                if (logger.isDebugEnabled()) {
171                    logger.debug("ht(H) = {}", H.leadingExpVector());
172                }
173                if (H.isONE()) {
174                    G.clear();
175                    G.add(H);
176                    return G; // since no threads are activated
177                }
178                logger.debug("H = {}", H);
179                if (!H.isZERO()) {
180                    logger.info("Sred = {}", H);
181                    //len = G.size();
182                    //l++;
183                    G.add(H);
184                    pairlist.put(H);
185                }
186            }
187        }
188        logger.debug("#sequential list = {}", G.size());
189        G = minimalGB(G);
190        logger.info("{}", pairlist);
191        return G;
192    }
193
194
195    /**
196     * Extended Groebner base using pairlist class.
197     * @param modv module variable number.
198     * @param F polynomial list.
199     * @return a container for an extended Groebner base of F.
200     */
201    @Override
202    public ExtendedGB<C> extGB(int modv, List<GenPolynomial<C>> F) {
203        if (F == null || F.isEmpty()) {
204            throw new IllegalArgumentException("null or empty F not allowed");
205        }
206        List<GenPolynomial<C>> G = new ArrayList<GenPolynomial<C>>();
207        List<List<GenPolynomial<C>>> F2G = new ArrayList<List<GenPolynomial<C>>>();
208        List<List<GenPolynomial<C>>> G2F = new ArrayList<List<GenPolynomial<C>>>();
209        OrderedDPairlist<C> pairlist = null;
210        boolean oneInGB = false;
211        int len = F.size();
212
213        List<GenPolynomial<C>> row = null;
214        List<GenPolynomial<C>> rows = null;
215        List<GenPolynomial<C>> rowh = null;
216        GenPolynomialRing<C> ring = null;
217        GenPolynomial<C> H;
218        GenPolynomial<C> p;
219
220        int nzlen = 0;
221        for (GenPolynomial<C> f : F) {
222            if (f.length() > 0) {
223                nzlen++;
224            }
225            if (ring == null) {
226                ring = f.ring;
227            }
228        }
229        GenPolynomial<C> one = ring.getONE();
230        int k = 0;
231        ListIterator<GenPolynomial<C>> it = F.listIterator();
232        while (it.hasNext()) {
233            p = it.next();
234            if (p.length() > 0) {
235                row = blas.genVector(nzlen, null);
236                row.set(k, one);
237                k++;
238                if (p.isUnit()) {
239                    G.clear();
240                    G.add(p);
241                    G2F.clear();
242                    G2F.add(row);
243                    oneInGB = true;
244                    break;
245                }
246                G.add(p);
247                //logger.info("p row = {}", row);
248                G2F.add(row);
249                if (pairlist == null) {
250                    //pairlist = strategy.create(modv, p.ring);
251                    pairlist = new OrderedDPairlist<C>(modv, p.ring);
252                    if (p.ring.coFac.isField()) {
253                        //throw new RuntimeException("coefficients from a field");
254                        logger.warn("coefficients from a field " + p.ring.coFac);
255                    }
256                }
257                // putOne not required
258                pairlist.put(p);
259            } else {
260                len--;
261            }
262        }
263        ExtendedGB<C> exgb;
264        if (len <= 1 || oneInGB) {
265            // adjust F2G
266            for (GenPolynomial<C> f : F) {
267                row = blas.genVector(G.size(), null);
268                H = dred.normalform(row, G, f);
269                if (!H.isZERO()) {
270                    logger.error("nonzero H = {}", H);
271                }
272                logger.info("f row = {}", row);
273                F2G.add(row);
274            }
275            exgb = new ExtendedGB<C>(F, G, F2G, G2F);
276            //System.out.println("exgb #1 = " + exgb);
277            return exgb;
278        }
279
280        Pair<C> pair;
281        int i, j;
282        GenPolynomial<C> pi, pj;
283        GenPolynomial<C> S, D;
284        k = 0;
285        while (pairlist.hasNext() && !oneInGB) {
286            pair = pairlist.removeNext();
287            if (pair == null) {
288                continue;
289            }
290            i = pair.i;
291            j = pair.j;
292            pi = pair.pi;
293            pj = pair.pj;
294            if (debug) {
295                logger.info("i, pi    = {}, {}", i, pi);
296                logger.info("j, pj    = {}, {}", j, pj);
297            }
298            H = null;
299
300            // D-polynomial case ----------------------
301            rows = blas.genVector(G.size() + 1, null);
302            D = dred.GPolynomial(rows, i, pi, j, pj);
303            logger.info("Gpol = {}", D);
304            if (debug) {
305                logger.info("is reduction D = " + dred.isReductionNF(rows, G, D, ring.getZERO()));
306            }
307            if (!D.isZERO()) { //&& !dred.isTopReducible(G, D)
308                //continue;
309                rowh = blas.genVector(G.size() + 1, null);
310                H = dred.normalform(rowh, G, D);
311                if (debug) {
312                    logger.info("is reduction H = " + dred.isReductionNF(rowh, G, D, H));
313                }
314                //H = H.monic();
315                int s = H.leadingBaseCoefficient().signum();
316                if (s < 0) {
317                    logger.info("negate: H_D rowd, rowh = {}, {}", rows, rowh);
318                    H = H.negate();
319                    rows = blas.vectorNegate(rows);
320                    rowh = blas.vectorNegate(rowh);
321                }
322                if (H.isONE()) {
323                    // G.clear();
324                    G.add(H);
325                    oneInGB = true;
326                } else if (!H.isZERO()) {
327                    logger.info("H_G red = {}", H);
328                    G.add(H);
329                    pairlist.put(H);
330                }
331                //System.out.println("rowd = " + rows);
332                //System.out.println("rowh = " + rowh);
333                row = blas.vectorCombineRep(rows, rowh);
334                logger.debug("H_G row = {}", row);
335                if (!H.isZERO()) {
336                    G2F.add(row);
337                }
338                if (debug) {
339                    logger.debug("ht(H) = {}", H.leadingExpVector());
340                    logger.info("is reduction D,H = " + dred.isReductionNF(row, G, H, ring.getZERO()));
341                }
342            }
343
344            // S-polynomial case ----------------------
345            rows = blas.genVector(G.size() + 1, null);
346            S = dred.SPolynomial(rows, i, pi, j, pj);
347            logger.info("Spol = {}", S);
348            if (debug) {
349                logger.info("is reduction S = " + dred.isReductionNF(rows, G, S, ring.getZERO()));
350            }
351            rowh = blas.genVector(G.size() + 1, null);
352            if (!S.isZERO()) { //&& !dred.isTopReducible(G, S)
353                //continue;
354                H = dred.normalform(rowh, G, S);
355                if (debug) {
356                    logger.info("Spol_red = {}", H);
357                    logger.info("is reduction H = " + dred.isReductionNF(rowh, G, S, H));
358                }
359                //H = H.monic();
360                int s = H.leadingBaseCoefficient().signum();
361                if (s < 0) {
362                    logger.info("negate: H_S rows, rowh = {}, {}", rows, rowh);
363                    H = H.negate();
364                    //rowh = rowh.negate(); //rowh.set(G.size(), one.negate());
365                    rows = blas.vectorNegate(rows);
366                    rowh = blas.vectorNegate(rowh);
367                }
368                //logger.info("Spol_red_norm = {}", H);
369                if (H.isONE()) {
370                    //G.clear();
371                    G.add(H);
372                    oneInGB = true;
373                } else if (!H.isZERO()) {
374                    logger.info("H_S red = {}", H);
375                    G.add(H);
376                    pairlist.put(H);
377                }
378                //System.out.println("rows = " + rows);
379                //System.out.println("rowh = " + rowh);
380                row = blas.vectorCombineRep(rows, rowh);
381                logger.debug("H_S row = {}", row);
382                if (!H.isZERO()) {
383                    G2F.add(row);
384                }
385                if (debug) {
386                    logger.debug("ht(H) = {}", H.leadingExpVector());
387                    logger.info("is reduction S,H = " + dred.isReductionNF(row, G, H, ring.getZERO()));
388                }
389            }
390        }
391        if (true || debug) {
392            exgb = new ExtendedGB<C>(F, G, F2G, G2F);
393            boolean t = isReductionMatrix(exgb);
394            if (!t) {
395                logger.info("exgb unnorm = {}", exgb);
396                logger.info("exgb t_1 = {}", t);
397            }
398        }
399        G2F = normalizeMatrix(F.size(), G2F);
400        if (debug) {
401            exgb = new ExtendedGB<C>(F, G, F2G, G2F);
402            boolean t = isReductionMatrix(exgb);
403            if (!t) {
404                logger.info("exgb norm nonmin = {}", exgb);
405                logger.info("exgb t_2 = {}", t);
406            }
407        }
408        exgb = minimalExtendedGB(F.size(), G, G2F);
409        G = exgb.G;
410        G2F = exgb.G2F;
411        if (true || debug) {
412            exgb = new ExtendedGB<C>(F, G, F2G, G2F);
413            boolean t = isMinReductionMatrix(exgb);
414            if (!t) {
415                logger.info("exgb minGB = {}", exgb);
416                logger.info("exgb t_3 = {}", t);
417            }
418        }
419        exgb = new ExtendedGB<C>(F, G, F2G, G2F);
420        // setup matrices F and F2G
421        for (GenPolynomial<C> f : F) {
422            row = blas.genVector(G.size() + 1, null);
423            H = dred.normalform(row, G, f);
424            if (!H.isZERO()) {
425                logger.error("nonzero H, G = {}, {}", H, G);
426                throw new RuntimeException("H != 0");
427            }
428            F2G.add(row);
429        }
430        exgb = new ExtendedGB<C>(F, G, F2G, G2F);
431        if (debug) {
432            boolean t = isMinReductionMatrix(exgb);
433            if (!t) {
434                logger.info("exgb +F+F2G = {}", exgb);
435                logger.info("exgb t_4 = {}", t);
436            }
437        }
438        return exgb;
439    }
440
441
442    /**
443     * Minimal extended groebner basis.
444     * @param flen length of rows.
445     * @param Gp a Groebner base.
446     * @param M a reduction matrix.
447     * @return a (partially) reduced Groebner base of Gp in a (fake) container.
448     */
449    @Override
450    public ExtendedGB<C> minimalExtendedGB(int flen, List<GenPolynomial<C>> Gp,
451                    List<List<GenPolynomial<C>>> M) {
452        if (Gp == null) {
453            return null; //new ExtendedGB<C>(null,Gp,null,M);
454        }
455        List<GenPolynomial<C>> G, F, T;
456        G = new ArrayList<GenPolynomial<C>>(Gp);
457        F = new ArrayList<GenPolynomial<C>>(Gp.size());
458
459        List<List<GenPolynomial<C>>> Mg, Mf;
460        Mg = new ArrayList<List<GenPolynomial<C>>>(M.size());
461        Mf = new ArrayList<List<GenPolynomial<C>>>(M.size());
462        for (List<GenPolynomial<C>> r : M) {
463            // must be copied also
464            List<GenPolynomial<C>> row = new ArrayList<GenPolynomial<C>>(r);
465            Mg.add(row);
466        }
467
468        boolean mt;
469        ListIterator<GenPolynomial<C>> it;
470        ArrayList<Integer> ix = new ArrayList<Integer>();
471        ArrayList<Integer> jx = new ArrayList<Integer>();
472        int k = 0;
473        //System.out.println("flen, Gp, M = " + flen + ", " + Gp.size() + ", " + M.size() );
474        GenPolynomialRing<C> pfac = null;
475        GenPolynomial<C> pi, pj;
476        ExpVector ei, ej;
477        C ai, aj, r;
478        while (G.size() > 0) {
479            pi = G.remove(0);
480            ei = pi.leadingExpVector();
481            ai = pi.leadingBaseCoefficient();
482            if (pfac == null) {
483                pfac = pi.ring;
484            }
485
486            it = G.listIterator();
487            mt = false;
488            while (it.hasNext() && !mt) {
489                pj = it.next();
490                ej = pj.leadingExpVector();
491                mt = ei.multipleOf(ej);
492                if (mt) {
493                    aj = pj.leadingBaseCoefficient();
494                    r = ai.remainder(aj);
495                    mt = r.isZERO(); // && mt
496                }
497            }
498            it = F.listIterator();
499            while (it.hasNext() && !mt) {
500                pj = it.next();
501                ej = pj.leadingExpVector();
502                mt = ei.multipleOf(ej);
503                if (mt) {
504                    aj = pj.leadingBaseCoefficient();
505                    r = ai.remainder(aj);
506                    mt = r.isZERO(); // && mt
507                }
508            }
509            T = new ArrayList<GenPolynomial<C>>(G);
510            T.addAll(F);
511            //GenPolynomial<C> t = dred.normalform(T, pi);
512            //System.out.println("t, mt, t==0: " + t + ", " + mt + ", " + t.isZERO());
513            if (!mt) {
514                F.add(pi);
515                ix.add(k);
516                //System.out.println("ix: " + ix);
517            } else { // drop polynomial and corresponding row and column ??
518                jx.add(k);
519                //System.out.println("jx: " + jx);
520            }
521            k++;
522        }
523        logger.info("ix, jx, #M = {}, {}, {}", ix, jx, Mg.size());
524        int fix = -1; // copied polys
525        // copy Mg to Mf as indicated by ix
526        for (int i = 0; i < ix.size(); i++) {
527            int u = ix.get(i);
528            if (u >= flen && fix == -1) {
529                fix = Mf.size();
530            }
531            //System.out.println("copy u_ix, fix = " + u + ", " + fix);
532            if (u >= 0) {
533                List<GenPolynomial<C>> row = Mg.get(u);
534                Mf.add(row);
535            }
536        }
537        // for (int i = 0; i < jx.size(); i++) {
538        //     int u = jx.get(i);
539        //     //System.out.println("copy u_jx = " + u);
540        //     if (u >= 0) {
541        //         List<GenPolynomial<C>> row = blas.genVector(flen,null);
542        //         Mf.add(u, row);
543        //         F.add(u, pfac.getZERO());
544        //     }
545        // }
546        if (F.size() <= 1 || fix == -1) {
547            return new ExtendedGB<C>(null, F, null, Mf);
548        }
549        // must return, since extended normalform has not correct order of polys
550        return new ExtendedGB<C>(null, F, null, Mf);
551        /*
552        G = F;
553        F = new ArrayList<GenPolynomial<C>>( G.size() );
554        List<GenPolynomial<C>> temp;
555        k = 0;
556        final int len = G.size();
557        while ( G.size() > 0 ) {
558            a = G.remove(0);
559            if ( k >= fix ) { // dont touch copied polys
560               row = Mf.get( k );
561               //System.out.println("doing k = " + k + ", " + a);
562               // must keep order, but removed polys missing
563               temp = new ArrayList<GenPolynomial<C>>( len );
564               temp.addAll( F );
565               temp.add( a.ring.getZERO() ); // ??
566               temp.addAll( G );
567               //System.out.println("row before = " + row);
568               a = red.normalform( row, temp, a );
569               //System.out.println("row after  = " + row);
570            }
571            F.add( a );
572            k++;
573        }
574        // does Mf need renormalization?
575        */
576    }
577
578
579    /**
580     * Inverse for element modulo ideal.
581     * @param h polynomial
582     * @param F polynomial list
583     * @return inverse of h with respect to ideal(F), if defined
584     */
585    public GenPolynomial<C> inverse(GenPolynomial<C> h, List<GenPolynomial<C>> F) {
586        if (h == null || h.isZERO()) {
587            throw new NotInvertibleException("zero not invertible");
588        }
589        if (F == null || F.size() == 0) {
590            throw new NotInvertibleException("zero ideal");
591        }
592        if (h.isUnit()) {
593            return h.inverse();
594        }
595        // prepare by GB precomputation
596        List<GenPolynomial<C>> G = GB(F);
597        List<GenPolynomial<C>> I = new ArrayList<GenPolynomial<C>>(1 + G.size());
598        I.add(h);
599        I.addAll(G);
600        // now compute extended gcd(h,G)
601        ExtendedGB<C> X = extGB(I);
602        List<GenPolynomial<C>> hG = X.G;
603        //System.out.println("hG = " + hG);
604        GenPolynomial<C> one = null;
605        int i = -1;
606        for (GenPolynomial<C> p : hG) {
607            i++;
608            if (p == null) {
609                continue;
610            }
611            if (p.isUnit()) {
612                one = p;
613                break;
614            }
615        }
616        if (one == null) {
617            throw new NotInvertibleException("h = " + h);
618        }
619        List<GenPolynomial<C>> row = X.G2F.get(i); // != -1
620        GenPolynomial<C> g = row.get(0);
621        if (g == null || g.isZERO()) {
622            throw new NotInvertibleException("h = " + h);
623        }
624        // adjust g to get g*h == 1 mod ideal(G)
625        GenPolynomial<C> f = g.multiply(h);
626        GenPolynomial<C> k = red.normalform(G, f);
627        //System.out.println("g = " + g + ", h = " + h + ", f = " + f + ", k = " + k);
628        if (k.signum() < 0) { // then is -1 or inv-G(0)
629            //System.out.println("k < 0: " + G);
630            g = g.sum(G.get(0)); //.negate();
631        }
632        return g;
633    }
634
635}