001 /*
002 * $Id: GreatestCommonDivisorModular.java 3355 2010-10-23 16:01:52Z kredel $
003 */
004
005 package edu.jas.ufd;
006
007
008 import org.apache.log4j.Logger;
009
010 import edu.jas.structure.GcdRingElem;
011 import edu.jas.arith.BigInteger;
012 import edu.jas.arith.Modular;
013 import edu.jas.arith.ModLongRing;
014 import edu.jas.arith.ModIntegerRing;
015 import edu.jas.arith.ModularRingFactory;
016 import edu.jas.arith.PrimeList;
017 import edu.jas.poly.ExpVector;
018 import edu.jas.poly.GenPolynomial;
019 import edu.jas.poly.GenPolynomialRing;
020 import edu.jas.poly.PolyUtil;
021
022
023 /**
024 * Greatest common divisor algorithms with modular computation and chinese
025 * remainder algorithm.
026 * @author Heinz Kredel
027 */
028
029 public class GreatestCommonDivisorModular<MOD extends GcdRingElem<MOD> & Modular>
030 extends GreatestCommonDivisorAbstract<BigInteger> {
031
032
033 private static final Logger logger = Logger.getLogger(GreatestCommonDivisorModular.class);
034
035
036 private final boolean debug = logger.isDebugEnabled(); //logger.isInfoEnabled();
037
038
039 /*
040 * Modular gcd algorithm to use.
041 */
042 protected final GreatestCommonDivisorAbstract<MOD> mufd;
043
044
045 /*
046 * Integer gcd algorithm for fall back.
047 */
048 protected final GreatestCommonDivisorAbstract<BigInteger> iufd = new GreatestCommonDivisorSubres<BigInteger>();
049
050
051 /**
052 * Constructor to set recursive algorithm. Use modular evaluation GCD
053 * algorithm.
054 */
055 public GreatestCommonDivisorModular() {
056 this(false);
057 }
058
059
060 /**
061 * Constructor to set recursive algorithm.
062 * @param simple , true if the simple PRS should be used.
063 */
064 public GreatestCommonDivisorModular(boolean simple) {
065 if (simple) {
066 mufd = new GreatestCommonDivisorSimple<MOD>();
067 } else {
068 mufd = new GreatestCommonDivisorModEval<MOD>();
069 }
070 }
071
072
073 /**
074 * Univariate GenPolynomial greatest comon divisor. Delegate to subresultant
075 * baseGcd, should not be needed.
076 * @param P univariate GenPolynomial.
077 * @param S univariate GenPolynomial.
078 * @return gcd(P,S).
079 */
080 @Override
081 public GenPolynomial<BigInteger> baseGcd(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) {
082 return iufd.baseGcd(P, S);
083 }
084
085
086 /**
087 * Univariate GenPolynomial recursive greatest comon divisor. Delegate to
088 * subresultant recursiveGcd, should not be needed.
089 * @param P univariate recursive GenPolynomial.
090 * @param S univariate recursive GenPolynomial.
091 * @return gcd(P,S).
092 */
093 @Override
094 public GenPolynomial<GenPolynomial<BigInteger>> recursiveUnivariateGcd(
095 GenPolynomial<GenPolynomial<BigInteger>> P, GenPolynomial<GenPolynomial<BigInteger>> S) {
096 return iufd.recursiveUnivariateGcd(P, S);
097 }
098
099
100 /**
101 * GenPolynomial greatest comon divisor, modular algorithm.
102 * @param P GenPolynomial.
103 * @param S GenPolynomial.
104 * @return gcd(P,S).
105 */
106 @Override
107 public GenPolynomial<BigInteger> gcd(GenPolynomial<BigInteger> P, GenPolynomial<BigInteger> S) {
108 if (S == null || S.isZERO()) {
109 return P;
110 }
111 if (P == null || P.isZERO()) {
112 return S;
113 }
114 GenPolynomialRing<BigInteger> fac = P.ring;
115 // special case for univariate polynomials
116 if (fac.nvar <= 1) {
117 GenPolynomial<BigInteger> T = baseGcd(P, S);
118 return T;
119 }
120 long e = P.degree(0);
121 long f = S.degree(0);
122 GenPolynomial<BigInteger> q;
123 GenPolynomial<BigInteger> r;
124 if (f > e) {
125 r = P;
126 q = S;
127 long g = f;
128 f = e;
129 e = g;
130 } else {
131 q = P;
132 r = S;
133 }
134 r = r.abs();
135 q = q.abs();
136 // compute contents and primitive parts
137 BigInteger a = baseContent(r);
138 BigInteger b = baseContent(q);
139 // gcd of coefficient contents
140 BigInteger c = gcd(a, b); // indirection
141 r = divide(r, a); // indirection
142 q = divide(q, b); // indirection
143 if (r.isONE()) {
144 return r.multiply(c);
145 }
146 if (q.isONE()) {
147 return q.multiply(c);
148 }
149 // compute normalization factor
150 BigInteger ac = r.leadingBaseCoefficient();
151 BigInteger bc = q.leadingBaseCoefficient();
152 BigInteger cc = gcd(ac, bc); // indirection
153 // compute norms
154 BigInteger an = r.maxNorm();
155 BigInteger bn = q.maxNorm();
156 BigInteger n = (an.compareTo(bn) < 0 ? bn : an);
157 n = n.multiply(cc).multiply(n.fromInteger(2));
158 // compute degree vectors
159 ExpVector rdegv = r.degreeVector();
160 ExpVector qdegv = q.degreeVector();
161 //compute factor coefficient bounds
162 BigInteger af = an.multiply(PolyUtil.factorBound(rdegv));
163 BigInteger bf = bn.multiply(PolyUtil.factorBound(qdegv));
164 BigInteger cf = (af.compareTo(bf) < 0 ? bf : af);
165 cf = cf.multiply(cc.multiply(cc.fromInteger(8)));
166 //initialize prime list and degree vector
167 PrimeList primes = new PrimeList();
168 int pn = 10; //primes.size();
169 ExpVector wdegv = rdegv.subst(0, rdegv.getVal(0) + 1);
170 // +1 seems to be a hack for the unlucky prime test
171 ModularRingFactory<MOD> cofac;
172 ModularRingFactory<MOD> cofacM = null;
173 GenPolynomial<MOD> qm;
174 GenPolynomial<MOD> rm;
175 GenPolynomialRing<MOD> mfac;
176 GenPolynomialRing<MOD> rfac = null;
177 int i = 0;
178 BigInteger M = null;
179 BigInteger cfe = null;
180 GenPolynomial<MOD> cp = null;
181 GenPolynomial<MOD> cm = null;
182 GenPolynomial<BigInteger> cpi = null;
183 if (debug) {
184 logger.debug("c = " + c);
185 logger.debug("cc = " + cc);
186 logger.debug("n = " + n);
187 logger.debug("cf = " + cf);
188 logger.info("wdegv = " + wdegv);
189 }
190 for (java.math.BigInteger p : primes) {
191 //System.out.println("next run ++++++++++++++++++++++++++++++++++");
192 if ( p.longValue() == 2L ) { // skip 2
193 continue;
194 }
195 if (++i >= pn) {
196 logger.error("prime list exhausted, pn = " + pn);
197 return iufd.gcd(P, S);
198 //throw new ArithmeticException("prime list exhausted");
199 }
200 // initialize coefficient factory and map normalization factor
201 if ( ModLongRing.MAX_LONG.compareTo( p ) > 0 ) {
202 cofac = (ModularRingFactory) new ModLongRing(p, true);
203 } else {
204 cofac = (ModularRingFactory) new ModIntegerRing(p, true);
205 }
206 MOD nf = cofac.fromInteger(cc.getVal());
207 if (nf.isZERO()) {
208 continue;
209 }
210 // initialize polynomial factory and map polynomials
211 mfac = new GenPolynomialRing<MOD>(cofac, fac.nvar, fac.tord, fac.getVars());
212 qm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, q);
213 if (!qm.degreeVector().equals(qdegv)) {
214 continue;
215 }
216 rm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, r);
217 if (!rm.degreeVector().equals(rdegv)) {
218 continue;
219 }
220 if (debug) {
221 logger.info("cofac = " + cofac.getIntegerModul());
222 }
223 // compute modular gcd
224 cm = mufd.gcd(rm, qm);
225 // test for constant g.c.d
226 if (cm.isConstant()) {
227 logger.debug("cm, constant = " + cm);
228 return fac.getONE().multiply(c);
229 //return cm.abs().multiply( c );
230 }
231 // test for unlucky prime
232 ExpVector mdegv = cm.degreeVector();
233 if (wdegv.equals(mdegv)) { // TL = 0
234 // prime ok, next round
235 if (M != null) {
236 if (M.compareTo(cfe) > 0) {
237 System.out.println("M > cfe: " + M + " > " + cfe);
238 // continue; // why should this be required?
239 }
240 }
241 } else { // TL = 3
242 boolean ok = false;
243 if (wdegv.multipleOf(mdegv)) { // TL = 2 // EVMT(wdegv,mdegv)
244 M = null; // init chinese remainder
245 ok = true; // prime ok
246 }
247 if (mdegv.multipleOf(wdegv)) { // TL = 1 // EVMT(mdegv,wdegv)
248 continue; // skip this prime
249 }
250 if (!ok) {
251 M = null; // discard chinese remainder and previous work
252 continue; // prime not ok
253 }
254 }
255 //--wdegv = mdegv;
256 // prepare chinese remainder algorithm
257 cm = cm.multiply(nf);
258 if (M == null) {
259 // initialize chinese remainder algorithm
260 M = new BigInteger(p);
261 cofacM = cofac;
262 rfac = mfac;
263 cp = cm;
264 wdegv = wdegv.gcd(mdegv); //EVGCD(wdegv,mdegv);
265 cfe = cf;
266 for (int k = 0; k < wdegv.length(); k++) {
267 cfe = cfe.multiply(new BigInteger(wdegv.getVal(k) + 1));
268 }
269 } else {
270 // apply chinese remainder algorithm
271 MOD mi = cofac.fromInteger(M.getVal());
272 mi = mi.inverse(); // mod p
273 M = M.multiply(new BigInteger(p));
274 if ( ModLongRing.MAX_LONG.compareTo( M.getVal() ) > 0 ) {
275 cofacM = (ModularRingFactory) new ModLongRing(M.getVal());
276 } else {
277 cofacM = (ModularRingFactory) new ModIntegerRing(M.getVal());
278 }
279 rfac = new GenPolynomialRing<MOD>(cofacM, fac.nvar, fac.tord, fac.getVars());
280 cp = PolyUtil.<MOD> chineseRemainder(rfac, cp, mi, cm);
281 }
282 // test for completion
283 if (n.compareTo(M) <= 0) {
284 break;
285 }
286 // must use integer.sumNorm
287 cpi = PolyUtil.<MOD>integerFromModularCoefficients(fac, cp);
288 BigInteger cmn = cpi.sumNorm();
289 cmn = cmn.multiply(cmn.fromInteger(4));
290 //if ( cmn.compareTo( M ) <= 0 ) {
291 // does not work: only if cofactors are also considered?
292 // break;
293 //}
294 if (i % 2 != 0 && !cp.isZERO()) {
295 // check if done on every second prime
296 GenPolynomial<BigInteger> x;
297 x = PolyUtil.<MOD>integerFromModularCoefficients(fac, cp);
298 x = basePrimitivePart(x);
299 if (!PolyUtil.<BigInteger> basePseudoRemainder(q, x).isZERO()) {
300 continue;
301 }
302 if (!PolyUtil.<BigInteger> basePseudoRemainder(r, x).isZERO()) {
303 continue;
304 }
305 logger.info("done on exact division, #primes = " + i);
306 break;
307 }
308 }
309 if (debug) {
310 logger.info("done on M = " + M + ", #primes = " + i);
311 }
312 // remove normalization
313 q = PolyUtil.<MOD>integerFromModularCoefficients(fac, cp);
314 q = basePrimitivePart(q);
315 return q.abs().multiply(c);
316 }
317
318 }