001/*
002 * $Id$
003 */
004
005package edu.jas.gbufd;
006
007
008import java.util.ArrayList;
009import java.util.List;
010import java.util.Map;
011
012import org.apache.logging.log4j.LogManager;
013import org.apache.logging.log4j.Logger;
014
015import edu.jas.gb.GroebnerBaseAbstract;
016import edu.jas.gb.SolvableGroebnerBaseAbstract;
017import edu.jas.gb.SolvableReductionAbstract;
018import edu.jas.gb.SolvableReductionSeq;
019import edu.jas.gb.WordGroebnerBaseAbstract;
020import edu.jas.gb.WordGroebnerBaseSeq;
021import edu.jas.poly.ExpVector;
022import edu.jas.poly.GenPolynomial;
023import edu.jas.poly.GenPolynomialRing;
024import edu.jas.poly.GenSolvablePolynomial;
025import edu.jas.poly.GenSolvablePolynomialRing;
026import edu.jas.poly.GenWordPolynomial;
027import edu.jas.poly.GenWordPolynomialRing;
028import edu.jas.poly.PolyUtil;
029import edu.jas.structure.GcdRingElem;
030import edu.jas.structure.RingElem;
031
032
033/**
034 * Package gbufd utilities.
035 * @author Heinz Kredel
036 */
037
038public class PolyGBUtil {
039
040
041    private static final Logger logger = LogManager.getLogger(PolyGBUtil.class);
042
043
044    private static final boolean debug = logger.isDebugEnabled();
045
046
047    /**
048     * Test for resultant.
049     * @param A generic polynomial.
050     * @param B generic polynomial.
051     * @param r generic polynomial.
052     * @return true if res(A,B) isContained in ideal(A,B), else false.
053     */
054    public static <C extends GcdRingElem<C>> boolean isResultant(GenPolynomial<C> A, GenPolynomial<C> B,
055                    GenPolynomial<C> r) {
056        if (r == null || r.isZERO()) {
057            return true;
058        }
059        GroebnerBaseAbstract<C> bb = GBFactory.<C> getImplementation(r.ring.coFac);
060        List<GenPolynomial<C>> F = new ArrayList<GenPolynomial<C>>(2);
061        F.add(A);
062        F.add(B);
063        List<GenPolynomial<C>> G = bb.GB(F);
064        //System.out.println("G = " + G);
065        GenPolynomial<C> n = bb.red.normalform(G, r);
066        //System.out.println("n = " + n);
067        return n.isZERO();
068    }
069
070
071    /**
072     * Top pseudo reduction wrt the main variables.
073     * @param P generic polynomial.
074     * @param A list of generic polynomials sorted according to appearing main
075     *            variables.
076     * @return top pseudo remainder of P wrt. A for the appearing variables.
077     */
078    public static <C extends RingElem<C>> GenPolynomial<C> topPseudoRemainder(List<GenPolynomial<C>> A,
079                    GenPolynomial<C> P) {
080        if (A == null || A.isEmpty()) {
081            return P.monic();
082        }
083        if (P.isZERO()) {
084            return P;
085        }
086        //System.out.println("remainder, P = " + P);
087        GenPolynomialRing<C> pfac = A.get(0).ring;
088        if (pfac.nvar <= 1) { // recursion base 
089            GenPolynomial<C> R = PolyUtil.<C> baseSparsePseudoRemainder(P, A.get(0));
090            return R.monic();
091        }
092        // select polynomials according to the main variable
093        GenPolynomialRing<GenPolynomial<C>> rfac = pfac.recursive(1);
094        GenPolynomial<C> Q = A.get(0); // wrong, must eventually search polynomial
095        GenPolynomial<GenPolynomial<C>> qr = PolyUtil.<C> recursive(rfac, Q);
096        GenPolynomial<GenPolynomial<C>> pr = PolyUtil.<C> recursive(rfac, P);
097        GenPolynomial<GenPolynomial<C>> rr;
098        if (qr.isONE()) {
099            return P.ring.getZERO();
100        }
101        if (qr.degree(0) > 0) {
102            rr = PolyUtil.<C> recursiveSparsePseudoRemainder(pr, qr);
103            //System.out.println("remainder, pr = " + pr);
104            //System.out.println("remainder, qr = " + qr);
105            //System.out.println("remainder, rr = " + rr);
106        } else {
107            rr = pr;
108        }
109        if (rr.degree(0) > 0) {
110            GenPolynomial<C> R = PolyUtil.<C> distribute(pfac, rr);
111            return R.monic();
112            // not further reduced wrt. other variables = top-reduction only
113        }
114        List<GenPolynomial<C>> zeroDeg = zeroDegrees(A);
115        GenPolynomial<C> R = topPseudoRemainder(zeroDeg, rr.leadingBaseCoefficient());
116        R = R.extend(pfac, 0, 0L);
117        return R.monic();
118    }
119
120
121    /**
122     * Top coefficient pseudo remainder of the leading coefficient of P wrt A in
123     * the main variables.
124     * @param P generic polynomial in n+1 variables.
125     * @param A list of generic polynomials in n variables sorted according to
126     *            appearing main variables.
127     * @return pseudo remainder of the leading coefficient of P wrt A.
128     */
129    public static <C extends RingElem<C>> GenPolynomial<C> topCoefficientPseudoRemainder(
130                    List<GenPolynomial<C>> A, GenPolynomial<C> P) {
131        if (A == null || A.isEmpty()) {
132            return P.monic();
133        }
134        if (P.isZERO()) {
135            return P;
136        }
137        GenPolynomialRing<C> pfac = P.ring;
138        GenPolynomialRing<C> pfac1 = A.get(0).ring;
139        if (pfac1.nvar <= 1) { // recursion base 
140            GenPolynomial<C> a = A.get(0);
141            GenPolynomialRing<GenPolynomial<C>> rfac = pfac.recursive(pfac.nvar - 1);
142            GenPolynomial<GenPolynomial<C>> pr = PolyUtil.<C> recursive(rfac, P);
143            // ldcf(P,x_m) = q a + r 
144            GenPolynomial<GenPolynomial<C>> rr = PolyGBUtil.<C> coefficientPseudoRemainderBase(pr, a);
145            GenPolynomial<C> R = PolyUtil.<C> distribute(pfac, rr);
146            return R.monic();
147        }
148        // select polynomials according to the main variable
149        GenPolynomialRing<GenPolynomial<C>> rfac1 = pfac1.recursive(1);
150        int nv = pfac.nvar - pfac1.nvar;
151        GenPolynomialRing<GenPolynomial<C>> rfac = pfac.recursive(1 + nv);
152        GenPolynomialRing<GenPolynomial<GenPolynomial<C>>> rfac2 = rfac.recursive(nv);
153        if (debug) {
154            logger.info("rfac = {}", rfac);
155        }
156        GenPolynomial<GenPolynomial<C>> pr = PolyUtil.<C> recursive(rfac, P);
157        GenPolynomial<GenPolynomial<GenPolynomial<C>>> pr2 = PolyUtil.<GenPolynomial<C>> recursive(rfac2, pr);
158        //System.out.println("recursion, pr2 = " + pr2);
159        GenPolynomial<C> Q = A.get(0);
160        GenPolynomial<GenPolynomial<C>> qr = PolyUtil.<C> recursive(rfac1, Q);
161        GenPolynomial<GenPolynomial<GenPolynomial<C>>> rr;
162        if (qr.isONE()) {
163            return P.ring.getZERO();
164        }
165        if (qr.degree(0) > 0) {
166            // pseudo remainder:  ldcf(P,x_m) = a q + r 
167            rr = PolyGBUtil.<C> coefficientPseudoRemainder(pr2, qr);
168            //System.out.println("recursion, qr  = " + qr);
169            //System.out.println("recursion, pr  = " + pr2);
170            //System.out.println("recursion, rr  = " + rr);
171        } else {
172            rr = pr2;
173        }
174        // reduction wrt. the other variables
175        List<GenPolynomial<C>> zeroDeg = zeroDegrees(A);
176        GenPolynomial<GenPolynomial<C>> Rr = PolyUtil.<GenPolynomial<C>> distribute(rfac, rr);
177        GenPolynomial<C> R = PolyUtil.<C> distribute(pfac, Rr);
178        R = topCoefficientPseudoRemainder(zeroDeg, R);
179        return R.monic();
180    }
181
182
183    /**
184     * Polynomial leading coefficient pseudo remainder.
185     * @param P generic polynomial in n+1 variables.
186     * @param A generic polynomial in n variables.
187     * @return pseudo remainder of the leading coefficient of P wrt A, with
188     *         ldcf(A)<sup>m'</sup> P = quotient * A + remainder.
189     */
190    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<GenPolynomial<C>>> coefficientPseudoRemainder(
191                    GenPolynomial<GenPolynomial<GenPolynomial<C>>> P, GenPolynomial<GenPolynomial<C>> A) {
192        if (A == null || A.isZERO()) { // findbugs
193            throw new ArithmeticException(P + " division by zero " + A);
194        }
195        if (A.isONE()) {
196            return P.ring.getZERO();
197        }
198        if (P.isZERO() || P.isONE()) {
199            return P;
200        }
201        GenPolynomialRing<GenPolynomial<GenPolynomial<C>>> pfac = P.ring;
202        GenPolynomialRing<GenPolynomial<C>> afac = A.ring; // == pfac.coFac
203        GenPolynomial<GenPolynomial<GenPolynomial<C>>> r = P;
204        GenPolynomial<GenPolynomial<C>> h;
205        GenPolynomial<GenPolynomial<GenPolynomial<C>>> hr;
206        GenPolynomial<GenPolynomial<C>> ldcf = P.leadingBaseCoefficient();
207        long m = ldcf.degree(0);
208        long n = A.degree(0);
209        GenPolynomial<C> c = A.leadingBaseCoefficient();
210        GenPolynomial<GenPolynomial<C>> cc = afac.getZERO().sum(c);
211        //System.out.println("cc = " + cc);
212        ExpVector e = A.leadingExpVector();
213        for (long i = m; i >= n; i--) {
214            if (r.isZERO()) {
215                return r;
216            }
217            GenPolynomial<GenPolynomial<C>> p = r.leadingBaseCoefficient();
218            ExpVector g = r.leadingExpVector();
219            long k = p.degree(0);
220            if (i == k) {
221                GenPolynomial<C> pl = p.leadingBaseCoefficient();
222                ExpVector f = p.leadingExpVector();
223                f = f.subtract(e);
224                r = r.multiply(cc); // coeff cc
225                h = A.multiply(pl, f); // coeff ac
226                hr = new GenPolynomial<GenPolynomial<GenPolynomial<C>>>(pfac, h, g);
227                r = r.subtract(hr);
228            } else {
229                r = r.multiply(cc);
230            }
231            //System.out.println("r = " + r);
232        }
233        if (r.degree(0) < P.degree(0)) { // recursion for degree
234            r = coefficientPseudoRemainder(r, A);
235        }
236        return r;
237    }
238
239
240    /**
241     * Polynomial leading coefficient pseudo remainder, base case.
242     * @param P generic polynomial in 1+1 variables.
243     * @param A generic polynomial in 1 variable.
244     * @return pseudo remainder of the leading coefficient of P wrt. A, with
245     *         ldcf(A)<sup>m'</sup> P = quotient * A + remainder.
246     */
247    public static <C extends RingElem<C>> GenPolynomial<GenPolynomial<C>> coefficientPseudoRemainderBase(
248                    GenPolynomial<GenPolynomial<C>> P, GenPolynomial<C> A) {
249        if (A == null || A.isZERO()) { // findbugs
250            throw new ArithmeticException(P + " division by zero " + A);
251        }
252        if (A.isONE()) {
253            return P.ring.getZERO();
254        }
255        if (P.isZERO() || P.isONE()) {
256            return P;
257        }
258        GenPolynomialRing<GenPolynomial<C>> pfac = P.ring;
259        GenPolynomialRing<C> afac = A.ring; // == pfac.coFac
260        GenPolynomial<GenPolynomial<C>> r = P;
261        GenPolynomial<C> h;
262        GenPolynomial<GenPolynomial<C>> hr;
263        GenPolynomial<C> ldcf = P.leadingBaseCoefficient();
264        long m = ldcf.degree(0);
265        long n = A.degree(0);
266        C c = A.leadingBaseCoefficient();
267        GenPolynomial<C> cc = afac.getZERO().sum(c);
268        //System.out.println("cc = " + cc);
269        ExpVector e = A.leadingExpVector();
270        for (long i = m; i >= n; i--) {
271            if (r.isZERO()) {
272                return r;
273            }
274            GenPolynomial<C> p = r.leadingBaseCoefficient();
275            ExpVector g = r.leadingExpVector();
276            long k = p.degree(0);
277            if (i == k) {
278                C pl = p.leadingBaseCoefficient();
279                ExpVector f = p.leadingExpVector();
280                f = f.subtract(e);
281                r = r.multiply(cc); // coeff cc
282                h = A.multiply(pl, f); // coeff ac
283                hr = new GenPolynomial<GenPolynomial<C>>(pfac, h, g);
284                r = r.subtract(hr);
285            } else {
286                r = r.multiply(cc);
287            }
288            //System.out.println("r = " + r);
289        }
290        if (r.degree(0) < P.degree(0)) { // recursion for degree
291            r = coefficientPseudoRemainderBase(r, A);
292        }
293        return r;
294    }
295
296
297    /**
298     * Extract polynomials with degree zero in the main variable.
299     * @param A list of generic polynomials in n variables.
300     * @return Z = [a_i] with deg(a_i,x_n) = 0 and in n-1 variables.
301     */
302    public static <C extends RingElem<C>> List<GenPolynomial<C>> zeroDegrees(List<GenPolynomial<C>> A) {
303        if (A == null || A.isEmpty()) {
304            return A;
305        }
306        GenPolynomialRing<C> pfac = A.get(0).ring;
307        GenPolynomialRing<GenPolynomial<C>> rfac = pfac.recursive(1);
308        List<GenPolynomial<C>> zeroDeg = new ArrayList<GenPolynomial<C>>(A.size());
309        for (int i = 0; i < A.size(); i++) {
310            GenPolynomial<C> q = A.get(i);
311            GenPolynomial<GenPolynomial<C>> fr = PolyUtil.<C> recursive(rfac, q);
312            if (fr.degree(0) == 0) {
313                zeroDeg.add(fr.leadingBaseCoefficient());
314            }
315        }
316        return zeroDeg;
317    }
318
319
320    /**
321     * Intersection. Generators for the intersection of ideals.
322     * @param pfac polynomial ring
323     * @param A list of polynomials
324     * @param B list of polynomials
325     * @return generators for (A \cap B)
326     */
327    public static <C extends GcdRingElem<C>> List<GenPolynomial<C>> intersect(GenPolynomialRing<C> pfac,
328                    List<GenPolynomial<C>> A, List<GenPolynomial<C>> B) {
329        if (A == null || A.isEmpty()) { // (0)
330            return B;
331        }
332        if (B == null || B.isEmpty()) { // (0)
333            return A;
334        }
335        int s = A.size() + B.size();
336        List<GenPolynomial<C>> c = new ArrayList<GenPolynomial<C>>(s);
337        GenPolynomialRing<C> tfac = pfac.extend(1);
338        // term order is also adjusted
339        for (GenPolynomial<C> p : A) {
340            p = p.extend(tfac, 0, 1L); // t*p
341            c.add(p);
342        }
343        for (GenPolynomial<C> p : B) {
344            GenPolynomial<C> q = p.extend(tfac, 0, 1L);
345            GenPolynomial<C> r = p.extend(tfac, 0, 0L);
346            p = r.subtract(q); // (1-t)*p
347            c.add(p);
348        }
349        GroebnerBaseAbstract<C> bb = GBFactory.<C> getImplementation(tfac.coFac);
350        logger.warn("intersect computing GB");
351        List<GenPolynomial<C>> G = bb.GB(c);
352        if (debug) {
353            logger.debug("intersect GB = {}", G);
354        }
355        List<GenPolynomial<C>> I = PolyUtil.<C> intersect(pfac, G);
356        return I;
357    }
358
359
360    /**
361     * Intersection. Generators for the intersection of ideals.
362     * @param pfac solvable polynomial ring
363     * @param A list of polynomials
364     * @param B list of polynomials
365     * @return generators for (A \cap B)
366     */
367    public static <C extends GcdRingElem<C>> List<GenSolvablePolynomial<C>> intersect(
368                    GenSolvablePolynomialRing<C> pfac, List<GenSolvablePolynomial<C>> A,
369                    List<GenSolvablePolynomial<C>> B) {
370        if (A == null || A.isEmpty()) { // (0)
371            return B;
372        }
373        if (B == null || B.isEmpty()) { // (0)
374            return A;
375        }
376        int s = A.size() + B.size();
377        List<GenSolvablePolynomial<C>> c = new ArrayList<GenSolvablePolynomial<C>>(s);
378        GenSolvablePolynomialRing<C> tfac = pfac.extend(1);
379        // term order is also adjusted
380        for (GenSolvablePolynomial<C> p : A) {
381            p = (GenSolvablePolynomial<C>) p.extend(tfac, 0, 1L); // t*p
382            c.add(p);
383        }
384        for (GenSolvablePolynomial<C> p : B) {
385            GenSolvablePolynomial<C> q = (GenSolvablePolynomial<C>) p.extend(tfac, 0, 1L);
386            GenSolvablePolynomial<C> r = (GenSolvablePolynomial<C>) p.extend(tfac, 0, 0L);
387            p = (GenSolvablePolynomial<C>) r.subtract(q); // (1-t)*p
388            c.add(p);
389        }
390        SolvableGroebnerBaseAbstract<C> sbb = SGBFactory.<C> getImplementation(tfac.coFac);
391        //new SolvableGroebnerBaseSeq<C>();
392        logger.warn("intersect computing GB");
393        List<GenSolvablePolynomial<C>> g = sbb.leftGB(c);
394        //List<GenSolvablePolynomial<C>> g = sbb.twosidedGB(c);
395        if (debug) {
396            logger.debug("intersect GB = {}", g);
397        }
398        List<GenSolvablePolynomial<C>> I = PolyUtil.<C> intersect(pfac, g);
399        return I;
400    }
401
402
403    /**
404     * Intersection. Generators for the intersection of word ideals.
405     * @param pfac word polynomial ring
406     * @param A list of word polynomials
407     * @param B list of word polynomials
408     * @return generators for (A \cap B) if it exists
409     */
410    public static <C extends GcdRingElem<C>> List<GenWordPolynomial<C>> intersect(
411                    GenWordPolynomialRing<C> pfac, List<GenWordPolynomial<C>> A,
412                    List<GenWordPolynomial<C>> B) {
413        logger.warn("new WordGroebnerBaseSeq generated");
414        return intersect(pfac, A, B, new WordGroebnerBaseSeq<C>());
415    }
416
417
418    /**
419     * Intersection. Generators for the intersection of word ideals.
420     * @param pfac word polynomial ring
421     * @param A list of word polynomials
422     * @param B list of word polynomials
423     * @param bb Groebner Base engine
424     * @return generators for (A \cap B) if it exists
425     */
426    public static <C extends GcdRingElem<C>> List<GenWordPolynomial<C>> intersect(
427                    GenWordPolynomialRing<C> pfac, List<GenWordPolynomial<C>> A,
428                    List<GenWordPolynomial<C>> B, WordGroebnerBaseAbstract<C> bb) {
429        if (A == null || A.isEmpty()) { // (0)
430            return B;
431        }
432        if (B == null || B.isEmpty()) { // (0)
433            return A;
434        }
435        int s = A.size() + B.size();
436        List<GenWordPolynomial<C>> L = new ArrayList<GenWordPolynomial<C>>(s);
437        GenWordPolynomialRing<C> tfac = pfac.extend(1);
438        List<GenWordPolynomial<C>> gens = tfac.univariateList();
439        //System.out.println("gens = " + gens);
440        GenWordPolynomial<C> t = gens.get(gens.size() - 1);
441        //System.out.println("t = " + t);
442        // make t commute with other variables
443        for (GenWordPolynomial<C> p : gens) {
444            if (t == p) {
445                continue;
446            }
447            GenWordPolynomial<C> c = t.multiply(p).subtract(p.multiply(t)); // t p - p t
448            L.add(c);
449        }
450        for (GenWordPolynomial<C> p : A) {
451            p = tfac.valueOf(p).multiply(t); // t p
452            L.add(p);
453        }
454        for (GenWordPolynomial<C> p : B) {
455            GenWordPolynomial<C> q = tfac.valueOf(p).multiply(t);
456            GenWordPolynomial<C> r = tfac.valueOf(p);
457            p = r.subtract(q); // (1-t) p
458            L.add(p);
459        }
460        //System.out.println("L = " + L);
461        //WordGroebnerBaseAbstract<C> bb = new WordGroebnerBaseSeq<C>(); // todo timestatus
462        logger.warn("intersect computing GB");
463        List<GenWordPolynomial<C>> G = bb.GB(L);
464        //System.out.println("G = " + G);
465        if (debug) {
466            logger.debug("intersect GB = {}", G);
467        }
468        List<GenWordPolynomial<C>> I = PolyUtil.<C> intersect(pfac, G);
469        return I;
470    }
471
472
473    /**
474     * Solvable quotient and remainder via reduction.
475     * @param n first solvable polynomial.
476     * @param d second solvable polynomial.
477     * @return [ n/d, n - (n/d)*d ]
478     */
479    @SuppressWarnings({"cast","unchecked"})
480    public static <C extends GcdRingElem<C>> GenSolvablePolynomial<C>[] quotientRemainder(
481                    GenSolvablePolynomial<C> n, GenSolvablePolynomial<C> d) {
482        GenSolvablePolynomial<C>[] res = (GenSolvablePolynomial<C>[]) new GenSolvablePolynomial[2];
483        if (d.isZERO()) {
484            throw new RuntimeException("division by zero: " + n + "/" + d);
485        }
486        if (n.isZERO()) {
487            res[0] = n;
488            res[1] = n;
489            return res;
490        }
491        GenSolvablePolynomialRing<C> r = n.ring;
492        if (d.isONE()) {
493            res[0] = n;
494            res[1] = r.getZERO();
495            return res;
496        }
497        // divide
498        List<GenSolvablePolynomial<C>> Q = new ArrayList<GenSolvablePolynomial<C>>(1);
499        Q.add(r.getZERO());
500        List<GenSolvablePolynomial<C>> D = new ArrayList<GenSolvablePolynomial<C>>(1);
501        D.add(d);
502        SolvableReductionAbstract<C> sred = new SolvableReductionSeq<C>();
503        res[1] = sred.rightNormalform(Q, D, n); // left
504        res[0] = Q.get(0);
505        return res;
506    }
507
508
509    /**
510     * Subring generators.
511     * @param A list of polynomials in n variables.
512     * @return a Groebner base of polynomials in m &gt; n variables generating
513     *         the subring of K[A].
514     */
515    public static <C extends GcdRingElem<C>> List<GenPolynomial<C>> subRing(List<GenPolynomial<C>> A) {
516        if (A == null || A.isEmpty()) {
517            return A;
518        }
519        GenPolynomialRing<C> pfac = A.get(0).ring;
520        logger.debug("pfac = {}", pfac); //.toScript()
521        int n = pfac.nvar;
522        List<GenPolynomial<C>> Ap = new ArrayList<GenPolynomial<C>>();
523        for (GenPolynomial<C> a : A) {
524            if (a == null || a.isZERO() || a.isONE()) {
525                continue;
526            }
527            Ap.add(a);
528        }
529        int k = Ap.size();
530        if (k == 0) {
531            return Ap;
532        }
533        GenPolynomialRing<C> rfac = pfac.extendLower(k);
534        logger.debug("rfac = {}", rfac); //.toScript()
535        assert rfac.nvar == n + k : "rfac.nvar == n+k";
536        List<GenPolynomial<C>> sr = new ArrayList<GenPolynomial<C>>();
537        int i = 0;
538        for (GenPolynomial<C> a : Ap) {
539            GenPolynomial<C> b = a.extendLower(rfac, 0, 0L);
540            b = b.subtract(pfac.getONE().extendLower(rfac, i++, 1L));
541            //System.out.println("a = " + a);
542            //System.out.println("b = " + b);
543            sr.add(b);
544        }
545        GroebnerBaseAbstract<C> bb = GBFactory.<C> getImplementation(pfac.coFac);
546        List<GenPolynomial<C>> srg = bb.GB(sr);
547        return srg;
548    }
549
550
551    /**
552     * Subring membership.
553     * @param A Groebner base of polynomials in m &gt; n variables generating
554     *            the subring of elements of K[A].
555     * @param g polynomial in n variables.
556     * @return true, if g \in K[A], else false.
557     */
558    public static <C extends GcdRingElem<C>> boolean subRingMember(List<GenPolynomial<C>> A,
559                    GenPolynomial<C> g) {
560        if (A == null || A.isEmpty()) {
561            return true;
562        }
563        GenPolynomialRing<C> pfac = A.get(0).ring;
564        GenPolynomial<C> m = g;
565        if (pfac.nvar != g.ring.nvar) {
566            m = m.extendLower(pfac, 0, 0L);
567        } else {
568            throw new IllegalArgumentException("g must be extended: " + pfac.nvar + " == " + g.ring.nvar
569                            + " did you mean method subRingAndMember()?");
570        }
571        //ReductionAbstract<C> red = new ReductionSeq<C>();
572        GroebnerBaseAbstract<C> bb = GBFactory.<C> getImplementation(pfac.coFac);
573        GenPolynomial<C> r = bb.red.normalform(A, m);
574        //System.out.println("r = " + r);
575        GenPolynomialRing<C> cfac = pfac.contract(g.ring.nvar);
576        logger.debug("cfac = {}", cfac); //.toScript()
577        Map<ExpVector, GenPolynomial<C>> map = r.contract(cfac);
578        //System.out.println("map = " + map);
579        boolean t = map.size() == 1 && map.keySet().contains(g.ring.evzero);
580        if (!t) {
581            System.out.println("false: map = " + map);
582        }
583        return t;
584    }
585
586
587    /**
588     * Subring and membership test.
589     * @param A list of polynomials in n variables.
590     * @param g polynomial in n variables.
591     * @return true, if g \in K[A], else false.
592     */
593    public static <C extends GcdRingElem<C>> boolean subRingAndMember(List<GenPolynomial<C>> A,
594                    GenPolynomial<C> g) {
595        if (A == null || A.isEmpty()) {
596            return true;
597        }
598        List<GenPolynomial<C>> srg = PolyGBUtil.<C> subRing(A);
599        return PolyGBUtil.<C> subRingMember(srg, g);
600    }
601
602
603    /**
604     * Chinese remainder theorem.
605     * @param F = ( F_i ) list of list of polynomials in n variables.
606     * @param A = ( f_i ) list of polynomials in n variables.
607     * @return p \in \Cap_i (f_i + ideal(F_i)) if it exists, else null.
608     */
609    public static <C extends GcdRingElem<C>> GenPolynomial<C> chineseRemainderTheorem(
610                    List<List<GenPolynomial<C>>> F, List<GenPolynomial<C>> A) {
611        if (F == null || F.isEmpty() || A == null || A.isEmpty()) {
612            throw new IllegalArgumentException("F and A may not be empty or null");
613        }
614        int m = F.size();
615        if (m != A.size()) {
616            throw new IllegalArgumentException("size(F) and size(A) must be equal");
617        }
618        GenPolynomialRing<C> pfac = A.get(0).ring;
619        logger.debug("pfac = {}", pfac); //.toScript()
620        GenPolynomialRing<C> rfac = pfac.extend(m);
621        logger.debug("rfac = {}", rfac); //.toScript()
622        GenPolynomial<C> y = rfac.getONE();
623        GenPolynomial<C> f = rfac.getZERO();
624        int i = 0;
625        List<GenPolynomial<C>> Fp = new ArrayList<GenPolynomial<C>>(); //m**2?
626        //System.out.println("A = " + A);
627        for (List<GenPolynomial<C>> Fi : F) {
628            GenPolynomial<C> Yi = pfac.getONE().extend(rfac, i, 1L);
629            y = y.subtract(Yi);
630            List<GenPolynomial<C>> Fip = new ArrayList<GenPolynomial<C>>(Fi.size());
631            //System.out.println("Fi = " + Fi);
632            for (GenPolynomial<C> a : Fi) {
633                GenPolynomial<C> b = a.extend(rfac, i, 1L);
634                Fip.add(b);
635            }
636            //System.out.println("Fip = " + Fip);
637            Fp.addAll(Fip);
638            GenPolynomial<C> a = A.get(i);
639            //System.out.println("a = " + a);
640            f = f.sum(a.extend(rfac, i, 1L));
641            i++;
642        }
643        Fp.add(y);
644        //System.out.println("f = " + f);
645        //System.out.println("Fp = " + Fp);
646        GroebnerBaseAbstract<C> bb = GBFactory.<C> getImplementation(pfac.coFac);
647        List<GenPolynomial<C>> Fpp = bb.GB(Fp);
648        //System.out.println("Fpp = " + Fpp);
649        GenPolynomial<C> h = bb.red.normalform(Fpp, f);
650        //System.out.println("h = " + h);
651        ////PseudoReduction<C> pr = new PseudoReductionSeq<C>();
652        ////PseudoReductionEntry<C> fz = pr.normalformFactor(Fpp, f);
653        ////System.out.println("fz = " + fz);
654        List<GenPolynomial<C>> H = new ArrayList<GenPolynomial<C>>();
655        H.add(h);
656        H = PolyUtil.<C> intersect(pfac, H);
657        //System.out.println("H != (): " + (! H.isEmpty()));
658        if (H.isEmpty()) {
659            return null;
660        }
661        return H.get(0);
662    }
663
664
665    /**
666     * Is Chinese remainder.
667     * @param F = ( F_i ) list of list of polynomials in n variables.
668     * @param A = ( f_i ) list of polynomials in n variables.
669     * @param h polynomial in n variables.
670     * @return true if h \in \Cap_i (f_i + ideal(F_i)), else false.
671     */
672    public static <C extends GcdRingElem<C>> boolean isChineseRemainder(List<List<GenPolynomial<C>>> F,
673                    List<GenPolynomial<C>> A, GenPolynomial<C> h) {
674        if (h == null) {
675            return false;
676        }
677        if (F == null || F.isEmpty() || A == null || A.isEmpty()) {
678            return false;
679        }
680        if (F.size() != A.size()) {
681            return false;
682        }
683        GenPolynomialRing<C> pfac = h.ring;
684        if (!pfac.coFac.isField()) {
685            //System.out.println("pfac = " + pfac.toScript());
686            logger.error("only for field coefficients: {}", pfac); //.toScript()
687            //throw new IllegalArgumentException("only for field coefficients: " + pfac.toScript());
688        }
689        GroebnerBaseAbstract<C> bb = GBFactory.<C> getImplementation(pfac.coFac);
690        int i = 0;
691        for (List<GenPolynomial<C>> Fi : F) {
692            List<GenPolynomial<C>> Fp = bb.GB(Fi);
693            //System.out.println("Fp = " + Fp);
694            GenPolynomial<C> a = A.get(i);
695            GenPolynomial<C> fi = bb.red.normalform(Fp, h.subtract(a));
696            if (!fi.isZERO()) {
697                //System.out.println("Fp = " + Fp + ", Fi = " + Fi);
698                //System.out.println("h  = " + h  + ", a  = " + a  + ", fi  = " + fi);
699                logger.info("h-a = {}, fi  = {}", h.subtract(a), fi);
700                return false;
701            }
702            i++;
703        }
704        return true;
705    }
706
707
708    /**
709     * Chinese remainder theorem, interpolation.
710     * @param fac polynomial ring over K in n variables.
711     * @param E = ( E_i ), E_i = ( e_ij ) list of list of elements of K, the evaluation points.
712     * @param V = ( f_i ) list of elements of K, the evaluation values.
713     * @return p \in K[X1,...,Xn], with p(E_i) = f_i, if it exists, else null.
714     */
715    @SuppressWarnings({"cast","unchecked"})
716    public static <C extends GcdRingElem<C>> GenPolynomial<C> CRTInterpolation(
717                  GenPolynomialRing<C> fac, List<List<C>> E, List<C> V) {
718        if (E == null || E.isEmpty() || V == null || V.isEmpty()) {
719            throw new IllegalArgumentException("E and V may not be empty or null");
720        }
721        int m = E.size();
722        if (m != V.size()) {
723            throw new IllegalArgumentException("size(E) and size(V) must be equal");
724        }
725        //System.out.println("fac = " + fac.toScript());
726        List<List<GenPolynomial<C>>> F = new ArrayList<List<GenPolynomial<C>>>(E.size());
727        List<GenPolynomial<C>> A = new ArrayList<GenPolynomial<C>>(V.size());
728        List<GenPolynomial<C>> gen = (List<GenPolynomial<C>>) fac.univariateList();
729        //System.out.println("gen = " + gen);
730        int i = 0;
731        for (List<C> Ei : E) {
732            List<GenPolynomial<C>> Fi = new ArrayList<GenPolynomial<C>>();
733            int j = 0;
734            for (C eij : Ei) {
735                GenPolynomial<C> ep = gen.get(j);
736                ep = ep.subtract( fac.valueOf(eij) );
737                Fi.add(ep);
738                j++;
739            }
740            F.add(Fi);
741            String ai = " " + V.get(i); // (sic) .toString() not possible
742            //System.out.println("ai = " + ai);
743            GenPolynomial<C> ap = fac.valueOf(ai);
744            A.add(ap);
745            i++;
746        }
747        //System.out.println("F = " + F);
748        //System.out.println("A = " + A);
749        GenPolynomial<C> p = PolyGBUtil. <C>chineseRemainderTheorem(F,A);
750        //System.out.println("p = " + p);
751        //System.out.println("t = " + PolyGBUtil. <C>isChineseRemainder(F,A,p));
752        return p;
753    }
754
755}