001 /*
002 * $Id: GreatestCommonDivisorModEval.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.arith.Modular;
011 import edu.jas.arith.ModularRingFactory;
012 import edu.jas.poly.ExpVector;
013 import edu.jas.poly.GenPolynomial;
014 import edu.jas.poly.GenPolynomialRing;
015 import edu.jas.poly.PolyUtil;
016 import edu.jas.structure.GcdRingElem;
017
018
019 /**
020 * Greatest common divisor algorithms with modular evaluation algorithm for
021 * recursion.
022 * @author Heinz Kredel
023 */
024
025 public class GreatestCommonDivisorModEval <MOD extends GcdRingElem<MOD> & Modular>
026 extends GreatestCommonDivisorAbstract<MOD> {
027
028
029 private static final Logger logger = Logger.getLogger(GreatestCommonDivisorModEval.class);
030
031
032 private final boolean debug = logger.isDebugEnabled();
033
034
035 /**
036 * Modular gcd algorithm to use.
037 */
038 protected final GreatestCommonDivisorAbstract<MOD> mufd
039 // == new GreatestCommonDivisorModular();
040 = new GreatestCommonDivisorSimple<MOD>();
041 // = new GreatestCommonDivisorPrimitive<MOD>();
042 // = new GreatestCommonDivisorSubres<MOD>();
043
044
045 /**
046 * Univariate GenPolynomial greatest comon divisor. Delegate to subresultant
047 * baseGcd, should not be needed.
048 * @param P univariate GenPolynomial.
049 * @param S univariate GenPolynomial.
050 * @return gcd(P,S).
051 */
052 @Override
053 public GenPolynomial<MOD> baseGcd(GenPolynomial<MOD> P, GenPolynomial<MOD> S) {
054 return mufd.baseGcd(P, S);
055 }
056
057
058 /**
059 * Univariate GenPolynomial recursive greatest comon divisor. Delegate to
060 * subresultant recursiveGcd, should not be needed.
061 * @param P univariate recursive GenPolynomial.
062 * @param S univariate recursive GenPolynomial.
063 * @return gcd(P,S).
064 */
065 @Override
066 public GenPolynomial<GenPolynomial<MOD>> recursiveUnivariateGcd(
067 GenPolynomial<GenPolynomial<MOD>> P, GenPolynomial<GenPolynomial<MOD>> S) {
068 return mufd.recursiveUnivariateGcd(P, S);
069 }
070
071
072 /**
073 * GenPolynomial greatest comon divisor, modular evaluation algorithm.
074 * Method name must be different because of parameter type erasure.
075 * @param P GenPolynomial.
076 * @param S GenPolynomial.
077 * @return gcd(P,S).
078 */
079 @Override
080 public GenPolynomial<MOD> gcd(GenPolynomial<MOD> P, GenPolynomial<MOD> S) {
081 if (S == null || S.isZERO()) {
082 return P;
083 }
084 if (P == null || P.isZERO()) {
085 return S;
086 }
087 GenPolynomialRing<MOD> fac = P.ring;
088 // special case for univariate polynomials
089 if (fac.nvar <= 1) {
090 GenPolynomial<MOD> T = baseGcd(P, S);
091 return T;
092 }
093 long e = P.degree(0);
094 long f = S.degree(0);
095 GenPolynomial<MOD> q;
096 GenPolynomial<MOD> 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 // setup factories
110 ModularRingFactory<MOD> cofac = (ModularRingFactory<MOD>) P.ring.coFac;
111 if (!cofac.isField()) {
112 logger.warn("cofac is not a field: " + cofac);
113 }
114 GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(cofac, fac.nvar - 1, fac.tord);
115 GenPolynomialRing<MOD> ufac = new GenPolynomialRing<MOD>(cofac, 1, fac.tord);
116 GenPolynomialRing<GenPolynomial<MOD>> rfac = new GenPolynomialRing<GenPolynomial<MOD>>(
117 ufac, fac.nvar - 1, fac.tord);
118 // convert polynomials
119 GenPolynomial<GenPolynomial<MOD>> qr;
120 GenPolynomial<GenPolynomial<MOD>> rr;
121 qr = PolyUtil.<MOD> recursive(rfac, q);
122 rr = PolyUtil.<MOD> recursive(rfac, r);
123
124 // compute univariate contents and primitive parts
125 GenPolynomial<MOD> a = recursiveContent(rr);
126 GenPolynomial<MOD> b = recursiveContent(qr);
127 // gcd of univariate coefficient contents
128 GenPolynomial<MOD> c = gcd(a, b);
129 rr = PolyUtil.<MOD> recursiveDivide(rr, a);
130 qr = PolyUtil.<MOD> recursiveDivide(qr, b);
131 if (rr.isONE()) {
132 rr = rr.multiply(c);
133 r = PolyUtil.<MOD> distribute(fac, rr);
134 return r;
135 }
136 if (qr.isONE()) {
137 qr = qr.multiply(c);
138 q = PolyUtil.<MOD> distribute(fac, qr);
139 return q;
140 }
141 // compute normalization factor
142 GenPolynomial<MOD> ac = rr.leadingBaseCoefficient();
143 GenPolynomial<MOD> bc = qr.leadingBaseCoefficient();
144 GenPolynomial<MOD> cc = gcd(ac, bc);
145 // compute degrees and degree vectors
146 ExpVector rdegv = rr.degreeVector();
147 ExpVector qdegv = qr.degreeVector();
148 long rd0 = PolyUtil.<MOD> coeffMaxDegree(rr);
149 long qd0 = PolyUtil.<MOD> coeffMaxDegree(qr);
150 long cd0 = cc.degree(0);
151 long G = (rd0 >= qd0 ? rd0 : qd0) + cd0;
152
153 // initialize element and degree vector
154 ExpVector wdegv = rdegv.subst(0, rdegv.getVal(0) + 1);
155 // +1 seems to be a hack for the unlucky prime test
156 MOD inc = cofac.getONE();
157 long i = 0;
158 long en = cofac.getIntegerModul().longValue() - 1; // just a stopper
159 MOD end = cofac.fromInteger(en);
160 MOD mi;
161 GenPolynomial<MOD> M = null;
162 GenPolynomial<MOD> mn;
163 GenPolynomial<MOD> qm;
164 GenPolynomial<MOD> rm;
165 GenPolynomial<MOD> cm;
166 GenPolynomial<GenPolynomial<MOD>> cp = null;
167 if (debug) {
168 logger.debug("c = " + c);
169 logger.debug("cc = " + cc);
170 logger.debug("G = " + G);
171 logger.info("wdegv = " + wdegv);
172 }
173 for (MOD d = cofac.getZERO(); d.compareTo(end) <= 0; d = d.sum(inc)) {
174 if (++i >= en) {
175 logger.error("elements of Z_p exhausted, en = " + en);
176 return mufd.gcd(P, S);
177 //throw new ArithmeticException("prime list exhausted");
178 }
179 // map normalization factor
180 MOD nf = PolyUtil.<MOD> evaluateMain(cofac, cc, d);
181 if (nf.isZERO()) {
182 continue;
183 }
184 // map polynomials
185 qm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, qr, d);
186 if (!qm.degreeVector().equals(qdegv)) {
187 continue;
188 }
189 rm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, rr, d);
190 if (!rm.degreeVector().equals(rdegv)) {
191 continue;
192 }
193 if (debug) {
194 logger.debug("eval d = " + d);
195 }
196 // compute modular gcd in recursion
197 //System.out.println("recursion +++++++++++++++++++++++++++++++++++");
198 cm = gcd(rm, qm);
199 //System.out.println("recursion -----------------------------------");
200 //System.out.println("cm = " + cm);
201 // test for constant g.c.d
202 if (cm.isConstant()) {
203 logger.debug("cm.isConstant = " + cm + ", c = " + c);
204 if (c.ring.nvar < cm.ring.nvar) {
205 c = c.extend(mfac, 0, 0);
206 }
207 cm = cm.abs().multiply(c);
208 q = cm.extend(fac, 0, 0);
209 logger.debug("q = " + q + ", c = " + c);
210 return q;
211 }
212 // test for unlucky prime
213 ExpVector mdegv = cm.degreeVector();
214 if (wdegv.equals(mdegv)) { // TL = 0
215 // prime ok, next round
216 if (M != null) {
217 if (M.degree(0) > G) {
218 logger.info("deg(M) > G: " + M.degree(0) + " > " + G);
219 // continue; // why should this be required?
220 }
221 }
222 } else { // TL = 3
223 boolean ok = false;
224 if (wdegv.multipleOf(mdegv)) { // TL = 2 // EVMT(wdegv,mdegv)
225 M = null; // init chinese remainder
226 ok = true; // prime ok
227 }
228 if (mdegv.multipleOf(wdegv)) { // TL = 1 // EVMT(mdegv,wdegv)
229 continue; // skip this prime
230 }
231 if (!ok) {
232 M = null; // discard chinese remainder and previous work
233 continue; // prime not ok
234 }
235 }
236 // prepare interpolation algorithm
237 cm = cm.multiply(nf);
238 if (M == null) {
239 // initialize interpolation
240 M = ufac.getONE();
241 cp = rfac.getZERO();
242 wdegv = wdegv.gcd(mdegv); //EVGCD(wdegv,mdegv);
243 }
244 // interpolate
245 mi = PolyUtil.<MOD> evaluateMain(cofac, M, d);
246 mi = mi.inverse(); // mod p
247 cp = PolyUtil.interpolate(rfac, cp, M, mi, cm, d);
248 mn = ufac.getONE().multiply(d);
249 mn = ufac.univariate(0).subtract(mn);
250 M = M.multiply(mn);
251 // test for completion
252 if (M.degree(0) > G) {
253 break;
254 }
255 //long cmn = PolyUtil.<MOD>coeffMaxDegree(cp);
256 //if ( M.degree(0) > cmn ) {
257 // does not work: only if cofactors are also considered?
258 // break;
259 //}
260 }
261 // remove normalization
262 cp = recursivePrimitivePart(cp).abs();
263 cp = cp.multiply(c);
264 q = PolyUtil.<MOD> distribute(fac, cp);
265 return q;
266 }
267
268 }