001/* 002 * $Id: GreatestCommonDivisorHensel.java 5047 2014-12-30 17:44:11Z kredel $ 003 */ 004 005package edu.jas.ufd; 006 007 008import java.util.ArrayList; 009import java.util.Iterator; 010import java.util.List; 011 012import org.apache.log4j.Logger; 013 014import edu.jas.arith.BigInteger; 015import edu.jas.arith.ModIntegerRing; 016import edu.jas.arith.ModLongRing; 017import edu.jas.arith.Modular; 018import edu.jas.arith.ModularRingFactory; 019import edu.jas.arith.PrimeList; 020import edu.jas.poly.ExpVector; 021import edu.jas.poly.GenPolynomial; 022import edu.jas.poly.GenPolynomialRing; 023import edu.jas.poly.PolyUtil; 024import edu.jas.structure.GcdRingElem; 025import edu.jas.structure.NotInvertibleException; 026import edu.jas.structure.Power; 027import edu.jas.structure.RingFactory; 028 029 030// import edu.jas.application.Ideal; 031 032 033/** 034 * Greatest common divisor algorithms with subresultant polynomial remainder 035 * sequence and univariate Hensel lifting. 036 * @author Heinz Kredel 037 */ 038 039public class GreatestCommonDivisorHensel<MOD extends GcdRingElem<MOD> & Modular> extends 040 GreatestCommonDivisorAbstract<BigInteger> { 041 042 043 private static final Logger logger = Logger.getLogger(GreatestCommonDivisorHensel.class); 044 045 046 private final boolean debug = logger.isDebugEnabled(); 047 048 049 /** 050 * Flag for linear or quadratic Hensel lift. 051 */ 052 public final boolean quadratic; 053 054 055 /** 056 * Fall back gcd algorithm. 057 */ 058 public final GreatestCommonDivisorAbstract<BigInteger> iufd; 059 060 061 /* 062 * Internal dispatcher. 063 */ 064 private final GreatestCommonDivisorAbstract<BigInteger> ufd; 065 066 067 /** 068 * Constructor. 069 */ 070 public GreatestCommonDivisorHensel() { 071 this(true); 072 } 073 074 075 /** 076 * Constructor. 077 * @param quadratic use quadratic Hensel lift. 078 */ 079 public GreatestCommonDivisorHensel(boolean quadratic) { 080 this.quadratic = quadratic; 081 iufd = new GreatestCommonDivisorSubres<BigInteger>(); 082 ufd = this; //iufd; 083 } 084 085 086 /** 087 * Univariate GenPolynomial greatest comon divisor. Uses univariate Hensel 088 * lifting. 089 * @param P univariate GenPolynomial. 090 * @param S univariate GenPolynomial. 091 * @return gcd(P,S). 092 */ 093 @Override 094 @SuppressWarnings("unchecked") 095 public GenPolynomial<BigInteger> baseGcd(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) { 096 if (S == null || S.isZERO()) { 097 return P; 098 } 099 if (P == null || P.isZERO()) { 100 return S; 101 } 102 if (P.ring.nvar > 1) { 103 throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial"); 104 } 105 GenPolynomialRing<BigInteger> fac = P.ring; 106 long e = P.degree(0); 107 long f = S.degree(0); 108 GenPolynomial<BigInteger> q; 109 GenPolynomial<BigInteger> r; 110 if (f > e) { 111 r = P; 112 q = S; 113 long g = f; 114 f = e; 115 e = g; 116 } else { 117 q = P; 118 r = S; 119 } 120 if (debug) { 121 logger.debug("degrees: e = " + e + ", f = " + f); 122 } 123 r = r.abs(); 124 q = q.abs(); 125 // compute contents and primitive parts 126 BigInteger a = baseContent(r); 127 BigInteger b = baseContent(q); 128 // gcd of coefficient contents 129 BigInteger c = gcd(a, b); // indirection 130 r = divide(r, a); // indirection 131 q = divide(q, b); // indirection 132 if (r.isONE()) { 133 return r.multiply(c); 134 } 135 if (q.isONE()) { 136 return q.multiply(c); 137 } 138 // compute normalization factor 139 BigInteger ac = r.leadingBaseCoefficient(); 140 BigInteger bc = q.leadingBaseCoefficient(); 141 BigInteger cc = gcd(ac, bc); // indirection 142 // compute degree vectors, only univeriate 143 ExpVector rdegv = r.degreeVector(); 144 ExpVector qdegv = q.degreeVector(); 145 //initialize prime list and degree vector 146 PrimeList primes = new PrimeList(PrimeList.Range.medium); 147 int pn = 50; //primes.size(); 148 149 ModularRingFactory<MOD> cofac; 150 GenPolynomial<MOD> qm; 151 GenPolynomial<MOD> qmf; 152 GenPolynomial<MOD> rm; 153 GenPolynomial<MOD> rmf; 154 GenPolynomial<MOD> cmf; 155 GenPolynomialRing<MOD> mfac; 156 GenPolynomial<MOD> cm = null; 157 GenPolynomial<MOD>[] ecm = null; 158 GenPolynomial<MOD> sm = null; 159 GenPolynomial<MOD> tm = null; 160 HenselApprox<MOD> lift = null; 161 if (debug) { 162 logger.debug("c = " + c); 163 logger.debug("cc = " + cc); 164 logger.debug("primes = " + primes); 165 } 166 167 int i = 0; 168 for (java.math.BigInteger p : primes) { 169 //System.out.println("next run ++++++++++++++++++++++++++++++++++"); 170 if (++i >= pn) { 171 logger.error("prime list exhausted, pn = " + pn); 172 //logger.info("primes = " + primes); 173 return iufd.baseGcd(P, S); 174 //throw new ArithmeticException("prime list exhausted"); 175 } 176 // initialize coefficient factory and map normalization factor 177 //cofac = new ModIntegerRing(p, true); 178 if (ModLongRing.MAX_LONG.compareTo(p) > 0) { 179 cofac = (ModularRingFactory) new ModLongRing(p, true); 180 } else { 181 cofac = (ModularRingFactory) new ModIntegerRing(p, true); 182 } 183 MOD nf = cofac.fromInteger(cc.getVal()); 184 if (nf.isZERO()) { 185 continue; 186 } 187 nf = cofac.fromInteger(q.leadingBaseCoefficient().getVal()); 188 if (nf.isZERO()) { 189 continue; 190 } 191 nf = cofac.fromInteger(r.leadingBaseCoefficient().getVal()); 192 if (nf.isZERO()) { 193 continue; 194 } 195 // initialize polynomial factory and map polynomials 196 mfac = new GenPolynomialRing<MOD>(cofac, fac.nvar, fac.tord, fac.getVars()); 197 qm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, q); 198 if (!qm.degreeVector().equals(qdegv)) { 199 continue; 200 } 201 rm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, r); 202 if (!rm.degreeVector().equals(rdegv)) { 203 continue; 204 } 205 if (debug) { 206 logger.info("cofac = " + cofac.getIntegerModul()); 207 } 208 209 // compute univariate modular gcd 210 cm = qm.gcd(rm); 211 212 // test for constant g.c.d 213 if (cm.isConstant()) { 214 logger.debug("cm, constant = " + cm); 215 return fac.getONE().multiply(c); 216 } 217 218 // compute factors and gcd with factor 219 GenPolynomial<BigInteger> crq; 220 rmf = rm.divide(cm); // rm = cm * rmf 221 ecm = cm.egcd(rmf); 222 if (ecm[0].isONE()) { 223 //logger.debug("gcd() first factor " + rmf); 224 crq = r; 225 cmf = rmf; 226 sm = ecm[1]; 227 tm = ecm[2]; 228 } else { 229 qmf = qm.divide(cm); // qm = cm * qmf 230 ecm = cm.egcd(qmf); 231 if (ecm[0].isONE()) { 232 //logger.debug("gcd() second factor " + qmf); 233 crq = q; 234 cmf = qmf; 235 sm = ecm[1]; 236 tm = ecm[2]; 237 } else { 238 logger.info("both gcd != 1: Hensel not applicable"); 239 return iufd.baseGcd(P, S); 240 } 241 } 242 BigInteger cn = crq.maxNorm(); 243 cn = cn.multiply(crq.leadingBaseCoefficient().abs()); 244 cn = cn.multiply(cn.fromInteger(2)); 245 if (debug) { 246 System.out.println("crq = " + crq); 247 System.out.println("cm = " + cm); 248 System.out.println("cmf = " + cmf); 249 System.out.println("sm = " + sm); 250 System.out.println("tm = " + tm); 251 System.out.println("cn = " + cn); 252 } 253 try { 254 if (quadratic) { 255 lift = HenselUtil.liftHenselQuadratic(crq, cn, cm, cmf, sm, tm); 256 } else { 257 lift = HenselUtil.liftHensel(crq, cn, cm, cmf, sm, tm); 258 } 259 } catch (NoLiftingException nle) { 260 logger.info("giving up on Hensel gcd reverting to Subres gcd " + nle); 261 return iufd.baseGcd(P, S); 262 } 263 q = lift.A; 264 if (debug) { 265 System.out.println("q = " + q); 266 System.out.println("qf = " + lift.B); 267 } 268 q = basePrimitivePart(q); 269 q = q.multiply(c).abs(); 270 if (PolyUtil.<BigInteger> baseSparsePseudoRemainder(P, q).isZERO() 271 && PolyUtil.<BigInteger> baseSparsePseudoRemainder(S, q).isZERO()) { 272 break; 273 } 274 logger.info("final devision not successfull"); 275 //System.out.println("P rem q = " + PolyUtil.<BigInteger>basePseudoRemainder(P,q)); 276 //System.out.println("S rem q = " + PolyUtil.<BigInteger>basePseudoRemainder(S,q)); 277 //break; 278 } 279 return q; 280 } 281 282 283 /** 284 * Univariate GenPolynomial recursive greatest comon divisor. Uses 285 * multivariate Hensel list. 286 * @param P univariate recursive GenPolynomial. 287 * @param S univariate recursive GenPolynomial. 288 * @return gcd(P,S). 289 */ 290 @Override 291 @SuppressWarnings("unchecked") 292 public GenPolynomial<GenPolynomial<BigInteger>> recursiveUnivariateGcd( 293 GenPolynomial<GenPolynomial<BigInteger>> P, GenPolynomial<GenPolynomial<BigInteger>> S) { 294 if (S == null || S.isZERO()) { 295 return P; 296 } 297 if (P == null || P.isZERO()) { 298 return S; 299 } 300 if (P.ring.nvar > 1) { 301 throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial"); 302 } 303 long e = P.degree(0); 304 long f = S.degree(0); 305 GenPolynomial<GenPolynomial<BigInteger>> q, r, s; 306 if (f > e) { 307 r = P; 308 q = S; 309 long g = f; 310 f = e; 311 e = g; 312 } else { 313 q = P; 314 r = S; 315 } 316 if (debug) { 317 logger.debug("degrees: e = " + e + ", f = " + f); 318 } 319 r = r.abs(); 320 q = q.abs(); 321 //logger.info("r: " + r + ", q: " + q); 322 323 GenPolynomial<BigInteger> a = ufd.recursiveContent(r); 324 GenPolynomial<BigInteger> b = ufd.recursiveContent(q); 325 326 GenPolynomial<BigInteger> c = ufd.gcd(a, b); // go to recursion 327 //System.out.println("rgcd c = " + c); 328 r = PolyUtil.<BigInteger> recursiveDivide(r, a); 329 q = PolyUtil.<BigInteger> recursiveDivide(q, b); 330 a = PolyUtil.<BigInteger> basePseudoDivide(a, c); // unused ? 331 b = PolyUtil.<BigInteger> basePseudoDivide(b, c); // unused ? 332 if (r.isONE()) { 333 return r.multiply(c); 334 } 335 if (q.isONE()) { 336 return q.multiply(c); 337 } 338 // check constant ldcf, TODO general case 339 GenPolynomial<BigInteger> la, lb, lc, lh; 340 la = r.leadingBaseCoefficient(); 341 lb = q.leadingBaseCoefficient(); 342 lc = ufd.gcd(la, lb); 343 //logger.info("la = " + la + ", lb = " + lb + ", lc = " + lc); 344 if (!lc.isConstant()) { 345 //continue; // easy way out 346 GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(r, q); 347 T = T.abs().multiply(c); 348 logger.info("non monic ldcf (" + lc + ") not implemented: " + T + "= gcd(" + r + "," + q + ") * " 349 + c); 350 return T; 351 } 352 353 // convert from Z[y1,...,yr][x] to Z[x][y1,...,yr] to Z[x,y1,...,yr] 354 GenPolynomial<GenPolynomial<BigInteger>> qs = PolyUtil.<BigInteger> switchVariables(q); 355 GenPolynomial<GenPolynomial<BigInteger>> rs = PolyUtil.<BigInteger> switchVariables(r); 356 357 GenPolynomialRing<GenPolynomial<BigInteger>> rfac = qs.ring; 358 RingFactory<GenPolynomial<BigInteger>> rrfac = rfac.coFac; 359 GenPolynomialRing<BigInteger> cfac = (GenPolynomialRing<BigInteger>) rrfac; 360 GenPolynomialRing<BigInteger> dfac = cfac.extend(rfac.getVars()); 361 //System.out.println("pfac = " + P.ring.toScript()); 362 //System.out.println("rfac = " + rfac.toScript()); 363 //System.out.println("dfac = " + dfac.toScript()); 364 GenPolynomial<BigInteger> qd = PolyUtil.<BigInteger> distribute(dfac, qs); 365 GenPolynomial<BigInteger> rd = PolyUtil.<BigInteger> distribute(dfac, rs); 366 367 // compute normalization factor 368 BigInteger ac = rd.leadingBaseCoefficient(); 369 BigInteger bc = qd.leadingBaseCoefficient(); 370 BigInteger cc = gcd(ac, bc); // indirection 371 372 //initialize prime list 373 PrimeList primes = new PrimeList(PrimeList.Range.medium); 374 Iterator<java.math.BigInteger> primeIter = primes.iterator(); 375 int pn = 50; //primes.size(); 376 377 // double check variables 378 // need qe,re,qd,rd,a,b 379 GenPolynomial<BigInteger> qe0, re0, ce0 = null; 380 381 for (int i = 0; i < 11; i++) { // meta loop 382 //System.out.println("======== run " + dfac.nvar + ", " + i); 383 java.math.BigInteger p = null; //new java.math.BigInteger("19"); //primes.next(); 384 // 5 small, 4 medium and 2 large size primes 385 if (i == 0) { // medium size 386 primes = new PrimeList(PrimeList.Range.medium); 387 primeIter = primes.iterator(); 388 } 389 if (i == 4) { // small size 390 primes = new PrimeList(PrimeList.Range.small); 391 primeIter = primes.iterator(); 392 p = primeIter.next(); // 2 393 p = primeIter.next(); // 3 394 p = primeIter.next(); // 5 395 p = primeIter.next(); // 7 396 } 397 if (i == 9) { // large size 398 primes = new PrimeList(PrimeList.Range.large); 399 primeIter = primes.iterator(); 400 } 401 ModularRingFactory<MOD> cofac = null; 402 int pi = 0; 403 while (pi < pn && primeIter.hasNext()) { 404 p = primeIter.next(); 405 //p = new java.math.BigInteger("19"); 406 logger.info("prime = " + p); 407 // initialize coefficient factory and map normalization factor and polynomials 408 ModularRingFactory<MOD> cf = null; 409 if (ModLongRing.MAX_LONG.compareTo(p) > 0) { 410 cf = (ModularRingFactory) new ModLongRing(p, true); 411 } else { 412 cf = (ModularRingFactory) new ModIntegerRing(p, true); 413 } 414 MOD nf = cf.fromInteger(cc.getVal()); 415 if (nf.isZERO()) { 416 continue; 417 } 418 nf = cf.fromInteger(q.leadingBaseCoefficient().leadingBaseCoefficient().getVal()); 419 if (nf.isZERO()) { 420 continue; 421 } 422 nf = cf.fromInteger(r.leadingBaseCoefficient().leadingBaseCoefficient().getVal()); 423 if (nf.isZERO()) { 424 continue; 425 } 426 cofac = cf; 427 break; 428 } 429 if (cofac == null) { // no lucky prime found 430 GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(q, r); 431 logger.info("no lucky prime, gave up on Hensel: " + T + "= gcd(" + r + "," + q + ")"); 432 return T.abs().multiply(c); //.abs(); 433 } 434 //System.out.println("cofac = " + cofac); 435 436 // search evaluation points and evaluate 437 List<BigInteger> V = new ArrayList<BigInteger>(P.ring.nvar); 438 GenPolynomialRing<BigInteger> ckfac = dfac; 439 GenPolynomial<BigInteger> qe = qd; 440 GenPolynomial<BigInteger> re = rd; 441 GenPolynomial<BigInteger> qei; 442 GenPolynomial<BigInteger> rei; 443 for (int j = dfac.nvar; j > 1; j--) { 444 // evaluation to univariate case 445 long degq = qe.degree(ckfac.nvar - 2); 446 long degr = re.degree(ckfac.nvar - 2); 447 ckfac = ckfac.contract(1); 448 long vi = 1L; //(long)(dfac.nvar-j); // 1L; 0 not so good for small p 449 if (p.longValue() > 1000L) { 450 //vi = (long)j+1L; 451 vi = 0L; 452 } 453 // search small evaluation point 454 while (true) { 455 MOD vp = cofac.fromInteger(vi++); 456 //System.out.println("vp = " + vp); 457 if (vp.isZERO() && vi != 1L) { // all elements of Z_p exhausted 458 qe = null; 459 re = null; 460 break; 461 } 462 BigInteger vii = new BigInteger(vi - 1); 463 qei = PolyUtil.<BigInteger> evaluateMain(ckfac, qe, vii); 464 rei = PolyUtil.<BigInteger> evaluateMain(ckfac, re, vii); 465 //System.out.println("qei = " + qei); 466 //System.out.println("rei = " + rei); 467 468 // check lucky evaluation point 469 if (degq != qei.degree(ckfac.nvar - 1)) { 470 //System.out.println("degv(qe) = " + qe.degreeVector()); 471 //System.out.println("deg(qe) = " + degq + ", deg(qe) = " + qei.degree(ckfac.nvar-1)); 472 continue; 473 } 474 if (degr != rei.degree(ckfac.nvar - 1)) { 475 //System.out.println("degv(re) = " + re.degreeVector()); 476 //System.out.println("deg(re) = " + degr + ", deg(re) = " + rei.degree(ckfac.nvar-1)); 477 continue; 478 } 479 V.add(vii); 480 qe = qei; 481 re = rei; 482 break; 483 } 484 if (qe == null && re == null) { 485 break; 486 } 487 } 488 if (qe == null && re == null) { 489 continue; 490 } 491 logger.info("evaluation points = " + V); 492 493 // recursion base: 494 GenPolynomial<BigInteger> ce = ufd.baseGcd(qe, re); 495 if (ce.isConstant()) { 496 return P.ring.getONE().multiply(c); 497 } 498 logger.info("base gcd = " + ce); 499 500 // double check 501 // need qe,re,qd,rd,a,b 502 if (i == 0) { 503 qe0 = qe; 504 re0 = re; 505 ce0 = ce; 506 continue; 507 } 508 long d0 = ce0.degree(0); 509 long d1 = ce.degree(0); 510 //System.out.println("d0, d1 = " + d0 + ", " + d1); 511 if (d1 < d0) { 512 qe0 = qe; 513 re0 = re; 514 ce0 = ce; 515 continue; 516 } else if (d1 > d0) { 517 continue; 518 } 519 // d0 == d1 is ok 520 long dx = r.degree(0); 521 //System.out.println("d0, dx = " + d0 + ", " + dx); 522 if (d0 == dx) { // gcd == r ? 523 if (PolyUtil.<BigInteger> recursivePseudoRemainder(q, r).isZERO()) { 524 r = r.abs().multiply(c); //.abs(); 525 logger.info("exit with r | q : " + r); 526 return r; 527 } 528 continue; 529 } 530 // norm 531 BigInteger mn = null; //mn = mn.multiply(cc).multiply(mn.fromInteger(2)); 532 // prepare lifting, chose factor polynomials 533 GenPolynomial<BigInteger> re1 = PolyUtil.<BigInteger> basePseudoDivide(re, ce); 534 GenPolynomial<BigInteger> qe1 = PolyUtil.<BigInteger> basePseudoDivide(qe, ce); 535 GenPolynomial<BigInteger> ui, he, pe; 536 GenPolynomial<BigInteger> g, gi, lui; 537 GenPolynomial<BigInteger> gcr, gcq; 538 gcr = ufd.baseGcd(re1, ce); 539 gcq = ufd.baseGcd(qe1, ce); 540 if (gcr.isONE() && gcq.isONE()) { // both gcds == 1: chose smaller ldcf 541 if (la.totalDegree() > lb.totalDegree()) { 542 ui = qd; 543 s = q; 544 he = qe1; 545 pe = qe; 546 BigInteger bn = qd.maxNorm(); 547 mn = bn.multiply(cc).multiply(new BigInteger(2L)); 548 g = lb; 549 logger.debug("select deg: ui = qd, g = b"); //, qe1 = " + qe1); // + ", qe = " + qe); 550 } else { 551 ui = rd; 552 s = r; 553 he = re1; 554 pe = re; 555 BigInteger an = rd.maxNorm(); 556 mn = an.multiply(cc).multiply(new BigInteger(2L)); 557 g = la; 558 logger.debug("select deg: ui = rd, g = a"); //, re1 = " + re1); // + ", re = " + re); 559 } 560 } else if (gcr.isONE()) { 561 ui = rd; 562 s = r; 563 he = re1; 564 pe = re; 565 BigInteger an = rd.maxNorm(); 566 mn = an.multiply(cc).multiply(new BigInteger(2L)); 567 g = la; 568 logger.debug("select: ui = rd, g = a"); //, re1 = " + re1); // + ", re = " + re); 569 } else if (gcq.isONE()) { 570 ui = qd; 571 s = q; 572 he = qe1; 573 pe = qe; 574 BigInteger bn = qd.maxNorm(); 575 mn = bn.multiply(cc).multiply(new BigInteger(2L)); 576 g = lb; 577 logger.debug("select: ui = qd, g = b"); //, qe1 = " + qe1); // + ", qe = " + qe); 578 } else { // both gcds != 1: method not applicable 579 logger.info("both gcds != 1: method not applicable"); 580 break; 581 } 582 lui = lc; //s.leadingBaseCoefficient(); 583 lh = PolyUtil.<BigInteger> basePseudoDivide(g, lui); 584 BigInteger ge = PolyUtil.<BigInteger> evaluateAll(g.ring.coFac, lui, V); 585 if (ge.isZERO()) { 586 continue; 587 } 588 BigInteger geh = PolyUtil.<BigInteger> evaluateAll(g.ring.coFac, lh, V); 589 if (geh.isZERO()) { 590 continue; 591 } 592 BigInteger gg = PolyUtil.<BigInteger> evaluateAll(g.ring.coFac, g, V); 593 if (gg.isZERO()) { 594 continue; 595 } 596 //System.out.println("ge = " + ge + ", geh = " + geh + ", gg = " + gg + ", pe = " + pe); 597 // 598 //ce = ce.multiply(geh); //ge); 599 // 600 he = he.multiply(ge); //gg); //ge); //geh); 601 // 602 gi = lui.extendLower(dfac, 0, 0L); //lui. // g. 603 // 604 ui = ui.multiply(gi); // gi !.multiply(gi) 605 //System.out.println("ui = " + ui + ", deg(ui) = " + ui.degreeVector()); 606 //System.out.println("ce = " + ce + ", he = " + he + ", ge = " + ge); 607 logger.info("gcd(ldcf): " + lui + ", ldcf cofactor: " + lh + ", base cofactor: " + he); 608 609 long k = Power.logarithm(new BigInteger(p), mn); 610 //System.out.println("mn = " + mn); 611 //System.out.println("k = " + k); 612 613 BigInteger qp = Power.positivePower(cofac.getIntegerModul(), k); // == p 614 ModularRingFactory<MOD> muqfac; 615 if (ModLongRing.MAX_LONG.compareTo(qp.getVal()) > 0) { 616 muqfac = (ModularRingFactory) new ModLongRing(qp.getVal(), true); // nearly a field 617 } else { 618 muqfac = (ModularRingFactory) new ModIntegerRing(qp.getVal(), true); // nearly a field 619 } 620 GenPolynomialRing<MOD> mucpfac = new GenPolynomialRing<MOD>(muqfac, ckfac); 621 //System.out.println("mucpfac = " + mucpfac.toScript()); 622 if (muqfac.fromInteger(ge.getVal()).isZERO()) { 623 continue; 624 } 625 //GenPolynomial<BigInteger> xxx = invertPoly(muqfac,lui,V); 626 //System.out.println("inv(lui) = " + xxx + ", muqfac = " + muqfac + ", lui = " + lui); 627 //ce = ce.multiply(xxx); //.leadingBaseCoefficient()); 628 //xxx = invertPoly(muqfac,lh,V); 629 //System.out.println("inv(lh) = " + xxx + ", muqfac = " + muqfac + ", lh = " + lh); 630 //he = he.multiply(xxx); //.leadingBaseCoefficient()); 631 632 GenPolynomial<MOD> cm = PolyUtil.<MOD> fromIntegerCoefficients(mucpfac, ce); 633 GenPolynomial<MOD> hm = PolyUtil.<MOD> fromIntegerCoefficients(mucpfac, he); 634 if (cm.degree(0) != ce.degree(0) || hm.degree(0) != he.degree(0)) { 635 continue; 636 } 637 if (cm.isZERO() || hm.isZERO()) { 638 continue; 639 } 640 logger.info("univariate modulo p^k: " + cm + ", " + hm); 641 642 // convert C from Z[...] to Z_q[...] 643 GenPolynomialRing<MOD> qcfac = new GenPolynomialRing<MOD>(muqfac, dfac); 644 GenPolynomial<MOD> uq = PolyUtil.<MOD> fromIntegerCoefficients(qcfac, ui); 645 if (!ui.leadingExpVector().equals(uq.leadingExpVector())) { 646 logger.info("ev(ui) = " + ui.leadingExpVector() + ", ev(uq) = " + uq.leadingExpVector()); 647 continue; 648 } 649 logger.info("multivariate modulo p^k: " + uq); 650 651 List<GenPolynomial<MOD>> F = new ArrayList<GenPolynomial<MOD>>(2); 652 F.add(cm); 653 F.add(hm); 654 List<GenPolynomial<BigInteger>> G = new ArrayList<GenPolynomial<BigInteger>>(2); 655 G.add(lui.ring.getONE()); //lui: lui.ring.getONE()); // TODO 656 G.add(lui.ring.getONE()); //lh: lui); 657 List<GenPolynomial<MOD>> lift; 658 try { 659 //lift = HenselMultUtil.<MOD> liftHenselFull(ui, F, V, k, G); 660 lift = HenselMultUtil.<MOD> liftHensel(ui, uq, F, V, k, G); 661 logger.info("lift = " + lift); 662 } catch (NoLiftingException nle) { 663 logger.info("NoLiftingException"); 664 //System.out.println("exception : " + nle); 665 continue; 666 } catch (ArithmeticException ae) { 667 logger.info("ArithmeticException"); 668 //System.out.println("exception : " + ae); 669 continue; 670 } catch (NotInvertibleException ni) { 671 logger.info("NotInvertibleException"); 672 //System.out.println("exception : " + ni); 673 continue; 674 } 675 //if (!HenselMultUtil.<MOD> isHenselLift(ui, uq, F, k, lift)) { // not meaningfull test 676 // logger.info("isHenselLift: false"); 677 // //continue; 678 //} 679 680 // convert Ci from Z_{p^k}[x,y1,...,yr] to Z[x,y1,...,yr] to Z[x][y1,...,yr] to Z[y1,...,yr][x] 681 GenPolynomial<BigInteger> ci = PolyUtil.integerFromModularCoefficients(dfac, lift.get(0)); 682 ci = basePrimitivePart(ci); 683 GenPolynomial<GenPolynomial<BigInteger>> Cr = PolyUtil.<BigInteger> recursive(rfac, ci); 684 GenPolynomial<GenPolynomial<BigInteger>> Cs = PolyUtil.<BigInteger> switchVariables(Cr); 685 if (!Cs.ring.equals(P.ring)) { 686 System.out.println("Cs.ring = " + Cs.ring + ", P.ring = " + P.ring); 687 } 688 GenPolynomial<GenPolynomial<BigInteger>> Q = ufd.recursivePrimitivePart(Cs); 689 Q = ufd.baseRecursivePrimitivePart(Q); 690 Q = Q.abs().multiply(c); //.abs(); 691 GenPolynomial<GenPolynomial<BigInteger>> Pq, Sq; 692 Pq = PolyUtil.<BigInteger> recursivePseudoRemainder(P, Q); 693 Sq = PolyUtil.<BigInteger> recursivePseudoRemainder(S, Q); 694 if (Pq.isZERO() && Sq.isZERO()) { 695 logger.info("gcd normal exit: " + Q); 696 return Q; 697 } 698 logger.info("bad Q = " + Q); // + ", Pq = " + Pq + ", Sq = " + Sq); 699 } // end for meta loop 700 // Hensel gcd failed 701 GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(r, q); 702 T = T.abs().multiply(c); 703 logger.info("no lucky prime or evaluation points, gave up on Hensel: " + T + "= gcd(" + r + "," + q 704 + ")"); 705 return T; 706 } 707 708 709 GenPolynomial<BigInteger> invertPoly(ModularRingFactory<MOD> mfac, GenPolynomial<BigInteger> li, 710 List<BigInteger> V) { 711 if (li == null || li.isZERO()) { 712 throw new RuntimeException("li not invertible: " + li); 713 } 714 if (li.isONE()) { 715 return li; 716 } 717 //System.out.println("mfac = " + mfac + ", V = " + V +", li = " + li); 718 GenPolynomialRing<BigInteger> pfac = li.ring; 719 GenPolynomialRing<MOD> mpfac = new GenPolynomialRing<MOD>(mfac, pfac); 720 GenPolynomial<MOD> lm = PolyUtil.<MOD> fromIntegerCoefficients(mpfac, li); 721 //System.out.println("pfac = " + pfac + ", lm = " + lm); 722 List<GenPolynomial<MOD>> lid = new ArrayList<GenPolynomial<MOD>>(V.size()); 723 int i = 0; 724 for (BigInteger bi : V) { 725 MOD m = mfac.fromInteger(bi.getVal()); 726 GenPolynomial<MOD> mp = mpfac.univariate(i); 727 mp = mp.subtract(m); // X_i - v_i 728 lid.add(mp); 729 i++; 730 } 731 //System.out.println("lid = " + lid); 732 //Ideal<MOD> id = new Ideal<MOD>(mpfac,lid,true); // is a GB 733 //System.out.println("id = " + id); 734 GenPolynomial<MOD> mi = lm; //id.inverse(lm); 735 //System.out.println("mi = " + mi); 736 GenPolynomial<BigInteger> inv = PolyUtil.integerFromModularCoefficients(pfac, mi); 737 return inv; 738 } 739 740}