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