001/* 002 * $Id$ 003 */ 004 005package edu.jas.gb; 006 007 008import java.util.ArrayList; 009import java.util.List; 010import java.util.ListIterator; 011 012import org.apache.logging.log4j.LogManager; 013import org.apache.logging.log4j.Logger; 014 015import edu.jas.poly.ExpVector; 016import edu.jas.poly.GenPolynomial; 017import edu.jas.poly.GenPolynomialRing; 018import edu.jas.structure.RingElem; 019import edu.jas.structure.NotInvertibleException; 020 021 022/** 023 * D-Groebner Base sequential algorithm. Implements D-Groebner bases and GB 024 * test. <b>Note:</b> Minimal reduced GBs are not unique. see BWK, section 10.1. 025 * @param <C> coefficient type 026 * @author Heinz Kredel 027 */ 028 029public class DGroebnerBaseSeq<C extends RingElem<C>> extends GroebnerBaseAbstract<C> { 030 031 032 private static final Logger logger = LogManager.getLogger(DGroebnerBaseSeq.class); 033 034 035 private static final boolean debug = logger.isDebugEnabled(); 036 037 038 /** 039 * Reduction engine. 040 */ 041 protected DReduction<C> dred; // shadow super.red ?? 042 043 044 /** 045 * Constructor. 046 */ 047 public DGroebnerBaseSeq() { 048 this(new DReductionSeq<C>()); 049 } 050 051 052 /** 053 * Constructor. 054 * @param dred D-Reduction engine 055 */ 056 public DGroebnerBaseSeq(DReduction<C> dred) { 057 super(dred); 058 this.dred = dred; 059 assert this.dred == super.red; 060 } 061 062 063 /** 064 * D-Groebner base test. 065 * @param modv module variable number. 066 * @param F polynomial list. 067 * @return true, if F is a D-Groebner base, else false. 068 */ 069 @Override 070 public boolean isGB(int modv, List<GenPolynomial<C>> F) { 071 GenPolynomial<C> pi, pj, s, d; 072 for (int i = 0; i < F.size(); i++) { 073 pi = F.get(i); 074 for (int j = i + 1; j < F.size(); j++) { 075 pj = F.get(j); 076 if (!dred.moduleCriterion(modv, pi, pj)) { 077 continue; 078 } 079 d = dred.GPolynomial(pi, pj); 080 if (!d.isZERO()) { 081 // better check for top reduction only 082 d = dred.normalform(F, d); 083 } 084 if (!d.isZERO()) { 085 System.out.println("d-pol(" + pi + "," + pj + ") != 0: " + d); 086 return false; 087 } 088 // works ok 089 if (!dred.criterion4(pi, pj)) { 090 continue; 091 } 092 s = dred.SPolynomial(pi, pj); 093 if (!s.isZERO()) { 094 s = dred.normalform(F, s); 095 } 096 if (!s.isZERO()) { 097 System.out.println("s-pol(" + i + "," + j + ") != 0: " + s); 098 return false; 099 } 100 } 101 } 102 return true; 103 } 104 105 106 /** 107 * D-Groebner base using pairlist class. 108 * @param modv module variable number. 109 * @param F polynomial list. 110 * @return GB(F) a D-Groebner base of F. 111 */ 112 public List<GenPolynomial<C>> GB(int modv, List<GenPolynomial<C>> F) { 113 List<GenPolynomial<C>> G = normalizeZerosOnes(F); 114 if (G.size() <= 1) { 115 return G; 116 } 117 GenPolynomialRing<C> ring = G.get(0).ring; 118 OrderedDPairlist<C> pairlist = new OrderedDPairlist<C>(modv, ring); 119 pairlist.put(G); 120 121 Pair<C> pair; 122 GenPolynomial<C> pi, pj, S, D, H; 123 while (pairlist.hasNext()) { 124 pair = pairlist.removeNext(); 125 //System.out.println("pair = " + pair); 126 if (pair == null) 127 continue; 128 pi = pair.pi; 129 pj = pair.pj; 130 if (debug) { 131 logger.debug("pi = {}", pi); 132 logger.debug("pj = {}", pj); 133 } 134 H = null; 135 136 // D-polynomial case ---------------------- 137 D = dred.GPolynomial(pi, pj); 138 //System.out.println("D_d = " + D); 139 if (!D.isZERO() && !dred.isTopReducible(G, D)) { 140 H = dred.normalform(G, D); 141 if (H.isONE()) { 142 G.clear(); 143 G.add(H); 144 return G; // since no threads are activated 145 } 146 if (!H.isZERO()) { 147 logger.info("Dred = {}", H); 148 //l++; 149 G.add(H); 150 pairlist.put(H); 151 } 152 } 153 154 // S-polynomial case ----------------------- 155 if (pair.getUseCriterion3() && pair.getUseCriterion4()) { 156 S = dred.SPolynomial(pi, pj); 157 //System.out.println("S_d = " + S); 158 if (S.isZERO()) { 159 pair.setZero(); 160 //continue; 161 } 162 if (logger.isDebugEnabled()) { 163 logger.debug("ht(S) = {}", S.leadingExpVector()); 164 } 165 H = dred.normalform(G, S); 166 if (H.isZERO()) { 167 pair.setZero(); 168 //continue; 169 } 170 if (logger.isDebugEnabled()) { 171 logger.debug("ht(H) = {}", H.leadingExpVector()); 172 } 173 if (H.isONE()) { 174 G.clear(); 175 G.add(H); 176 return G; // since no threads are activated 177 } 178 logger.debug("H = {}", H); 179 if (!H.isZERO()) { 180 logger.info("Sred = {}", H); 181 //len = G.size(); 182 //l++; 183 G.add(H); 184 pairlist.put(H); 185 } 186 } 187 } 188 logger.debug("#sequential list = {}", G.size()); 189 G = minimalGB(G); 190 logger.info("{}", pairlist); 191 return G; 192 } 193 194 195 /** 196 * Extended Groebner base using pairlist class. 197 * @param modv module variable number. 198 * @param F polynomial list. 199 * @return a container for an extended Groebner base of F. 200 */ 201 @Override 202 public ExtendedGB<C> extGB(int modv, List<GenPolynomial<C>> F) { 203 if (F == null || F.isEmpty()) { 204 throw new IllegalArgumentException("null or empty F not allowed"); 205 } 206 List<GenPolynomial<C>> G = new ArrayList<GenPolynomial<C>>(); 207 List<List<GenPolynomial<C>>> F2G = new ArrayList<List<GenPolynomial<C>>>(); 208 List<List<GenPolynomial<C>>> G2F = new ArrayList<List<GenPolynomial<C>>>(); 209 OrderedDPairlist<C> pairlist = null; 210 boolean oneInGB = false; 211 int len = F.size(); 212 213 List<GenPolynomial<C>> row = null; 214 List<GenPolynomial<C>> rows = null; 215 List<GenPolynomial<C>> rowh = null; 216 GenPolynomialRing<C> ring = null; 217 GenPolynomial<C> H; 218 GenPolynomial<C> p; 219 220 int nzlen = 0; 221 for (GenPolynomial<C> f : F) { 222 if (f.length() > 0) { 223 nzlen++; 224 } 225 if (ring == null) { 226 ring = f.ring; 227 } 228 } 229 GenPolynomial<C> one = ring.getONE(); 230 int k = 0; 231 ListIterator<GenPolynomial<C>> it = F.listIterator(); 232 while (it.hasNext()) { 233 p = it.next(); 234 if (p.length() > 0) { 235 row = blas.genVector(nzlen, null); 236 row.set(k, one); 237 k++; 238 if (p.isUnit()) { 239 G.clear(); 240 G.add(p); 241 G2F.clear(); 242 G2F.add(row); 243 oneInGB = true; 244 break; 245 } 246 G.add(p); 247 //logger.info("p row = {}", row); 248 G2F.add(row); 249 if (pairlist == null) { 250 //pairlist = strategy.create(modv, p.ring); 251 pairlist = new OrderedDPairlist<C>(modv, p.ring); 252 if (p.ring.coFac.isField()) { 253 //throw new RuntimeException("coefficients from a field"); 254 logger.warn("coefficients from a field " + p.ring.coFac); 255 } 256 } 257 // putOne not required 258 pairlist.put(p); 259 } else { 260 len--; 261 } 262 } 263 ExtendedGB<C> exgb; 264 if (len <= 1 || oneInGB) { 265 // adjust F2G 266 for (GenPolynomial<C> f : F) { 267 row = blas.genVector(G.size(), null); 268 H = dred.normalform(row, G, f); 269 if (!H.isZERO()) { 270 logger.error("nonzero H = {}", H); 271 } 272 logger.info("f row = {}", row); 273 F2G.add(row); 274 } 275 exgb = new ExtendedGB<C>(F, G, F2G, G2F); 276 //System.out.println("exgb #1 = " + exgb); 277 return exgb; 278 } 279 280 Pair<C> pair; 281 int i, j; 282 GenPolynomial<C> pi, pj; 283 GenPolynomial<C> S, D; 284 k = 0; 285 while (pairlist.hasNext() && !oneInGB) { 286 pair = pairlist.removeNext(); 287 if (pair == null) { 288 continue; 289 } 290 i = pair.i; 291 j = pair.j; 292 pi = pair.pi; 293 pj = pair.pj; 294 if (debug) { 295 logger.info("i, pi = {}, {}", i, pi); 296 logger.info("j, pj = {}, {}", j, pj); 297 } 298 H = null; 299 300 // D-polynomial case ---------------------- 301 rows = blas.genVector(G.size() + 1, null); 302 D = dred.GPolynomial(rows, i, pi, j, pj); 303 logger.info("Gpol = {}", D); 304 if (debug) { 305 logger.info("is reduction D = " + dred.isReductionNF(rows, G, D, ring.getZERO())); 306 } 307 if (!D.isZERO()) { //&& !dred.isTopReducible(G, D) 308 //continue; 309 rowh = blas.genVector(G.size() + 1, null); 310 H = dred.normalform(rowh, G, D); 311 if (debug) { 312 logger.info("is reduction H = " + dred.isReductionNF(rowh, G, D, H)); 313 } 314 //H = H.monic(); 315 int s = H.leadingBaseCoefficient().signum(); 316 if (s < 0) { 317 logger.info("negate: H_D rowd, rowh = {}, {}", rows, rowh); 318 H = H.negate(); 319 rows = blas.vectorNegate(rows); 320 rowh = blas.vectorNegate(rowh); 321 } 322 if (H.isONE()) { 323 // G.clear(); 324 G.add(H); 325 oneInGB = true; 326 } else if (!H.isZERO()) { 327 logger.info("H_G red = {}", H); 328 G.add(H); 329 pairlist.put(H); 330 } 331 //System.out.println("rowd = " + rows); 332 //System.out.println("rowh = " + rowh); 333 row = blas.vectorCombineRep(rows, rowh); 334 logger.debug("H_G row = {}", row); 335 if (!H.isZERO()) { 336 G2F.add(row); 337 } 338 if (debug) { 339 logger.debug("ht(H) = {}", H.leadingExpVector()); 340 logger.info("is reduction D,H = " + dred.isReductionNF(row, G, H, ring.getZERO())); 341 } 342 } 343 344 // S-polynomial case ---------------------- 345 rows = blas.genVector(G.size() + 1, null); 346 S = dred.SPolynomial(rows, i, pi, j, pj); 347 logger.info("Spol = {}", S); 348 if (debug) { 349 logger.info("is reduction S = " + dred.isReductionNF(rows, G, S, ring.getZERO())); 350 } 351 rowh = blas.genVector(G.size() + 1, null); 352 if (!S.isZERO()) { //&& !dred.isTopReducible(G, S) 353 //continue; 354 H = dred.normalform(rowh, G, S); 355 if (debug) { 356 logger.info("Spol_red = {}", H); 357 logger.info("is reduction H = " + dred.isReductionNF(rowh, G, S, H)); 358 } 359 //H = H.monic(); 360 int s = H.leadingBaseCoefficient().signum(); 361 if (s < 0) { 362 logger.info("negate: H_S rows, rowh = {}, {}", rows, rowh); 363 H = H.negate(); 364 //rowh = rowh.negate(); //rowh.set(G.size(), one.negate()); 365 rows = blas.vectorNegate(rows); 366 rowh = blas.vectorNegate(rowh); 367 } 368 //logger.info("Spol_red_norm = {}", H); 369 if (H.isONE()) { 370 //G.clear(); 371 G.add(H); 372 oneInGB = true; 373 } else if (!H.isZERO()) { 374 logger.info("H_S red = {}", H); 375 G.add(H); 376 pairlist.put(H); 377 } 378 //System.out.println("rows = " + rows); 379 //System.out.println("rowh = " + rowh); 380 row = blas.vectorCombineRep(rows, rowh); 381 logger.debug("H_S row = {}", row); 382 if (!H.isZERO()) { 383 G2F.add(row); 384 } 385 if (debug) { 386 logger.debug("ht(H) = {}", H.leadingExpVector()); 387 logger.info("is reduction S,H = " + dred.isReductionNF(row, G, H, ring.getZERO())); 388 } 389 } 390 } 391 if (true || debug) { 392 exgb = new ExtendedGB<C>(F, G, F2G, G2F); 393 boolean t = isReductionMatrix(exgb); 394 if (!t) { 395 logger.info("exgb unnorm = {}", exgb); 396 logger.info("exgb t_1 = {}", t); 397 } 398 } 399 G2F = normalizeMatrix(F.size(), G2F); 400 if (debug) { 401 exgb = new ExtendedGB<C>(F, G, F2G, G2F); 402 boolean t = isReductionMatrix(exgb); 403 if (!t) { 404 logger.info("exgb norm nonmin = {}", exgb); 405 logger.info("exgb t_2 = {}", t); 406 } 407 } 408 exgb = minimalExtendedGB(F.size(), G, G2F); 409 G = exgb.G; 410 G2F = exgb.G2F; 411 if (true || debug) { 412 exgb = new ExtendedGB<C>(F, G, F2G, G2F); 413 boolean t = isMinReductionMatrix(exgb); 414 if (!t) { 415 logger.info("exgb minGB = {}", exgb); 416 logger.info("exgb t_3 = {}", t); 417 } 418 } 419 exgb = new ExtendedGB<C>(F, G, F2G, G2F); 420 // setup matrices F and F2G 421 for (GenPolynomial<C> f : F) { 422 row = blas.genVector(G.size() + 1, null); 423 H = dred.normalform(row, G, f); 424 if (!H.isZERO()) { 425 logger.error("nonzero H, G = {}, {}", H, G); 426 throw new RuntimeException("H != 0"); 427 } 428 F2G.add(row); 429 } 430 exgb = new ExtendedGB<C>(F, G, F2G, G2F); 431 if (debug) { 432 boolean t = isMinReductionMatrix(exgb); 433 if (!t) { 434 logger.info("exgb +F+F2G = {}", exgb); 435 logger.info("exgb t_4 = {}", t); 436 } 437 } 438 return exgb; 439 } 440 441 442 /** 443 * Minimal extended groebner basis. 444 * @param flen length of rows. 445 * @param Gp a Groebner base. 446 * @param M a reduction matrix. 447 * @return a (partially) reduced Groebner base of Gp in a (fake) container. 448 */ 449 @Override 450 public ExtendedGB<C> minimalExtendedGB(int flen, List<GenPolynomial<C>> Gp, 451 List<List<GenPolynomial<C>>> M) { 452 if (Gp == null) { 453 return null; //new ExtendedGB<C>(null,Gp,null,M); 454 } 455 List<GenPolynomial<C>> G, F, T; 456 G = new ArrayList<GenPolynomial<C>>(Gp); 457 F = new ArrayList<GenPolynomial<C>>(Gp.size()); 458 459 List<List<GenPolynomial<C>>> Mg, Mf; 460 Mg = new ArrayList<List<GenPolynomial<C>>>(M.size()); 461 Mf = new ArrayList<List<GenPolynomial<C>>>(M.size()); 462 for (List<GenPolynomial<C>> r : M) { 463 // must be copied also 464 List<GenPolynomial<C>> row = new ArrayList<GenPolynomial<C>>(r); 465 Mg.add(row); 466 } 467 468 boolean mt; 469 ListIterator<GenPolynomial<C>> it; 470 ArrayList<Integer> ix = new ArrayList<Integer>(); 471 ArrayList<Integer> jx = new ArrayList<Integer>(); 472 int k = 0; 473 //System.out.println("flen, Gp, M = " + flen + ", " + Gp.size() + ", " + M.size() ); 474 GenPolynomialRing<C> pfac = null; 475 GenPolynomial<C> pi, pj; 476 ExpVector ei, ej; 477 C ai, aj, r; 478 while (G.size() > 0) { 479 pi = G.remove(0); 480 ei = pi.leadingExpVector(); 481 ai = pi.leadingBaseCoefficient(); 482 if (pfac == null) { 483 pfac = pi.ring; 484 } 485 486 it = G.listIterator(); 487 mt = false; 488 while (it.hasNext() && !mt) { 489 pj = it.next(); 490 ej = pj.leadingExpVector(); 491 mt = ei.multipleOf(ej); 492 if (mt) { 493 aj = pj.leadingBaseCoefficient(); 494 r = ai.remainder(aj); 495 mt = r.isZERO(); // && mt 496 } 497 } 498 it = F.listIterator(); 499 while (it.hasNext() && !mt) { 500 pj = it.next(); 501 ej = pj.leadingExpVector(); 502 mt = ei.multipleOf(ej); 503 if (mt) { 504 aj = pj.leadingBaseCoefficient(); 505 r = ai.remainder(aj); 506 mt = r.isZERO(); // && mt 507 } 508 } 509 T = new ArrayList<GenPolynomial<C>>(G); 510 T.addAll(F); 511 //GenPolynomial<C> t = dred.normalform(T, pi); 512 //System.out.println("t, mt, t==0: " + t + ", " + mt + ", " + t.isZERO()); 513 if (!mt) { 514 F.add(pi); 515 ix.add(k); 516 //System.out.println("ix: " + ix); 517 } else { // drop polynomial and corresponding row and column ?? 518 jx.add(k); 519 //System.out.println("jx: " + jx); 520 } 521 k++; 522 } 523 logger.info("ix, jx, #M = {}, {}, {}", ix, jx, Mg.size()); 524 int fix = -1; // copied polys 525 // copy Mg to Mf as indicated by ix 526 for (int i = 0; i < ix.size(); i++) { 527 int u = ix.get(i); 528 if (u >= flen && fix == -1) { 529 fix = Mf.size(); 530 } 531 //System.out.println("copy u_ix, fix = " + u + ", " + fix); 532 if (u >= 0) { 533 List<GenPolynomial<C>> row = Mg.get(u); 534 Mf.add(row); 535 } 536 } 537 // for (int i = 0; i < jx.size(); i++) { 538 // int u = jx.get(i); 539 // //System.out.println("copy u_jx = " + u); 540 // if (u >= 0) { 541 // List<GenPolynomial<C>> row = blas.genVector(flen,null); 542 // Mf.add(u, row); 543 // F.add(u, pfac.getZERO()); 544 // } 545 // } 546 if (F.size() <= 1 || fix == -1) { 547 return new ExtendedGB<C>(null, F, null, Mf); 548 } 549 // must return, since extended normalform has not correct order of polys 550 return new ExtendedGB<C>(null, F, null, Mf); 551 /* 552 G = F; 553 F = new ArrayList<GenPolynomial<C>>( G.size() ); 554 List<GenPolynomial<C>> temp; 555 k = 0; 556 final int len = G.size(); 557 while ( G.size() > 0 ) { 558 a = G.remove(0); 559 if ( k >= fix ) { // dont touch copied polys 560 row = Mf.get( k ); 561 //System.out.println("doing k = " + k + ", " + a); 562 // must keep order, but removed polys missing 563 temp = new ArrayList<GenPolynomial<C>>( len ); 564 temp.addAll( F ); 565 temp.add( a.ring.getZERO() ); // ?? 566 temp.addAll( G ); 567 //System.out.println("row before = " + row); 568 a = red.normalform( row, temp, a ); 569 //System.out.println("row after = " + row); 570 } 571 F.add( a ); 572 k++; 573 } 574 // does Mf need renormalization? 575 */ 576 } 577 578 579 /** 580 * Inverse for element modulo ideal. 581 * @param h polynomial 582 * @param F polynomial list 583 * @return inverse of h with respect to ideal(F), if defined 584 */ 585 public GenPolynomial<C> inverse(GenPolynomial<C> h, List<GenPolynomial<C>> F) { 586 if (h == null || h.isZERO()) { 587 throw new NotInvertibleException("zero not invertible"); 588 } 589 if (F == null || F.size() == 0) { 590 throw new NotInvertibleException("zero ideal"); 591 } 592 if (h.isUnit()) { 593 return h.inverse(); 594 } 595 // prepare by GB precomputation 596 List<GenPolynomial<C>> G = GB(F); 597 List<GenPolynomial<C>> I = new ArrayList<GenPolynomial<C>>(1 + G.size()); 598 I.add(h); 599 I.addAll(G); 600 // now compute extended gcd(h,G) 601 ExtendedGB<C> X = extGB(I); 602 List<GenPolynomial<C>> hG = X.G; 603 //System.out.println("hG = " + hG); 604 GenPolynomial<C> one = null; 605 int i = -1; 606 for (GenPolynomial<C> p : hG) { 607 i++; 608 if (p == null) { 609 continue; 610 } 611 if (p.isUnit()) { 612 one = p; 613 break; 614 } 615 } 616 if (one == null) { 617 throw new NotInvertibleException("h = " + h); 618 } 619 List<GenPolynomial<C>> row = X.G2F.get(i); // != -1 620 GenPolynomial<C> g = row.get(0); 621 if (g == null || g.isZERO()) { 622 throw new NotInvertibleException("h = " + h); 623 } 624 // adjust g to get g*h == 1 mod ideal(G) 625 GenPolynomial<C> f = g.multiply(h); 626 GenPolynomial<C> k = red.normalform(G, f); 627 //System.out.println("g = " + g + ", h = " + h + ", f = " + f + ", k = " + k); 628 if (k.signum() < 0) { // then is -1 or inv-G(0) 629 //System.out.println("k < 0: " + G); 630 g = g.sum(G.get(0)); //.negate(); 631 } 632 return g; 633 } 634 635}