001/* 002 * $Id$ 003 */ 004 005package edu.jas.gbufd; 006 007 008import java.util.ArrayList; 009import java.util.List; 010import java.util.Set; 011import java.util.SortedSet; 012 013import org.apache.logging.log4j.Logger; 014import org.apache.logging.log4j.LogManager; 015 016import edu.jas.gb.GroebnerBaseAbstract; 017import edu.jas.gb.ReductionAbstract; 018import edu.jas.poly.ExpVector; 019import edu.jas.poly.GenPolynomial; 020import edu.jas.poly.GenPolynomialRing; 021import edu.jas.poly.Monomial; 022import edu.jas.poly.PolyUtil; 023import edu.jas.poly.PolynomialList; 024import edu.jas.poly.TermOrder; 025import edu.jas.poly.TermOrderByName; 026import edu.jas.structure.GcdRingElem; 027import edu.jas.structure.RingFactory; 028 029 030/** 031 * Groebner Base sequential Groebner Walk algorithm. Implements Groebner base 032 * computation via Groebner Walk algorithm. See "The generic Groebner walk" by 033 * Fukuda, Jensen, Lauritzen, Thomas, 2005. 034 * 035 * The start term order t1 can be set by a constructor. The target term order t2 036 * is taken from the input polynomials term order. 037 * 038 * @param <C> coefficient type 039 * @author Heinz Kredel 040 * 041 * @see edu.jas.application.GBAlgorithmBuilder 042 * @see edu.jas.gbufd.GBFactory 043 */ 044public class GroebnerBaseWalk<C extends GcdRingElem<C>> extends GroebnerBaseAbstract<C> { 045 046 047 private static final Logger logger = LogManager.getLogger(GroebnerBaseWalk.class); 048 049 050 private static final boolean debug = logger.isDebugEnabled(); 051 052 053 /** 054 * The backing GB algorithm implementation. 055 */ 056 protected GroebnerBaseAbstract<C> sgb; 057 058 059 /** 060 * The start term order t1. 061 */ 062 //protected TermOrder startTO = TermOrderByName.IGRLEX.blockOrder(2); 063 protected TermOrder startTO = TermOrderByName.IGRLEX; 064 065 066 /** 067 * Print intermediate GB after this number of iterations. 068 */ 069 int iterPrint = 100; 070 071 072 /** 073 * Constructor. 074 */ 075 public GroebnerBaseWalk() { 076 super(); 077 sgb = null; 078 } 079 080 081 /** 082 * Constructor. 083 * @param coFac coefficient ring of polynomial ring. 084 */ 085 public GroebnerBaseWalk(RingFactory<C> coFac) { 086 this(GBFactory.<C> getImplementation(coFac)); 087 } 088 089 090 /** 091 * Constructor. 092 * @param coFac coefficient ring of polynomial ring. 093 * @param t1 start term order. 094 */ 095 public GroebnerBaseWalk(RingFactory<C> coFac, TermOrder t1) { 096 this(GBFactory.<C> getImplementation(coFac), t1); 097 } 098 099 100 /** 101 * Constructor. 102 * @param gb backing GB algorithm. 103 */ 104 public GroebnerBaseWalk(GroebnerBaseAbstract<C> gb) { 105 super(gb.red, gb.strategy); 106 sgb = gb; 107 } 108 109 110 /** 111 * Constructor. 112 * @param gb backing GB algorithm. 113 * @param t1 start term order. 114 */ 115 public GroebnerBaseWalk(GroebnerBaseAbstract<C> gb, TermOrder t1) { 116 this(gb); 117 startTO = t1; 118 } 119 120 121 /** 122 * Get the String representation with GB engine. 123 * @see java.lang.Object#toString() 124 */ 125 @Override 126 public String toString() { 127 if (sgb == null) { 128 return "GroebnerBaseWalk(" + startTO.toScript() + ")"; 129 } 130 return "GroebnerBaseWalk( " + sgb.toString() + ", " + startTO.toScript() + " )"; 131 } 132 133 134 /** 135 * Groebner base using Groebner Walk algorithm. 136 * @param modv module variable number. 137 * @param F polynomial list in target term order. 138 * @return GB(F) a INVLEX / target term order Groebner base of F. 139 */ 140 public List<GenPolynomial<C>> GB(int modv, List<GenPolynomial<C>> F) { 141 List<GenPolynomial<C>> G = normalizeZerosOnes(F); 142 G = PolyUtil.<C> monic(G); 143 if (G == null || G.size() == 0) { 144 return G; 145 } 146 GenPolynomialRing<C> pfac = G.get(0).ring; 147 if (!pfac.coFac.isField()) { 148 throw new IllegalArgumentException("coefficients not from a field: " + pfac.coFac); 149 } 150 if (G.size() <= 1) { 151 GenPolynomial<C> p = pfac.copy(G.get(0)); // change term order 152 G.clear(); 153 G.add(p); 154 return G; 155 } 156 // compute in graded / start term order 157 TermOrder grord = startTO; //new TermOrder(TermOrder.IGRLEX); 158 GenPolynomialRing<C> gfac = new GenPolynomialRing<C>(pfac, grord); 159 if (debug) { 160 logger.info("gfac = {}", gfac); //.toScript() 161 } 162 List<GenPolynomial<C>> Fp = gfac.copy(F); 163 logger.info("Term order: graded = {}, target = {}", grord, pfac.tord); 164 165 // compute graded / start term order Groebner base 166 if (sgb == null) { 167 sgb = GBFactory.<C> getImplementation(pfac.coFac, strategy); 168 } 169 List<GenPolynomial<C>> Gp = sgb.GB(modv, Fp); 170 logger.info("graded / start GB = {}", Gp); 171 if (grord.equals(pfac.tord)) { 172 return Gp; 173 } 174 if (Gp.size() == 1) { // also dimension -1 175 GenPolynomial<C> p = pfac.copy(Gp.get(0)); // change term order 176 G.clear(); 177 G.add(p); 178 return G; 179 } 180 // info on FGLM 181 int z = commonZeroTest(Gp); 182 if (z == 0) { 183 logger.info("ideal zero dimensional, can use also FGLM algorithm"); 184 } 185 // compute target term order Groebner base via Groebner Walk 186 G = walkGroebnerToTarget(modv, Gp, pfac); 187 return G; 188 } 189 190 191 /** 192 * Converts Groebner bases w.r.t. total degree / start term order to 193 * Groebner base w.r.t to inverse lexicographical / target term order. 194 * @param modv module variable number. 195 * @param Gl Groebner base with respect to graded / start term order. 196 * @param ufac target polynomial ring and term order. 197 * @return Groebner base w.r.t inverse lexicographical / target term order 198 */ 199 public List<GenPolynomial<C>> walkGroebnerToTarget(int modv, List<GenPolynomial<C>> Gl, 200 GenPolynomialRing<C> ufac) { 201 if (Gl == null || Gl.size() == 0) { 202 throw new IllegalArgumentException("G may not be null or empty"); 203 } 204 //Polynomial ring of input Groebner basis Gl 205 GenPolynomialRing<C> ring = Gl.get(0).ring; 206 if (debug) { 207 logger.info("ring = {}", ring); //.toScript() 208 } 209 //Polynomial ring of newGB with INVLEX / target term order 210 //TermOrder lexi = ufac.tord; //new TermOrder(TermOrder.INVLEX); 211 logger.info("G walk from ev1 = {} to ev2 = {}", ring.tord.toScript(), ufac.tord.toScript()); 212 List<GenPolynomial<C>> Giter = Gl; 213 // determine initial marks 214 List<ExpVector> marks = new ArrayList<ExpVector>(); 215 List<Monomial<C>> M = new ArrayList<Monomial<C>>(); 216 for (GenPolynomial<C> p : Giter) { 217 marks.add(p.leadingExpVector()); 218 M.add(new Monomial<C>(p.leadingMonomial())); 219 } 220 logger.info("marks = {}", marks); 221 List<Monomial<C>> Mp = new ArrayList<Monomial<C>>(M); 222 // weight matrix for target term order 223 long[][] ufweight = TermOrderByName.weightForOrder(ufac.tord, ring.nvar); 224 //TermOrder word = TermOrder.reverseWeight(ufweight); 225 TermOrder word = new TermOrder(ufweight); 226 logger.info("weight order: {}", word); 227 ufweight = word.getWeight(); // because of weightDeg usage 228 229 // loop through term orders 230 ExpVector w = null; 231 int iter = 0; // count #loops 232 boolean done = false; 233 while (!done) { 234 iter++; 235 // determine V and w 236 PolynomialList<C> Pl = new PolynomialList<C>(ring, Giter); 237 SortedSet<ExpVector> delta = Pl.deltaExpVectors(marks); 238 //logger.info("delta(marks) = {}", delta); 239 logger.info("w_old = {}", w); 240 ExpVector v = facetNormal(ring.tord, ufac.tord, delta, ring.evzero, ufweight); 241 logger.info("minimal v = {}", v); 242 if (v == null) { 243 done = true; 244 break; // finished 245 } 246 w = v; 247 // determine facet polynomials for w 248 List<GenPolynomial<C>> iG = new ArrayList<GenPolynomial<C>>(); 249 int i = 0; 250 for (GenPolynomial<C> f : Giter) { 251 ExpVector h = marks.get(i++); 252 GenPolynomial<C> ing = f.leadingFacetPolynomial(h, w); 253 if (debug) { 254 logger.info("ing_g = [{}], lt(ing) = {}, f = {}", ing, 255 ing.ring.toScript(ing.leadingExpVector()), 256 f.ring.toScript(f.leadingExpVector())); 257 } 258 iG.add(ing); 259 } 260 List<GenPolynomial<C>> inOmega = ufac.copy(iG); 261 if (debug) { 262 logger.info("inOmega = {}", inOmega); 263 logger.info("inOmega.ring: {}", inOmega.get(0).ring.toScript()); 264 } 265 266 // INVLEX / target term order GB of inOmega 267 List<GenPolynomial<C>> inOG = sgb.GB(modv, inOmega); 268 if (debug) { 269 logger.info("GB(inOmega) = {}", inOG); 270 } 271 // remark polynomials 272 marks.clear(); 273 M.clear(); 274 for (GenPolynomial<C> p : inOG) { 275 marks.add(p.leadingExpVector()); 276 M.add(new Monomial<C>(p.leadingMonomial())); 277 } 278 logger.info("new marks/M = {}", marks); 279 // lift and reduce 280 List<GenPolynomial<C>> G = liftReductas(M, Mp, Giter, inOG); 281 if (debug || (iter % iterPrint == 0)) { 282 logger.info("lift({}) inOG, new GB: {}", iter, G); 283 } 284 if (G.size() == 1) { // will not happen 285 GenPolynomial<C> p = ufac.copy(G.get(0)); // change term order 286 G.clear(); 287 G.add(p); 288 return G; 289 } 290 // iterate 291 Giter = G; 292 Mp.clear(); 293 Mp.addAll(M); 294 } 295 return Giter; 296 } 297 298 299 /** 300 * Determine new facet normal. 301 * @param t1 old term order. 302 * @param t2 new term order. 303 * @param delta exponent vectors deltas. 304 * @param zero exponent vector. 305 * @param t2weight weight representation of t2. 306 * @return new facet normal v or null if no new facet normal exists. 307 */ 308 public ExpVector facetNormal(TermOrder t1, TermOrder t2, Set<ExpVector> delta, ExpVector zero, 309 long[][] t2weight) { 310 TermOrder.EVComparator ev1 = t1.getAscendComparator(); 311 TermOrder.EVComparator ev2 = t2.getAscendComparator(); 312 ExpVector v = null; 313 long d = 0; // = Long.MIN_VALUE; 314 for (ExpVector e : delta) { 315 if (ev1.compare(zero, e) >= 0 || ev2.compare(zero, e) <= 0) { //ring.evzero 316 //logger.info("skip e = {}", e); 317 continue; 318 } 319 int s = 0; 320 long d2 = 0; 321 ExpVector vt = null; 322 ExpVector et = null; 323 if (v == null) { 324 v = e; 325 logger.info("init v = {}", v); 326 continue; 327 } 328 for (long[] tau : t2weight) { 329 //logger.info("current tau = {}", Arrays.toString(tau)); 330 //compare 331 d = v.weightDeg(tau); 332 d2 = e.weightDeg(tau); 333 vt = v.scalarMultiply(d2); 334 et = e.scalarMultiply(d); 335 s = ev1.compare(et, vt); 336 if (s == 0) { 337 continue; // next tau 338 } else if (s > 0) { // <, (> by example) 339 v = e; 340 break; 341 } else { 342 break; 343 } 344 } 345 //logger.info("step s = {}, d = {}, d2 = {}, e = {}, v = {}, et = {}, vt = {}", s, d, d2, e, 346 // v, et, vt); 347 } 348 return v; 349 } 350 351 352 /** 353 * Cleanup and terminate ThreadPool. 354 */ 355 @Override 356 public void terminate() { 357 if (sgb == null) { 358 return; 359 } 360 sgb.terminate(); 361 } 362 363 364 /** 365 * Cancel ThreadPool. 366 */ 367 @Override 368 public int cancel() { 369 if (sgb == null) { 370 return 0; 371 } 372 return sgb.cancel(); 373 } 374 375 376 /** 377 * Lift leading polynomials to full Groebner base with respect to term 378 * order. 379 * @param M new leading monomial list of polynomials as marks. 380 * @param Mp old leading monomial list of polynomials as marks. 381 * @param G Groebner base polynomials for lift. 382 * @param A polynomial list of leading omega polynomials to lift. 383 * @return lift(A) a Groebner base wrt M of ideal(G). 384 */ 385 public List<GenPolynomial<C>> liftReductas(List<Monomial<C>> M, List<Monomial<C>> Mp, 386 List<GenPolynomial<C>> G, List<GenPolynomial<C>> A) { 387 if (G == null || M == null || Mp == null || G.isEmpty()) { 388 throw new IllegalArgumentException("null or empty lists not allowed"); 389 } 390 if (A == null || A.isEmpty()) { 391 return A; 392 } 393 if (G.size() != Mp.size() || A.size() != M.size()) { 394 throw new IllegalArgumentException("equal sized lists required"); 395 } 396 GenPolynomial<C> pol = G.get(0); 397 if (pol == null) { 398 throw new IllegalArgumentException("null polynomial not allowed"); 399 } 400 // remove mark monomials 401 List<GenPolynomial<C>> Gp = new ArrayList<GenPolynomial<C>>(G.size()); 402 int i = 0; 403 int len = G.size(); 404 for (i = 0; i < len; i++) { 405 Monomial<C> mon = Mp.get(i); 406 GenPolynomial<C> s = G.get(i).subtract(mon); 407 Gp.add(s); 408 } 409 if (debug) { 410 logger.info("lifter GB: Gp = {}, Mp = {}", Gp, Mp); 411 } 412 // compute f^Gp for f in A 413 //GenPolynomialRing<C> oring = pol.ring; 414 logger.info("liftReductas: G = {}, Mp = {}", G.size(), Mp.size()); // ", old = {}", oring.toScript()); 415 List<GenPolynomial<C>> Ap = A; //oring.copy(A); 416 //logger.info("to lift Ap = {}", Ap); 417 ReductionAbstract<C> sred = (ReductionAbstract<C>) sgb.red; //new ReductionSeq<C>(); 418 List<GenPolynomial<C>> red = new ArrayList<GenPolynomial<C>>(); 419 GenPolynomialRing<C> tring = A.get(0).ring; 420 len = Ap.size(); 421 for (i = 0; i < len; i++) { 422 GenPolynomial<C> a = Ap.get(i); 423 GenPolynomial<C> r = sred.normalformMarked(Mp, Gp, a); 424 red.add(r); 425 } 426 logger.info("liftReductas: red(A) = {}", red.size()); 427 // combine f - f^Gp in tring 428 if (debug) { 429 logger.info("tring = {}", tring); //.toScript() + ", M = {}", M); 430 } 431 List<GenPolynomial<C>> nb = new ArrayList<GenPolynomial<C>>(red.size()); 432 for (i = 0; i < A.size(); i++) { 433 GenPolynomial<C> x = tring.copy(A.get(i)); // Ap? A! 434 GenPolynomial<C> r = tring.copy(red.get(i)); 435 GenPolynomial<C> s = x.subtract(r); 436 Monomial<C> m = M.get(i); 437 s.doAddTo(m.coefficient().negate(), m.exponent()); // remove marked term 438 if (!s.coefficient(m.exponent()).isZERO()) { 439 System.out.println("L-M: x = " + x + ", r = " + r); 440 throw new IllegalArgumentException("mark not removed: " + s + ", m = " + m); 441 } 442 nb.add(s); 443 } 444 if (debug) { 445 logger.info("lifted-M, nb = {}", nb.size()); 446 } 447 // minimal GB with preserved marks 448 //Collections.reverse(nb); // important for lex GB 449 len = nb.size(); 450 i = 0; 451 while (i < len) { 452 GenPolynomial<C> a = nb.remove(0); 453 Monomial<C> m = M.remove(0); // in step with element from nb 454 if (debug) { 455 logger.info("doing {}, lt = {}", a, tring.toScript(m.exponent())); 456 } 457 //a = sgb.red.normalform(nb, a); 458 a = sred.normalformMarked(M, nb, a); 459 if (debug) { 460 logger.info("done, a = {}, lt = {}", a, tring.toScript(a.leadingExpVector())); 461 } 462 nb.add(a); // adds as last 463 M.add(m); 464 i++; 465 } 466 // re-add mark after minimal 467 for (i = 0; i < len; i++) { 468 GenPolynomial<C> a = nb.get(i); 469 Monomial<C> m = M.get(i); 470 a.doAddTo(m.coefficient(), m.exponent()); // re-add marked term 471 nb.set(i, a); 472 } 473 logger.info("liftReductas: nb = {}, M = {}", nb.size(), M.size()); 474 //Collections.reverse(nb); // undo reverse 475 return nb; 476 } 477 478}