001/*
002 * $Id: SolvableSyzygySeq.java 5310 2015-10-25 18:35:31Z kredel $
003 */
004
005package edu.jas.gbufd;
006
007
008import java.util.ArrayList;
009import java.util.Arrays;
010import java.util.Iterator;
011import java.util.List;
012
013import org.apache.log4j.Logger;
014
015import edu.jas.gb.SolvableExtendedGB;
016import edu.jas.gb.SolvableGroebnerBase;
017import edu.jas.poly.GenPolynomial;
018import edu.jas.poly.GenSolvablePolynomial;
019import edu.jas.poly.GenSolvablePolynomialRing;
020import edu.jas.poly.ModuleList;
021import edu.jas.poly.PolynomialList;
022import edu.jas.structure.GcdRingElem;
023import edu.jas.structure.RingFactory;
024
025
026/**
027 * Syzygy sequential class for solvable polynomials. Implements Syzygy
028 * computations and tests with Groebner bases.
029 * @param <C> coefficient type
030 * @author Heinz Kredel
031 */
032
033public class SolvableSyzygySeq<C extends GcdRingElem<C>> extends SolvableSyzygyAbstract<C> {
034
035
036    private static final Logger logger = Logger.getLogger(SolvableSyzygySeq.class);
037
038
039    private static final boolean debug = logger.isDebugEnabled();
040
041
042    private static boolean assertEnabled = false;
043
044
045    static {
046        assert assertEnabled = true; // official hack to check assertions enabled
047    }
048
049
050    /**
051     * Groebner basis engine.
052     */
053    protected SolvableGroebnerBase<C> sbb;
054
055
056    /*
057     * Module Groebner basis engine.
058    //protected ModSolvableGroebnerBase<C> msbb;
059     */
060
061
062    /**
063     * Constructor.
064     * @param cf coefficient ring.
065     */
066    public SolvableSyzygySeq(RingFactory<C> cf) {
067        sbb = SGBFactory.getImplementation(cf);
068        //msbb = new ModSolvableGroebnerBaseSeq<C>(cf);
069    }
070
071
072    /**
073     * Resolution of a module. Only with direct GBs.
074     * @param M a module list of a Groebner basis.
075     * @return a resolution of M.
076     */
077    public List<SolvResPart<C>> resolution(ModuleList<C> M) {
078        List<SolvResPart<C>> R = new ArrayList<SolvResPart<C>>();
079        ModuleList<C> MM = M;
080        ModuleList<C> GM;
081        ModuleList<C> Z;
082        while (true) {
083            GM = sbb.leftGB(MM);
084            Z = leftZeroRelations(GM);
085            R.add(new SolvResPart<C>(MM, GM, Z));
086            if (Z == null || Z.list == null || Z.list.size() == 0) {
087                break;
088            }
089            MM = Z;
090        }
091        return R;
092    }
093
094
095    /**
096     * Resolution of a polynomial list. Only with direct GBs.
097     * @param F a polynomial list of a Groebner basis.
098     * @return a resolution of F.
099     */
100    @SuppressWarnings("unchecked")
101    public List/*<SolvResPart<C>|SolvResPolPart<C>>*/resolution(PolynomialList<C> F) {
102        List<List<GenSolvablePolynomial<C>>> Z;
103        ModuleList<C> Zm;
104        List<GenSolvablePolynomial<C>> G;
105        PolynomialList<C> Gl;
106
107        G = sbb.leftGB(F.castToSolvableList());
108        Z = leftZeroRelations(G);
109        Gl = new PolynomialList<C>((GenSolvablePolynomialRing<C>) F.ring, G);
110        Zm = new ModuleList<C>((GenSolvablePolynomialRing<C>) F.ring, Z);
111
112        List R = resolution(Zm);
113        R.add(0, new SolvResPolPart<C>(F, Gl, Zm));
114        return R;
115    }
116
117
118    /**
119     * Resolution of a module.
120     * @param M a module list of an arbitrary basis.
121     * @return a resolution of M.
122     */
123    public List<SolvResPart<C>> resolutionArbitrary(ModuleList<C> M) {
124        List<SolvResPart<C>> R = new ArrayList<SolvResPart<C>>();
125        ModuleList<C> MM = M;
126        ModuleList<C> GM = null;
127        ModuleList<C> Z;
128        while (true) {
129            //GM = sbb.leftGB(MM);
130            Z = leftZeroRelationsArbitrary(MM);
131            R.add(new SolvResPart<C>(MM, GM, Z));
132            if (Z == null || Z.list == null || Z.list.size() == 0) {
133                break;
134            }
135            MM = Z;
136        }
137        return R;
138    }
139
140
141    /**
142     * Resolution of a polynomial list.
143     * @param F a polynomial list of an arbitrary basis.
144     * @return a resolution of F.
145     */
146    @SuppressWarnings("unchecked")
147    public List/*<SolvResPart<C>|SolvResPolPart<C>>*/resolutionArbitrary(PolynomialList<C> F) {
148        List<List<GenSolvablePolynomial<C>>> Z;
149        ModuleList<C> Zm;
150        PolynomialList<C> Gl = null;
151        //List<GenSolvablePolynomial<C>> G = sbb.leftGB( F.castToSolvableList() );
152        //Gl = new PolynomialList<C>((GenSolvablePolynomialRing<C>)F.ring, G);
153        Z = leftZeroRelationsArbitrary(F.castToSolvableList());
154        Zm = new ModuleList<C>((GenSolvablePolynomialRing<C>) F.ring, Z);
155
156        List R = resolutionArbitrary(Zm);
157        R.add(0, new SolvResPolPart<C>(F, Gl, Zm));
158        return R;
159    }
160
161
162    /**
163     * Left syzygy module from arbitrary base.
164     * @param modv number of module variables.
165     * @param F a solvable polynomial list.
166     * @return syz(F), a basis for the module of left syzygies for F.
167     */
168    @SuppressWarnings("unchecked")
169    public List<List<GenSolvablePolynomial<C>>> leftZeroRelationsArbitrary(int modv,
170                    List<GenSolvablePolynomial<C>> F) {
171        if (F == null) {
172            return null; //leftZeroRelations( modv, F );
173        }
174        if (F.size() <= 1) {
175            return leftZeroRelations(modv, F);
176        }
177        final int lenf = F.size();
178        SolvableExtendedGB<C> exgb = sbb.extLeftGB(modv, F);
179        if (debug) {
180            logger.info("exgb = " + exgb);
181        }
182        if (assertEnabled) {
183            logger.info("check1 exgb start");
184            if (!sbb.isLeftReductionMatrix(exgb)) {
185                logger.error("is reduction matrix ? false");
186            }
187            logger.info("check1 exgb end");
188        }
189        List<GenSolvablePolynomial<C>> G = exgb.G;
190        List<List<GenSolvablePolynomial<C>>> G2F = exgb.G2F;
191        List<List<GenSolvablePolynomial<C>>> F2G = exgb.F2G;
192
193        List<List<GenSolvablePolynomial<C>>> sg = leftZeroRelations(modv, G);
194        GenSolvablePolynomialRing<C> ring = G.get(0).ring;
195        ModuleList<C> S = new ModuleList<C>(ring, sg);
196        if (debug) {
197            logger.info("syz = " + S);
198        }
199        if (assertEnabled) {
200            logger.info("check2 left syz start");
201            if (!isLeftZeroRelation(sg, G)) {
202                logger.error("is syzygy ? false");
203            }
204            logger.info("check2 left syz end");
205        }
206
207        List<List<GenSolvablePolynomial<C>>> sf;
208        sf = new ArrayList<List<GenSolvablePolynomial<C>>>(sg.size());
209        //List<GenPolynomial<C>> row;
210
211        for (List<GenSolvablePolynomial<C>> r : sg) {
212            Iterator<GenSolvablePolynomial<C>> it = r.iterator();
213            Iterator<List<GenSolvablePolynomial<C>>> jt = G2F.iterator();
214
215            List<GenSolvablePolynomial<C>> rf;
216            rf = new ArrayList<GenSolvablePolynomial<C>>(lenf);
217            for (int m = 0; m < lenf; m++) {
218                rf.add(ring.getZERO());
219            }
220            while (it.hasNext() && jt.hasNext()) {
221                GenSolvablePolynomial<C> si = it.next();
222                List<GenSolvablePolynomial<C>> ai = jt.next();
223                //System.out.println("si = " + si);
224                //System.out.println("ai = " + ai);
225                if (si == null || ai == null) {
226                    continue;
227                }
228                // pi has wrong type:
229                List<GenPolynomial<C>> pi = blas.scalarProduct(si, PolynomialList.<C> castToList(ai));
230                //System.out.println("pi = " + pi);
231                rf = PolynomialList.<C> castToSolvableList(blas.vectorAdd(PolynomialList.<C> castToList(rf),
232                                pi));
233            }
234            if (it.hasNext() || jt.hasNext()) {
235                logger.error("leftZeroRelationsArbitrary wrong sizes");
236            }
237            //System.out.println("\nrf = " + rf + "\n");
238            sf.add(rf);
239        }
240        if (assertEnabled) {
241            logger.info("check3 left syz start");
242            if (!isLeftZeroRelation(sf, F)) {
243                logger.error("is partial syz sf ? false");
244            }
245            logger.info("check3 left syz end");
246        }
247
248        List<List<GenSolvablePolynomial<C>>> M;
249        M = new ArrayList<List<GenSolvablePolynomial<C>>>(lenf);
250        for (List<GenSolvablePolynomial<C>> r : F2G) {
251            Iterator<GenSolvablePolynomial<C>> it = r.iterator();
252            Iterator<List<GenSolvablePolynomial<C>>> jt = G2F.iterator();
253
254            List<GenSolvablePolynomial<C>> rf;
255            rf = new ArrayList<GenSolvablePolynomial<C>>(lenf);
256            for (int m = 0; m < lenf; m++) {
257                rf.add(ring.getZERO());
258            }
259            while (it.hasNext() && jt.hasNext()) {
260                GenSolvablePolynomial<C> si = it.next();
261                List<GenSolvablePolynomial<C>> ai = jt.next();
262                //System.out.println("si = " + si);
263                //System.out.println("ai = " + ai);
264                if (si == null || ai == null) {
265                    continue;
266                }
267                //pi has wrong type, should be: List<GenSolvablePolynomial<C>>
268                List<GenPolynomial<C>> pi = blas.scalarProduct(si, PolynomialList.<C> castToList(ai));
269                //System.out.println("pi = " + pi);
270                rf = PolynomialList.<C> castToSolvableList(blas.vectorAdd(PolynomialList.<C> castToList(rf),
271                                pi));
272            }
273            if (it.hasNext() || jt.hasNext()) {
274                logger.error("zeroRelationsArbitrary wrong sizes");
275            }
276            //System.out.println("\nMg Mf = " + rf + "\n");
277            M.add(rf);
278        }
279        //ModuleList<C> ML = new ModuleList<C>( ring, M );
280        //System.out.println("syz ML = " + ML);
281        // debug only:
282        //List<GenSolvablePolynomial<C>> F2 = new ArrayList<GenSolvablePolynomial<C>>( F.size() );
283        /* not true in general
284        List<GenPolynomial<C>> Fp = PolynomialList.<C>castToList(F);
285        for ( List<GenSolvablePolynomial<C>> rr: M ) {
286            GenSolvablePolynomial<C> rrg = PolynomialList.<C>castToSolvableList(blas.scalarProduct(Fp,PolynomialList.<C>castToList(rr)));
287            F2.add( rrg );
288        }
289        PolynomialList<C> pF = new PolynomialList<C>( ring, F );
290        PolynomialList<C> pF2 = new PolynomialList<C>( ring, F2 );
291        if ( ! pF.equals( pF2 ) ) {
292           logger.error("is FAB = F ? false");
293           //System.out.println("pF  = " + pF.list.size());
294           //System.out.println("pF2 = " + pF2.list.size());
295        }
296        */
297        int sflen = sf.size();
298        List<List<GenSolvablePolynomial<C>>> M2;
299        M2 = new ArrayList<List<GenSolvablePolynomial<C>>>(lenf);
300        int i = 0;
301        for (List<GenSolvablePolynomial<C>> ri : M) {
302            List<GenSolvablePolynomial<C>> r2i;
303            r2i = new ArrayList<GenSolvablePolynomial<C>>(ri.size());
304            int j = 0;
305            for (GenSolvablePolynomial<C> rij : ri) {
306                GenSolvablePolynomial<C> p = null;
307                if (i == j) {
308                    p = (GenSolvablePolynomial<C>) ring.getONE().subtract(rij);
309                } else {
310                    if (rij != null) {
311                        p = (GenSolvablePolynomial<C>) rij.negate();
312                    }
313                }
314                r2i.add(p);
315                j++;
316            }
317            M2.add(r2i);
318            if (!blas.isZero(PolynomialList.<C> castToList(r2i))) {
319                sf.add(r2i);
320            }
321            i++;
322        }
323        ModuleList<C> M2L = new ModuleList<C>(ring, M2);
324        if (debug) {
325            logger.debug("syz M2L = " + M2L);
326        }
327
328        if (debug) {
329            ModuleList<C> SF = new ModuleList<C>(ring, sf);
330            logger.debug("syz sf = " + SF);
331            logger.debug("#syz " + sflen + ", " + sf.size());
332        }
333        if (assertEnabled) {
334            logger.info("check4 left syz start");
335            if (!isLeftZeroRelation(sf, F)) {
336                logger.error("is syz sf ? false");
337            }
338            logger.info("check4 left syz end");
339        }
340        return sf;
341    }
342
343
344    /**
345     * Left Ore condition. Generators for the left Ore condition of two solvable
346     * polynomials.
347     * @param a solvable polynomial
348     * @param b solvable polynomial
349     * @return [p,q] with p*a = q*b
350     */
351    @SuppressWarnings({ "cast", "unchecked" })
352    public GenSolvablePolynomial<C>[] leftOreCond(GenSolvablePolynomial<C> a, GenSolvablePolynomial<C> b) {
353        if (a == null || a.isZERO() || b == null || b.isZERO()) {
354            throw new IllegalArgumentException("a and b must be non zero");
355        }
356        GenSolvablePolynomialRing<C> pfac = a.ring;
357        GenSolvablePolynomial<C>[] oc = (GenSolvablePolynomial<C>[]) new GenSolvablePolynomial[2];
358        if (a.equals(b)) {
359            oc[0] = pfac.getONE();
360            oc[1] = pfac.getONE();
361            return oc;
362        }
363        if (a.isConstant()) {
364            if (pfac.coFac.isCommutative()) { // ??
365                oc[0] = b;
366                oc[1] = a;
367                return oc;
368            }
369            oc[1] = pfac.getONE();
370            C c = a.leadingBaseCoefficient().inverse();
371            oc[0] = b.multiply(c);
372            return oc;
373        }
374        if (b.isConstant()) {
375            if (pfac.coFac.isCommutative()) { // ??
376                oc[0] = b;
377                oc[1] = a;
378                return oc;
379            }
380            oc[0] = pfac.getONE();
381            C c = b.leadingBaseCoefficient().inverse();
382            oc[1] = a.multiply(c);
383            return oc;
384        }
385        logger.info("computing left Ore condition: " + a + ", " + b);
386        List<GenSolvablePolynomial<C>> F = new ArrayList<GenSolvablePolynomial<C>>(2);
387        F.add(a);
388        F.add(b);
389        List<List<GenSolvablePolynomial<C>>> Gz = leftZeroRelationsArbitrary(F);
390        /*
391        if (Gz.size() < 0) { // always false
392            //System.out.println("Gz = " + Gz);
393            ModuleList<C> M = new ModuleList<C>(pfac, Gz);
394            ModuleList<C> GM = sbb.leftGB(M);
395            //System.out.println("GM = " + GM);
396            Gz = GM.castToSolvableList();
397        }
398        */
399        List<GenSolvablePolynomial<C>> G1 = null;
400        GenSolvablePolynomial<C> g1 = null;
401        for (List<GenSolvablePolynomial<C>> Gi : Gz) {
402            //System.out.println("Gi = " + Gi);
403            if (Gi.get(0).isZERO()) {
404                continue;
405            }
406            if (G1 == null) {
407                G1 = Gi;
408            }
409            if (g1 == null) {
410                g1 = G1.get(0);
411            } else if (g1.compareTo(Gi.get(0)) > 0) { //g1.degree() > Gi.get(0).degree() 
412                G1 = Gi;
413                g1 = G1.get(0);
414            }
415        }
416        oc[0] = g1; //G1.get(0);
417        oc[1] = (GenSolvablePolynomial<C>) G1.get(1).negate();
418        //logger.info("Ore multiple: " + oc[0].multiply(a) + ", " + Arrays.toString(oc));
419        return oc;
420    }
421
422
423    /**
424     * Right Ore condition. Generators for the right Ore condition of two
425     * solvable polynomials.
426     * @param a solvable polynomial
427     * @param b solvable polynomial
428     * @return [p,q] with a*p = b*q
429     */
430    @SuppressWarnings({ "cast", "unchecked" })
431    public GenSolvablePolynomial<C>[] rightOreCond(GenSolvablePolynomial<C> a, GenSolvablePolynomial<C> b) {
432        if (a == null || a.isZERO() || b == null || b.isZERO()) {
433            throw new IllegalArgumentException("a and b must be non zero");
434        }
435        GenSolvablePolynomialRing<C> pfac = a.ring;
436        GenSolvablePolynomial<C>[] oc = (GenSolvablePolynomial<C>[]) new GenSolvablePolynomial[2];
437        if (a.equals(b)) {
438            oc[0] = pfac.getONE();
439            oc[1] = pfac.getONE();
440            return oc;
441        }
442        if (a.isConstant()) {
443            if (pfac.coFac.isCommutative()) { // ??
444                oc[0] = b;
445                oc[1] = a;
446                return oc;
447            }
448            oc[1] = pfac.getONE();
449            C c = a.leadingBaseCoefficient().inverse();
450            oc[0] = b.multiply(c);
451            return oc;
452        }
453        if (b.isConstant()) {
454            if (pfac.coFac.isCommutative()) { // ??
455                oc[0] = b;
456                oc[1] = a;
457                return oc;
458            }
459            oc[0] = pfac.getONE();
460            C c = b.leadingBaseCoefficient().inverse();
461            oc[1] = a.multiply(c);
462            return oc;
463        }
464        logger.info("computing right Ore condition: " + a + ", " + b);
465        List<GenSolvablePolynomial<C>> F = new ArrayList<GenSolvablePolynomial<C>>(2);
466        F.add(a);
467        F.add(b);
468        List<List<GenSolvablePolynomial<C>>> Gz = rightZeroRelationsArbitrary(F);
469        List<GenSolvablePolynomial<C>> G1 = null;
470        GenSolvablePolynomial<C> g1 = null;
471        for (List<GenSolvablePolynomial<C>> Gi : Gz) {
472            if (Gi.get(0).isZERO()) {
473                continue;
474            }
475            if (G1 == null) {
476                G1 = Gi;
477            }
478            if (g1 == null) {
479                g1 = G1.get(0);
480            } else if (g1.compareTo(Gi.get(0)) > 0) {
481                G1 = Gi;
482                g1 = G1.get(0);
483            }
484        }
485        oc[0] = G1.get(0);
486        oc[1] = (GenSolvablePolynomial<C>) G1.get(1).negate();
487        //logger.info("Ore multiple: " + a.multiply(oc[0]) + ", " + Arrays.toString(oc));
488        return oc;
489    }
490
491
492    /**
493     * Left simplifier. Method of Apel &amp; Lassner (1987).
494     * @param a solvable polynomial
495     * @param b solvable polynomial
496     * @return [p,q] with a/b = p/q and q is minimal and monic
497     */
498    @Override
499    @SuppressWarnings({ "unchecked" })
500    public GenSolvablePolynomial<C>[] leftSimplifier(GenSolvablePolynomial<C> a, GenSolvablePolynomial<C> b) {
501        if (a == null || a.isZERO() || b == null || b.isZERO()) {
502            throw new IllegalArgumentException("a and b must be non zero");
503        }
504        GenSolvablePolynomial<C>[] oc = null;
505        if (a.isConstant() || b.isConstant()) {
506            oc = new GenSolvablePolynomial[] { a, b };
507            return oc;
508        }
509        if (a.totalDegree() > 3 || b.totalDegree() > 3) { // how avoid too long running GBs ?
510            //if (a.totalDegree() + b.totalDegree() > 6) { // how avoid too long running GBs ?
511            // && a.length() < 10 && b.length() < 10
512            logger.warn("skipping GB computation: degs = " + a.totalDegree() + ", " + b.totalDegree());
513            oc = new GenSolvablePolynomial[] { a, b };
514            return oc;
515        }
516        //GenSolvablePolynomialRing<C> pfac = a.ring;
517        oc = rightOreCond(a, b);
518        logger.info("oc = " + Arrays.toString(oc)); // + ", a = " + a + ", b = " + b);
519        List<GenSolvablePolynomial<C>> F = new ArrayList<GenSolvablePolynomial<C>>(oc.length);
520        // opposite order and undo negation
521        F.add((GenSolvablePolynomial<C>) oc[1].negate());
522        F.add(oc[0]);
523        //logger.info("F = " + F);
524        List<List<GenSolvablePolynomial<C>>> Gz = leftZeroRelationsArbitrary(F);
525        //logger.info("Gz: " + Gz);
526        List<GenSolvablePolynomial<C>> G1 = new ArrayList<GenSolvablePolynomial<C>>(Gz.size());
527        List<GenSolvablePolynomial<C>> G2 = new ArrayList<GenSolvablePolynomial<C>>(Gz.size());
528        for (List<GenSolvablePolynomial<C>> ll : Gz) {
529            if (!ll.get(0).isZERO()) { // && !ll.get(1).isZERO()
530                G1.add(ll.get(0)); // denominators
531                G2.add(ll.get(1)); // numerators
532            }
533        }
534        logger.info("G1(den): " + G1 + ", G2(num): " + G2);
535        SolvableExtendedGB<C> exgb = sbb.extLeftGB(G1);
536        logger.info("exgb.F: " + exgb.F + ", exgb.G: " + exgb.G);
537        List<GenSolvablePolynomial<C>> G = exgb.G;
538        int m = 0;
539        GenSolvablePolynomial<C> min = null;
540        for (int i = 0; i < G.size(); i++) {
541            if (min == null) {
542                min = G.get(i);
543                m = i;
544            } else if (min.compareTo(G.get(i)) > 0) {
545                min = G.get(i);
546                m = i;
547            }
548        }
549        //wrong: blas.scalarProduct(G2,exgb.G2F.get(m));
550        GenSolvablePolynomial<C> min2 = (GenSolvablePolynomial<C>) blas.scalarProduct(
551                        PolynomialList.<C> castToList(exgb.G2F.get(m)), PolynomialList.<C> castToList(G2));
552        logger.info("min(den): " + min + ", min(num): " + min2 + ", m = " + m + ", " + exgb.G2F.get(m));
553        // opposite order
554        GenSolvablePolynomial<C> n = min2; // nominator
555        GenSolvablePolynomial<C> d = min; // denominator
556        // normalize
557        if (d.signum() < 0) {
558            n = (GenSolvablePolynomial<C>) n.negate();
559            d = (GenSolvablePolynomial<C>) d.negate();
560        }
561        C lc = d.leadingBaseCoefficient();
562        if (!lc.isONE() && lc.isUnit()) {
563            lc = lc.inverse();
564            n = n.multiplyLeft(lc);
565            d = d.multiplyLeft(lc);
566        }
567        if (debug) {
568            int t = compare(a, b, n, d);
569            if (t != 0) {
570                oc[0] = a; // nominator
571                oc[1] = b; // denominator
572                throw new RuntimeException("simp wrong, giving up: t = " + t);
573                //logger.error("simp wrong, giving up: t = " + t);
574                //return oc;
575            }
576        }
577        oc[0] = n; // nominator
578        oc[1] = d; // denominator
579        return oc;
580    }
581
582}