001/* 002 * $Id: GreatestCommonDivisorModEval.java 4025 2012-07-23 16:41:43Z kredel $ 003 */ 004 005package edu.jas.ufd; 006 007 008import org.apache.log4j.Logger; 009 010import edu.jas.arith.Modular; 011import edu.jas.arith.ModularRingFactory; 012import edu.jas.poly.ExpVector; 013import edu.jas.poly.GenPolynomial; 014import edu.jas.poly.GenPolynomialRing; 015import edu.jas.poly.PolyUtil; 016import 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 025public 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 GreatestCommonDivisorSimple<MOD>(); 040 // = new GreatestCommonDivisorPrimitive<MOD>(); 041 // not okay: = new GreatestCommonDivisorSubres<MOD>(); 042 043 044 /** 045 * Univariate GenPolynomial greatest common divisor. 046 * @param P univariate GenPolynomial. 047 * @param S univariate GenPolynomial. 048 * @return gcd(P,S). 049 */ 050 @Override 051 public GenPolynomial<MOD> baseGcd(GenPolynomial<MOD> P, GenPolynomial<MOD> S) { 052 // required as recursion base 053 return mufd.baseGcd(P, S); 054 } 055 056 057 /** 058 * Recursive univariate GenPolynomial greatest common divisor. 059 * @param P univariate recursive GenPolynomial. 060 * @param S univariate recursive GenPolynomial. 061 * @return gcd(P,S). 062 */ 063 @Override 064 public GenPolynomial<GenPolynomial<MOD>> recursiveUnivariateGcd( 065 GenPolynomial<GenPolynomial<MOD>> P, GenPolynomial<GenPolynomial<MOD>> S) { 066 return mufd.recursiveUnivariateGcd(P, S); 067 } 068 069 070 /** 071 * GenPolynomial greatest common divisor, modular evaluation algorithm. 072 * @param P GenPolynomial. 073 * @param S GenPolynomial. 074 * @return gcd(P,S). 075 */ 076 @Override 077 public GenPolynomial<MOD> gcd(GenPolynomial<MOD> P, GenPolynomial<MOD> S) { 078 if (S == null || S.isZERO()) { 079 return P; 080 } 081 if (P == null || P.isZERO()) { 082 return S; 083 } 084 GenPolynomialRing<MOD> fac = P.ring; 085 // recusion base case for univariate polynomials 086 if (fac.nvar <= 1) { 087 GenPolynomial<MOD> T = baseGcd(P, S); 088 return T; 089 } 090 long e = P.degree(fac.nvar-1); 091 long f = S.degree(fac.nvar-1); 092 if ( e == 0 && f == 0 ) { 093 GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(1); 094 GenPolynomial<MOD> Pc = PolyUtil.<MOD> recursive(rfac, P).leadingBaseCoefficient(); 095 GenPolynomial<MOD> Sc = PolyUtil.<MOD> recursive(rfac, S).leadingBaseCoefficient(); 096 GenPolynomial<MOD> r = gcd(Pc,Sc); 097 return r.extend(fac,0,0L); 098 } 099 GenPolynomial<MOD> q; 100 GenPolynomial<MOD> r; 101 if (f > e) { 102 r = P; 103 q = S; 104 long g = f; 105 f = e; 106 e = g; 107 } else { 108 q = P; 109 r = S; 110 } 111 if (debug) { 112 logger.debug("degrees: e = " + e + ", f = " + f); 113 } 114 r = r.abs(); 115 q = q.abs(); 116 // setup factories 117 ModularRingFactory<MOD> cofac = (ModularRingFactory<MOD>) P.ring.coFac; 118 if (!cofac.isField()) { 119 logger.warn("cofac is not a field: " + cofac); 120 } 121 GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(fac.nvar - 1); 122 GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(cofac, rfac); 123 GenPolynomialRing<MOD> ufac = (GenPolynomialRing<MOD>) rfac.coFac; 124 //GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(cofac, fac.nvar - 1, fac.tord); 125 //GenPolynomialRing<MOD> ufac = new GenPolynomialRing<MOD>(cofac, 1, fac.tord); 126 //GenPolynomialRing<GenPolynomial<MOD>> rfac = new GenPolynomialRing<GenPolynomial<MOD>>(ufac, fac.nvar - 1, fac.tord); 127 // convert polynomials 128 GenPolynomial<GenPolynomial<MOD>> qr; 129 GenPolynomial<GenPolynomial<MOD>> rr; 130 qr = PolyUtil.<MOD> recursive(rfac, q); 131 rr = PolyUtil.<MOD> recursive(rfac, r); 132 133 // compute univariate contents and primitive parts 134 GenPolynomial<MOD> a = recursiveContent(rr); 135 GenPolynomial<MOD> b = recursiveContent(qr); 136 // gcd of univariate coefficient contents 137 GenPolynomial<MOD> c = gcd(a, b); 138 rr = PolyUtil.<MOD> recursiveDivide(rr, a); 139 qr = PolyUtil.<MOD> recursiveDivide(qr, b); 140 if (rr.isONE()) { 141 rr = rr.multiply(c); 142 r = PolyUtil.<MOD> distribute(fac, rr); 143 return r; 144 } 145 if (qr.isONE()) { 146 qr = qr.multiply(c); 147 q = PolyUtil.<MOD> distribute(fac, qr); 148 return q; 149 } 150 // compute normalization factor 151 GenPolynomial<MOD> ac = rr.leadingBaseCoefficient(); 152 GenPolynomial<MOD> bc = qr.leadingBaseCoefficient(); 153 GenPolynomial<MOD> cc = gcd(ac, bc); 154 // compute degrees and degree vectors 155 ExpVector rdegv = rr.degreeVector(); 156 ExpVector qdegv = qr.degreeVector(); 157 long rd0 = PolyUtil.<MOD> coeffMaxDegree(rr); 158 long qd0 = PolyUtil.<MOD> coeffMaxDegree(qr); 159 long cd0 = cc.degree(0); 160 long G = (rd0 >= qd0 ? rd0 : qd0) + cd0; 161 162 // initialize element and degree vector 163 ExpVector wdegv = rdegv.subst(0, rdegv.getVal(0) + 1); 164 // +1 seems to be a hack for the unlucky prime test 165 MOD inc = cofac.getONE(); 166 long i = 0; 167 long en = cofac.getIntegerModul().longValue() - 1; // just a stopper 168 MOD end = cofac.fromInteger(en); 169 MOD mi; 170 GenPolynomial<MOD> M = null; 171 GenPolynomial<MOD> mn; 172 GenPolynomial<MOD> qm; 173 GenPolynomial<MOD> rm; 174 GenPolynomial<MOD> cm; 175 GenPolynomial<GenPolynomial<MOD>> cp = null; 176 if (debug) { 177 logger.debug("c = " + c); 178 logger.debug("cc = " + cc); 179 logger.debug("G = " + G); 180 logger.info("wdegv = " + wdegv); 181 } 182 for (MOD d = cofac.getZERO(); d.compareTo(end) <= 0; d = d.sum(inc)) { 183 if (++i >= en) { 184 logger.warn("elements of Z_p exhausted, en = " + en); 185 return mufd.gcd(P, S); 186 //throw new ArithmeticException("prime list exhausted"); 187 } 188 // map normalization factor 189 MOD nf = PolyUtil.<MOD> evaluateMain(cofac, cc, d); 190 if (nf.isZERO()) { 191 continue; 192 } 193 // map polynomials 194 qm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, qr, d); 195 if (qm.isZERO() || !qm.degreeVector().equals(qdegv)) { 196 continue; 197 } 198 rm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, rr, d); 199 if (rm.isZERO() || !rm.degreeVector().equals(rdegv)) { 200 continue; 201 } 202 if (debug) { 203 logger.debug("eval d = " + d); 204 } 205 // compute modular gcd in recursion 206 cm = gcd(rm, qm); 207 //System.out.println("cm = " + cm); 208 // test for constant g.c.d 209 if (cm.isConstant()) { 210 logger.debug("cm.isConstant = " + cm + ", c = " + c); 211 if (c.ring.nvar < cm.ring.nvar) { 212 c = c.extend(mfac, 0, 0); 213 } 214 cm = cm.abs().multiply(c); 215 q = cm.extend(fac, 0, 0); 216 logger.debug("q = " + q + ", c = " + c); 217 return q; 218 } 219 // test for unlucky prime 220 ExpVector mdegv = cm.degreeVector(); 221 if (wdegv.equals(mdegv)) { // TL = 0 222 // prime ok, next round 223 if (M != null) { 224 if (M.degree(0) > G) { 225 logger.info("deg(M) > G: " + M.degree(0) + " > " + G); 226 // continue; // why should this be required? 227 } 228 } 229 } else { // TL = 3 230 boolean ok = false; 231 if (wdegv.multipleOf(mdegv)) { // TL = 2 232 M = null; // init chinese remainder 233 ok = true; // prime ok 234 } 235 if (mdegv.multipleOf(wdegv)) { // TL = 1 236 continue; // skip this prime 237 } 238 if (!ok) { 239 M = null; // discard chinese remainder and previous work 240 continue; // prime not ok 241 } 242 } 243 // prepare interpolation algorithm 244 cm = cm.multiply(nf); 245 if (M == null) { 246 // initialize interpolation 247 M = ufac.getONE(); 248 cp = rfac.getZERO(); 249 wdegv = wdegv.gcd(mdegv); //EVGCD(wdegv,mdegv); 250 } 251 // interpolate 252 mi = PolyUtil.<MOD> evaluateMain(cofac, M, d); 253 mi = mi.inverse(); // mod p 254 cp = PolyUtil.interpolate(rfac, cp, M, mi, cm, d); 255 mn = ufac.getONE().multiply(d); 256 mn = ufac.univariate(0).subtract(mn); 257 M = M.multiply(mn); 258 // test for completion 259 if (M.degree(0) > G) { 260 break; 261 } 262 //long cmn = PolyUtil.<MOD>coeffMaxDegree(cp); 263 //if ( M.degree(0) > cmn ) { 264 // does not work: only if cofactors are also considered? 265 // break; 266 //} 267 } 268 // remove normalization 269 cp = recursivePrimitivePart(cp).abs(); 270 cp = cp.multiply(c); 271 q = PolyUtil.<MOD> distribute(fac, cp); 272 return q; 273 } 274 275 276 /** 277 * Univariate GenPolynomial resultant. 278 * @param P univariate GenPolynomial. 279 * @param S univariate GenPolynomial. 280 * @return res(P,S). 281 */ 282 @Override 283 public GenPolynomial<MOD> baseResultant(GenPolynomial<MOD> P, GenPolynomial<MOD> S) { 284 // required as recursion base 285 return mufd.baseResultant(P, S); 286 } 287 288 289 /** 290 * Univariate GenPolynomial recursive resultant. 291 * @param P univariate recursive GenPolynomial. 292 * @param S univariate recursive GenPolynomial. 293 * @return res(P,S). 294 */ 295 @Override 296 public GenPolynomial<GenPolynomial<MOD>> recursiveUnivariateResultant(GenPolynomial<GenPolynomial<MOD>> P, 297 GenPolynomial<GenPolynomial<MOD>> S) { 298 // only in this class 299 return recursiveResultant(P,S); 300 } 301 302 303 /** 304 * GenPolynomial resultant, modular evaluation algorithm. 305 * @param P GenPolynomial. 306 * @param S GenPolynomial. 307 * @return res(P,S). 308 */ 309 @Override 310 public GenPolynomial<MOD> resultant(GenPolynomial<MOD> P, GenPolynomial<MOD> S) { 311 if (S == null || S.isZERO()) { 312 return S; 313 } 314 if (P == null || P.isZERO()) { 315 return P; 316 } 317 GenPolynomialRing<MOD> fac = P.ring; 318 // recusion base case for univariate polynomials 319 if (fac.nvar <= 1) { 320 GenPolynomial<MOD> T = mufd.baseResultant(P, S); 321 return T; 322 } 323 long e = P.degree(fac.nvar-1); 324 long f = S.degree(fac.nvar-1); 325 if ( e == 0 && f == 0 ) { 326 GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(1); 327 GenPolynomial<MOD> Pc = PolyUtil.<MOD> recursive(rfac, P).leadingBaseCoefficient(); 328 GenPolynomial<MOD> Sc = PolyUtil.<MOD> recursive(rfac, S).leadingBaseCoefficient(); 329 GenPolynomial<MOD> r = resultant(Pc,Sc); 330 return r.extend(fac,0,0L); 331 } 332 GenPolynomial<MOD> q; 333 GenPolynomial<MOD> r; 334 if (f > e) { 335 r = P; 336 q = S; 337 long g = f; 338 f = e; 339 e = g; 340 } else { 341 q = P; 342 r = S; 343 } 344 if (debug) { 345 logger.debug("degrees: e = " + e + ", f = " + f); 346 } 347 // setup factories 348 ModularRingFactory<MOD> cofac = (ModularRingFactory<MOD>) P.ring.coFac; 349 if (!cofac.isField()) { 350 logger.warn("cofac is not a field: " + cofac); 351 } 352 GenPolynomialRing<GenPolynomial<MOD>> rfac = fac.recursive(fac.nvar - 1); 353 GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(cofac, rfac); 354 GenPolynomialRing<MOD> ufac = (GenPolynomialRing<MOD>) rfac.coFac; 355 356 // convert polynomials 357 GenPolynomial<GenPolynomial<MOD>> qr = PolyUtil.<MOD> recursive(rfac, q); 358 GenPolynomial<GenPolynomial<MOD>> rr = PolyUtil.<MOD> recursive(rfac, r); 359 360 // compute degrees and degree vectors 361 ExpVector qdegv = qr.degreeVector(); 362 ExpVector rdegv = rr.degreeVector(); 363 364 long qd0 = PolyUtil.<MOD> coeffMaxDegree(qr); 365 long rd0 = PolyUtil.<MOD> coeffMaxDegree(rr); 366 qd0 = ( qd0 == 0 ? 1 : qd0 ); 367 rd0 = ( rd0 == 0 ? 1 : rd0 ); 368 long qd1 = qr.degree(); 369 long rd1 = rr.degree(); 370 qd1 = ( qd1 == 0 ? 1 : qd1 ); 371 rd1 = ( rd1 == 0 ? 1 : rd1 ); 372 long G = qd0 * rd1 + rd0 * qd1 + 1; 373 374 // initialize element 375 MOD inc = cofac.getONE(); 376 long i = 0; 377 long en = cofac.getIntegerModul().longValue() - 1; // just a stopper 378 MOD end = cofac.fromInteger(en); 379 MOD mi; 380 GenPolynomial<MOD> M = null; 381 GenPolynomial<MOD> mn; 382 GenPolynomial<MOD> qm; 383 GenPolynomial<MOD> rm; 384 GenPolynomial<MOD> cm; 385 GenPolynomial<GenPolynomial<MOD>> cp = null; 386 if (debug) { 387 //logger.info("qr = " + qr + ", q = " + q); 388 //logger.info("rr = " + rr + ", r = " + r); 389 //logger.info("qd0 = " + qd0); 390 //logger.info("rd0 = " + rd0); 391 logger.info("G = " + G); 392 //logger.info("rdegv = " + rdegv); // + ", rr.degree(0) = " + rr.degree(0)); 393 //logger.info("qdegv = " + qdegv); // + ", qr.degree(0) = " + qr.degree(0)); 394 } 395 for (MOD d = cofac.getZERO(); d.compareTo(end) <= 0; d = d.sum(inc)) { 396 if (++i >= en) { 397 logger.warn("elements of Z_p exhausted, en = " + en + ", p = " + cofac.getIntegerModul()); 398 return mufd.resultant(P, S); 399 //throw new ArithmeticException("prime list exhausted"); 400 } 401 // map polynomials 402 qm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, qr, d); 403 //logger.info("qr(" + d + ") = " + qm + ", qr = " + qr); 404 if (qm.isZERO() || !qm.degreeVector().equals(qdegv)) { 405 if (debug) { 406 logger.info("un-lucky evaluation point " + d + ", qm = " + qm.degreeVector() + " < " + qdegv); 407 } 408 continue; 409 } 410 rm = PolyUtil.<MOD> evaluateFirstRec(ufac, mfac, rr, d); 411 //logger.info("rr(" + d + ") = " + rm + ", rr = " + rr); 412 if (rm.isZERO() || !rm.degreeVector().equals(rdegv)) { 413 if (debug) { 414 logger.info("un-lucky evaluation point " + d + ", rm = " + rm.degreeVector() + " < " + rdegv); 415 } 416 continue; 417 } 418 // compute modular resultant in recursion 419 cm = resultant(rm, qm); 420 //System.out.println("cm = " + cm); 421 422 // prepare interpolation algorithm 423 if (M == null) { 424 // initialize interpolation 425 M = ufac.getONE(); 426 cp = rfac.getZERO(); 427 } 428 // interpolate 429 mi = PolyUtil.<MOD> evaluateMain(cofac, M, d); 430 mi = mi.inverse(); // mod p 431 cp = PolyUtil.interpolate(rfac, cp, M, mi, cm, d); 432 //logger.info("cp = " + cp); 433 mn = ufac.getONE().multiply(d); 434 mn = ufac.univariate(0).subtract(mn); 435 M = M.multiply(mn); 436 // test for completion 437 if (M.degree(0) > G) { 438 if (debug) { 439 logger.info("last lucky evaluation point " + d); 440 } 441 break; 442 } 443 //logger.info("M = " + M); 444 } 445 // distribute 446 q = PolyUtil.<MOD> distribute(fac, cp); 447 return q; 448 } 449 450}