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