001 /*
002 * $Id: Ideal.java 3671 2011-06-25 21:11:34Z kredel $
003 */
004
005 package edu.jas.application;
006
007
008 import java.io.Serializable;
009 import java.util.ArrayList;
010 import java.util.Arrays;
011 import java.util.HashSet;
012 import java.util.Iterator;
013 import java.util.List;
014 import java.util.Map;
015 import java.util.Set;
016 import java.util.SortedMap;
017 import java.util.TreeMap;
018
019 import org.apache.log4j.Logger;
020
021 import edu.jas.gb.ExtendedGB;
022 import edu.jas.gb.GroebnerBaseAbstract;
023 import edu.jas.gb.Reduction;
024 import edu.jas.gbufd.GBFactory;
025 import edu.jas.gbufd.GroebnerBasePartial;
026 import edu.jas.poly.AlgebraicNumber;
027 import edu.jas.poly.AlgebraicNumberRing;
028 import edu.jas.poly.ExpVector;
029 import edu.jas.poly.GenPolynomial;
030 import edu.jas.poly.GenPolynomialRing;
031 import edu.jas.poly.OptimizedPolynomialList;
032 import edu.jas.poly.PolyUtil;
033 import edu.jas.poly.PolynomialList;
034 import edu.jas.poly.TermOrder;
035 import edu.jas.poly.TermOrderOptimization;
036 import edu.jas.structure.GcdRingElem;
037 import edu.jas.structure.NotInvertibleException;
038 import edu.jas.structure.Power;
039 import edu.jas.structure.RingFactory;
040 import edu.jas.ufd.FactorAbstract;
041 //import edu.jas.ufd.FactorFactory;
042 import edu.jas.ufd.GCDFactory;
043 import edu.jas.ufd.GreatestCommonDivisor;
044 import edu.jas.ufd.PolyUfdUtil;
045 import edu.jas.ufd.Quotient;
046 import edu.jas.ufd.QuotientRing;
047 import edu.jas.ufd.SquarefreeAbstract;
048 import edu.jas.ufd.SquarefreeFactory;
049
050
051 /**
052 * Ideal implements some methods for ideal arithmetic, for example intersection,
053 * quotient and zero and positive dimensional ideal decomposition.
054 * @author Heinz Kredel
055 */
056 public class Ideal<C extends GcdRingElem<C>> implements Comparable<Ideal<C>>, Serializable, Cloneable {
057
058
059 private static final Logger logger = Logger.getLogger(Ideal.class);
060
061
062 private final boolean debug = true || logger.isDebugEnabled();
063
064
065 /**
066 * The data structure is a PolynomialList.
067 */
068 protected PolynomialList<C> list;
069
070
071 /**
072 * Indicator if list is a Groebner Base.
073 */
074 protected boolean isGB;
075
076
077 /**
078 * Indicator if test has been performed if this is a Groebner Base.
079 */
080 protected boolean testGB;
081
082
083 /**
084 * Indicator if list has optimized term order.
085 */
086 protected boolean isTopt;
087
088
089 /**
090 * Groebner base engine.
091 */
092 protected final GroebnerBaseAbstract<C> bb;
093
094
095 /**
096 * Reduction engine.
097 */
098 protected final Reduction<C> red;
099
100
101 /**
102 * Squarefree decomposition engine.
103 */
104 protected final SquarefreeAbstract<C> engine;
105
106
107 /**
108 * Constructor.
109 * @param ring polynomial ring
110 */
111 public Ideal(GenPolynomialRing<C> ring) {
112 this(ring, new ArrayList<GenPolynomial<C>>());
113 }
114
115
116 /**
117 * Constructor.
118 * @param ring polynomial ring
119 * @param F list of polynomials
120 */
121 public Ideal(GenPolynomialRing<C> ring, List<GenPolynomial<C>> F) {
122 this(new PolynomialList<C>(ring, F));
123 }
124
125
126 /**
127 * Constructor.
128 * @param ring polynomial ring
129 * @param F list of polynomials
130 * @param gb true if F is known to be a Groebner Base, else false
131 */
132 public Ideal(GenPolynomialRing<C> ring, List<GenPolynomial<C>> F, boolean gb) {
133 this(new PolynomialList<C>(ring, F), gb);
134 }
135
136
137 /**
138 * Constructor.
139 * @param ring polynomial ring
140 * @param F list of polynomials
141 * @param gb true if F is known to be a Groebner Base, else false
142 * @param topt true if term order is optimized, else false
143 */
144 public Ideal(GenPolynomialRing<C> ring, List<GenPolynomial<C>> F, boolean gb, boolean topt) {
145 this(new PolynomialList<C>(ring, F), gb, topt);
146 }
147
148
149 /**
150 * Constructor.
151 * @param list polynomial list
152 */
153 public Ideal(PolynomialList<C> list) {
154 this(list, false);
155 }
156
157
158 /**
159 * Constructor.
160 * @param list polynomial list
161 * @param bb Groebner Base engine
162 * @param red Reduction engine
163 */
164 public Ideal(PolynomialList<C> list, GroebnerBaseAbstract<C> bb, Reduction<C> red) {
165 this(list, false, bb, red);
166 }
167
168
169 /**
170 * Constructor.
171 * @param list polynomial list
172 * @param gb true if list is known to be a Groebner Base, else false
173 */
174 public Ideal(PolynomialList<C> list, boolean gb) {
175 //this(list, gb, new GroebnerBaseSeqPairSeq<C>(), new ReductionSeq<C>());
176 this(list, gb, GBFactory.getImplementation(list.ring.coFac));
177 }
178
179
180 /**
181 * Constructor.
182 * @param list polynomial list
183 * @param gb true if list is known to be a Groebner Base, else false
184 * @param topt true if term order is optimized, else false
185 */
186 public Ideal(PolynomialList<C> list, boolean gb, boolean topt) {
187 //this(list, gb, topt, new GroebnerBaseSeqPairSeq<C>(), new ReductionSeq<C>());
188 this(list, gb, topt, GBFactory.getImplementation(list.ring.coFac));
189 }
190
191
192 /**
193 * Constructor.
194 * @param list polynomial list
195 * @param gb true if list is known to be a Groebner Base, else false
196 * @param bb Groebner Base engine
197 * @param red Reduction engine
198 */
199 public Ideal(PolynomialList<C> list, boolean gb, GroebnerBaseAbstract<C> bb, Reduction<C> red) {
200 this(list, gb, false, bb, red);
201 }
202
203
204 /**
205 * Constructor.
206 * @param list polynomial list
207 * @param gb true if list is known to be a Groebner Base, else false
208 * @param bb Groebner Base engine
209 */
210 public Ideal(PolynomialList<C> list, boolean gb, GroebnerBaseAbstract<C> bb) {
211 this(list, gb, false, bb, bb.red);
212 }
213
214
215 /**
216 * Constructor.
217 * @param list polynomial list
218 * @param gb true if list is known to be a Groebner Base, else false
219 * @param topt true if term order is optimized, else false
220 * @param bb Groebner Base engine
221 */
222 public Ideal(PolynomialList<C> list, boolean gb, boolean topt, GroebnerBaseAbstract<C> bb) {
223 this(list, gb, topt, bb, bb.red);
224 }
225
226
227 /**
228 * Constructor.
229 * @param list polynomial list
230 * @param gb true if list is known to be a Groebner Base, else false
231 * @param topt true if term order is optimized, else false
232 * @param bb Groebner Base engine
233 * @param red Reduction engine
234 */
235 public Ideal(PolynomialList<C> list, boolean gb, boolean topt, GroebnerBaseAbstract<C> bb,
236 Reduction<C> red) {
237 if (list == null || list.list == null) {
238 throw new IllegalArgumentException("list and list.list may not be null");
239 }
240 this.list = list;
241 this.isGB = gb;
242 this.isTopt = topt;
243 this.testGB = (gb ? true : false); // ??
244 this.bb = bb;
245 this.red = red;
246 this.engine = SquarefreeFactory.<C> getImplementation(list.ring.coFac);
247 }
248
249
250 /**
251 * Clone this.
252 * @return a copy of this.
253 */
254 @Override
255 public Ideal<C> clone() {
256 return new Ideal<C>(list.clone(), isGB, isTopt, bb, red);
257 }
258
259
260 /**
261 * Get the List of GenPolynomials.
262 * @return list.list
263 */
264 public List<GenPolynomial<C>> getList() {
265 return list.list;
266 }
267
268
269 /**
270 * Get the GenPolynomialRing.
271 * @return list.ring
272 */
273 public GenPolynomialRing<C> getRing() {
274 return list.ring;
275 }
276
277
278 /**
279 * Get the zero ideal.
280 * @return ideal(0)
281 */
282 public Ideal<C> getZERO() {
283 List<GenPolynomial<C>> z = new ArrayList<GenPolynomial<C>>(0);
284 PolynomialList<C> pl = new PolynomialList<C>(getRing(), z);
285 return new Ideal<C>(pl, true, isTopt, bb, red);
286 }
287
288
289 /**
290 * Get the one ideal.
291 * @return ideal(1)
292 */
293 public Ideal<C> getONE() {
294 List<GenPolynomial<C>> one = new ArrayList<GenPolynomial<C>>(1);
295 one.add(list.ring.getONE());
296 PolynomialList<C> pl = new PolynomialList<C>(getRing(), one);
297 return new Ideal<C>(pl, true, isTopt, bb, red);
298 }
299
300
301 /**
302 * String representation of the ideal.
303 * @see java.lang.Object#toString()
304 */
305 @Override
306 public String toString() {
307 return list.toString();
308 }
309
310
311 /**
312 * Get a scripting compatible string representation.
313 * @return script compatible representation for this Element.
314 * @see edu.jas.structure.Element#toScript()
315 */
316 public String toScript() {
317 // Python case
318 return list.toScript();
319 }
320
321
322 /**
323 * Comparison with any other object. Note: If both ideals are not Groebner
324 * Bases, then false may be returned even the ideals are equal.
325 * @see java.lang.Object#equals(java.lang.Object)
326 */
327 @Override
328 @SuppressWarnings("unchecked")
329 public boolean equals(Object b) {
330 if (!(b instanceof Ideal)) {
331 logger.warn("equals no Ideal");
332 return false;
333 }
334 Ideal<C> B = null;
335 try {
336 B = (Ideal<C>) b;
337 } catch (ClassCastException ignored) {
338 return false;
339 }
340 //if ( isGB && B.isGB ) {
341 // return list.equals( B.list ); requires also monic polys
342 //} else { // compute GBs ?
343 return this.contains(B) && B.contains(this);
344 //}
345 }
346
347
348 /**
349 * Ideal list comparison.
350 * @param L other Ideal.
351 * @return compareTo() of polynomial lists.
352 */
353 public int compareTo(Ideal<C> L) {
354 return list.compareTo(L.list);
355 }
356
357
358 /**
359 * Hash code for this ideal.
360 * @see java.lang.Object#hashCode()
361 */
362 @Override
363 public int hashCode() {
364 int h;
365 h = list.hashCode();
366 if (isGB) {
367 h = h << 1;
368 }
369 if (testGB) {
370 h += 1;
371 }
372 return h;
373 }
374
375
376 /**
377 * Test if ZERO ideal.
378 * @return true, if this is the 0 ideal, else false
379 */
380 public boolean isZERO() {
381 return list.isZERO();
382 }
383
384
385 /**
386 * Test if ONE is contained in the ideal. To test for a proper ideal use
387 * <code>! id.isONE()</code>.
388 * @return true, if this is the 1 ideal, else false
389 */
390 public boolean isONE() {
391 return list.isONE();
392 }
393
394
395 /**
396 * Optimize the term order.
397 */
398 public void doToptimize() {
399 if (isTopt) {
400 return;
401 }
402 list = TermOrderOptimization.<C> optimizeTermOrder(list);
403 isTopt = true;
404 if (isGB) {
405 isGB = false;
406 doGB();
407 }
408 return;
409 }
410
411
412 /**
413 * Test if this is a Groebner base.
414 * @return true, if this is a Groebner base, else false
415 */
416 public boolean isGB() {
417 if (testGB) {
418 return isGB;
419 }
420 logger.warn("isGB computing");
421 isGB = bb.isGB(getList());
422 testGB = true;
423 return isGB;
424 }
425
426
427 /**
428 * Do Groebner Base. compute the Groebner Base for this ideal.
429 */
430 public void doGB() {
431 if (isGB && testGB) {
432 return;
433 }
434 //logger.warn("GB computing");
435 List<GenPolynomial<C>> G = getList();
436 logger.info("GB computing = " + G);
437 G = bb.GB(G);
438 if (isTopt) {
439 List<Integer> perm = ((OptimizedPolynomialList<C>) list).perm;
440 list = new OptimizedPolynomialList<C>(perm, getRing(), G);
441 } else {
442 list = new PolynomialList<C>(getRing(), G);
443 }
444 isGB = true;
445 testGB = true;
446 return;
447 }
448
449
450 /**
451 * Groebner Base. Get a Groebner Base for this ideal.
452 * @return GB(this)
453 */
454 public Ideal<C> GB() {
455 if (isGB) {
456 return this;
457 }
458 doGB();
459 return this;
460 }
461
462
463 /**
464 * Ideal containment. Test if B is contained in this ideal. Note: this is
465 * eventually modified to become a Groebner Base.
466 * @param B ideal
467 * @return true, if B is contained in this, else false
468 */
469 public boolean contains(Ideal<C> B) {
470 if (B == null || B.isZERO()) {
471 return true;
472 }
473 return contains(B.getList());
474 }
475
476
477 /**
478 * Ideal containment. Test if b is contained in this ideal. Note: this is
479 * eventually modified to become a Groebner Base.
480 * @param b polynomial
481 * @return true, if b is contained in this, else false
482 */
483 public boolean contains(GenPolynomial<C> b) {
484 if (b == null || b.isZERO()) {
485 return true;
486 }
487 if (this.isONE()) {
488 return true;
489 }
490 if (this.isZERO()) {
491 return false;
492 }
493 if (!isGB) {
494 doGB();
495 }
496 GenPolynomial<C> z;
497 z = red.normalform(getList(), b);
498 if (z == null || z.isZERO()) {
499 return true;
500 }
501 return false;
502 }
503
504
505 /**
506 * Ideal containment. Test if each b in B is contained in this ideal. Note:
507 * this is eventually modified to become a Groebner Base.
508 * @param B list of polynomials
509 * @return true, if each b in B is contained in this, else false
510 */
511 public boolean contains(List<GenPolynomial<C>> B) {
512 if (B == null || B.size() == 0) {
513 return true;
514 }
515 if (this.isONE()) {
516 return true;
517 }
518 if (!isGB) {
519 doGB();
520 }
521 for (GenPolynomial<C> b : B) {
522 if (b == null) {
523 continue;
524 }
525 GenPolynomial<C> z = red.normalform(getList(), b);
526 if (!z.isZERO()) {
527 //System.out.println("contains nf(b) != 0: " + b);
528 return false;
529 }
530 }
531 return true;
532 }
533
534
535 /**
536 * Summation. Generators for the sum of ideals. Note: if both ideals are
537 * Groebner bases, a Groebner base is returned.
538 * @param B ideal
539 * @return ideal(this+B)
540 */
541 public Ideal<C> sum(Ideal<C> B) {
542 if (B == null || B.isZERO()) {
543 return this;
544 }
545 if (this.isZERO()) {
546 return B;
547 }
548 int s = getList().size() + B.getList().size();
549 List<GenPolynomial<C>> c;
550 c = new ArrayList<GenPolynomial<C>>(s);
551 c.addAll(getList());
552 c.addAll(B.getList());
553 Ideal<C> I = new Ideal<C>(getRing(), c, false);
554 if (isGB && B.isGB) {
555 I.doGB();
556 }
557 return I;
558 }
559
560
561 /**
562 * Summation. Generators for the sum of ideal and a polynomial. Note: if
563 * this ideal is a Groebner base, a Groebner base is returned.
564 * @param b polynomial
565 * @return ideal(this+{b})
566 */
567 public Ideal<C> sum(GenPolynomial<C> b) {
568 if (b == null || b.isZERO()) {
569 return this;
570 }
571 int s = getList().size() + 1;
572 List<GenPolynomial<C>> c;
573 c = new ArrayList<GenPolynomial<C>>(s);
574 c.addAll(getList());
575 c.add(b);
576 Ideal<C> I = new Ideal<C>(getRing(), c, false);
577 if (isGB) {
578 I.doGB();
579 }
580 return I;
581 }
582
583
584 /**
585 * Summation. Generators for the sum of this ideal and a list of
586 * polynomials. Note: if this ideal is a Groebner base, a Groebner base is
587 * returned.
588 * @param L list of polynomials
589 * @return ideal(this+L)
590 */
591 public Ideal<C> sum(List<GenPolynomial<C>> L) {
592 if (L == null || L.isEmpty()) {
593 return this;
594 }
595 int s = getList().size() + L.size();
596 List<GenPolynomial<C>> c = new ArrayList<GenPolynomial<C>>(s);
597 c.addAll(getList());
598 c.addAll(L);
599 Ideal<C> I = new Ideal<C>(getRing(), c, false);
600 if (isGB) {
601 I.doGB();
602 }
603 return I;
604 }
605
606
607 /**
608 * Product. Generators for the product of ideals. Note: if both ideals are
609 * Groebner bases, a Groebner base is returned.
610 * @param B ideal
611 * @return ideal(this*B)
612 */
613 public Ideal<C> product(Ideal<C> B) {
614 if (B == null || B.isZERO()) {
615 return B;
616 }
617 if (this.isZERO()) {
618 return this;
619 }
620 int s = getList().size() * B.getList().size();
621 List<GenPolynomial<C>> c;
622 c = new ArrayList<GenPolynomial<C>>(s);
623 for (GenPolynomial<C> p : getList()) {
624 for (GenPolynomial<C> q : B.getList()) {
625 q = p.multiply(q);
626 c.add(q);
627 }
628 }
629 Ideal<C> I = new Ideal<C>(getRing(), c, false);
630 if (isGB && B.isGB) {
631 I.doGB();
632 }
633 return I;
634 }
635
636
637 /**
638 * Intersection. Generators for the intersection of ideals. Using an
639 * iterative algorithm.
640 * @param Bl list of ideals
641 * @return ideal(cap_i B_i), a Groebner base
642 */
643 public Ideal<C> intersect(List<Ideal<C>> Bl) {
644 if (Bl == null || Bl.size() == 0) {
645 return getZERO();
646 }
647 Ideal<C> I = null;
648 for (Ideal<C> B : Bl) {
649 if (I == null) {
650 I = B;
651 continue;
652 }
653 if (I.isONE()) {
654 return I;
655 }
656 I = I.intersect(B);
657 }
658 return I;
659 }
660
661
662 /**
663 * Intersection. Generators for the intersection of ideals.
664 * @param B ideal
665 * @return ideal(this \cap B), a Groebner base
666 */
667 public Ideal<C> intersect(Ideal<C> B) {
668 if (B == null || B.isZERO()) { // (0)
669 return B;
670 }
671 if (this.isZERO()) {
672 return this;
673 }
674 int s = getList().size() + B.getList().size();
675 List<GenPolynomial<C>> c;
676 c = new ArrayList<GenPolynomial<C>>(s);
677 List<GenPolynomial<C>> a = getList();
678 List<GenPolynomial<C>> b = B.getList();
679
680 GenPolynomialRing<C> tfac = getRing().extend(1);
681 // term order is also adjusted
682 for (GenPolynomial<C> p : a) {
683 p = p.extend(tfac, 0, 1L); // t*p
684 c.add(p);
685 }
686 for (GenPolynomial<C> p : b) {
687 GenPolynomial<C> q = p.extend(tfac, 0, 1L);
688 GenPolynomial<C> r = p.extend(tfac, 0, 0L);
689 p = r.subtract(q); // (1-t)*p
690 c.add(p);
691 }
692 logger.warn("intersect computing GB");
693 List<GenPolynomial<C>> g = bb.GB(c);
694 if (debug) {
695 logger.debug("intersect GB = " + g);
696 }
697 Ideal<C> E = new Ideal<C>(tfac, g, true);
698 Ideal<C> I = E.intersect(getRing());
699 return I;
700 }
701
702
703 /**
704 * Intersection. Generators for the intersection of a ideal with a
705 * polynomial ring. The polynomial ring of this ideal must be a contraction
706 * of R and the TermOrder must be an elimination order.
707 * @param R polynomial ring
708 * @return ideal(this \cap R)
709 */
710 public Ideal<C> intersect(GenPolynomialRing<C> R) {
711 if (R == null) {
712 throw new IllegalArgumentException("R may not be null");
713 }
714 int d = getRing().nvar - R.nvar;
715 if (d <= 0) {
716 return this;
717 }
718 List<GenPolynomial<C>> H = new ArrayList<GenPolynomial<C>>(getList().size());
719 for (GenPolynomial<C> p : getList()) {
720 Map<ExpVector, GenPolynomial<C>> m = null;
721 m = p.contract(R);
722 if (debug) {
723 logger.debug("intersect contract m = " + m);
724 }
725 if (m.size() == 1) { // contains one power of variables
726 for (ExpVector e : m.keySet()) {
727 if (e.isZERO()) {
728 H.add(m.get(e));
729 }
730 }
731 }
732 }
733 GenPolynomialRing<C> tfac = getRing().contract(d);
734 if (tfac.equals(R)) { // check
735 return new Ideal<C>(R, H, isGB, isTopt);
736 }
737 logger.info("tfac, R = " + tfac + ", " + R);
738 // throw new RuntimeException("contract(this) != R");
739 return new Ideal<C>(R, H); // compute GB
740 }
741
742
743 /**
744 * Eliminate. Generators for the intersection of a ideal with a polynomial
745 * ring. The polynomial rings must have variable names.
746 * @param R polynomial ring
747 * @return ideal(this \cap R)
748 */
749 public Ideal<C> eliminate(GenPolynomialRing<C> R) {
750 if (R == null) {
751 throw new IllegalArgumentException("R may not be null");
752 }
753 if (list.ring.equals(R)) {
754 return this;
755 }
756 String[] ename = R.getVars();
757 Ideal<C> I = eliminate(ename);
758 return I.intersect(R);
759 }
760
761
762 /**
763 * Eliminate. Preparation of generators for the intersection of a ideal with
764 * a polynomial ring.
765 * @param ename variables for the elimination ring.
766 * @return ideal(this) in K[ename,{vars \ ename}])
767 */
768 public Ideal<C> eliminate(String... ename) {
769 //System.out.println("ename = " + Arrays.toString(ename));
770 if (ename == null) {
771 throw new IllegalArgumentException("ename may not be null");
772 }
773 String[] aname = getRing().getVars();
774 //System.out.println("aname = " + Arrays.toString(aname));
775 if (aname == null) {
776 throw new IllegalArgumentException("aname may not be null");
777 }
778
779 GroebnerBasePartial<C> bbp = new GroebnerBasePartial<C>(bb, null);
780 String[] rname = GroebnerBasePartial.remainingVars(aname, ename);
781 //System.out.println("rname = " + Arrays.toString(rname));
782 PolynomialList<C> Pl = null;
783 if (rname.length == 0) {
784 if (Arrays.equals(aname, ename)) {
785 return this;
786 } else {
787 Pl = bbp.partialGB(getList(), ename); // normal GB
788 }
789 } else {
790 Pl = bbp.elimPartialGB(getList(), rname, ename); // reversed!
791 }
792 //System.out.println("Pl = " + Pl);
793 if (debug) {
794 logger.debug("elimination GB = " + Pl);
795 }
796 Ideal<C> I = new Ideal<C>(Pl, true);
797 return I;
798 }
799
800
801 /**
802 * Quotient. Generators for the ideal quotient.
803 * @param h polynomial
804 * @return ideal(this : h), a Groebner base
805 */
806 public Ideal<C> quotient(GenPolynomial<C> h) {
807 if (h == null) { // == (0)
808 return this;
809 }
810 if (h.isZERO()) {
811 return this;
812 }
813 if (this.isZERO()) {
814 return this;
815 }
816 List<GenPolynomial<C>> H;
817 H = new ArrayList<GenPolynomial<C>>(1);
818 H.add(h);
819 Ideal<C> Hi = new Ideal<C>(getRing(), H, true);
820
821 Ideal<C> I = this.intersect(Hi);
822
823 List<GenPolynomial<C>> Q;
824 Q = new ArrayList<GenPolynomial<C>>(I.getList().size());
825 for (GenPolynomial<C> q : I.getList()) {
826 q = q.divide(h); // remainder == 0
827 Q.add(q);
828 }
829 return new Ideal<C>(getRing(), Q, true /*false?*/);
830 }
831
832
833 /**
834 * Quotient. Generators for the ideal quotient.
835 * @param H ideal
836 * @return ideal(this : H), a Groebner base
837 */
838 public Ideal<C> quotient(Ideal<C> H) {
839 if (H == null) { // == (0)
840 return this;
841 }
842 if (H.isZERO()) {
843 return this;
844 }
845 if (this.isZERO()) {
846 return this;
847 }
848 Ideal<C> Q = null;
849 for (GenPolynomial<C> h : H.getList()) {
850 Ideal<C> Hi = this.quotient(h);
851 if (Q == null) {
852 Q = Hi;
853 } else {
854 Q = Q.intersect(Hi);
855 }
856 }
857 return Q;
858 }
859
860
861 /**
862 * Infinite quotient. Generators for the infinite ideal quotient.
863 * @param h polynomial
864 * @return ideal(this : h<sup>s</sup>), a Groebner base
865 */
866 public Ideal<C> infiniteQuotientRab(GenPolynomial<C> h) {
867 if (h == null || h.isZERO()) { // == (0)
868 return getONE();
869 }
870 if (h.isONE()) {
871 return this;
872 }
873 if (this.isZERO()) {
874 return this;
875 }
876 Ideal<C> I = this.GB(); // should be already
877 List<GenPolynomial<C>> a = I.getList();
878 List<GenPolynomial<C>> c;
879 c = new ArrayList<GenPolynomial<C>>(a.size() + 1);
880
881 GenPolynomialRing<C> tfac = getRing().extend(1);
882 // term order is also adjusted
883 for (GenPolynomial<C> p : a) {
884 p = p.extend(tfac, 0, 0L); // p
885 c.add(p);
886 }
887 GenPolynomial<C> q = h.extend(tfac, 0, 1L);
888 GenPolynomial<C> r = tfac.getONE(); // h.extend( tfac, 0, 0L );
889 GenPolynomial<C> hs = q.subtract(r); // 1 - t*h // (1-t)*h
890 c.add(hs);
891 logger.warn("infiniteQuotientRab computing GB ");
892 List<GenPolynomial<C>> g = bb.GB(c);
893 if (debug) {
894 logger.info("infiniteQuotientRab = " + tfac + ", c = " + c);
895 logger.info("infiniteQuotientRab GB = " + g);
896 }
897 Ideal<C> E = new Ideal<C>(tfac, g, true);
898 Ideal<C> Is = E.intersect(getRing());
899 return Is;
900 }
901
902
903 /**
904 * Infinite quotient exponent.
905 * @param h polynomial
906 * @param Q quotient this : h^\infinity
907 * @return s with Q = this : h<sup>s</sup>
908 */
909 public int infiniteQuotientExponent(GenPolynomial<C> h, Ideal<C> Q) {
910 int s = 0;
911 if (h == null) { // == 0
912 return s;
913 }
914 if (h.isZERO() || h.isONE()) {
915 return s;
916 }
917 if (this.isZERO() || this.isONE()) {
918 return s;
919 }
920 //see below: if (this.contains(Q)) {
921 // return s;
922 //}
923 GenPolynomial<C> p = getRing().getONE();
924 for (GenPolynomial<C> q : Q.getList()) {
925 if (this.contains(q)) {
926 continue;
927 }
928 //System.out.println("q = " + q + ", p = " + p + ", s = " + s);
929 GenPolynomial<C> qp = q.multiply(p);
930 while (!this.contains(qp)) {
931 p = p.multiply(h);
932 s++;
933 qp = q.multiply(p);
934 }
935 }
936 return s;
937 }
938
939
940 /**
941 * Infinite quotient. Generators for the infinite ideal quotient.
942 * @param h polynomial
943 * @return ideal(this : h<sup>s</sup>), a Groebner base
944 */
945 public Ideal<C> infiniteQuotient(GenPolynomial<C> h) {
946 if (h == null) { // == (0)
947 return this;
948 }
949 if (h.isZERO()) {
950 return this;
951 }
952 if (this.isZERO()) {
953 return this;
954 }
955 int s = 0;
956 Ideal<C> I = this.GB(); // should be already
957 GenPolynomial<C> hs = h;
958 Ideal<C> Is = I;
959
960 boolean eq = false;
961 while (!eq) {
962 Is = I.quotient(hs);
963 Is = Is.GB(); // should be already
964 logger.info("infiniteQuotient s = " + s);
965 eq = Is.contains(I); // I.contains(Is) always
966 if (!eq) {
967 I = Is;
968 s++;
969 // hs = hs.multiply( h );
970 }
971 }
972 return Is;
973 }
974
975
976 /**
977 * Radical membership test.
978 * @param h polynomial
979 * @return true if h is contained in the radical of ideal(this), else false.
980 */
981 public boolean isRadicalMember(GenPolynomial<C> h) {
982 if (h == null) { // == (0)
983 return true;
984 }
985 if (h.isZERO()) {
986 return true;
987 }
988 if (this.isZERO()) {
989 return true;
990 }
991 Ideal<C> x = infiniteQuotientRab(h);
992 if (debug) {
993 logger.debug("infiniteQuotientRab = " + x);
994 }
995 return x.isONE();
996 }
997
998
999 /**
1000 * Infinite quotient. Generators for the infinite ideal quotient.
1001 * @param h polynomial
1002 * @return ideal(this : h<sup>s</sup>), a Groebner base
1003 */
1004 public Ideal<C> infiniteQuotientOld(GenPolynomial<C> h) {
1005 if (h == null) { // == (0)
1006 return this;
1007 }
1008 if (h.isZERO()) {
1009 return this;
1010 }
1011 if (this.isZERO()) {
1012 return this;
1013 }
1014 int s = 0;
1015 Ideal<C> I = this.GB(); // should be already
1016 GenPolynomial<C> hs = h;
1017
1018 boolean eq = false;
1019 while (!eq) {
1020 Ideal<C> Is = I.quotient(hs);
1021 Is = Is.GB(); // should be already
1022 logger.debug("infiniteQuotient s = " + s);
1023 eq = Is.contains(I); // I.contains(Is) always
1024 if (!eq) {
1025 I = Is;
1026 s++;
1027 hs = hs.multiply(h);
1028 }
1029 }
1030 return I;
1031 }
1032
1033
1034 /**
1035 * Infinite Quotient. Generators for the ideal infinite quotient.
1036 * @param H ideal
1037 * @return ideal(this : H<sup>s</sup>), a Groebner base
1038 */
1039 public Ideal<C> infiniteQuotient(Ideal<C> H) {
1040 if (H == null) { // == (0)
1041 return this;
1042 }
1043 if (H.isZERO()) {
1044 return this;
1045 }
1046 if (this.isZERO()) {
1047 return this;
1048 }
1049 Ideal<C> Q = null;
1050 for (GenPolynomial<C> h : H.getList()) {
1051 Ideal<C> Hi = this.infiniteQuotient(h);
1052 if (Q == null) {
1053 Q = Hi;
1054 } else {
1055 Q = Q.intersect(Hi);
1056 }
1057 }
1058 return Q;
1059 }
1060
1061
1062 /**
1063 * Infinite Quotient. Generators for the ideal infinite quotient.
1064 * @param H ideal
1065 * @return ideal(this : H<sup>s</sup>), a Groebner base
1066 */
1067 public Ideal<C> infiniteQuotientRab(Ideal<C> H) {
1068 if (H == null) { // == (0)
1069 return this;
1070 }
1071 if (H.isZERO()) {
1072 return this;
1073 }
1074 if (this.isZERO()) {
1075 return this;
1076 }
1077 Ideal<C> Q = null;
1078 for (GenPolynomial<C> h : H.getList()) {
1079 Ideal<C> Hi = this.infiniteQuotientRab(h);
1080 if (Q == null) {
1081 Q = Hi;
1082 } else {
1083 Q = Q.intersect(Hi);
1084 }
1085 }
1086 return Q;
1087 }
1088
1089
1090 /**
1091 * Power. Generators for the power of this ideal. Note: if this ideal is
1092 * a Groebner base, a Groebner base is returned.
1093 * @param d integer
1094 * @return ideal(this^d)
1095 */
1096 public Ideal<C> power(int d) {
1097 if ( d <= 0 ) {
1098 return getONE();
1099 }
1100 if (this.isZERO() || this.isONE()) {
1101 return this;
1102 }
1103 Ideal<C> c = this;
1104 for (int i = 1; i < d; i++ ) {
1105 c = c.product(this);
1106 }
1107 return c;
1108 }
1109
1110
1111 /**
1112 * Normalform for element.
1113 * @param h polynomial
1114 * @return normalform of h with respect to this
1115 */
1116 public GenPolynomial<C> normalform(GenPolynomial<C> h) {
1117 if (h == null) {
1118 return h;
1119 }
1120 if (h.isZERO()) {
1121 return h;
1122 }
1123 if (this.isZERO()) {
1124 return h;
1125 }
1126 GenPolynomial<C> r;
1127 r = red.normalform(list.list, h);
1128 return r;
1129 }
1130
1131
1132 /**
1133 * Normalform for list of elements.
1134 * @param L polynomial list
1135 * @return list of normalforms of the elements of L with respect to this
1136 */
1137 public List<GenPolynomial<C>> normalform(List<GenPolynomial<C>> L) {
1138 if (L == null) {
1139 return L;
1140 }
1141 if (L.size() == 0) {
1142 return L;
1143 }
1144 if (this.isZERO()) {
1145 return L;
1146 }
1147 List<GenPolynomial<C>> M = new ArrayList<GenPolynomial<C>>(L.size());
1148 for (GenPolynomial<C> h : L) {
1149 GenPolynomial<C> r = normalform(h);
1150 if (r != null && !r.isZERO()) {
1151 M.add(r);
1152 }
1153 }
1154 return M;
1155 }
1156
1157
1158 /**
1159 * Inverse for element modulo this ideal.
1160 * @param h polynomial
1161 * @return inverse of h with respect to this, if defined
1162 */
1163 public GenPolynomial<C> inverse(GenPolynomial<C> h) {
1164 if (h == null || h.isZERO()) {
1165 throw new NotInvertibleException("zero not invertible");
1166 }
1167 if (this.isZERO()) {
1168 throw new NotInvertibleException("zero ideal");
1169 }
1170 if ( h.isUnit() ) {
1171 return h.inverse();
1172 }
1173 doGB();
1174 List<GenPolynomial<C>> F = new ArrayList<GenPolynomial<C>>(1 + list.list.size());
1175 F.add(h);
1176 F.addAll(list.list);
1177 //System.out.println("F = " + F);
1178 ExtendedGB<C> x = bb.extGB(F);
1179 List<GenPolynomial<C>> G = x.G;
1180 //System.out.println("G = " + G);
1181 GenPolynomial<C> one = null;
1182 int i = -1;
1183 for (GenPolynomial<C> p : G) {
1184 i++;
1185 if (p == null) {
1186 continue;
1187 }
1188 if (p.isUnit()) {
1189 one = p;
1190 break;
1191 }
1192 }
1193 if (one == null) {
1194 throw new NotInvertibleException(" h = " + h);
1195 }
1196 List<GenPolynomial<C>> row = x.G2F.get(i); // != -1
1197 GenPolynomial<C> g = row.get(0);
1198 if (g == null || g.isZERO()) {
1199 throw new NotInvertibleException(" h = " + h);
1200 }
1201 // adjust g to get g*h == 1
1202 GenPolynomial<C> f = g.multiply(h);
1203 GenPolynomial<C> k = red.normalform(list.list, f);
1204 if (!k.isONE()) {
1205 C lbc = k.leadingBaseCoefficient();
1206 lbc = lbc.inverse();
1207 g = g.multiply(lbc);
1208 }
1209 if (debug) {
1210 //logger.info("inv G = " + G);
1211 //logger.info("inv G2F = " + x.G2F);
1212 //logger.info("inv row "+i+" = " + row);
1213 //logger.info("inv h = " + h);
1214 //logger.info("inv g = " + g);
1215 //logger.info("inv f = " + f);
1216 f = g.multiply(h);
1217 k = red.normalform(list.list, f);
1218 logger.debug("inv k = " + k);
1219 if (!k.isUnit()) {
1220 throw new NotInvertibleException(" k = " + k);
1221 }
1222 }
1223 return g;
1224 }
1225
1226
1227 /**
1228 * Test if element is a unit modulo this ideal.
1229 * @param h polynomial
1230 * @return true if h is a unit with respect to this, else false
1231 */
1232 public boolean isUnit(GenPolynomial<C> h) {
1233 if (h == null || h.isZERO()) {
1234 return false;
1235 }
1236 if (this.isZERO()) {
1237 return false;
1238 }
1239 List<GenPolynomial<C>> F = new ArrayList<GenPolynomial<C>>(1 + list.list.size());
1240 F.add(h);
1241 F.addAll(list.list);
1242 List<GenPolynomial<C>> G = bb.GB(F);
1243 for (GenPolynomial<C> p : G) {
1244 if (p == null) {
1245 continue;
1246 }
1247 if (p.isUnit()) {
1248 return true;
1249 }
1250 }
1251 return false;
1252 }
1253
1254
1255 /**
1256 * Radical approximation. Squarefree generators for the ideal.
1257 * @return squarefree(this), a Groebner base
1258 */
1259 public Ideal<C> squarefree() {
1260 if (this.isZERO()) {
1261 return this;
1262 }
1263 Ideal<C> R = this;
1264 Ideal<C> Rp = null;
1265 List<GenPolynomial<C>> li, ri;
1266 while (true) {
1267 li = R.getList();
1268 ri = new ArrayList<GenPolynomial<C>>(li); //.size() );
1269 for (GenPolynomial<C> h : li) {
1270 GenPolynomial<C> r = engine.squarefreePart(h);
1271 ri.add(r);
1272 }
1273 Rp = new Ideal<C>(R.getRing(), ri, false);
1274 Rp.doGB();
1275 if (R.equals(Rp)) {
1276 break;
1277 }
1278 R = Rp;
1279 }
1280 return R;
1281 }
1282
1283
1284 /**
1285 * Ideal common zero test.
1286 * @return -1, 0 or 1 if dimension(this) &eq; -1, 0 or ≥ 1.
1287 */
1288 public int commonZeroTest() {
1289 if (this.isZERO()) {
1290 return 1;
1291 }
1292 if (!isGB) {
1293 doGB();
1294 }
1295 if (this.isONE()) {
1296 return -1;
1297 }
1298 return bb.commonZeroTest(getList());
1299 }
1300
1301
1302 /**
1303 * Test if this ideal is maximal.
1304 * @return true, if this is maximal and not one, else false.
1305 */
1306 public boolean isMaximal() {
1307 if ( commonZeroTest() != 0 ) {
1308 return false;
1309 }
1310 for (Long d : univariateDegrees()) {
1311 if (d > 1L) {
1312 // todo: test if irreducible
1313 return false;
1314 }
1315 }
1316 return true;
1317 }
1318
1319
1320 /**
1321 * Univariate head term degrees.
1322 * @return a list of the degrees of univariate head terms.
1323 */
1324 public List<Long> univariateDegrees() {
1325 List<Long> ud = new ArrayList<Long>();
1326 if (this.isZERO()) {
1327 return ud;
1328 }
1329 if (!isGB) {
1330 doGB();
1331 }
1332 if (this.isONE()) {
1333 return ud;
1334 }
1335 if (this.list.ring.nvar <= 0) {
1336 return ud;
1337 }
1338 //int uht = 0;
1339 Map<Integer, Long> v = new TreeMap<Integer, Long>(); // for non reduced GBs
1340 for (GenPolynomial<C> p : getList()) {
1341 ExpVector e = p.leadingExpVector();
1342 if (e == null) {
1343 continue;
1344 }
1345 int[] u = e.dependencyOnVariables();
1346 if (u == null) {
1347 continue;
1348 }
1349 if (u.length == 1) {
1350 //uht++;
1351 Long d = v.get(u[0]);
1352 if (d == null) {
1353 v.put(u[0], e.getVal(u[0]));
1354 }
1355 }
1356 }
1357 for (int i = 0; i < this.list.ring.nvar; i++) {
1358 Long d = v.get(i);
1359 ud.add(d);
1360 }
1361 //Collections.reverse(ud);
1362 return ud;
1363 }
1364
1365
1366 /**
1367 * Ideal dimension.
1368 * @return a dimension container (dim,maxIndep,list(maxIndep),vars).
1369 */
1370 public Dimension dimension() {
1371 int t = commonZeroTest();
1372 Set<Integer> S = new HashSet<Integer>();
1373 Set<Set<Integer>> M = new HashSet<Set<Integer>>();
1374 if (t <= 0) {
1375 return new Dimension(t, S, M, this.list.ring.getVars());
1376 }
1377 int d = 0;
1378 Set<Integer> U = new HashSet<Integer>();
1379 for (int i = 0; i < this.list.ring.nvar; i++) {
1380 U.add(i);
1381 }
1382 M = dimension(S, U, M);
1383 for (Set<Integer> m : M) {
1384 int dp = m.size();
1385 if (dp > d) {
1386 d = dp;
1387 S = m;
1388 }
1389 }
1390 return new Dimension(d, S, M, this.list.ring.getVars());
1391 }
1392
1393
1394 /**
1395 * Ideal dimension.
1396 * @param S is a set of independent variables.
1397 * @param U is a set of variables of unknown status.
1398 * @param M is a list of maximal sets of independent variables.
1399 * @return a list of maximal sets of independent variables, eventually
1400 * containing S.
1401 */
1402 protected Set<Set<Integer>> dimension(Set<Integer> S, Set<Integer> U, Set<Set<Integer>> M) {
1403 Set<Set<Integer>> MP = M;
1404 Set<Integer> UP = new HashSet<Integer>(U);
1405 for (Integer j : U) {
1406 UP.remove(j);
1407 Set<Integer> SP = new HashSet<Integer>(S);
1408 SP.add(j);
1409 if (!containsHT(SP, getList())) {
1410 MP = dimension(SP, UP, MP);
1411 }
1412 }
1413 boolean contained = false;
1414 for (Set<Integer> m : MP) {
1415 if (m.containsAll(S)) {
1416 contained = true;
1417 break;
1418 }
1419 }
1420 if (!contained) {
1421 MP.add(S);
1422 }
1423 return MP;
1424 }
1425
1426
1427 /**
1428 * Ideal head term containment test.
1429 * @param G list of polynomials.
1430 * @param H index set.
1431 * @return true, if the vaiables of the head terms of each polynomial in G
1432 * are contained in H, else false.
1433 */
1434 protected boolean containsHT(Set<Integer> H, List<GenPolynomial<C>> G) {
1435 Set<Integer> S = null;
1436 for (GenPolynomial<C> p : G) {
1437 if (p == null) {
1438 continue;
1439 }
1440 ExpVector e = p.leadingExpVector();
1441 if (e == null) {
1442 continue;
1443 }
1444 int[] v = e.dependencyOnVariables();
1445 if (v == null) {
1446 continue;
1447 }
1448 //System.out.println("v = " + Arrays.toString(v));
1449 if (S == null) { // revert indices
1450 S = new HashSet<Integer>(H.size());
1451 int r = e.length() - 1;
1452 for (Integer i : H) {
1453 S.add(r - i);
1454 }
1455 }
1456 if (contains(v, S)) { // v \subset S
1457 return true;
1458 }
1459 }
1460 return false;
1461 }
1462
1463
1464 /**
1465 * Set containment. is v \subset H.
1466 * @param v index array.
1467 * @param H index set.
1468 * @return true, if each element of v is contained in H, else false .
1469 */
1470 protected boolean contains(int[] v, Set<Integer> H) {
1471 for (int i = 0; i < v.length; i++) {
1472 if (!H.contains(v[i])) {
1473 return false;
1474 }
1475 }
1476 return true;
1477 }
1478
1479
1480 /**
1481 * Construct univariate polynomials of minimal degree in all variables in
1482 * zero dimensional ideal(G).
1483 * @return list of univariate polynomial of minimal degree in each variable
1484 * in ideal(G)
1485 */
1486 public List<GenPolynomial<C>> constructUnivariate() {
1487 List<GenPolynomial<C>> univs = new ArrayList<GenPolynomial<C>>();
1488 for (int i = list.ring.nvar - 1; i >= 0; i--) {
1489 GenPolynomial<C> u = constructUnivariate(i);
1490 univs.add(u);
1491 }
1492 return univs;
1493 }
1494
1495
1496 /**
1497 * Construct univariate polynomial of minimal degree in variable i in zero
1498 * dimensional ideal(G).
1499 * @param i variable index.
1500 * @return univariate polynomial of minimal degree in variable i in ideal(G)
1501 */
1502 public GenPolynomial<C> constructUnivariate(int i) {
1503 doGB();
1504 return constructUnivariate(i, getList());
1505 }
1506
1507
1508 /**
1509 * Construct univariate polynomial of minimal degree in variable i of a zero
1510 * dimensional ideal(G).
1511 * @param i variable index.
1512 * @param G list of polynomials, a monic reduced Gröbner base of a zero
1513 * dimensional ideal.
1514 * @return univariate polynomial of minimal degree in variable i in ideal(G)
1515 */
1516 public GenPolynomial<C> constructUnivariate(int i, List<GenPolynomial<C>> G) {
1517 if (G == null || G.size() == 0) {
1518 throw new IllegalArgumentException("G may not be null or empty");
1519 }
1520 List<Long> ud = univariateDegrees();
1521 if (ud == null || ud.size() <= i) {
1522 //logger.info("univ pol, ud = " + ud);
1523 throw new IllegalArgumentException("ideal(G) not zero dimensional " + ud);
1524 }
1525 int ll = 0;
1526 Long di = ud.get(i);
1527 if (di != null) {
1528 ll = (int) (long) di;
1529 } else {
1530 throw new IllegalArgumentException("ideal(G) not zero dimensional");
1531 }
1532 long vsdim = 1;
1533 for ( Long d : ud ) {
1534 if ( d != null ) {
1535 vsdim *= d;
1536 }
1537 }
1538 logger.info("univariate construction, deg = " + ll + ", vsdim = " + vsdim);
1539 GenPolynomialRing<C> pfac = G.get(0).ring;
1540 RingFactory<C> cfac = pfac.coFac;
1541 String var = pfac.getVars()[pfac.nvar - 1 - i];
1542 GenPolynomialRing<C> ufac = new GenPolynomialRing<C>(cfac, 1, new TermOrder(TermOrder.INVLEX),
1543 new String[] { var });
1544 GenPolynomial<C> pol = ufac.getZERO();
1545
1546 GenPolynomialRing<C> cpfac = new GenPolynomialRing<C>(cfac, ll, new TermOrder(TermOrder.INVLEX));
1547 GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(cpfac, pfac);
1548 GenPolynomial<GenPolynomial<C>> P = rfac.getZERO();
1549 for (int k = 0; k < ll; k++) {
1550 GenPolynomial<GenPolynomial<C>> Pp = rfac.univariate(i, k);
1551 GenPolynomial<C> cp = cpfac.univariate(cpfac.nvar - 1 - k);
1552 Pp = Pp.multiply(cp);
1553 P = P.sum(Pp);
1554 }
1555 if ( debug ) {
1556 logger.info("univariate construction, P = " + P);
1557 logger.info("univariate construction, G = " + G);
1558 //throw new RuntimeException("check");
1559 }
1560 GenPolynomial<C> X;
1561 GenPolynomial<C> XP;
1562 // solve system of linear equations for the coefficients of the univariate polynomial
1563 List<GenPolynomial<C>> ls;
1564 int z = -1;
1565 do {
1566 //System.out.println("ll = " + ll);
1567 GenPolynomial<GenPolynomial<C>> Pp = rfac.univariate(i, ll);
1568 GenPolynomial<C> cp = cpfac.univariate(cpfac.nvar - 1 - ll);
1569 Pp = Pp.multiply(cp);
1570 P = P.sum(Pp);
1571 X = pfac.univariate(i, ll);
1572 XP = red.normalform(G, X);
1573 //System.out.println("XP = " + XP);
1574 GenPolynomial<GenPolynomial<C>> XPp = PolyUtil.<C> toRecursive(rfac, XP);
1575 GenPolynomial<GenPolynomial<C>> XPs = XPp.sum(P);
1576 ls = new ArrayList<GenPolynomial<C>>(XPs.getMap().values());
1577 //System.out.println("ls,1 = " + ls);
1578 ls = red.irreducibleSet(ls);
1579 Ideal<C> L = new Ideal<C>(cpfac, ls, true);
1580 z = L.commonZeroTest();
1581 if (z != 0) {
1582 ll++;
1583 if ( ll > vsdim ) {
1584 logger.info("univariate construction, P = " + P);
1585 logger.info("univariate construction, nf(P) = " + XP);
1586 logger.info("G = " + G);
1587 throw new ArithmeticException("univariate polynomial degree greater than vector space dimansion");
1588 }
1589 cpfac = cpfac.extend(1);
1590 rfac = new GenPolynomialRing<GenPolynomial<C>>(cpfac, pfac);
1591 P = PolyUtil.<C> extendCoefficients(rfac, P, 0, 0L);
1592 XPp = PolyUtil.<C> extendCoefficients(rfac, XPp, 0, 1L);
1593 P = P.sum(XPp);
1594 }
1595 } while (z != 0); // && ll <= 5 && !XP.isZERO()
1596 // construct result polynomial
1597 pol = ufac.univariate(0, ll);
1598 for (GenPolynomial<C> pc : ls) {
1599 ExpVector e = pc.leadingExpVector();
1600 if (e == null) {
1601 continue;
1602 }
1603 int[] v = e.dependencyOnVariables();
1604 if (v == null || v.length == 0) {
1605 continue;
1606 }
1607 int vi = v[0];
1608 C tc = pc.trailingBaseCoefficient();
1609 tc = tc.negate();
1610 GenPolynomial<C> pi = ufac.univariate(0, ll - 1 - vi);
1611 pi = pi.multiply(tc);
1612 pol = pol.sum(pi);
1613 }
1614 if (logger.isInfoEnabled()) {
1615 logger.info("univariate construction, pol = " + pol);
1616 }
1617 return pol;
1618 }
1619
1620
1621 /**
1622 * Zero dimensional radical decompostition. See Seidenbergs lemma 92, and
1623 * BWK lemma 8.13.
1624 * @return intersection of radical ideals G_i with ideal(this) subseteq
1625 * cap_i( ideal(G_i) )
1626 */
1627 public List<IdealWithUniv<C>> zeroDimRadicalDecomposition() {
1628 List<IdealWithUniv<C>> dec = new ArrayList<IdealWithUniv<C>>();
1629 if (this.isZERO()) {
1630 return dec;
1631 }
1632 IdealWithUniv<C> iwu = new IdealWithUniv<C>(this, new ArrayList<GenPolynomial<C>>());
1633 dec.add(iwu);
1634 if (this.isONE()) {
1635 return dec;
1636 }
1637 if (list.ring.coFac.characteristic().signum() > 0 && !list.ring.coFac.isFinite()) {
1638 logger.warn("must use prime decomposition for char p and infinite coefficient rings, found "
1639 + list.ring.coFac.toScript());
1640 return zeroDimPrimeDecomposition();
1641 }
1642 for (int i = list.ring.nvar - 1; i >= 0; i--) {
1643 List<IdealWithUniv<C>> part = new ArrayList<IdealWithUniv<C>>();
1644 for (IdealWithUniv<C> id : dec) {
1645 //System.out.println("id = " + id + ", i = " + i);
1646 GenPolynomial<C> u = id.ideal.constructUnivariate(i);
1647 SortedMap<GenPolynomial<C>, Long> facs = engine.baseSquarefreeFactors(u);
1648 if (facs == null || facs.size() == 0 || (facs.size() == 1 && facs.get(facs.firstKey()) == 1L)) {
1649 List<GenPolynomial<C>> iup = new ArrayList<GenPolynomial<C>>();
1650 iup.addAll(id.upolys);
1651 iup.add(u);
1652 IdealWithUniv<C> Ipu = new IdealWithUniv<C>(id.ideal, iup);
1653 part.add(Ipu);
1654 continue; // irreducible
1655 }
1656 if (logger.isInfoEnabled()) {
1657 logger.info("squarefree facs = " + facs);
1658 }
1659 GenPolynomialRing<C> mfac = id.ideal.list.ring;
1660 int j = mfac.nvar - 1 - i;
1661 for (GenPolynomial<C> p : facs.keySet()) {
1662 // make p multivariate
1663 GenPolynomial<C> pm = p.extendUnivariate(mfac, j);
1664 // mfac.parse( p.toString() );
1665 //stem.out.println("pm = " + pm);
1666 Ideal<C> Ip = id.ideal.sum(pm);
1667 List<GenPolynomial<C>> iup = new ArrayList<GenPolynomial<C>>();
1668 iup.addAll(id.upolys);
1669 iup.add(p);
1670 IdealWithUniv<C> Ipu = new IdealWithUniv<C>(Ip, iup);
1671 if (debug) {
1672 logger.info("ideal with squarefree facs = " + Ipu);
1673 }
1674 part.add(Ipu);
1675 }
1676 }
1677 dec = part;
1678 part = new ArrayList<IdealWithUniv<C>>();
1679 }
1680 return dec;
1681 }
1682
1683
1684 /**
1685 * Test for Zero dimensional radical. See Seidenbergs lemma 92, and BWK
1686 * lemma 8.13.
1687 * @return true if this is an zero dimensional radical ideal, else false
1688 */
1689 public boolean isZeroDimRadical() {
1690 if (this.isZERO()) {
1691 return false;
1692 }
1693 if (this.isONE()) {
1694 return false; // not 0-dim
1695 }
1696 if (list.ring.coFac.characteristic().signum() > 0 && !list.ring.coFac.isFinite()) {
1697 logger.warn("radical only for char 0 or finite coefficient rings, but found "
1698 + list.ring.coFac.toScript());
1699 }
1700 for (int i = list.ring.nvar - 1; i >= 0; i--) {
1701 GenPolynomial<C> u = constructUnivariate(i);
1702 boolean t = engine.isSquarefree(u);
1703 if (!t) {
1704 System.out.println("not squarefree " + engine.squarefreePart(u) + ", " + u);
1705 return false;
1706 }
1707 }
1708 return true;
1709 }
1710
1711
1712 /**
1713 * Zero dimensional ideal irreducible decompostition. See algorithm DIRGZD
1714 * of BGK 1986 and also PREDEC of the Gröbner bases book 1993.
1715 * @return intersection H, of ideals G_i with ideal(this) subseteq cap_i(
1716 * ideal(G_i) ) and each ideal G_i has only irreducible minimal
1717 * univariate polynomials and the G_i are pairwise co-prime.
1718 */
1719 public List<IdealWithUniv<C>> zeroDimDecomposition() {
1720 List<IdealWithUniv<C>> dec = new ArrayList<IdealWithUniv<C>>();
1721 if (this.isZERO()) {
1722 return dec;
1723 }
1724 IdealWithUniv<C> iwu = new IdealWithUniv<C>(this, new ArrayList<GenPolynomial<C>>());
1725 dec.add(iwu);
1726 if (this.isONE()) {
1727 return dec;
1728 }
1729 FactorAbstract<C> ufd = FactorFactory.<C> getImplementation(list.ring.coFac);
1730 for (int i = list.ring.nvar - 1; i >= 0; i--) {
1731 List<IdealWithUniv<C>> part = new ArrayList<IdealWithUniv<C>>();
1732 for (IdealWithUniv<C> id : dec) {
1733 //System.out.println("id.ideal = " + id.ideal);
1734 GenPolynomial<C> u = id.ideal.constructUnivariate(i);
1735 SortedMap<GenPolynomial<C>, Long> facs = ufd.baseFactors(u);
1736 if (facs == null || facs.size() == 0 || (facs.size() == 1 && facs.get(facs.firstKey()) == 1L)) {
1737 List<GenPolynomial<C>> iup = new ArrayList<GenPolynomial<C>>();
1738 iup.addAll(id.upolys);
1739 iup.add(u);
1740 IdealWithUniv<C> Ipu = new IdealWithUniv<C>(id.ideal, iup);
1741 part.add(Ipu);
1742 continue; // irreducible
1743 }
1744 if (debug) {
1745 logger.info("irreducible facs = " + facs);
1746 }
1747 GenPolynomialRing<C> mfac = id.ideal.list.ring;
1748 int j = mfac.nvar - 1 - i;
1749 for (GenPolynomial<C> p : facs.keySet()) {
1750 // make p multivariate
1751 GenPolynomial<C> pm = p.extendUnivariate(mfac, j);
1752 // mfac.parse( p.toString() );
1753 //System.out.println("pm = " + pm);
1754 Ideal<C> Ip = id.ideal.sum(pm);
1755 List<GenPolynomial<C>> iup = new ArrayList<GenPolynomial<C>>();
1756 iup.addAll(id.upolys);
1757 iup.add(p);
1758 IdealWithUniv<C> Ipu = new IdealWithUniv<C>(Ip, iup);
1759 part.add(Ipu);
1760 }
1761 }
1762 dec = part;
1763 part = new ArrayList<IdealWithUniv<C>>();
1764 }
1765 return dec;
1766 }
1767
1768
1769 /**
1770 * Zero dimensional ideal irreducible decompostition extension. One step
1771 * decomposition via a minimal univariate polynomial in the lowest variable,
1772 * used after each normalPosition step.
1773 * @param upol list of univariate polynomials
1774 * @param og list of other generators for the ideal
1775 * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
1776 * ideal(G_i) ) and all minimal univariate polynomials of all G_i
1777 * are irreducible
1778 */
1779 public List<IdealWithUniv<C>> zeroDimDecompositionExtension(List<GenPolynomial<C>> upol,
1780 List<GenPolynomial<C>> og) {
1781 if (upol == null || upol.size() + 1 != list.ring.nvar) {
1782 throw new IllegalArgumentException("univariate polynomial list not correct " + upol);
1783 }
1784 List<IdealWithUniv<C>> dec = new ArrayList<IdealWithUniv<C>>();
1785 if (this.isZERO()) {
1786 return dec;
1787 }
1788 IdealWithUniv<C> iwu = new IdealWithUniv<C>(this, upol);
1789 if (this.isONE()) {
1790 dec.add(iwu);
1791 return dec;
1792 }
1793 FactorAbstract<C> ufd = FactorFactory.<C> getImplementation(list.ring.coFac);
1794 int i = list.ring.nvar - 1;
1795 //IdealWithUniv<C> id = new IdealWithUniv<C>(this,upol);
1796 GenPolynomial<C> u = this.constructUnivariate(i);
1797 SortedMap<GenPolynomial<C>, Long> facs = ufd.baseFactors(u);
1798 if (facs.size() == 1 && facs.get(facs.firstKey()) == 1L) {
1799 List<GenPolynomial<C>> iup = new ArrayList<GenPolynomial<C>>();
1800 iup.add(u); // new polynomial first
1801 iup.addAll(upol);
1802 IdealWithUniv<C> Ipu = new IdealWithUniv<C>(this, iup, og);
1803 dec.add(Ipu);
1804 return dec;
1805 }
1806 if (true) {
1807 logger.info("irreducible facs = " + facs);
1808 }
1809 GenPolynomialRing<C> mfac = list.ring;
1810 int j = mfac.nvar - 1 - i;
1811 for (GenPolynomial<C> p : facs.keySet()) {
1812 // make p multivariate
1813 GenPolynomial<C> pm = p.extendUnivariate(mfac, j);
1814 //System.out.println("pm = " + pm);
1815 Ideal<C> Ip = this.sum(pm);
1816 List<GenPolynomial<C>> iup = new ArrayList<GenPolynomial<C>>();
1817 iup.add(p); // new polynomial first
1818 iup.addAll(upol);
1819 IdealWithUniv<C> Ipu = new IdealWithUniv<C>(Ip, iup, og);
1820 dec.add(Ipu);
1821 }
1822 return dec;
1823 }
1824
1825
1826 /**
1827 * Test for zero dimensional ideal decompostition.
1828 * @param L intersection of ideals G_i with ideal(G) subseteq cap_i(
1829 * ideal(G_i) ) and all minimal univariate polynomials of all G_i
1830 * are irreducible
1831 * @return true if L is a zero dimensional irreducible decomposition of
1832 * this, else false
1833 */
1834 public boolean isZeroDimDecomposition(List<IdealWithUniv<C>> L) {
1835 if (L == null || L.size() == 0) {
1836 if (this.isZERO()) {
1837 return true;
1838 } else {
1839 return false;
1840 }
1841 }
1842 // add lower variables if L contains ideals from field extensions
1843 GenPolynomialRing<C> ofac = list.ring;
1844 int r = ofac.nvar;
1845 int rp = L.get(0).ideal.list.ring.nvar;
1846 int d = rp - r;
1847 //System.out.println("d = " + d);
1848 Ideal<C> Id = this;
1849 if (d > 0) {
1850 GenPolynomialRing<C> nfac = ofac.extendLower(d);
1851 //System.out.println("nfac = " + nfac);
1852 List<GenPolynomial<C>> elist = new ArrayList<GenPolynomial<C>>(list.list.size());
1853 for (GenPolynomial<C> p : getList()) {
1854 //System.out.println("p = " + p);
1855 GenPolynomial<C> q = p.extendLower(nfac, 0, 0L);
1856 //System.out.println("q = " + q);
1857 elist.add(q);
1858 }
1859 Id = new Ideal<C>(nfac, elist, isGB, isTopt);
1860 }
1861 // test if this is contained in the intersection
1862 for (IdealWithUniv<C> I : L) {
1863 boolean t = I.ideal.contains(Id);
1864 if (!t) {
1865 System.out.println("not contained " + this + " in " + I.ideal);
1866 return false;
1867 }
1868 }
1869 // test if all univariate polynomials are contained in the respective ideal
1870 //List<GenPolynomial<C>> upprod = new ArrayList<GenPolynomial<C>>(rp);
1871 for (IdealWithUniv<C> I : L) {
1872 GenPolynomialRing<C> mfac = I.ideal.list.ring;
1873 int i = 0;
1874 for (GenPolynomial<C> p : I.upolys) {
1875 GenPolynomial<C> pm = p.extendUnivariate(mfac, i++);
1876 //System.out.println("pm = " + pm + ", p = " + p);
1877 boolean t = I.ideal.contains(pm);
1878 if (!t) {
1879 System.out.println("not contained " + pm + " in " + I.ideal);
1880 return false;
1881 }
1882 }
1883 }
1884 return true;
1885 }
1886
1887
1888 /**
1889 * Compute normal position for variables i and j.
1890 * @param i first variable index
1891 * @param j second variable index
1892 * @param og other generators for the ideal
1893 * @return this + (z - x_j - t x_i) in the ring C[z, x_1, ..., x_r]
1894 */
1895 public IdealWithUniv<C> normalPositionFor(int i, int j, List<GenPolynomial<C>> og) {
1896 // extend variables by one
1897 GenPolynomialRing<C> ofac = list.ring;
1898 if (ofac.tord.getEvord() != TermOrder.INVLEX) {
1899 throw new IllegalArgumentException("invalid term order for normalPosition " + ofac.tord);
1900 }
1901 if ( ofac.characteristic().signum() == 0 ) {
1902 return normalPositionForChar0(i,j,og);
1903 }
1904 return normalPositionForCharP(i,j,og);
1905 }
1906
1907
1908 /**
1909 * Compute normal position for variables i and j, characteristic zero.
1910 * @param i first variable index
1911 * @param j second variable index
1912 * @param og other generators for the ideal
1913 * @return this + (z - x_j - t x_i) in the ring C[z, x_1, ..., x_r]
1914 */
1915 IdealWithUniv<C> normalPositionForChar0(int i, int j, List<GenPolynomial<C>> og) {
1916 // extend variables by one
1917 GenPolynomialRing<C> ofac = list.ring;
1918 GenPolynomialRing<C> nfac = ofac.extendLower(1);
1919 List<GenPolynomial<C>> elist = new ArrayList<GenPolynomial<C>>(list.list.size() + 1);
1920 for (GenPolynomial<C> p : getList()) {
1921 GenPolynomial<C> q = p.extendLower(nfac, 0, 0L);
1922 //System.out.println("q = " + q);
1923 elist.add(q);
1924 }
1925 List<GenPolynomial<C>> ogen = new ArrayList<GenPolynomial<C>>();
1926 if (og != null && og.size() > 0) {
1927 for (GenPolynomial<C> p : og) {
1928 GenPolynomial<C> q = p.extendLower(nfac, 0, 0L);
1929 //System.out.println("q = " + q);
1930 ogen.add(q);
1931 }
1932 }
1933 Ideal<C> I = new Ideal<C>(nfac, elist, true);
1934 //System.out.println("I = " + I);
1935 int ip = list.ring.nvar - 1 - i;
1936 int jp = list.ring.nvar - 1 - j;
1937 GenPolynomial<C> xi = nfac.univariate(ip);
1938 GenPolynomial<C> xj = nfac.univariate(jp);
1939 GenPolynomial<C> z = nfac.univariate(nfac.nvar - 1);
1940 // compute GBs until value of t is OK
1941 Ideal<C> Ip;
1942 GenPolynomial<C> zp;
1943 int t = 0;
1944 do {
1945 t--;
1946 // zp = z - ( xj - xi * t )
1947 zp = z.subtract(xj.subtract(xi.multiply( nfac.fromInteger(t) )));
1948 zp = zp.monic();
1949 Ip = I.sum(zp);
1950 //System.out.println("Ip = " + Ip);
1951 if (-t % 5 == 0) {
1952 logger.info("normal position, t = " + t);
1953 }
1954 } while (!Ip.isNormalPositionFor(i + 1, j + 1));
1955 if (debug) {
1956 logger.info("normal position = " + Ip);
1957 }
1958 ogen.add(zp);
1959 IdealWithUniv<C> Ips = new IdealWithUniv<C>(Ip, null, ogen);
1960 return Ips;
1961 }
1962
1963
1964 /**
1965 * Compute normal position for variables i and j, positive characteristic.
1966 * @param i first variable index
1967 * @param j second variable index
1968 * @param og other generators for the ideal
1969 * @return this + (z - x_j - t x_i) in the ring C[z, x_1, ..., x_r]
1970 */
1971 IdealWithUniv<C> normalPositionForCharP(int i, int j, List<GenPolynomial<C>> og) {
1972 // extend variables by one
1973 GenPolynomialRing<C> ofac = list.ring;
1974 GenPolynomialRing<C> nfac = ofac.extendLower(1);
1975 List<GenPolynomial<C>> elist = new ArrayList<GenPolynomial<C>>(list.list.size() + 1);
1976 for (GenPolynomial<C> p : getList()) {
1977 GenPolynomial<C> q = p.extendLower(nfac, 0, 0L);
1978 //System.out.println("q = " + q);
1979 elist.add(q);
1980 }
1981 List<GenPolynomial<C>> ogen = new ArrayList<GenPolynomial<C>>();
1982 if (og != null && og.size() > 0) {
1983 for (GenPolynomial<C> p : og) {
1984 GenPolynomial<C> q = p.extendLower(nfac, 0, 0L);
1985 //System.out.println("q = " + q);
1986 ogen.add(q);
1987 }
1988 }
1989 Ideal<C> I = new Ideal<C>(nfac, elist, true);
1990 //System.out.println("I = " + I);
1991 int ip = list.ring.nvar - 1 - i;
1992 int jp = list.ring.nvar - 1 - j;
1993 GenPolynomial<C> xi = nfac.univariate(ip);
1994 GenPolynomial<C> xj = nfac.univariate(jp);
1995 GenPolynomial<C> z = nfac.univariate(nfac.nvar - 1);
1996 // compute GBs until value of t is OK
1997 Ideal<C> Ip;
1998 GenPolynomial<C> zp;
1999 AlgebraicNumberRing<C> afac = null;
2000 Iterator<AlgebraicNumber<C>> aiter = null;
2001 //String obr = "";
2002 //String cbr = "";
2003 int t = 0;
2004 do {
2005 t--;
2006 // zp = z - ( xj - xi * t )
2007 GenPolynomial<C> tn;
2008 if (afac == null) {
2009 tn = nfac.fromInteger(t);
2010 if (tn.isZERO()) {
2011 RingFactory<C> fac = nfac.coFac;
2012 //int braces = 2;
2013 while (!(fac instanceof AlgebraicNumberRing)) {
2014 if (fac instanceof GenPolynomialRing) {
2015 GenPolynomialRing<C> pfac = (GenPolynomialRing) (Object) fac;
2016 fac = pfac.coFac;
2017 } else if (fac instanceof QuotientRing) {
2018 QuotientRing<C> pfac = (QuotientRing) (Object) fac;
2019 fac = pfac.ring.coFac;
2020 } else {
2021 throw new ArithmeticException("field elements exhausted, need algebraic extension of base ring");
2022 }
2023 //braces++;
2024 }
2025 //for (int ii = 0; ii < braces; ii++) {
2026 // obr += "{ ";
2027 // cbr += " }";
2028 //}
2029 afac = (AlgebraicNumberRing) (Object) fac;
2030 logger.info("afac = " + afac.toScript());
2031 aiter = afac.iterator();
2032 AlgebraicNumber<C> an = aiter.next();
2033 for (int kk = 0; kk < afac.characteristic().intValue(); kk++) {
2034 an = aiter.next();
2035 }
2036 //System.out.println("an,iter = " + an);
2037 //tn = nfac.parse(obr + an.toString() + cbr);
2038 tn = nfac.parse(an.toString());
2039 //System.out.println("tn = " + tn);
2040 if (false) {
2041 throw new RuntimeException("probe");
2042 }
2043 }
2044 } else {
2045 if (!aiter.hasNext()) {
2046 throw new ArithmeticException("field elements exhausted, normal position not reachable: !aiter.hasNext(): " + t);
2047 }
2048 AlgebraicNumber<C> an = aiter.next();
2049 //System.out.println("an,iter = " + an);
2050 //tn = nfac.parse(obr + an.toString() + cbr);
2051 tn = nfac.parse(an.toString());
2052 //System.out.println("tn = " + tn);
2053 }
2054 if (tn.isZERO()) {
2055 throw new ArithmeticException("field elements exhausted, normal position not reachable: tn == 0: " + t);
2056 }
2057 zp = z.subtract(xj.subtract(xi.multiply(tn)));
2058 zp = zp.monic();
2059 Ip = I.sum(zp);
2060 //System.out.println("Ip = " + Ip);
2061 if (-t % 4 == 0) {
2062 logger.info("normal position, t = " + t);
2063 logger.info("normal position, GB = " + Ip);
2064 if (t < -550) {
2065 throw new ArithmeticException("normal position not reached in " + t + " steps");
2066 }
2067 }
2068 } while (!Ip.isNormalPositionFor(i + 1, j + 1));
2069 if (debug) {
2070 logger.info("normal position = " + Ip);
2071 }
2072 ogen.add(zp);
2073 IdealWithUniv<C> Ips = new IdealWithUniv<C>(Ip, null, ogen);
2074 return Ips;
2075 }
2076
2077
2078 /**
2079 * Test if this ideal is in normal position for variables i and j.
2080 * @param i first variable index
2081 * @param j second variable index
2082 * @return true if this is in normal position with respect to i and j
2083 */
2084 public boolean isNormalPositionFor(int i, int j) {
2085 // called in extended ring!
2086 int ip = list.ring.nvar - 1 - i;
2087 int jp = list.ring.nvar - 1 - j;
2088 boolean iOK = false;
2089 boolean jOK = false;
2090 for (GenPolynomial<C> p : getList()) {
2091 ExpVector e = p.leadingExpVector();
2092 int[] dov = e.dependencyOnVariables();
2093 //System.out.println("dov = " + Arrays.toString(dov));
2094 if (dov.length == 0) {
2095 throw new IllegalArgumentException("ideal dimension is not zero");
2096 }
2097 if (dov[0] == ip) {
2098 if (e.totalDeg() != 1) {
2099 return false;
2100 } else {
2101 iOK = true;
2102 }
2103 } else if (dov[0] == jp) {
2104 if (e.totalDeg() != 1) {
2105 return false;
2106 } else {
2107 jOK = true;
2108 }
2109 }
2110 if (iOK && jOK) {
2111 return true;
2112 }
2113 }
2114 return iOK && jOK;
2115 }
2116
2117
2118 /**
2119 * Normal position index, separate for polynomials with more than 2
2120 * variables. See also <a
2121 * href="http://krum.rz.uni-mannheim.de/mas/src/masring/DIPDEC0.mi.html"
2122 * >mas.masring.DIPDEC0#DIGISR</a>
2123 * @return (i,j) for non-normal variables
2124 */
2125 public int[] normalPositionIndex2Vars() {
2126 int[] np = null;
2127 int i = -1;
2128 int j = -1;
2129 for (GenPolynomial<C> p : getList()) {
2130 ExpVector e = p.leadingExpVector();
2131 int[] dov = e.dependencyOnVariables();
2132 //System.out.println("dov_head = " + Arrays.toString(dov));
2133 if (dov.length == 0) {
2134 throw new IllegalArgumentException("ideal dimension is not zero " + p);
2135 }
2136 // search bi-variate head terms
2137 if (dov.length >= 2) {
2138 i = dov[0];
2139 j = dov[1];
2140 break;
2141 }
2142 int n = dov[0];
2143 GenPolynomial<C> q = p.reductum();
2144 e = q.degreeVector();
2145 dov = e.dependencyOnVariables();
2146 //System.out.println("dov_red = " + Arrays.toString(dov));
2147 int k = Arrays.binarySearch(dov, n);
2148 int len = 2;
2149 if (k >= 0) {
2150 len = 3;
2151 }
2152 // search bi-variate reductas
2153 if (dov.length >= len) {
2154 switch (k) {
2155 case 0:
2156 i = dov[1];
2157 j = dov[2];
2158 break;
2159 case 1:
2160 i = dov[0];
2161 j = dov[2];
2162 break;
2163 case 2:
2164 i = dov[0];
2165 j = dov[1];
2166 break;
2167 default:
2168 i = dov[0];
2169 j = dov[1];
2170 break;
2171 }
2172 break;
2173 }
2174 }
2175 if (i < 0 || j < 0) {
2176 return np;
2177 }
2178 // adjust index
2179 i = list.ring.nvar - 1 - i;
2180 j = list.ring.nvar - 1 - j;
2181 np = new int[] { j, i }; // reverse
2182 logger.info("normalPositionIndex2Vars, np = " + Arrays.toString(np));
2183 return np;
2184 }
2185
2186
2187 /**
2188 * Normal position index, separate multiple univariate polynomials. See also
2189 * <a href="http://krum.rz.uni-mannheim.de/mas/src/masring/DIPDEC0.mi.html">
2190 * mas.masring.DIPDEC0#DIGISM</a>
2191 * @return (i,j) for non-normal variables
2192 */
2193 public int[] normalPositionIndexUnivars() {
2194 int[] np = null; //new int[] { -1, -1 };
2195 int i = -1;
2196 int j = -1;
2197 // search multiple univariate polynomials with degree >= 2
2198 for (GenPolynomial<C> p : getList()) {
2199 ExpVector e = p.degreeVector();
2200 int[] dov = e.dependencyOnVariables();
2201 long t = e.totalDeg(); // lt(p) would be enough
2202 //System.out.println("dov_univ = " + Arrays.toString(dov) + ", e = " + e);
2203 if (dov.length == 0) {
2204 throw new IllegalArgumentException("ideal dimension is not zero");
2205 }
2206 if (dov.length == 1 && t >= 2L) {
2207 if (i == -1) {
2208 i = dov[0];
2209 } else if (j == -1) {
2210 j = dov[0];
2211 if (i > j) {
2212 int x = i;
2213 i = j;
2214 j = x;
2215 }
2216 }
2217 }
2218 if (i >= 0 && j >= 0) {
2219 break;
2220 }
2221 }
2222 if (i < 0 || j < 0) {
2223 // search polynomials with univariate head term and degree >= 2
2224 for (GenPolynomial<C> p : getList()) {
2225 ExpVector e = p.leadingExpVector();
2226 long t = e.totalDeg();
2227 if (t >= 2) {
2228 e = p.degreeVector();
2229 int[] dov = e.dependencyOnVariables();
2230 //System.out.println("dov_univ2 = " + Arrays.toString(dov) + " e = " + e);
2231 if (dov.length == 0) {
2232 throw new IllegalArgumentException("ideal dimension is not zero");
2233 }
2234 if (dov.length >= 2) {
2235 i = dov[0];
2236 j = dov[1];
2237 }
2238 }
2239 if (i >= 0 && j >= 0) {
2240 break;
2241 }
2242 }
2243 }
2244 if (i < 0 || j < 0) {
2245 return np;
2246 }
2247 // adjust index
2248 i = list.ring.nvar - 1 - i;
2249 j = list.ring.nvar - 1 - j;
2250 np = new int[] { j, i }; // reverse
2251 logger.info("normalPositionIndexUnivars, np = " + Arrays.toString(np));
2252 return np;
2253 }
2254
2255
2256 /**
2257 * Zero dimensional ideal decompostition for real roots. See algorithm
2258 * mas.masring.DIPDEC0#DINTSR.
2259 * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
2260 * ideal(G_i) ) and each G_i contains at most bi-variate polynomials
2261 * and all univariate minimal polynomials are irreducible
2262 */
2263 public List<IdealWithUniv<C>> zeroDimRootDecomposition() {
2264 List<IdealWithUniv<C>> dec = zeroDimDecomposition();
2265 if (this.isZERO()) {
2266 return dec;
2267 }
2268 if (this.isONE()) {
2269 return dec;
2270 }
2271 List<IdealWithUniv<C>> rdec = new ArrayList<IdealWithUniv<C>>();
2272 while (dec.size() > 0) {
2273 IdealWithUniv<C> id = dec.remove(0);
2274 int[] ri = id.ideal.normalPositionIndex2Vars();
2275 if (ri == null || ri.length != 2) {
2276 rdec.add(id);
2277 } else {
2278 IdealWithUniv<C> I = id.ideal.normalPositionFor(ri[0], ri[1], id.others);
2279 List<IdealWithUniv<C>> rd = I.ideal.zeroDimDecompositionExtension(id.upolys, I.others);
2280 //System.out.println("r_rd = " + rd);
2281 dec.addAll(rd);
2282 }
2283 }
2284 return rdec;
2285 }
2286
2287
2288 /**
2289 * Zero dimensional ideal prime decompostition. See algorithm
2290 * mas.masring.DIPDEC0#DINTSS.
2291 * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
2292 * ideal(G_i) ) and each G_i is a prime ideal
2293 */
2294 public List<IdealWithUniv<C>> zeroDimPrimeDecomposition() {
2295 List<IdealWithUniv<C>> pdec = zeroDimPrimeDecompositionFE();
2296 List<IdealWithUniv<C>> dec = new ArrayList<IdealWithUniv<C>>();
2297 if ( pdec.size() == 1 ) { // already prime
2298 IdealWithUniv<C> Ip = pdec.get(0);
2299 int s = Ip.upolys.size() - getRing().nvar; // skip field ext univariate polys
2300 List<GenPolynomial<C>> upol = Ip.upolys.subList(s, Ip.upolys.size());
2301 Ip = new IdealWithUniv<C>(this,upol);
2302 dec.add(Ip);
2303 return dec;
2304 }
2305 for (IdealWithUniv<C> Ip : pdec) {
2306 if (Ip.ideal.getRing().nvar == getRing().nvar) { // no field extension
2307 dec.add(Ip);
2308 } else { // remove field extension
2309 // add other generators for performance
2310 Ideal<C> Id = Ip.ideal;
2311 if (Ip.others != null) {
2312 //System.out.println("adding Ip.others = " + Ip.others);
2313 List<GenPolynomial<C>> pp = new ArrayList<GenPolynomial<C>>();
2314 pp.addAll(Id.getList());
2315 pp.addAll(Ip.others);
2316 Id = new Ideal<C>(Id.getRing(), pp);
2317 }
2318 Ideal<C> Is = Id.eliminate(getRing());
2319 //System.out.println("Is = " + Is);
2320 int s = Ip.upolys.size() - getRing().nvar; // skip field ext univariate polys
2321 List<GenPolynomial<C>> upol = Ip.upolys.subList(s, Ip.upolys.size());
2322 IdealWithUniv<C> Iu = new IdealWithUniv<C>(Is, upol);
2323 //,Ip.others); used above and must be ignored here
2324 dec.add(Iu);
2325 }
2326 }
2327 return dec;
2328 }
2329
2330
2331 /**
2332 * Zero dimensional ideal prime decompostition, with field extension. See
2333 * algorithm mas.masring.DIPDEC0#DINTSS.
2334 * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
2335 * ideal(G_i) ) and each G_i is a prime ideal with eventually
2336 * containing field extension variables
2337 */
2338 public List<IdealWithUniv<C>> zeroDimPrimeDecompositionFE() {
2339 List<IdealWithUniv<C>> dec = zeroDimRootDecomposition();
2340 if (this.isZERO()) {
2341 return dec;
2342 }
2343 if (this.isONE()) {
2344 return dec;
2345 }
2346 List<IdealWithUniv<C>> rdec = new ArrayList<IdealWithUniv<C>>();
2347 while (dec.size() > 0) {
2348 IdealWithUniv<C> id = dec.remove(0);
2349 int[] ri = id.ideal.normalPositionIndexUnivars();
2350 if (ri == null || ri.length != 2) {
2351 rdec.add(id);
2352 } else {
2353 IdealWithUniv<C> I = id.ideal.normalPositionFor(ri[0], ri[1], id.others);
2354 List<IdealWithUniv<C>> rd = I.ideal.zeroDimDecompositionExtension(id.upolys, I.others);
2355 //System.out.println("rd = " + rd);
2356 dec.addAll(rd);
2357 }
2358 }
2359 return rdec;
2360 }
2361
2362
2363 /**
2364 * Zero dimensional ideal associated primary ideal. See algorithm
2365 * mas.masring.DIPIDEAL#DIRLPI.
2366 * @param P prime ideal associated to this
2367 * @return primary ideal of this with respect to the associated pime ideal P
2368 */
2369 public Ideal<C> primaryIdeal(Ideal<C> P) {
2370 Ideal<C> Qs = P;
2371 Ideal<C> Q;
2372 int e = 0;
2373 do {
2374 Q = Qs;
2375 e++;
2376 Qs = Q.product(P);
2377 } while (Qs.contains(this));
2378 boolean t;
2379 Ideal<C> As;
2380 do {
2381 As = this.sum(Qs);
2382 t = As.contains(Q);
2383 if (!t) {
2384 Q = Qs;
2385 e++;
2386 Qs = Q.product(P);
2387 }
2388 } while (!t);
2389 logger.info("exponent = " + e);
2390 return As;
2391 }
2392
2393
2394 /**
2395 * Zero dimensional ideal primary decompostition.
2396 * @return list of primary components of primary ideals G_i (pairwise
2397 * co-prime) with ideal(this) = cap_i( ideal(G_i) ) together with
2398 * the associated primes
2399 */
2400 public List<PrimaryComponent<C>> zeroDimPrimaryDecomposition() {
2401 List<IdealWithUniv<C>> pdec = zeroDimPrimeDecomposition();
2402 if (logger.isInfoEnabled()) {
2403 logger.info("prim decomp = " + pdec);
2404 }
2405 return zeroDimPrimaryDecomposition(pdec);
2406 }
2407
2408
2409 /**
2410 * Zero dimensional ideal elimination to original ring.
2411 * @param pdec list of prime ideals G_i
2412 * @return intersection of pairwise co-prime prime ideals G_i in the ring of
2413 * this with ideal(this) = cap_i( ideal(G_i) )
2414 */
2415 public List<IdealWithUniv<C>> zeroDimElimination(List<IdealWithUniv<C>> pdec) {
2416 List<IdealWithUniv<C>> dec = new ArrayList<IdealWithUniv<C>>();
2417 if (this.isZERO()) {
2418 return dec;
2419 }
2420 if (this.isONE()) {
2421 dec.add(pdec.get(0));
2422 return dec;
2423 }
2424 List<IdealWithUniv<C>> qdec = new ArrayList<IdealWithUniv<C>>();
2425 for (IdealWithUniv<C> Ip : pdec) {
2426 //System.out.println("Ip = " + Ip);
2427 List<GenPolynomial<C>> epol = new ArrayList<GenPolynomial<C>>();
2428 epol.addAll(Ip.ideal.getList());
2429 GenPolynomialRing<C> mfac = Ip.ideal.list.ring;
2430 int j = 0;
2431 // add univariate polynomials for performance
2432 for (GenPolynomial<C> p : Ip.upolys) {
2433 GenPolynomial<C> pm = p.extendUnivariate(mfac, j++);
2434 if (j != 1) { // skip double
2435 epol.add(pm);
2436 }
2437 }
2438 // add other generators for performance
2439 if (Ip.others != null) {
2440 epol.addAll(Ip.others);
2441 }
2442 Ideal<C> Ipp = new Ideal<C>(mfac, epol);
2443 // logger.info("eliminate_1 = " + Ipp);
2444 TermOrder to = null;
2445 if (mfac.tord.getEvord() != TermOrder.IGRLEX) {
2446 List<GenPolynomial<C>> epols = new ArrayList<GenPolynomial<C>>();
2447 to = new TermOrder(TermOrder.IGRLEX);
2448 GenPolynomialRing<C> smfac = new GenPolynomialRing<C>(mfac.coFac, mfac.nvar, to, mfac
2449 .getVars());
2450 for (GenPolynomial<C> p : epol) {
2451 GenPolynomial<C> pm = smfac.copy(p);
2452 epols.add(pm.monic());
2453 }
2454 //epol = epols;
2455 Ipp = new Ideal<C>(smfac, epols);
2456 }
2457 epol = red.irreducibleSet(Ipp.getList());
2458 Ipp = new Ideal<C>(Ipp.getRing(), epol);
2459 if (logger.isInfoEnabled()) {
2460 logger.info("eliminate = " + Ipp);
2461 }
2462 Ideal<C> Is = Ipp.eliminate(list.ring);
2463 //System.out.println("Is = " + Is);
2464 if (to != null && !Is.list.ring.equals(list.ring)) {
2465 List<GenPolynomial<C>> epols = new ArrayList<GenPolynomial<C>>();
2466 for (GenPolynomial<C> p : Is.getList()) {
2467 GenPolynomial<C> pm = list.ring.copy(p);
2468 epols.add(pm);
2469 }
2470 Is = new Ideal<C>(list.ring, epols);
2471 //System.out.println("Is = " + Is);
2472 }
2473 int k = Ip.upolys.size() - list.ring.nvar;
2474 List<GenPolynomial<C>> up = new ArrayList<GenPolynomial<C>>();
2475 for (int i = 0; i < list.ring.nvar; i++) {
2476 up.add(Ip.upolys.get(i + k));
2477 }
2478 IdealWithUniv<C> Ie = new IdealWithUniv<C>(Is, up);
2479 qdec.add(Ie);
2480 }
2481 return qdec;
2482 }
2483
2484
2485 /**
2486 * Zero dimensional ideal primary decompostition.
2487 * @param pdec list of prime ideals G_i with no field extensions
2488 * @return list of primary components of primary ideals G_i (pairwise
2489 * co-prime) with ideal(this) = cap_i( ideal(G_i) ) together with
2490 * the associated primes
2491 */
2492 public List<PrimaryComponent<C>> zeroDimPrimaryDecomposition(List<IdealWithUniv<C>> pdec) {
2493 List<PrimaryComponent<C>> dec = new ArrayList<PrimaryComponent<C>>();
2494 if (this.isZERO()) {
2495 return dec;
2496 }
2497 if (this.isONE()) {
2498 PrimaryComponent<C> pc = new PrimaryComponent<C>(pdec.get(0).ideal, pdec.get(0));
2499 dec.add(pc);
2500 return dec;
2501 }
2502 for (IdealWithUniv<C> Ip : pdec) {
2503 Ideal<C> Qs = this.primaryIdeal(Ip.ideal);
2504 PrimaryComponent<C> pc = new PrimaryComponent<C>(Qs, Ip);
2505 dec.add(pc);
2506 }
2507 return dec;
2508 }
2509
2510
2511 /**
2512 * Test for primary ideal decompostition.
2513 * @param L list of primary components G_i
2514 * @return true if ideal(this) == cap_i( ideal(G_i) )
2515 */
2516 public boolean isPrimaryDecomposition(List<PrimaryComponent<C>> L) {
2517 // test if this is contained in the intersection
2518 for (PrimaryComponent<C> I : L) {
2519 boolean t = I.primary.contains(this);
2520 if (!t) {
2521 System.out.println("not contained " + this + " in " + I);
2522 return false;
2523 }
2524 }
2525 Ideal<C> isec = null;
2526 for (PrimaryComponent<C> I : L) {
2527 if (isec == null) {
2528 isec = I.primary;
2529 } else {
2530 isec = isec.intersect(I.primary);
2531 }
2532 }
2533 return this.contains(isec);
2534 }
2535
2536
2537 /**
2538 * Ideal extension.
2539 * @param vars list of variables for a polynomial ring for extension
2540 * @return ideal G, with coefficients in
2541 * QuotientRing(GenPolynomialRing<C>(vars))
2542 */
2543 public IdealWithUniv<Quotient<C>> extension(String... vars) {
2544 GenPolynomialRing<C> fac = getRing();
2545 GenPolynomialRing<C> efac = new GenPolynomialRing<C>(fac.coFac, vars.length, fac.tord, vars);
2546 IdealWithUniv<Quotient<C>> ext = extension(efac);
2547 return ext;
2548 }
2549
2550
2551 /**
2552 * Ideal extension.
2553 * @param efac polynomial ring for extension
2554 * @return ideal G, with coefficients in QuotientRing(efac)
2555 */
2556 public IdealWithUniv<Quotient<C>> extension(GenPolynomialRing<C> efac) {
2557 QuotientRing<C> qfac = new QuotientRing<C>(efac);
2558 IdealWithUniv<Quotient<C>> ext = extension(qfac);
2559 return ext;
2560 }
2561
2562
2563 /**
2564 * Ideal extension.
2565 * @param qfac quotient polynomial ring for extension
2566 * @return ideal G, with coefficients in qfac
2567 */
2568 public IdealWithUniv<Quotient<C>> extension(QuotientRing<C> qfac) {
2569 GenPolynomialRing<C> fac = getRing();
2570 GenPolynomialRing<C> efac = qfac.ring;
2571 String[] rvars = GroebnerBasePartial.remainingVars(fac.getVars(), efac.getVars());
2572 //System.out.println("rvars = " + Arrays.toString(rvars));
2573
2574 GroebnerBasePartial<C> bbp = new GroebnerBasePartial<C>();
2575 //wrong: OptimizedPolynomialList<C> pgb = bbp.partialGB(getList(),rvars);
2576 OptimizedPolynomialList<C> pgb = bbp.elimPartialGB(getList(), rvars, efac.getVars());
2577 if (logger.isInfoEnabled()) {
2578 logger.info("rvars = " + Arrays.toString(rvars));
2579 logger.info("partialGB = " + pgb);
2580 }
2581
2582 GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(efac,
2583 rvars.length, fac.tord, rvars);
2584 List<GenPolynomial<C>> plist = pgb.list;
2585 List<GenPolynomial<GenPolynomial<C>>> rpgb = PolyUtil.<C> recursive(rfac, plist);
2586 //System.out.println("rfac = " + rfac);
2587 GenPolynomialRing<Quotient<C>> qpfac = new GenPolynomialRing<Quotient<C>>(qfac, rfac);
2588 List<GenPolynomial<Quotient<C>>> qpgb = PolyUfdUtil.<C> quotientFromIntegralCoefficients(qpfac, rpgb);
2589 //System.out.println("qpfac = " + qpfac);
2590
2591 // compute f
2592 GreatestCommonDivisor<C> ufd = GCDFactory.getImplementation(fac.coFac);
2593 GenPolynomial<C> f = null; // qfac.ring.getONE();
2594 for (GenPolynomial<GenPolynomial<C>> p : rpgb) {
2595 if (f == null) {
2596 f = p.leadingBaseCoefficient();
2597 } else {
2598 f = ufd.lcm(f, p.leadingBaseCoefficient());
2599 }
2600 }
2601 //SquarefreeAbstract<C> sqf = SquarefreeFactory.getImplementation(fac.coFac);
2602 //not required: f = sqf.squarefreePart(f);
2603 GenPolynomial<GenPolynomial<C>> fp = rfac.getONE().multiply(f);
2604 GenPolynomial<Quotient<C>> fq = PolyUfdUtil.<C> quotientFromIntegralCoefficients(qpfac, fp);
2605 if (logger.isInfoEnabled()) {
2606 logger.info("extension f = " + f);
2607 logger.info("ext = " + qpgb);
2608 }
2609 List<GenPolynomial<Quotient<C>>> upols = new ArrayList<GenPolynomial<Quotient<C>>>(0);
2610 List<GenPolynomial<Quotient<C>>> opols = new ArrayList<GenPolynomial<Quotient<C>>>(1);
2611 opols.add(fq);
2612
2613 qpgb = PolyUtil.<Quotient<C>> monic(qpgb);
2614 Ideal<Quotient<C>> ext = new Ideal<Quotient<C>>(qpfac, qpgb);
2615 IdealWithUniv<Quotient<C>> extu = new IdealWithUniv<Quotient<C>>(ext, upols, opols);
2616 return extu;
2617 }
2618
2619
2620 /**
2621 * Ideal contraction and permutation.
2622 * @param eideal extension ideal of this.
2623 * @return contraction ideal of eideal in this polynomial ring
2624 */
2625 public IdealWithUniv<C> permContraction(IdealWithUniv<Quotient<C>> eideal) {
2626 return Ideal.<C> permutation(getRing(), Ideal.<C> contraction(eideal));
2627 }
2628
2629
2630 /**
2631 * Ideal contraction.
2632 * @param eid extension ideal of this.
2633 * @return contraction ideal of eid in distributed polynomial ring
2634 */
2635 public static <C extends GcdRingElem<C>> IdealWithUniv<C> contraction(IdealWithUniv<Quotient<C>> eid) {
2636 Ideal<Quotient<C>> eideal = eid.ideal;
2637 List<GenPolynomial<Quotient<C>>> qgb = eideal.getList();
2638 QuotientRing<C> qfac = (QuotientRing<C>) eideal.getRing().coFac;
2639 GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(qfac.ring, eideal
2640 .getRing());
2641 GenPolynomialRing<C> dfac = qfac.ring.extend(eideal.getRing().getVars());
2642 TermOrder to = new TermOrder(qfac.ring.tord.getEvord());
2643 dfac = new GenPolynomialRing<C>(dfac.coFac, dfac.nvar, to, dfac.getVars());
2644 //System.out.println("qfac = " + qfac);
2645 //System.out.println("rfac = " + rfac);
2646 //System.out.println("dfac = " + dfac);
2647 // convert polynomials
2648 List<GenPolynomial<GenPolynomial<C>>> cgb = PolyUfdUtil.<C> integralFromQuotientCoefficients(rfac,
2649 qgb);
2650 List<GenPolynomial<C>> dgb = PolyUtil.<C> distribute(dfac, cgb);
2651 Ideal<C> cont = new Ideal<C>(dfac, dgb);
2652 // convert other polynomials
2653 List<GenPolynomial<C>> opols = new ArrayList<GenPolynomial<C>>();
2654 if (eid.others != null && eid.others.size() > 0) {
2655 List<GenPolynomial<GenPolynomial<C>>> orpol = PolyUfdUtil.<C> integralFromQuotientCoefficients(
2656 rfac, eid.others);
2657 List<GenPolynomial<C>> opol = PolyUtil.<C> distribute(dfac, orpol);
2658 opols.addAll(opol);
2659 }
2660 // convert univariate polynomials
2661 List<GenPolynomial<C>> upols = new ArrayList<GenPolynomial<C>>(0);
2662 int i = 0;
2663 for (GenPolynomial<Quotient<C>> p : eid.upolys) {
2664 GenPolynomial<Quotient<C>> pm = p.extendUnivariate(eideal.getRing(), i++);
2665 //System.out.println("pm = " + pm + ", p = " + p);
2666 GenPolynomial<GenPolynomial<C>> urpol = PolyUfdUtil
2667 .<C> integralFromQuotientCoefficients(rfac, pm);
2668 GenPolynomial<C> upol = PolyUtil.<C> distribute(dfac, urpol);
2669 upols.add(upol);
2670 //System.out.println("upol = " + upol);
2671 }
2672 // compute f
2673 GreatestCommonDivisor<C> ufd = GCDFactory.getImplementation(qfac.ring.coFac);
2674 GenPolynomial<C> f = null; // qfac.ring.getONE();
2675 for (GenPolynomial<GenPolynomial<C>> p : cgb) {
2676 if (f == null) {
2677 f = p.leadingBaseCoefficient();
2678 } else {
2679 f = ufd.lcm(f, p.leadingBaseCoefficient());
2680 }
2681 }
2682 GenPolynomial<GenPolynomial<C>> fp = rfac.getONE().multiply(f);
2683 f = PolyUtil.<C> distribute(dfac, fp);
2684 if (logger.isInfoEnabled()) {
2685 logger.info("contraction f = " + f);
2686 logger.info("cont = " + cont);
2687 }
2688 opols.add(f);
2689 if (f.isONE()) {
2690 IdealWithUniv<C> cf = new IdealWithUniv<C>(cont, upols, opols);
2691 return cf;
2692 }
2693 // compute ideal quotient by f
2694 Ideal<C> Q = cont.infiniteQuotientRab(f);
2695 IdealWithUniv<C> Qu = new IdealWithUniv<C>(Q, upols, opols);
2696 return Qu;
2697 }
2698
2699
2700 /**
2701 * Ideal permutation.
2702 * @param oring polynomial ring to which variables are back permuted.
2703 * @param Cont ideal to be permuted
2704 * @return permutation of cont in polynomial ring oring
2705 */
2706 public static <C extends GcdRingElem<C>> IdealWithUniv<C> permutation(GenPolynomialRing<C> oring,
2707 IdealWithUniv<C> Cont) {
2708 Ideal<C> cont = Cont.ideal;
2709 GenPolynomialRing<C> dfac = cont.getRing();
2710 // (back) permutation of variables
2711 String[] ovars = oring.getVars();
2712 String[] dvars = dfac.getVars();
2713 //System.out.println("ovars = " + Arrays.toString(ovars));
2714 //System.out.println("dvars = " + Arrays.toString(dvars));
2715 if (Arrays.equals(ovars, dvars)) { // nothing to do
2716 return Cont;
2717 }
2718 List<Integer> perm = GroebnerBasePartial.getPermutation(dvars, ovars);
2719 //System.out.println("perm = " + perm);
2720 GenPolynomialRing<C> pfac = TermOrderOptimization.<C> permutation(perm, cont.getRing());
2721 if (logger.isInfoEnabled()) {
2722 logger.info("pfac = " + pfac);
2723 }
2724 List<GenPolynomial<C>> ppolys = TermOrderOptimization.<C> permutation(perm, pfac, cont.getList());
2725 //System.out.println("ppolys = " + ppolys);
2726 cont = new Ideal<C>(pfac, ppolys);
2727 if (logger.isDebugEnabled()) {
2728 logger.info("perm cont = " + cont);
2729 }
2730 List<GenPolynomial<C>> opolys = TermOrderOptimization.<C> permutation(perm, pfac, Cont.others);
2731 //System.out.println("opolys = " + opolys);
2732 List<GenPolynomial<C>> upolys = TermOrderOptimization.<C> permutation(perm, pfac, Cont.upolys);
2733 //System.out.println("opolys = " + opolys);
2734 IdealWithUniv<C> Cu = new IdealWithUniv<C>(cont, upolys, opolys);
2735 return Cu;
2736 }
2737
2738
2739 /**
2740 * Ideal radical.
2741 * @return the radical ideal of this
2742 */
2743 public Ideal<C> radical() {
2744 List<IdealWithUniv<C>> rdec = radicalDecomposition();
2745 List<Ideal<C>> dec = new ArrayList<Ideal<C>>(rdec.size());
2746 for (IdealWithUniv<C> ru : rdec) {
2747 dec.add(ru.ideal);
2748 }
2749 Ideal<C> R = intersect(dec);
2750 return R;
2751 }
2752
2753
2754 /**
2755 * Ideal radical decompostition.
2756 * @return intersection of ideals G_i with radical(this) eq cap_i(
2757 * ideal(G_i) ) and each G_i is a radical ideal and the G_i are
2758 * pairwise co-prime
2759 */
2760 public List<IdealWithUniv<C>> radicalDecomposition() {
2761 // check dimension
2762 int z = commonZeroTest();
2763 List<IdealWithUniv<C>> dec = new ArrayList<IdealWithUniv<C>>();
2764 List<GenPolynomial<C>> ups = new ArrayList<GenPolynomial<C>>();
2765 // dimension -1
2766 if (z < 0) {
2767 IdealWithUniv<C> id = new IdealWithUniv<C>(this, ups);
2768 dec.add(id); // see GB book
2769 return dec;
2770 }
2771 // dimension 0
2772 if (z == 0) {
2773 dec = zeroDimRadicalDecomposition();
2774 return dec;
2775 }
2776 // dimension > 0
2777 if (this.isZERO()) {
2778 return dec;
2779 }
2780 if (list.ring.coFac.characteristic().signum() > 0 && !list.ring.coFac.isFinite()) {
2781 // must not be the case at this point
2782 logger.warn("must use prime decomposition for char p and infinite coefficient rings, found "
2783 + list.ring.coFac.toScript());
2784 return primeDecomposition();
2785 }
2786 Dimension dim = dimension();
2787 if (logger.isInfoEnabled()) {
2788 logger.info("dimension = " + dim);
2789 }
2790
2791 // shortest maximal independent set
2792 Set<Set<Integer>> M = dim.M;
2793 Set<Integer> min = null;
2794 for (Set<Integer> m : M) {
2795 if (min == null) {
2796 min = m;
2797 continue;
2798 }
2799 if (m.size() < min.size()) {
2800 min = m;
2801 }
2802 }
2803 //System.out.println("min = " + min);
2804 String[] mvars = new String[min.size()];
2805 int j = 0;
2806 for (Integer i : min) {
2807 mvars[j++] = dim.v[i];
2808 }
2809 if (logger.isInfoEnabled()) {
2810 logger.info("extension for variables = " + Arrays.toString(mvars));
2811 }
2812 // reduce to dimension zero
2813 IdealWithUniv<Quotient<C>> Ext = extension(mvars);
2814 if (logger.isInfoEnabled()) {
2815 logger.info("extension = " + Ext);
2816 }
2817
2818 List<IdealWithUniv<Quotient<C>>> edec = Ext.ideal.zeroDimRadicalDecomposition();
2819 if (logger.isInfoEnabled()) {
2820 logger.info("0-dim radical decomp = " + edec);
2821 }
2822 // remove field extensions are not required
2823 // reconstruct dimension
2824 for (IdealWithUniv<Quotient<C>> ep : edec) {
2825 IdealWithUniv<C> cont = permContraction(ep);
2826 //System.out.println("cont = " + cont);
2827 dec.add(cont);
2828 }
2829 IdealWithUniv<C> extcont = permContraction(Ext);
2830 //System.out.println("extcont = " + extcont);
2831
2832 // get f
2833 List<GenPolynomial<C>> ql = extcont.others;
2834 if (ql.size() == 0) { // should not happen
2835 return dec;
2836 }
2837 GenPolynomial<C> fx = ql.get(0);
2838 //System.out.println("cont(Ext) fx = " + fx + ", " + fx.ring);
2839 if (fx.isONE()) {
2840 return dec;
2841 }
2842 Ideal<C> T = sum(fx);
2843 //System.out.println("T.rec = " + T.getList());
2844 if (T.isONE()) {
2845 logger.info("1 in ideal for " + fx);
2846 return dec;
2847 }
2848 if (logger.isInfoEnabled()) {
2849 logger.info("radical decomp ext-cont fx = " + fx);
2850 logger.info("recursion radical decomp T = " + T);
2851 }
2852 // recursion:
2853 List<IdealWithUniv<C>> Tdec = T.radicalDecomposition();
2854 if (logger.isInfoEnabled()) {
2855 logger.info("recursion radical decomp = " + Tdec);
2856 }
2857 dec.addAll(Tdec);
2858 return dec;
2859 }
2860
2861
2862 /**
2863 * Ideal irreducible decompostition.
2864 * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
2865 * ideal(G_i) ) and each G_i is an ideal with irreducible univariate
2866 * polynomials (after extension to a zero dimensional ideal) and the
2867 * G_i are pairwise co-prime
2868 */
2869 public List<IdealWithUniv<C>> decomposition() {
2870 // check dimension
2871 int z = commonZeroTest();
2872 List<IdealWithUniv<C>> dec = new ArrayList<IdealWithUniv<C>>();
2873 List<GenPolynomial<C>> ups = new ArrayList<GenPolynomial<C>>();
2874 // dimension -1
2875 if (z < 0) {
2876 //IdealWithUniv<C> id = new IdealWithUniv<C>(this, ups);
2877 //dec.add(id); see GB book
2878 return dec;
2879 }
2880 // dimension 0
2881 if (z == 0) {
2882 dec = zeroDimDecomposition();
2883 return dec;
2884 }
2885 // dimension > 0
2886 if (this.isZERO()) {
2887 return dec;
2888 }
2889 Dimension dim = dimension();
2890 if (logger.isInfoEnabled()) {
2891 logger.info("dimension = " + dim);
2892 }
2893
2894 // shortest maximal independent set
2895 Set<Set<Integer>> M = dim.M;
2896 Set<Integer> min = null;
2897 for (Set<Integer> m : M) {
2898 if (min == null) {
2899 min = m;
2900 continue;
2901 }
2902 if (m.size() < min.size()) {
2903 min = m;
2904 }
2905 }
2906 //System.out.println("min = " + min);
2907 String[] mvars = new String[min.size()];
2908 int j = 0;
2909 for (Integer i : min) {
2910 mvars[j++] = dim.v[i];
2911 }
2912 if (logger.isInfoEnabled()) {
2913 logger.info("extension for variables = " + Arrays.toString(mvars));
2914 }
2915 // reduce to dimension zero
2916 IdealWithUniv<Quotient<C>> Ext = extension(mvars);
2917 if (logger.isInfoEnabled()) {
2918 logger.info("extension = " + Ext);
2919 }
2920
2921 List<IdealWithUniv<Quotient<C>>> edec = Ext.ideal.zeroDimDecomposition();
2922 if (logger.isInfoEnabled()) {
2923 logger.info("0-dim irred decomp = " + edec);
2924 }
2925 // remove field extensions are not required
2926 // reconstruct dimension
2927 for (IdealWithUniv<Quotient<C>> ep : edec) {
2928 IdealWithUniv<C> cont = permContraction(ep);
2929 //System.out.println("cont = " + cont);
2930 dec.add(cont);
2931 }
2932 IdealWithUniv<C> extcont = permContraction(Ext);
2933 //System.out.println("extcont = " + extcont);
2934
2935 // get f
2936 List<GenPolynomial<C>> ql = extcont.others;
2937 if (ql.size() == 0) { // should not happen
2938 return dec;
2939 }
2940 GenPolynomial<C> fx = ql.get(0);
2941 //System.out.println("cont(Ext) fx = " + fx + ", " + fx.ring);
2942 if (fx.isONE()) {
2943 return dec;
2944 }
2945 Ideal<C> T = sum(fx);
2946 //System.out.println("T.rec = " + T.getList());
2947 if (T.isONE()) {
2948 logger.info("1 in ideal for " + fx);
2949 return dec;
2950 }
2951 if (logger.isInfoEnabled()) {
2952 logger.info("irred decomp ext-cont fx = " + fx);
2953 logger.info("recursion irred decomp T = " + T);
2954 }
2955 // recursion:
2956 List<IdealWithUniv<C>> Tdec = T.decomposition();
2957 if (logger.isInfoEnabled()) {
2958 logger.info("recursion irred decomposition = " + Tdec);
2959 }
2960 dec.addAll(Tdec);
2961 return dec;
2962 }
2963
2964
2965 /**
2966 * Ideal prime decompostition.
2967 * @return intersection of ideals G_i with ideal(this) subseteq cap_i(
2968 * ideal(G_i) ) and each G_i is a prime ideal and the G_i are
2969 * pairwise co-prime
2970 */
2971 public List<IdealWithUniv<C>> primeDecomposition() {
2972 // check dimension
2973 int z = commonZeroTest();
2974 List<IdealWithUniv<C>> dec = new ArrayList<IdealWithUniv<C>>();
2975 List<GenPolynomial<C>> ups = new ArrayList<GenPolynomial<C>>();
2976 // dimension -1
2977 if (z < 0) {
2978 //IdealWithUniv<C> id = new IdealWithUniv<C>(this, ups);
2979 //dec.add(id); see GB book
2980 return dec;
2981 }
2982 // dimension 0
2983 if (z == 0) {
2984 dec = zeroDimPrimeDecomposition();
2985 return dec;
2986 }
2987 // dimension > 0
2988 if (this.isZERO()) {
2989 return dec;
2990 }
2991 Dimension dim = dimension();
2992 if (logger.isInfoEnabled()) {
2993 logger.info("dimension = " + dim);
2994 }
2995
2996 // shortest maximal independent set
2997 Set<Set<Integer>> M = dim.M;
2998 Set<Integer> min = null;
2999 for (Set<Integer> m : M) {
3000 if (min == null) {
3001 min = m;
3002 continue;
3003 }
3004 if (m.size() < min.size()) {
3005 min = m;
3006 }
3007 }
3008 //System.out.println("min = " + min);
3009 String[] mvars = new String[min.size()];
3010 int j = 0;
3011 for (Integer i : min) {
3012 mvars[j++] = dim.v[i];
3013 }
3014 if (logger.isInfoEnabled()) {
3015 logger.info("extension for variables = " + Arrays.toString(mvars));
3016 }
3017 // reduce to dimension zero
3018 IdealWithUniv<Quotient<C>> Ext = extension(mvars);
3019 if (logger.isInfoEnabled()) {
3020 logger.info("extension = " + Ext);
3021 }
3022 List<IdealWithUniv<Quotient<C>>> edec = Ext.ideal.zeroDimPrimeDecomposition();
3023 if (logger.isInfoEnabled()) {
3024 logger.info("0-dim prime decomp = " + edec);
3025 }
3026 // remove field extensions, already done
3027 // reconstruct dimension
3028 for (IdealWithUniv<Quotient<C>> ep : edec) {
3029 IdealWithUniv<C> cont = permContraction(ep);
3030 //System.out.println("cont = " + cont);
3031 dec.add(cont);
3032 }
3033 // get f
3034 IdealWithUniv<C> extcont = permContraction(Ext);
3035 //System.out.println("extcont = " + extcont);
3036 List<GenPolynomial<C>> ql = extcont.others;
3037 if (ql.size() == 0) { // should not happen
3038 return dec;
3039 }
3040 GenPolynomial<C> fx = ql.get(0);
3041 //System.out.println("cont(Ext) fx = " + fx + ", " + fx.ring);
3042 if (fx.isONE()) {
3043 return dec;
3044 }
3045 // compute exponent not required
3046 Ideal<C> T = sum(fx);
3047 //System.out.println("T.rec = " + T.getList());
3048 if (T.isONE()) {
3049 logger.info("1 in ideal for " + fx);
3050 return dec;
3051 }
3052 if (logger.isInfoEnabled()) {
3053 logger.info("prime decomp ext-cont fx = " + fx);
3054 logger.info("recursion prime decomp T = " + T);
3055 }
3056 // recursion:
3057 List<IdealWithUniv<C>> Tdec = T.primeDecomposition();
3058 if (logger.isInfoEnabled()) {
3059 logger.info("recursion prime decomp = " + Tdec);
3060 }
3061 dec.addAll(Tdec);
3062 return dec;
3063 }
3064
3065
3066 /**
3067 * Test for ideal decompostition.
3068 * @param L intersection of ideals G_i with ideal(G) eq cap_i(ideal(G_i) )
3069 * @return true if L is a decomposition of this, else false
3070 */
3071 public boolean isDecomposition(List<IdealWithUniv<C>> L) {
3072 if (L == null || L.size() == 0) {
3073 if (this.isZERO()) {
3074 return true;
3075 } else {
3076 return false;
3077 }
3078 }
3079 GenPolynomialRing<C> ofac = list.ring;
3080 int r = ofac.nvar;
3081 int rp = L.get(0).ideal.list.ring.nvar;
3082 int d = rp - r;
3083 //System.out.println("d = " + d);
3084 Ideal<C> Id = this;
3085 if (d > 0) { // add lower variables
3086 GenPolynomialRing<C> nfac = ofac.extendLower(d);
3087 //System.out.println("nfac = " + nfac);
3088 List<GenPolynomial<C>> elist = new ArrayList<GenPolynomial<C>>(list.list.size());
3089 for (GenPolynomial<C> p : getList()) {
3090 //System.out.println("p = " + p);
3091 GenPolynomial<C> q = p.extendLower(nfac, 0, 0L);
3092 //System.out.println("q = " + q);
3093 elist.add(q);
3094 }
3095 Id = new Ideal<C>(nfac, elist, isGB, isTopt);
3096 }
3097
3098 // test if this is contained in the intersection
3099 for (IdealWithUniv<C> I : L) {
3100 boolean t = I.ideal.contains(Id);
3101 if (!t) {
3102 System.out.println("not contained " + this + " in " + I.ideal);
3103 return false;
3104 }
3105 }
3106 // // test if all univariate polynomials are contained in the respective ideal
3107 // for (IdealWithUniv<C> I : L) {
3108 // GenPolynomialRing<C> mfac = I.ideal.list.ring;
3109 // int i = 0;
3110 // for (GenPolynomial<C> p : I.upolys) {
3111 // GenPolynomial<C> pm = p.extendUnivariate(mfac, i++);
3112 // //System.out.println("pm = " + pm + ", p = " + p);
3113 // boolean t = I.ideal.contains(pm);
3114 // if (!t) {
3115 // System.out.println("not contained " + pm + " in " + I.ideal);
3116 // return false;
3117 // }
3118 // }
3119 // }
3120 return true;
3121 }
3122
3123
3124 /**
3125 * Ideal primary decompostition.
3126 * @return list of primary components of primary ideals G_i (pairwise
3127 * co-prime) with ideal(this) = cap_i( ideal(G_i) ) together with
3128 * the associated primes
3129 */
3130 public List<PrimaryComponent<C>> primaryDecomposition() {
3131 // check dimension
3132 int z = commonZeroTest();
3133 List<PrimaryComponent<C>> dec = new ArrayList<PrimaryComponent<C>>();
3134 List<GenPolynomial<C>> ups = new ArrayList<GenPolynomial<C>>();
3135 // dimension -1
3136 if (z < 0) {
3137 //IdealWithUniv<C> id = new IdealWithUniv<C>(this, ups);
3138 //PrimaryComponent<C> pc = new PrimaryComponent<C>(this, id);
3139 //dec.add(pc); see GB book
3140 return dec;
3141 }
3142 // dimension 0
3143 if (z == 0) {
3144 dec = zeroDimPrimaryDecomposition();
3145 return dec;
3146 }
3147 // dimension > 0
3148 if (this.isZERO()) {
3149 return dec;
3150 }
3151 Dimension dim = dimension();
3152 if (logger.isInfoEnabled()) {
3153 logger.info("dimension = " + dim);
3154 }
3155
3156 // shortest maximal independent set
3157 Set<Set<Integer>> M = dim.M;
3158 Set<Integer> min = null;
3159 for (Set<Integer> m : M) {
3160 if (min == null) {
3161 min = m;
3162 continue;
3163 }
3164 if (m.size() < min.size()) {
3165 min = m;
3166 }
3167 }
3168 //System.out.println("min = " + min);
3169 String[] mvars = new String[min.size()];
3170 int j = 0;
3171 for (Integer i : min) {
3172 mvars[j++] = dim.v[i];
3173 }
3174 if (logger.isInfoEnabled()) {
3175 logger.info("extension for variables = " + Arrays.toString(mvars));
3176 }
3177 // reduce to dimension zero
3178 IdealWithUniv<Quotient<C>> Ext = extension(mvars);
3179 if (logger.isInfoEnabled()) {
3180 logger.info("extension = " + Ext);
3181 }
3182
3183 List<PrimaryComponent<Quotient<C>>> edec = Ext.ideal.zeroDimPrimaryDecomposition();
3184 if (logger.isInfoEnabled()) {
3185 logger.info("0-dim primary decomp = " + edec);
3186 }
3187 // remove field extensions, already done
3188 // reconstruct dimension
3189 List<GenPolynomial<Quotient<C>>> upq = new ArrayList<GenPolynomial<Quotient<C>>>();
3190 for (PrimaryComponent<Quotient<C>> ep : edec) {
3191 IdealWithUniv<Quotient<C>> epu = new IdealWithUniv<Quotient<C>>(ep.primary, upq);
3192 IdealWithUniv<C> contq = permContraction(epu);
3193 IdealWithUniv<C> contp = permContraction(ep.prime);
3194 PrimaryComponent<C> pc = new PrimaryComponent<C>(contq.ideal, contp);
3195 //System.out.println("pc = " + pc);
3196 dec.add(pc);
3197 }
3198
3199 // get f
3200 IdealWithUniv<C> extcont = permContraction(Ext);
3201 if (debug) {
3202 logger.info("cont(Ext) = " + extcont);
3203 }
3204 List<GenPolynomial<C>> ql = extcont.others;
3205 if (ql.size() == 0) { // should not happen
3206 return dec;
3207 }
3208 GenPolynomial<C> fx = ql.get(0);
3209 //System.out.println("cont(Ext) fx = " + fx + ", " + fx.ring);
3210 if (fx.isONE()) {
3211 return dec;
3212 }
3213 // compute exponent
3214 int s = this.infiniteQuotientExponent(fx, extcont.ideal);
3215 if (s == 0) {
3216 logger.info("exponent is 0 ");
3217 return dec;
3218 }
3219 if (s > 1) {
3220 fx = Power.<GenPolynomial<C>> positivePower(fx, s);
3221 }
3222 if (debug) {
3223 logger.info("exponent fx = " + s + ", fx^s = " + fx);
3224 }
3225
3226 Ideal<C> T = sum(fx);
3227 //System.out.println("T.rec = " + T.getList());
3228 if (T.isONE()) {
3229 logger.info("1 in ideal for " + fx);
3230 return dec;
3231 }
3232 if (logger.isInfoEnabled()) {
3233 logger.info("primmary decomp ext-cont fx = " + fx);
3234 logger.info("recursion primary decomp T = " + T);
3235 }
3236 // recursion:
3237 List<PrimaryComponent<C>> Tdec = T.primaryDecomposition();
3238 if (logger.isInfoEnabled()) {
3239 logger.info("recursion primary decomp = " + Tdec);
3240 }
3241 dec.addAll(Tdec);
3242 return dec;
3243 }
3244
3245 }