001 /*
002 * $Id: GreatestCommonDivisorHensel.java 3718 2011-08-03 19:58:03Z kredel $
003 */
004
005 package edu.jas.ufd;
006
007 import java.util.ArrayList;
008 import java.util.List;
009 import java.util.Iterator;
010
011 import org.apache.log4j.Logger;
012
013 import edu.jas.arith.BigInteger;
014 import edu.jas.arith.ModIntegerRing;
015 import edu.jas.arith.ModLongRing;
016 import edu.jas.arith.Modular;
017 import edu.jas.arith.ModularRingFactory;
018 import edu.jas.arith.PrimeList;
019 import edu.jas.poly.ExpVector;
020 import edu.jas.poly.GenPolynomial;
021 import edu.jas.poly.GenPolynomialRing;
022 import edu.jas.poly.PolyUtil;
023 import edu.jas.structure.Power;
024 import edu.jas.structure.GcdRingElem;
025 import edu.jas.structure.RingFactory;
026
027
028 /**
029 * Greatest common divisor algorithms with subresultant polynomial remainder
030 * sequence and univariate Hensel lifting.
031 * @author Heinz Kredel
032 */
033
034 public class GreatestCommonDivisorHensel<MOD extends GcdRingElem<MOD> & Modular> extends
035 GreatestCommonDivisorAbstract<BigInteger> {
036
037
038 private static final Logger logger = Logger.getLogger(GreatestCommonDivisorHensel.class);
039
040
041 private final boolean debug = /*true ||*/logger.isDebugEnabled();
042
043
044 /**
045 * Flag for linear or quadratic Hensel lift.
046 */
047 public final boolean quadratic;
048
049
050 /**
051 * Fall back gcd algorithm.
052 */
053 public final GreatestCommonDivisorAbstract<BigInteger> iufd;
054
055
056 /**
057 * Constructor.
058 */
059 public GreatestCommonDivisorHensel() {
060 this(true);
061 }
062
063
064 /**
065 * Constructor.
066 * @param quadratic use quadratic Hensel lift.
067 */
068 public GreatestCommonDivisorHensel(boolean quadratic) {
069 this.quadratic = quadratic;
070 iufd = new GreatestCommonDivisorSubres<BigInteger>();
071 }
072
073
074 /**
075 * Univariate GenPolynomial greatest comon divisor. Uses univariate Hensel
076 * lifting.
077 * @param P univariate GenPolynomial.
078 * @param S univariate GenPolynomial.
079 * @return gcd(P,S).
080 */
081 @Override
082 public GenPolynomial<BigInteger> baseGcd(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) {
083 if (S == null || S.isZERO()) {
084 return P;
085 }
086 if (P == null || P.isZERO()) {
087 return S;
088 }
089 if (P.ring.nvar > 1) {
090 throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
091 }
092 GenPolynomialRing<BigInteger> fac = P.ring;
093 long e = P.degree(0);
094 long f = S.degree(0);
095 GenPolynomial<BigInteger> q;
096 GenPolynomial<BigInteger> r;
097 if (f > e) {
098 r = P;
099 q = S;
100 long g = f;
101 f = e;
102 e = g;
103 } else {
104 q = P;
105 r = S;
106 }
107 r = r.abs();
108 q = q.abs();
109 // compute contents and primitive parts
110 BigInteger a = baseContent(r);
111 BigInteger b = baseContent(q);
112 // gcd of coefficient contents
113 BigInteger c = gcd(a, b); // indirection
114 r = divide(r, a); // indirection
115 q = divide(q, b); // indirection
116 if (r.isONE()) {
117 return r.multiply(c);
118 }
119 if (q.isONE()) {
120 return q.multiply(c);
121 }
122 // compute normalization factor
123 BigInteger ac = r.leadingBaseCoefficient();
124 BigInteger bc = q.leadingBaseCoefficient();
125 BigInteger cc = gcd(ac, bc); // indirection
126 // compute degree vectors, only univeriate
127 ExpVector rdegv = r.degreeVector();
128 ExpVector qdegv = q.degreeVector();
129 //initialize prime list and degree vector
130 PrimeList primes = new PrimeList(PrimeList.Range.medium);
131 int pn = 50; //primes.size();
132
133 ModularRingFactory<MOD> cofac;
134 GenPolynomial<MOD> qm;
135 GenPolynomial<MOD> qmf;
136 GenPolynomial<MOD> rm;
137 GenPolynomial<MOD> rmf;
138 GenPolynomial<MOD> cmf;
139 GenPolynomialRing<MOD> mfac;
140 GenPolynomial<MOD> cm = null;
141 GenPolynomial<MOD>[] ecm = null;
142 GenPolynomial<MOD> sm = null;
143 GenPolynomial<MOD> tm = null;
144 HenselApprox<MOD> lift = null;
145 if (debug) {
146 logger.debug("c = " + c);
147 logger.debug("cc = " + cc);
148 logger.debug("primes = " + primes);
149 }
150
151 int i = 0;
152 for (java.math.BigInteger p : primes) {
153 //System.out.println("next run ++++++++++++++++++++++++++++++++++");
154 if (++i >= pn) {
155 logger.error("prime list exhausted, pn = " + pn);
156 logger.info("primes = " + primes);
157 return iufd.baseGcd(P, S);
158 //throw new ArithmeticException("prime list exhausted");
159 }
160 // initialize coefficient factory and map normalization factor
161 //cofac = new ModIntegerRing(p, true);
162 if (ModLongRing.MAX_LONG.compareTo(p) > 0) {
163 cofac = (ModularRingFactory) new ModLongRing(p, true);
164 } else {
165 cofac = (ModularRingFactory) new ModIntegerRing(p, true);
166 }
167 MOD nf = cofac.fromInteger(cc.getVal());
168 if (nf.isZERO()) {
169 continue;
170 }
171 nf = cofac.fromInteger(q.leadingBaseCoefficient().getVal());
172 if (nf.isZERO()) {
173 continue;
174 }
175 nf = cofac.fromInteger(r.leadingBaseCoefficient().getVal());
176 if (nf.isZERO()) {
177 continue;
178 }
179 // initialize polynomial factory and map polynomials
180 mfac = new GenPolynomialRing<MOD>(cofac, fac.nvar, fac.tord, fac.getVars());
181 qm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, q);
182 if (!qm.degreeVector().equals(qdegv)) {
183 continue;
184 }
185 rm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, r);
186 if (!rm.degreeVector().equals(rdegv)) {
187 continue;
188 }
189 if (debug) {
190 logger.info("cofac = " + cofac.getIntegerModul());
191 }
192
193 // compute univariate modular gcd
194 cm = qm.gcd(rm);
195
196 // test for constant g.c.d
197 if (cm.isConstant()) {
198 logger.debug("cm, constant = " + cm);
199 return fac.getONE().multiply(c);
200 }
201
202 // compute factors and gcd with factor
203 GenPolynomial<BigInteger> crq;
204 rmf = rm.divide(cm); // rm = cm * rmf
205 ecm = cm.egcd(rmf);
206 if (ecm[0].isONE()) {
207 //logger.debug("gcd() first factor " + rmf);
208 crq = r;
209 cmf = rmf;
210 sm = ecm[1];
211 tm = ecm[2];
212 } else {
213 qmf = qm.divide(cm); // qm = cm * qmf
214 ecm = cm.egcd(qmf);
215 if (ecm[0].isONE()) {
216 //logger.debug("gcd() second factor " + qmf);
217 crq = q;
218 cmf = qmf;
219 sm = ecm[1];
220 tm = ecm[2];
221 } else {
222 logger.info("giving up on Hensel gcd reverting to Subres gcd");
223 return iufd.baseGcd(P, S);
224 }
225 }
226 BigInteger cn = crq.maxNorm();
227 cn = cn.multiply(crq.leadingBaseCoefficient().abs());
228 cn = cn.multiply(cn.fromInteger(2));
229 if (debug) {
230 System.out.println("crq = " + crq);
231 System.out.println("cm = " + cm);
232 System.out.println("cmf = " + cmf);
233 System.out.println("sm = " + sm);
234 System.out.println("tm = " + tm);
235 System.out.println("cn = " + cn);
236 }
237 try {
238 if (quadratic) {
239 lift = HenselUtil.liftHenselQuadratic(crq, cn, cm, cmf, sm, tm);
240 } else {
241 lift = HenselUtil.liftHensel(crq, cn, cm, cmf, sm, tm);
242 }
243 } catch (NoLiftingException nle) {
244 logger.info("giving up on Hensel gcd reverting to Subres gcd " + nle);
245 return iufd.baseGcd(P, S);
246 }
247 q = lift.A;
248 if (debug) {
249 System.out.println("q = " + q);
250 System.out.println("qf = " + lift.B);
251 }
252 q = basePrimitivePart(q);
253 q = q.multiply(c).abs();
254 if (PolyUtil.<BigInteger> basePseudoRemainder(P, q).isZERO()
255 && PolyUtil.<BigInteger> basePseudoRemainder(S, q).isZERO()) {
256 break;
257 } else { // else should not happen at this point
258 logger.info("final devision not successfull");
259 //System.out.println("P rem q = " + PolyUtil.<BigInteger>basePseudoRemainder(P,q));
260 //System.out.println("S rem q = " + PolyUtil.<BigInteger>basePseudoRemainder(S,q));
261 //break;
262 }
263 }
264 return q;
265 }
266
267
268 /**
269 * Univariate GenPolynomial recursive greatest comon divisor. Uses
270 * multivariate Hensel list.
271 * @param P univariate recursive GenPolynomial.
272 * @param S univariate recursive GenPolynomial.
273 * @return gcd(P,S).
274 */
275 @Override
276 public GenPolynomial<GenPolynomial<BigInteger>> recursiveUnivariateGcd(
277 GenPolynomial<GenPolynomial<BigInteger>> P,
278 GenPolynomial<GenPolynomial<BigInteger>> S) {
279 if (S == null || S.isZERO()) {
280 return P;
281 }
282 if (P == null || P.isZERO()) {
283 return S;
284 }
285 if (P.ring.nvar > 1) {
286 throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
287 }
288 long e = P.degree(0);
289 long f = S.degree(0);
290 GenPolynomial<GenPolynomial<BigInteger>> q;
291 GenPolynomial<GenPolynomial<BigInteger>> r;
292 if (f > e) {
293 r = P;
294 q = S;
295 long g = f;
296 f = e;
297 e = g;
298 } else {
299 q = P;
300 r = S;
301 }
302 r = r.abs();
303 q = q.abs();
304 GenPolynomial<BigInteger> a = recursiveContent(r);
305 GenPolynomial<BigInteger> b = recursiveContent(q);
306
307 GenPolynomial<BigInteger> c = gcd(a, b); // go to recursion
308 //System.out.println("rgcd c = " + c);
309 r = PolyUtil.<BigInteger> recursiveDivide(r, a);
310 q = PolyUtil.<BigInteger> recursiveDivide(q, b);
311 a = PolyUtil.<BigInteger> basePseudoDivide(a, c);
312 b = PolyUtil.<BigInteger> basePseudoDivide(b, c);
313 if (r.isONE()) {
314 return r.multiply(c);
315 }
316 if (q.isONE()) {
317 return q.multiply(c);
318 }
319 // convert from Z[y1,...,yr][x] to Z[x][y1,...,yr] to Z[x,y1,...,yr]
320 GenPolynomial<GenPolynomial<BigInteger>> qs = PolyUtil.<BigInteger> switchVariables(q);
321 GenPolynomial<GenPolynomial<BigInteger>> rs = PolyUtil.<BigInteger> switchVariables(r);
322
323 GenPolynomialRing<GenPolynomial<BigInteger>> rfac = qs.ring;
324 RingFactory<GenPolynomial<BigInteger>> rrfac = rfac.coFac;
325 GenPolynomialRing<BigInteger> cfac = (GenPolynomialRing<BigInteger>) rrfac;
326 GenPolynomialRing<BigInteger> dfac = cfac.extend(rfac.getVars());
327 //System.out.println("pfac = " + P.ring.toScript());
328 //System.out.println("rfac = " + rfac.toScript());
329 //System.out.println("dfac = " + dfac.toScript());
330
331 GenPolynomial<BigInteger> qd = PolyUtil.<BigInteger> distribute(dfac, qs);
332 GenPolynomial<BigInteger> rd = PolyUtil.<BigInteger> distribute(dfac, rs);
333 //System.out.println("qd = " + qd);
334 //System.out.println("rd = " + rd);
335 //System.out.println("a = " + a);
336 //System.out.println("b = " + b);
337 //System.out.println("c = " + c);
338
339 // compute normalization factor
340 BigInteger ac = rd.leadingBaseCoefficient();
341 BigInteger bc = qd.leadingBaseCoefficient();
342 BigInteger cc = gcd(ac, bc); // indirection
343
344 //initialize prime list
345 PrimeList primes = new PrimeList(PrimeList.Range.medium); // PrimeList.Range.medium);
346 Iterator<java.math.BigInteger> primeIter = primes.iterator();
347 int pn = 50; //primes.size();
348
349 // double check variables
350 // need qe,re,qd,rd,a,b
351 GenPolynomial<MOD> qe0;
352 GenPolynomial<MOD> re0;
353 GenPolynomial<MOD> ce0 = null;
354
355 for ( int i = 0; i < 11; i++ ) { // meta loop
356 //System.out.println("======================================================= run "
357 // + dfac.nvar + ", " + i);
358 java.math.BigInteger p = null; //new java.math.BigInteger("19"); //primes.next();
359 // 5 small, 5 medium and 1 large size primes
360 if ( i == 0 ) { // medium size
361 primes = new PrimeList(PrimeList.Range.medium);
362 primeIter = primes.iterator();
363 }
364 if ( i == 5 ) { // smal size
365 primes = new PrimeList(PrimeList.Range.small);
366 primeIter = primes.iterator();
367 p = primeIter.next(); // 2
368 p = primeIter.next(); // 3
369 p = primeIter.next(); // 5
370 p = primeIter.next(); // 7
371 }
372 if ( i == 10 ) { // large size
373 primes = new PrimeList(PrimeList.Range.large);
374 primeIter = primes.iterator();
375 }
376 ModularRingFactory<MOD> cofac = null;
377 int pi = 0;
378 while ( pi < pn && primeIter.hasNext() ) {
379 p = primeIter.next();
380 //p = new java.math.BigInteger("19");
381 logger.info("prime = " + p);
382 // initialize coefficient factory and map normalization factor and polynomials
383 ModularRingFactory<MOD> cf = null;
384 if (ModLongRing.MAX_LONG.compareTo(p) > 0) {
385 cf = (ModularRingFactory) new ModLongRing(p, true);
386 } else {
387 cf = (ModularRingFactory) new ModIntegerRing(p, true);
388 }
389 MOD nf = cf.fromInteger(cc.getVal());
390 if (nf.isZERO()) {
391 continue;
392 }
393 nf = cf.fromInteger(q.leadingBaseCoefficient().leadingBaseCoefficient().getVal());
394 if (nf.isZERO()) {
395 continue;
396 }
397 nf = cf.fromInteger(r.leadingBaseCoefficient().leadingBaseCoefficient().getVal());
398 if (nf.isZERO()) {
399 continue;
400 }
401 cofac = cf;
402 break;
403 }
404 if ( cofac == null ) { // no lucky prime found
405 System.out.println("giving up on Hensel gcd reverting to Subres gcd");
406 GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(P, S);
407 return T.abs().multiply(c); //.abs();
408 }
409 //System.out.println("cofac = " + cofac);
410
411 List<MOD> V = new ArrayList<MOD>(P.ring.nvar);
412 GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(cofac, dfac);
413 //System.out.println("mfac = " + mfac.toScript());
414 GenPolynomial<MOD> qm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, qd);
415 GenPolynomial<MOD> rm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, rd);
416 //System.out.println("qm = " + qm);
417 //System.out.println("rm = " + rm);
418
419 // search evaluation point and evaluate
420 GenPolynomialRing<MOD> ckfac = mfac;
421 GenPolynomial<MOD> qe = qm;
422 GenPolynomial<MOD> re = rm;
423 GenPolynomial<MOD> qep;
424 GenPolynomial<MOD> rep;
425 for ( int j = dfac.nvar; j > 1; j-- ) {
426 // evaluation to univariate case
427 long degq = qe.degree(ckfac.nvar-2);
428 long degr = re.degree(ckfac.nvar-2);
429 ckfac = ckfac.contract(1);
430 long vi = 1L; //(long)(dfac.nvar-j); // 1L; 0 not so good for small p
431 if ( p.longValue() > 1000L ) {
432 //vi = (long)j+1L;
433 vi = 0L;
434 }
435 // search small evaluation point
436 while( true ) {
437 MOD vp = cofac.fromInteger(vi++);
438 //System.out.println("vp = " + vp);
439 if ( vp.isZERO() && vi != 1L ) { // all elements of Z_p exhausted
440 qe = null;
441 re = null;
442 break;
443 //GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(P, S);
444 //return T.abs().multiply(c); //.abs();
445 }
446 qep = PolyUtil.<MOD> evaluateMain(ckfac,qe,vp);
447 rep = PolyUtil.<MOD> evaluateMain(ckfac,re,vp);
448 //System.out.println("qep = " + qep);
449 //System.out.println("rep = " + rep);
450
451 // check lucky evaluation point
452 MOD ql = qep.leadingBaseCoefficient();
453 MOD rl = rep.leadingBaseCoefficient();
454 if (ql.isZERO()||rl.isZERO()) { // nearly non sense
455 continue;
456 }
457 if (degq != qep.degree(ckfac.nvar-1)) {
458 //System.out.println("degv(qe) = " + qe.degreeVector());
459 //System.out.println("deg(qe) = " + degq + ", deg(qe) = " + qep.degree(ckfac.nvar-1));
460 continue;
461 }
462 if (degr != rep.degree(ckfac.nvar-1)) {
463 //System.out.println("degv(re) = " + re.degreeVector());
464 //System.out.println("deg(re) = " + degr + ", deg(re) = " + rep.degree(ckfac.nvar-1));
465 continue;
466 }
467 V.add(vp);
468 qe = qep;
469 re = rep;
470 break;
471 }
472 if ( qe == null && re == null ) {
473 break;
474 }
475 }
476 if ( qe == null && re == null ) {
477 continue;
478 }
479 logger.info("evaluation points = " + V);
480 //System.out.println("qe = " + qe);
481 //System.out.println("re = " + re);
482
483 // recursion base:
484 GreatestCommonDivisorAbstract<MOD> mufd = GCDFactory.getImplementation(cofac);
485 GenPolynomial<MOD> ce = mufd.baseGcd(qe,re);
486 if ( ce.isConstant() ) {
487 return P.ring.getONE().multiply(c);
488 }
489 //System.out.println("ce = " + ce);
490 logger.info("base gcd = " + ce);
491 //System.out.println("c = " + c);
492
493 // double check
494 // need qe,re,qd,rd,a,b
495 if ( i == 0 ) {
496 qe0 = qe;
497 re0 = re;
498 ce0 = ce;
499 continue;
500 } else {
501 long d0 = ce0.degree(0);
502 long d1 = ce.degree(0);
503 //System.out.println("d0, d1 = " + d0 + ", " + d1);
504 if ( d1 < d0 ) {
505 qe0 = qe;
506 re0 = re;
507 ce0 = ce;
508 continue;
509 } else if ( d1 > d0 ) {
510 continue;
511 }
512 // d0 == d1 is ok
513 long dx = r.degree(0);
514 //System.out.println("d0, dx = " + d0 + ", " + dx);
515 if ( d0 == dx ) { // gcd == r ??
516 if (PolyUtil.<BigInteger> recursivePseudoRemainder(q,r).isZERO() ) {
517 r = r.abs().multiply(c); //.abs();
518 logger.info("exit with r | q : " + r);
519 return r;
520 }
521 continue;
522 }
523 }
524
525 // norm
526 BigInteger mn = null;
527 //mn = mn.multiply(cc).multiply(mn.fromInteger(2));
528
529 // prepare lifting
530 GenPolynomial<MOD> re1 = PolyUtil.<MOD> basePseudoDivide(re,ce);
531 GenPolynomial<MOD> qe1 = PolyUtil.<MOD> basePseudoDivide(qe,ce);
532 GenPolynomial<BigInteger> ui;
533 GenPolynomial<MOD> he;
534 GenPolynomial<BigInteger> g;
535 if ( mufd.baseGcd(re1,ce).isONE() ) {
536 ui = rd;
537 he = re1;
538 BigInteger an = rd.maxNorm();
539 mn = an.multiply(cc).multiply(new BigInteger(2L));
540 g = a;
541 } else if ( mufd.baseGcd(qe1,ce).isONE() ) {
542 ui = qd;
543 he = qe1;
544 BigInteger bn = qd.maxNorm();
545 mn = bn.multiply(cc).multiply(new BigInteger(2L));
546 g = b;
547 } else {
548 break;
549 //System.out.println("giving up on Hensel gcd reverting to Subres gcd");
550 }
551 //System.out.println("ui = " + ui);
552
553 long k = Power.logarithm(new BigInteger(p),mn);
554 //System.out.println("mn = " + mn);
555 //System.out.println("k = " + k);
556
557 List<GenPolynomial<MOD>> F = new ArrayList<GenPolynomial<MOD>>(2);
558 F.add(ce);
559 F.add(he);
560 List<GenPolynomial<BigInteger>> G = new ArrayList<GenPolynomial<BigInteger>>(2);
561 G.add(g.ring.getONE()); // TODO
562 G.add(g.ring.getONE());
563 List<GenPolynomial<MOD>> lift;
564 try {
565 lift = HenselMultUtil.<MOD> liftHenselFull(ui,F,V,k,G);
566 logger.info("lift = " + lift);
567 } catch ( NoLiftingException nle ) {
568 //System.out.println("exception : " + nle);
569 continue;
570 } catch ( ArithmeticException ae ) {
571 //System.out.println("exception : " + ae);
572 continue;
573 }
574
575 // convert Ci from Z_{p^k}[x,y1,...,yr] to Z[x,y1,...,yr] to Z[x][y1,...,yr] to Z[y1,...,yr][x]
576 GenPolynomial<BigInteger> ci = PolyUtil.integerFromModularCoefficients( dfac, lift.get(0) );
577 GenPolynomial<GenPolynomial<BigInteger>> Cr = PolyUtil.<BigInteger> recursive(rfac,ci);
578 GenPolynomial<GenPolynomial<BigInteger>> Cs = PolyUtil.<BigInteger> switchVariables(Cr);
579 if ( ! Cs.ring.equals(P.ring) ) {
580 System.out.println("Cs.ring = " + Cs.ring + ", P.ring = " + P.ring);
581 }
582 q = recursivePrimitivePart(Cs);
583 q = q.abs().multiply(c); //.abs();
584 if ( PolyUtil.<BigInteger> recursivePseudoRemainder(P,q).isZERO()
585 && PolyUtil.<BigInteger> recursivePseudoRemainder(S,q).isZERO() ) {
586 return q;
587 } else {
588 //System.out.println("bad q = " + q);
589 //System.out.println("bad P/q = " + PolyUtil.<BigInteger> recursivePseudoRemainder(P,q));
590 //System.out.println("bad S/q = " + PolyUtil.<BigInteger> recursivePseudoRemainder(S,q));
591 }
592 } // end for meta loop
593 // Hensel gcd failed
594 GenPolynomial<GenPolynomial<BigInteger>> T = iufd.recursiveUnivariateGcd(P, S);
595 T = T.abs().multiply(c); //.abs();
596 logger.info("giving up on Hensel gcd reverting to Subres gcd " + T);
597 //System.out.println("giving up on Hensel gcd reverted to Subres gcd: " + T);
598 return T;
599 }
600
601 }