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