1"""Euclidean algorithms, GCDs, LCMs and polynomial remainder sequences.""" 2 3import math 4import operator 5 6from ..config import query 7from ..ntheory import nextprime 8from ..ntheory.modular import crt, symmetric_residue 9from .polyerrors import HeuristicGCDFailed, HomomorphismFailed 10 11 12class _GCD: 13 """Mixin class for computing gcd.""" 14 15 def gcd(self, f, g): 16 """Returns GCD of ``f`` and ``g``.""" 17 if not f and not g: 18 return self.zero 19 elif not f: 20 return self._gcd_zero(g) 21 elif not g: 22 return self._gcd_zero(f) 23 elif f.is_term: 24 return self._gcd_term(f, g) 25 elif g.is_term: 26 return self._gcd_term(g, f) 27 28 J, (f, g) = self._deflate(f, g) 29 h = self._gcd(f, g) 30 31 return self._inflate(h, J) 32 33 def lcm(self, f, g): 34 """Returns LCM of ``f`` and ``g``.""" 35 domain = self.domain 36 37 if not domain.is_Field: 38 fc, f = f.primitive() 39 gc, g = g.primitive() 40 c = domain.lcm(fc, gc) 41 42 h = (f*g)//f.gcd(g) 43 44 if not domain.is_Field: 45 return h*c 46 else: 47 return h.monic() 48 49 def _deflate(self, *polys): 50 J = [0]*self.ngens 51 52 for p in polys: 53 for monom in p: 54 for i, m in enumerate(monom): 55 J[i] = math.gcd(J[i], m) 56 57 for i, b in enumerate(J): 58 if not b: 59 J[i] = 1 60 61 J = tuple(J) 62 63 if all(b == 1 for b in J): 64 return J, polys 65 66 H = [] 67 68 for p in polys: 69 h = self.zero 70 71 for I, coeff in p.items(): 72 N = [i//j for i, j in zip(I, J)] 73 h[N] = coeff 74 75 H.append(h) 76 77 return J, H 78 79 def _inflate(self, f, J): 80 poly = self.zero 81 82 for I, coeff in f.items(): 83 N = [i*j for i, j in zip(I, J)] 84 poly[N] = coeff 85 86 return poly 87 88 def _gcd_zero(self, f): 89 if self.domain.is_Field: 90 return f.monic() 91 else: 92 return f if self.is_normal(f) else -f 93 94 def _gcd_term(self, f, g): 95 domain = self.domain 96 ground_gcd = domain.gcd 97 _mgcd, _cgcd = f.LT 98 if domain.is_Field: 99 for mg, cg in g.items(): 100 _mgcd = _mgcd.gcd(mg) 101 _cgcd = domain.one 102 else: 103 for mg, cg in g.items(): 104 _mgcd = _mgcd.gcd(mg) 105 _cgcd = ground_gcd(_cgcd, cg) 106 return self.term_new(_mgcd, _cgcd) 107 108 def _gcd(self, f, g): 109 domain = self.domain 110 111 if domain.is_RationalField: 112 return self._gcd_QQ(f, g) 113 elif domain.is_IntegerRing: 114 return self._gcd_ZZ(f, g) 115 elif domain.is_AlgebraicField: 116 return self._gcd_AA(f, g) 117 elif not domain.is_Exact: 118 exact = domain.get_exact() 119 ring = self.clone(domain=exact) 120 f, g = map(operator.methodcaller('set_domain', exact), (f, g)) 121 h = ring._gcd(f, g) 122 return h.set_domain(domain) 123 elif domain.is_Field: 124 return self._ff_prs_gcd(f, g) 125 else: 126 return self._rr_prs_gcd(f, g) 127 128 def _gcd_ZZ(self, f, g): 129 from .modulargcd import modgcd 130 131 if query('USE_HEU_GCD'): 132 try: 133 return self._zz_heu_gcd(f, g) 134 except HeuristicGCDFailed: 135 pass 136 137 _gcd_zz_methods = {'modgcd': modgcd, 138 'prs': self._rr_prs_gcd} 139 140 method = _gcd_zz_methods[query('FALLBACK_GCD_ZZ_METHOD')] 141 return method(f, g) 142 143 def _gcd_QQ(self, f, g): 144 domain = self.domain 145 146 _, f = f.clear_denoms(convert=True) 147 _, g = g.clear_denoms(convert=True) 148 149 ring = self.clone(domain=domain.ring) 150 151 h = ring._gcd_ZZ(f, g) 152 h = h.set_ring(self) 153 154 return h.monic() 155 156 def _gcd_AA(self, f, g): 157 from .modulargcd import func_field_modgcd 158 159 _gcd_aa_methods = {'modgcd': func_field_modgcd, 160 'prs': self._ff_prs_gcd} 161 162 method = _gcd_aa_methods[query('GCD_AA_METHOD')] 163 return method(f, g) 164 165 def _zz_heu_gcd(self, f, g): 166 """ 167 Heuristic polynomial GCD in ``Z[X]``. 168 169 Given univariate polynomials ``f`` and ``g`` in ``Z[X]``, returns 170 their GCD, i.e. polynomial ``h``:: 171 172 h = gcd(f, g) 173 174 The algorithm is purely heuristic which means it may fail to compute 175 the GCD. This will be signaled by raising an exception. In this case 176 you will need to switch to another GCD method. 177 178 The algorithm computes the polynomial GCD by evaluating polynomials 179 ``f`` and ``g`` at certain points and computing (fast) integer GCD 180 of those evaluations. The polynomial GCD is recovered from the integer 181 image by interpolation. The evaluation proces reduces f and g variable 182 by variable into a large integer. The final step is to verify if the 183 interpolated polynomial is the correct GCD. 184 185 References 186 ========== 187 188 * :cite:`Liao1995heuristic` 189 190 """ 191 assert self == f.ring == g.ring and self.domain.is_IntegerRing 192 193 x0 = self.gens[0] 194 domain = self.domain 195 196 fc = f.content() 197 gc = g.content() 198 199 gcd = self.domain.gcd(fc, gc) 200 201 f = f.quo_ground(gcd) 202 g = g.quo_ground(gcd) 203 204 f_norm = f.max_norm() 205 g_norm = g.max_norm() 206 207 B = domain(2*min(f_norm, g_norm) + 29) 208 209 x = max(min(B, 99*domain.sqrt(B)), 210 2*min(f_norm // abs(f.LC), 211 g_norm // abs(g.LC)) + 4) 212 213 cofactors = domain.cofactors if self.is_univariate else self.drop(0).cofactors 214 215 for _ in range(query('HEU_GCD_MAX')): 216 ff = f.eval(x0, x) 217 gg = g.eval(x0, x) 218 219 if ff and gg: 220 h, cff, cfg = cofactors(ff, gg) 221 h = self._gcd_interpolate(h, x) 222 h = h.primitive()[1] 223 224 if not f % h and not g % h: 225 return h*gcd 226 227 cff = self._gcd_interpolate(cff, x) 228 229 h, r = divmod(f, cff) 230 231 if not r and not g % h: 232 return h*gcd 233 234 cfg = self._gcd_interpolate(cfg, x) 235 236 h, r = divmod(g, cfg) 237 238 if not r and not f % h: 239 return h*gcd 240 241 x = 73794*x * domain.sqrt(domain.sqrt(x)) // 27011 242 243 raise HeuristicGCDFailed('no luck') 244 245 def _gcd_interpolate(self, h, x): 246 """Interpolate polynomial GCD from integer GCD.""" 247 f, i = self.zero, 0 248 X = self.gens[0] 249 250 while h: 251 g = h % x 252 253 if self.is_univariate: 254 g = symmetric_residue(g, x) 255 256 h = (h - g) // x 257 258 f += X**i*g 259 i += 1 260 261 if not self.is_normal(f): 262 f = -f 263 264 return f 265 266 def _rr_prs_gcd(self, f, g): 267 """Computes polynomial GCD using subresultants over a ring.""" 268 if self.is_multivariate: 269 ring, f, g = map(operator.methodcaller('eject', *self.gens[1:]), 270 (self, f, g)) 271 h = ring._rr_prs_gcd(f, g) 272 return h.inject() 273 274 domain = self.domain 275 276 fc, ff = f.primitive() 277 gc, fg = g.primitive() 278 279 h = ff.subresultants(fg)[-1] 280 _, h = h.primitive() 281 c = domain.gcd(fc, gc) 282 283 return h*c 284 285 def _ff_prs_gcd(self, f, g): 286 """Computes polynomial GCD using subresultants over a field.""" 287 if self.is_multivariate: 288 ring, F, G = map(operator.methodcaller('eject', *self.gens[1:]), 289 (self, f, g)) 290 291 fc, F = F.primitive() 292 gc, G = G.primitive() 293 294 F, G = map(operator.methodcaller('inject'), (F, G)) 295 296 h = F.subresultants(G)[-1] 297 c = ring.domain._ff_prs_gcd(fc, gc) 298 h = h.eject(*self.gens[1:]) 299 _, h = h.primitive() 300 h = h.inject() 301 h *= c 302 303 return h.monic() 304 305 h = f.subresultants(g)[-1] 306 307 return h.monic() 308 309 def _primitive_prs(self, f, g): 310 """ 311 Subresultant PRS algorithm in `K[X]`. 312 313 Computes the last non-zero scalar subresultant of `f` and `g` 314 and subresultant polynomial remainder sequence (PRS). 315 316 The first subdeterminant is set to 1 by convention to match 317 the polynomial and the scalar subdeterminants. 318 If 'deg(f) < deg(g)', the subresultants of '(g,f)' are computed. 319 320 Examples 321 ======== 322 323 >>> _, x = ring('x', ZZ) 324 325 >>> (x**2 + 1).resultant(x**2 - 1, includePRS=True) 326 (4, [x**2 + 1, x**2 - 1, -2]) 327 328 References 329 ========== 330 331 * :cite:`Brown1978prs` 332 * :cite:`Geddes1992algorithms`, example 7.6 333 334 """ 335 domain = self.domain 336 337 if self.is_multivariate: 338 ring, f, g = map(operator.methodcaller('eject', *self.gens[1:]), 339 (self, f, g)) 340 res = ring._primitive_prs(f, g) 341 return res[0], [_.inject() for _ in res[1]] 342 343 n = f.degree() 344 m = g.degree() 345 346 if n < m: 347 res, sub = self._primitive_prs(g, f) 348 if res: 349 res *= (-domain.one)**(n*m) 350 return res, sub 351 352 c, r = domain.zero, [] 353 354 if not f: 355 return c, r 356 357 r.append(f) 358 359 if not g: 360 return c, r 361 362 r.append(g) 363 d = n - m 364 365 h = f.prem(g) 366 h *= (-domain.one)**(d + 1) 367 368 lc = g.LC 369 c = -lc**d 370 371 while h: 372 k = h.degree() 373 r.append(h) 374 375 f, g, m, d = g, h, k, m - k 376 377 h = f.prem(g) 378 h = h.quo_ground(-lc*c**d) 379 380 lc = g.LC 381 c = domain.quo((-lc)**d, c**(d - 1)) 382 383 if r[-1].degree() > 0: 384 c = domain.zero 385 else: 386 c = -c 387 388 return c, r 389 390 def _collins_resultant(self, f, g): 391 """ 392 Collins's modular resultant algorithm in `ZZ[X]` or `QQ[X]`. 393 394 References 395 ========== 396 397 * :cite:`Collins1971mod`, algorithm PRES 398 399 """ 400 domain = self.domain 401 402 if not f or not g: 403 return self.drop(0).zero 404 405 n = f.degree() 406 m = g.degree() 407 408 if domain.is_RationalField: 409 cf, f = f.clear_denoms(convert=True) 410 cg, g = g.clear_denoms(convert=True) 411 412 ring = self.clone(domain=domain.ring) 413 r = ring._collins_resultant(f, g) 414 r = r.set_domain(domain) 415 416 c = cf**n * cg**m 417 c = domain.convert(c, domain.ring) 418 419 return r.quo_ground(c) 420 421 assert domain.is_IntegerRing 422 423 A = f.max_norm() 424 B = g.max_norm() 425 426 a, b = f.LC, g.LC 427 428 B = domain(2)*domain.factorial(domain(n + m))*A**m*B**n 429 new_ring = self.drop(0) 430 r, p, P = new_ring.zero, domain.one, domain.one 431 432 while P <= B: 433 while True: 434 p = domain(nextprime(p)) 435 if a % p and b % p: 436 break 437 438 p_domain = domain.finite_field(p) 439 F, G = map(operator.methodcaller('set_domain', p_domain), (f, g)) 440 441 try: 442 R = self.clone(domain=p_domain)._modular_resultant(F, G) 443 except HomomorphismFailed: 444 continue 445 446 if P == 1: 447 r = R 448 else: 449 def _crt(r, R): 450 return domain(crt([P, p], map(domain.convert, [r, R]), 451 check=False, symmetric=True)[0]) 452 453 if new_ring.is_PolynomialRing: 454 r_new = new_ring.zero 455 456 for monom in r | R: 457 r_new[monom] = _crt(r.get(monom, 0), R.get(monom, 0)) 458 r = r_new 459 else: 460 r = _crt(r, R) 461 462 P *= p 463 464 return r 465 466 def _modular_resultant(self, f, g): 467 """ 468 Compute resultant of `f` and `g` in `GF(p)[X]`. 469 470 References 471 ========== 472 473 * :cite:`Collins1971mod`, algorithm CPRES 474 475 """ 476 domain = self.domain 477 478 assert domain.is_FiniteField 479 480 if self.is_univariate: 481 return self._primitive_prs(f, g)[0] 482 483 n = f.degree() 484 m = g.degree() 485 486 N = f.degree(1) 487 M = g.degree(1) 488 489 B = n*M + m*N 490 491 new_ring = self.drop(0) 492 r = new_ring.zero 493 D = self.eject(1).domain.one 494 domain_elts = iter(range(domain.order)) 495 496 while D.degree() <= B: 497 while True: 498 try: 499 a = next(domain_elts) 500 except StopIteration: 501 raise HomomorphismFailed('no luck') 502 503 F = f.eval(x=1, a=a) 504 505 if F.degree() == n: 506 G = g.eval(x=1, a=a) 507 508 if G.degree() == m: 509 break 510 511 R = self.drop(1)._modular_resultant(F, G) 512 e = r.eval(x=0, a=a) 513 514 if new_ring.is_univariate: 515 R = new_ring.ground_new(R) 516 e = new_ring.ground_new(e) 517 else: 518 R = R.set_ring(new_ring) 519 e = e.set_ring(new_ring) 520 521 d = D * D.eval(x=0, a=a)**-1 522 d = d.set_ring(new_ring) 523 r += d*(R - e) 524 D *= D.ring.gens[0] - a 525 526 return r 527