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