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