~dinko-metalac/calculus-app2/trunk

« back to all changes in this revision

Viewing changes to lib/py/sympy/polys/dispersion.py

  • Committer: dinko.metalac at gmail
  • Date: 2015-04-14 13:28:14 UTC
  • Revision ID: dinko.metalac@gmail.com-20150414132814-j25k3qd7sq3warup
new sympy

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
from __future__ import print_function, division
 
2
 
 
3
from sympy.core import S
 
4
from sympy.polys import Poly
 
5
from sympy.polys.polytools import factor_list
 
6
 
 
7
 
 
8
def dispersionset(p, q=None, *gens, **args):
 
9
    r"""Compute the *dispersion set* of two polynomials.
 
10
 
 
11
    For two polynomials `f(x)` and `g(x)` with `\deg f > 0`
 
12
    and `\deg g > 0` the dispersion set `\operatorname{J}(f, g)` is defined as:
 
13
 
 
14
    .. math::
 
15
        \operatorname{J}(f, g)
 
16
        & := \{a \in \mathbb{N}_0 | \gcd(f(x), g(x+a)) \neq 1\} \\
 
17
        &  = \{a \in \mathbb{N}_0 | \deg \gcd(f(x), g(x+a)) \geq 1\}
 
18
 
 
19
    For a single polynomial one defines `\operatorname{J}(f) := \operatorname{J}(f, f)`.
 
20
 
 
21
    Examples
 
22
    ========
 
23
 
 
24
    >>> from sympy import poly
 
25
    >>> from sympy.polys.dispersion import dispersion, dispersionset
 
26
    >>> from sympy.abc import x
 
27
 
 
28
    Dispersion set and dispersion of a simple polynomial:
 
29
 
 
30
    >>> fp = poly((x - 3)*(x + 3), x)
 
31
    >>> sorted(dispersionset(fp))
 
32
    [0, 6]
 
33
    >>> dispersion(fp)
 
34
    6
 
35
 
 
36
    Note that the definition of the dispersion is not symmetric:
 
37
 
 
38
    >>> fp = poly(x**4 - 3*x**2 + 1, x)
 
39
    >>> gp = fp.shift(-3)
 
40
    >>> sorted(dispersionset(fp, gp))
 
41
    [2, 3, 4]
 
42
    >>> dispersion(fp, gp)
 
43
    4
 
44
    >>> sorted(dispersionset(gp, fp))
 
45
    []
 
46
    >>> dispersion(gp, fp)
 
47
    -oo
 
48
 
 
49
    Computing the dispersion also works over field extensions:
 
50
 
 
51
    >>> from sympy import sqrt
 
52
    >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>')
 
53
    >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>')
 
54
    >>> sorted(dispersionset(fp, gp))
 
55
    [2]
 
56
    >>> sorted(dispersionset(gp, fp))
 
57
    [1, 4]
 
58
 
 
59
    We can even perform the computations for polynomials
 
60
    having symbolic coefficients:
 
61
 
 
62
    >>> from sympy.abc import a
 
63
    >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x)
 
64
    >>> sorted(dispersionset(fp))
 
65
    [0, 1]
 
66
 
 
67
    See Also
 
68
    ========
 
69
 
 
70
    dispersion
 
71
 
 
72
    References
 
73
    ==========
 
74
 
 
75
    1. [ManWright94]_
 
76
    2. [Koepf98]_
 
77
    3. [Abramov71]_
 
78
    4. [Man93]_
 
79
    """
 
80
    # Check for valid input
 
81
    same = False if q is not None else True
 
82
    if same:
 
83
        q = p
 
84
 
 
85
    p = Poly(p, *gens, **args)
 
86
    q = Poly(q, *gens, **args)
 
87
 
 
88
    if not p.is_univariate or not q.is_univariate:
 
89
        raise ValueError("Polynomials need to be univariate")
 
90
 
 
91
    # The generator
 
92
    if not p.gen == q.gen:
 
93
        raise ValueError("Polynomials must have the same generator")
 
94
    gen = p.gen
 
95
 
 
96
    # We define the dispersion of constant polynomials to be zero
 
97
    if p.degree() < 1 or q.degree() < 1:
 
98
        return set([0])
 
99
 
 
100
    # Factor p and q over the rationals
 
101
    fp = p.factor_list()
 
102
    fq = q.factor_list() if not same else fp
 
103
 
 
104
    # Iterate over all pairs of factors
 
105
    J = set([])
 
106
    for s, unused in fp[1]:
 
107
        for t, unused in fq[1]:
 
108
            m = s.degree()
 
109
            n = t.degree()
 
110
            if n != m:
 
111
                continue
 
112
            an = s.LC()
 
113
            bn = t.LC()
 
114
            if not (an - bn).is_zero:
 
115
                continue
 
116
            # Note that the roles of `s` and `t` below are switched
 
117
            # w.r.t. the original paper. This is for consistency
 
118
            # with the description in the book of W. Koepf.
 
119
            anm1 = s.coeff_monomial(gen**(m-1))
 
120
            bnm1 = t.coeff_monomial(gen**(n-1))
 
121
            alpha = (anm1 - bnm1) / S(n*bn)
 
122
            if not alpha.is_integer:
 
123
                continue
 
124
            if alpha < 0 or alpha in J:
 
125
                continue
 
126
            if n > 1 and not (s - t.shift(alpha)).is_zero:
 
127
                continue
 
128
            J.add(alpha)
 
129
 
 
130
    return J
 
131
 
 
132
 
 
133
def dispersion(p, q=None, *gens, **args):
 
134
    r"""Compute the *dispersion* of polynomials.
 
135
 
 
136
    For two polynomials `f(x)` and `g(x)` with `\deg f > 0`
 
137
    and `\deg g > 0` the dispersion `\operatorname{dis}(f, g)` is defined as:
 
138
 
 
139
    .. math::
 
140
        \operatorname{dis}(f, g)
 
141
        & := \max\{ J(f,g) \cup \{0\} \} \\
 
142
        &  = \max\{ \{a \in \mathbb{N} | \gcd(f(x), g(x+a)) \neq 1\} \cup \{0\} \}
 
143
 
 
144
    and for a single polynomial `\operatorname{dis}(f) := \operatorname{dis}(f, f)`.
 
145
    Note that we make the definition `\max\{\} := -\infty`.
 
146
 
 
147
    Examples
 
148
    ========
 
149
 
 
150
    >>> from sympy import poly
 
151
    >>> from sympy.polys.dispersion import dispersion, dispersionset
 
152
    >>> from sympy.abc import x
 
153
 
 
154
    Dispersion set and dispersion of a simple polynomial:
 
155
 
 
156
    >>> fp = poly((x - 3)*(x + 3), x)
 
157
    >>> sorted(dispersionset(fp))
 
158
    [0, 6]
 
159
    >>> dispersion(fp)
 
160
    6
 
161
 
 
162
    Note that the definition of the dispersion is not symmetric:
 
163
 
 
164
    >>> fp = poly(x**4 - 3*x**2 + 1, x)
 
165
    >>> gp = fp.shift(-3)
 
166
    >>> sorted(dispersionset(fp, gp))
 
167
    [2, 3, 4]
 
168
    >>> dispersion(fp, gp)
 
169
    4
 
170
    >>> sorted(dispersionset(gp, fp))
 
171
    []
 
172
    >>> dispersion(gp, fp)
 
173
    -oo
 
174
 
 
175
    The maximum of an empty set is defined to be `-\infty`
 
176
    as seen in this example.
 
177
 
 
178
    Computing the dispersion also works over field extensions:
 
179
 
 
180
    >>> from sympy import sqrt
 
181
    >>> fp = poly(x**2 + sqrt(5)*x - 1, x, domain='QQ<sqrt(5)>')
 
182
    >>> gp = poly(x**2 + (2 + sqrt(5))*x + sqrt(5), x, domain='QQ<sqrt(5)>')
 
183
    >>> sorted(dispersionset(fp, gp))
 
184
    [2]
 
185
    >>> sorted(dispersionset(gp, fp))
 
186
    [1, 4]
 
187
 
 
188
    We can even perform the computations for polynomials
 
189
    having symbolic coefficients:
 
190
 
 
191
    >>> from sympy.abc import a
 
192
    >>> fp = poly(4*x**4 + (4*a + 8)*x**3 + (a**2 + 6*a + 4)*x**2 + (a**2 + 2*a)*x, x)
 
193
    >>> sorted(dispersionset(fp))
 
194
    [0, 1]
 
195
 
 
196
    See Also
 
197
    ========
 
198
 
 
199
    dispersionset
 
200
 
 
201
    References
 
202
    ==========
 
203
 
 
204
    1. [ManWright94]_
 
205
    2. [Koepf98]_
 
206
    3. [Abramov71]_
 
207
    4. [Man93]_
 
208
    """
 
209
    J = dispersionset(p, q, *gens, **args)
 
210
    if not J:
 
211
        # Definition for maximum of empty set
 
212
        j = S.NegativeInfinity
 
213
    else:
 
214
        j = max(J)
 
215
    return j