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