001/*
002 * $Id: GroebnerBaseSeqPairSeq.java 3414 2010-12-19 12:47:23Z kredel $
003 */
004
005package edu.jas.gb;
006
007
008import java.util.ArrayList;
009import java.util.List;
010import java.util.ListIterator;
011
012import org.apache.log4j.Logger;
013
014import edu.jas.poly.GenPolynomial;
015import edu.jas.poly.GenPolynomialRing;
016import edu.jas.structure.RingElem;
017
018
019/**
020 * Groebner Base sequential algorithm. Implements Groebner bases and GB test.
021 * Uses sequential pair list class.
022 * @param <C> coefficient type
023 * @author Heinz Kredel
024 */
025
026public class GroebnerBaseSeqPairSeq<C extends RingElem<C>> extends GroebnerBaseAbstract<C> {
027
028
029    private static final Logger logger = Logger.getLogger(GroebnerBaseSeqPairSeq.class);
030
031
032    private final boolean debug = logger.isDebugEnabled();
033
034
035    /**
036     * Constructor.
037     */
038    public GroebnerBaseSeqPairSeq() {
039        super();
040    }
041
042
043    /**
044     * Constructor.
045     * @param red Reduction engine
046     */
047    public GroebnerBaseSeqPairSeq(Reduction<C> red) {
048        super(red);
049    }
050
051
052    /**
053     * Groebner base using pairlist class.
054     * @param modv module variable number.
055     * @param F polynomial list.
056     * @return GB(F) a Groebner base of F.
057     */
058    public List<GenPolynomial<C>> GB(int modv, List<GenPolynomial<C>> F) {
059        GenPolynomial<C> p;
060        List<GenPolynomial<C>> G = new ArrayList<GenPolynomial<C>>();
061        CriticalPairList<C> pairlist = null;
062        int len = F.size();
063        ListIterator<GenPolynomial<C>> it = F.listIterator();
064        while (it.hasNext()) {
065            p = it.next();
066            if (p.length() > 0) {
067                p = p.monic();
068                if (p.isONE()) {
069                    G.clear();
070                    G.add(p);
071                    return G; // since no threads are activated
072                }
073                G.add(p);
074                if (pairlist == null) {
075                    pairlist = new CriticalPairList<C>(modv, p.ring);
076                }
077                // putOne not required
078                pairlist.put(p);
079            } else {
080                len--;
081            }
082        }
083        if (len <= 1) {
084            return G; // since no threads are activated
085        }
086
087        CriticalPair<C> pair;
088        GenPolynomial<C> pi;
089        GenPolynomial<C> pj;
090        GenPolynomial<C> S;
091        GenPolynomial<C> H;
092        while (pairlist.hasNext()) {
093            pair = pairlist.getNext();
094            if (pair == null) {
095                pairlist.update(); // ?
096                continue;
097            }
098            pi = pair.pi;
099            pj = pair.pj;
100            if (debug) {
101                logger.debug("pi    = " + pi);
102                logger.debug("pj    = " + pj);
103            }
104
105            S = red.SPolynomial(pi, pj);
106            if (S.isZERO()) {
107                pairlist.update(pair, S);
108                continue;
109            }
110            if (debug) {
111                logger.debug("ht(S) = " + S.leadingExpVector());
112            }
113
114            H = red.normalform(G, S);
115            if (H.isZERO()) {
116                pairlist.update(pair, H);
117                continue;
118            }
119            if (debug) {
120                logger.debug("ht(H) = " + H.leadingExpVector());
121            }
122
123            H = H.monic();
124            if (H.isONE()) {
125                // pairlist.record( pair, H );
126                G.clear();
127                G.add(H);
128                return G; // since no threads are activated
129            }
130            if (debug) {
131                logger.debug("H = " + H);
132            }
133            G.add(H);
134            pairlist.update(pair, H);
135            //pairlist.update();
136        }
137        logger.debug("#sequential list = " + G.size());
138        G = minimalGB(G);
139        logger.info("" + pairlist);
140        return G;
141    }
142
143
144    /**
145     * Extended Groebner base using critical pair class.
146     * @param modv module variable number.
147     * @param F polynomial list.
148     * @return a container for an extended Groebner base of F.
149     */
150    @Override
151    public ExtendedGB<C> extGB(int modv, List<GenPolynomial<C>> F) {
152        List<GenPolynomial<C>> G = new ArrayList<GenPolynomial<C>>();
153        List<List<GenPolynomial<C>>> F2G = new ArrayList<List<GenPolynomial<C>>>();
154        List<List<GenPolynomial<C>>> G2F = new ArrayList<List<GenPolynomial<C>>>();
155        CriticalPairList<C> pairlist = null;
156        boolean oneInGB = false;
157        int len = F.size();
158
159        List<GenPolynomial<C>> row = null;
160        List<GenPolynomial<C>> rows = null;
161        List<GenPolynomial<C>> rowh = null;
162        GenPolynomialRing<C> ring = null;
163        GenPolynomial<C> H;
164        GenPolynomial<C> p;
165
166        int nzlen = 0;
167        for (GenPolynomial<C> f : F) {
168            if (f.length() > 0) {
169                nzlen++;
170            }
171            if (ring == null) {
172                ring = f.ring;
173            }
174        }
175        GenPolynomial<C> mone = ring.getONE(); //.negate();
176        int k = 0;
177        ListIterator<GenPolynomial<C>> it = F.listIterator();
178        while (it.hasNext()) {
179            p = it.next();
180            if (p.length() > 0) {
181                row = new ArrayList<GenPolynomial<C>>(nzlen);
182                for (int j = 0; j < nzlen; j++) {
183                    row.add(null);
184                }
185                //C c = p.leadingBaseCoefficient();
186                //c = c.inverse();
187                //p = p.multiply( c );
188                row.set(k, mone); //.multiply(c) );
189                k++;
190                if (p.isUnit()) {
191                    G.clear();
192                    G.add(p);
193                    G2F.clear();
194                    G2F.add(row);
195                    oneInGB = true;
196                    break;
197                }
198                G.add(p);
199                G2F.add(row);
200                if (pairlist == null) {
201                    pairlist = new CriticalPairList<C>(modv, p.ring);
202                }
203                // putOne not required
204                pairlist.put(p);
205            } else {
206                len--;
207            }
208        }
209        ExtendedGB<C> exgb;
210        if (len <= 1 || oneInGB) {
211            // adjust F2G
212            for (GenPolynomial<C> f : F) {
213                row = new ArrayList<GenPolynomial<C>>(G.size());
214                for (int j = 0; j < G.size(); j++) {
215                    row.add(null);
216                }
217                H = red.normalform(row, G, f);
218                if (!H.isZERO()) {
219                    logger.error("nonzero H = " + H);
220                }
221                F2G.add(row);
222            }
223            exgb = new ExtendedGB<C>(F, G, F2G, G2F);
224            //System.out.println("exgb 1 = " + exgb);
225            return exgb;
226        }
227
228        CriticalPair<C> pair;
229        int i, j;
230        GenPolynomial<C> pi;
231        GenPolynomial<C> pj;
232        GenPolynomial<C> S;
233        GenPolynomial<C> x;
234        GenPolynomial<C> y;
235        //GenPolynomial<C> z;
236        while (pairlist.hasNext() && !oneInGB) {
237            pair = pairlist.getNext();
238            if (pair == null) {
239                pairlist.update(); // ?
240                continue;
241            }
242            i = pair.i;
243            j = pair.j;
244            pi = pair.pi;
245            pj = pair.pj;
246            if (debug) {
247                logger.info("i, pi    = " + i + ", " + pi);
248                logger.info("j, pj    = " + j + ", " + pj);
249            }
250
251            rows = new ArrayList<GenPolynomial<C>>(G.size());
252            for (int m = 0; m < G.size(); m++) {
253                rows.add(null);
254            }
255            S = red.SPolynomial(rows, i, pi, j, pj);
256            if (debug) {
257                logger.debug("is reduction S = " + red.isReductionNF(rows, G, ring.getZERO(), S));
258            }
259            if (S.isZERO()) {
260                pairlist.update(pair, S);
261                // do not add to G2F
262                continue;
263            }
264            if (debug) {
265                logger.debug("ht(S) = " + S.leadingExpVector());
266            }
267
268            rowh = new ArrayList<GenPolynomial<C>>(G.size());
269            for (int m = 0; m < G.size(); m++) {
270                rowh.add(null);
271            }
272            H = red.normalform(rowh, G, S);
273            if (debug) {
274                logger.debug("is reduction H = " + red.isReductionNF(rowh, G, S, H));
275            }
276            if (H.isZERO()) {
277                pairlist.update(pair, H);
278                // do not add to G2F
279                continue;
280            }
281            if (debug) {
282                logger.debug("ht(H) = " + H.leadingExpVector());
283            }
284
285            row = new ArrayList<GenPolynomial<C>>(G.size() + 1);
286            for (int m = 0; m < G.size(); m++) {
287                x = rows.get(m);
288                if (x != null) {
289                    //System.out.println("ms = " + m + " " + x);
290                    x = x.negate();
291                }
292                y = rowh.get(m);
293                if (y != null) {
294                    y = y.negate();
295                    //System.out.println("mh = " + m + " " + y);
296                }
297                if (x == null) {
298                    x = y;
299                } else {
300                    x = x.sum(y);
301                }
302                //System.out.println("mx = " + m + " " + x);
303                row.add(x);
304            }
305            if (debug) {
306                logger.debug("is reduction 0+sum(row,G) == H : "
307                        + red.isReductionNF(row, G, H, ring.getZERO()));
308            }
309            row.add(null);
310
311
312            //  H = H.monic();
313            C c = H.leadingBaseCoefficient();
314            c = c.inverse();
315            H = H.multiply(c);
316            row = blas.scalarProduct(mone.multiply(c), row);
317            row.set(G.size(), mone);
318            if (H.isONE()) {
319                // pairlist.record( pair, H );
320                // G.clear(); 
321                G.add(H);
322                G2F.add(row);
323                oneInGB = true;
324                break;
325            }
326            if (debug) {
327                logger.debug("H = " + H);
328            }
329            G.add(H);
330            pairlist.update(pair, H);
331            G2F.add(row);
332        }
333        if (debug) {
334            exgb = new ExtendedGB<C>(F, G, F2G, G2F);
335            logger.info("exgb unnorm = " + exgb);
336        }
337        G2F = normalizeMatrix(F.size(), G2F);
338        if (debug) {
339            exgb = new ExtendedGB<C>(F, G, F2G, G2F);
340            logger.info("exgb nonmin = " + exgb);
341            boolean t2 = isReductionMatrix(exgb);
342            logger.info("exgb t2 = " + t2);
343        }
344        exgb = minimalExtendedGB(F.size(), G, G2F);
345        G = exgb.G;
346        G2F = exgb.G2F;
347        logger.debug("#sequential list = " + G.size());
348        logger.info("" + pairlist); 
349        // setup matrices F and F2G
350        for (GenPolynomial<C> f : F) {
351            row = new ArrayList<GenPolynomial<C>>(G.size());
352            for (int m = 0; m < G.size(); m++) {
353                row.add(null);
354            }
355            H = red.normalform(row, G, f);
356            if (!H.isZERO()) {
357                logger.error("nonzero H = " + H);
358            }
359            F2G.add(row);
360        }
361        exgb = new ExtendedGB<C>(F, G, F2G, G2F);
362        if (debug) {
363            logger.info("exgb nonmin = " + exgb);
364            boolean t2 = isReductionMatrix(exgb);
365            logger.info("exgb t2 = " + t2);
366        }
367        return exgb;
368    }
369
370}