sympy/polynomials/div_.py

   1 """Polynomial division algorithms for use with class Polynomial"""
   2
   3 from sympy.polynomials.base import *
   4 from sympy.core import sympify
   5
   6 def div(f, g, var=None, order=None, coeff=None):
   7     """Division with remainder.
   8
   9     Usage:
  10     ======
  11         The input consists of a polynomial f and either one or a list
  12         of polynomials g. When these are just SymPy expressions, you
  13         can additionally specify the variables and monomial orders
  14         with 'var' and 'order', respectively. Only f is checked for
  15         the input types, the rest is assumed to match.
  16
  17         If 'coeff' is set to 'int', only term divisions with proper
  18         coefficient divisions are allowed. That is, 3*x divides 6*x*y,
  19         but not 2*x**2.
  20
  21         The output's type is Polynomial, but there is a wrapper, see
  22         L{wrapper.div} With only one element in g given, the resulting
  23         list of quotients is unpacked, also. The second output is the
  24         remainder.
  25
  26     Notes:
  27     ======
  28         The loop then iterates over all terms of f, checking if any of
  29         the elements in g (or just g if it is the sole divisor) have a
  30         leading term dividing the current term of f. If yes, the
  31         proper multiple of the element of g is subtracted from f, so
  32         that the term is eliminated, otherwise added to the remainder,
  33         until f is 0.
  34
  35         This way, the algorithm doesn't stop, when the leading terms
  36         of g don't divide the leading term of f, but rather try to
  37         reduce all of f's other terms. Of course, the known univariate
  38         single-polynomial division with remainder is a special case of
  39         this function.
  40
  41         The result can depend on the order of the elements in g. For
  42         unicity, you need that their form a Groebner base of the ideal
  43         they generate, see L{groebner_.groebner}.
  44
  45     Examples:
  46     =========
  47         >>> x, y = symbols('xy')
  48         >>> q, r = div(x**2 + 6*x + 1, 3*x - 1)
  49         >>> print q
  50         19/9 + (1/3)*x
  51         >>> print r
  52         28/9
  53         >>> q, r = div(x**2 + 6*x + 1, 3*x - 1, coeff='int')
  54         >>> print q
  55         2
  56         >>> print r
  57         3 + x**2
  58         >>> q, r = div(2*x**3*y**2 - x*y + y**3, [x-y, y**2], [x,y], 'lex')
  59         >>> print q[0]
  60         -y + 2*y**4 + 2*x*y**3 + 2*x**2*y**2
  61         >>> print q[1]
  62         -1 + y + 2*y**3
  63         >>> print r
  64         0
  65
  66     References:
  67     ===========
  68         Cox, Little, O'Shea: Ideals, Varieties and Algorithms,
  69         Springer, 2. edition, p. 62
  70
  71     """
  72
  73     if not isinstance(g, list):
  74         g = [g]
  75     # Only f is checked, the rest is assumed to match.
  76     if not isinstance(f, Polynomial):
  77         f = sympify(f)
  78         g = map(lambda x: sympify(x), g)
  79         if isinstance(var, Symbol):
  80             var = [var]
  81         if var is None:
  82             var = merge_var(f.atoms(type=Symbol),
  83                             *[g_i.atoms(type=Symbol) for g_i in g])
  84         f = Polynomial(f, var=var, order=order)
  85         g = map(lambda x: Polynomial(x, var=var, order=order), g)
  86
  87     # Begin computation.
  88     r = Polynomial(S.Zero, var=f.var, order=f.order)
  89     q = []
  90     for i in range(0,len(g)):
  91         q.append(r)
  92
  93     while f.sympy_expr is not S.Zero:
  94         for g_i in g:
  95             if g_i.sympy_expr is S.Zero: # Avoid division by 0.
  96                 continue
  97             # Check if leading term of f is divisible by that of g_i.
  98             # When coeff equals 'int', check the coefficient's
  99             # divisibility, too.
 100             td = term_div(f.coeffs[0], g_i.coeffs[0])
 101             if (coeff != 'int' or isinstance(td[0], Integer)) \
 102                and all([e.is_nonnegative for e in td[1:]]):
 103                 quot = Polynomial(coeffs=(td,), var=f.var, order=f.order)
 104                 q[g.index(g_i)] += quot
 105                 f -= quot*g_i
 106                 break
 107         else: # No division occured, add the leading term to remainder.
 108             lt = f.leading_term()
 109             r += lt
 110             f -= lt
 111
 112     if len(q) == 1:
 113         return q[0], r
 114     else:
 115         return q, r
 116
 117
 118 def gcd(f, g, var=None, order=None, coeff=None):
 119     """Greatest common divisor.
 120
 121     Usage:
 122     ======
 123         Input are two polynomials, either as SymPy expressions or as
 124         instances of Polynomial. In the first case, the variables and
 125         monomial order can be specified with 'var' and 'order',
 126         respectively.
 127
 128         If 'coeff' is set to 'int', the content of each polynomial is
 129         checked and their gcd multiplied to the result. Otherwise it
 130         is monic, that is, of leading coefficient 1.
 131
 132         The output's type is Polynomial, but there is a wrapper, see
 133         L{wrapper.gcd}.
 134
 135     Notes:
 136     ======
 137         With only one variable, euclid's algorithm is used directly,
 138         which is reasonably fast. But in the multivariate case, we
 139         have to compute the gcd using the least common multiple, which
 140         relies on Groebner bases. This is based on the formula:
 141         f*g = gcd(f, g)*lcm(f, g)
 142
 143     Examples:
 144     =========
 145         >>> x, y = symbols('xy')
 146         >>> print gcd(4*x**2*y, 6*x*y**2)
 147         x*y
 148         >>> print gcd(4*x**2*y, 6*x*y**2, coeff='int')
 149         2*x*y
 150
 151     References:
 152     ===========
 153         Cox, Little, O'Shea: Ideals, Varieties and Algorithms,
 154         Springer, 2. edition, p. 41 & p. 187
 155
 156     See also L{div}, L{lcm}.
 157
 158     """
 159
 160     # Check if f is a Polynomial already, g is assumed to match.
 161     if not isinstance(f, Polynomial):
 162         f = sympify(f)
 163         g = sympify(g)
 164         if var is None:
 165             var = merge_var(f.atoms(type=Symbol), g.atoms(type=Symbol))
 166         f = Polynomial(f, var=var, order=order)
 167         g = Polynomial(g, var=var, order=order)
 168
 169     # Check if we need to keep an integer factor.
 170     if coeff == 'int':
 171         cf, f = f.as_primitive()
 172         cg, g = g.as_primitive()
 173         c = Integer(numbers.gcd(int(cf), int(cg)))
 174     else:
 175         c = S.One
 176
 177     if len(f.var) == 0: # Constant result.
 178         return Polynomial(c, var=f.var, order=f.order)
 179     elif len(f.var) == 1: # Use euclidean algorithm.
 180         while g.sympy_expr is not S.Zero:
 181             # Remove leading coefficient, to simplify computations.
 182             lc, g = g.as_monic()
 183             f, g = g, div(f, g)[-1]
 184     else: # Use lcm and product to get multivariate gcd.
 185         l = lcm(f, g)
 186         q, r = div(f*g, l)
 187         assert r.sympy_expr is S.Zero
 188         lc, f = q.as_monic()
 189
 190     return Polynomial(coeffs=tuple([(c*t[0],) + t[1:] for t in f.coeffs]),
 191                       var=f.var, order=f.order)
 192
 193
 194 def lcm(f, g, var=None, order=None, coeff=None):
 195     """Least common multiple.
 196
 197     Usage:
 198     ======
 199         Input are two polynomials, either as SymPy expressions or as
 200         instances of Polynomial. In the first case, the variables and
 201         monomial order can be specified with 'var' and 'order',
 202         respectively.
 203
 204         If 'coeff' is set to 'int', the content of each polynomial is
 205         checked and their lcm multiplied to the result. Otherwise it
 206         is monic, that is, of leading coefficient 1.
 207
 208         The output's type is Polynomial, but there is a wrapper, see
 209         L{wrapper.lcm}.
 210
 211     Notes:
 212     ======
 213         With only one variable, the gcd is used to get the lcm from
 214         the product, via f*g = gcd(f, g)*lcm(f, g).
 215
 216         In the univariate case, we compute the unique generator of the
 217         intersection of the two ideals, generated by f and g. This is
 218         done by computing a lexicographic Groebner base of
 219         [t*f. (t-1)*g], with t a dummy variable at the first place,
 220         then filtering trough the base elements not containing t.
 221
 222     Examples:
 223     =========
 224         >>> x, y = symbols('xy')
 225         >>> print lcm(4*x**2*y, 6*x*y**2)
 226         x**2*y**2
 227         >>> print lcm(4*x**2*y, 6*x*y**2, coeff='int')
 228         12*x**2*y**2
 229
 230     References:
 231     ===========
 232         Cox, Little, O'Shea: Ideals, Varieties and Algorithms,
 233         Springer, 2. edition, p. 187
 234
 235     See also L{div}, L{gcd}.
 236
 237     """
 238
 239     # Check if f is a Polynomial already, g is assumed to match.
 240     if not isinstance(f, Polynomial):
 241         f = sympify(f)
 242         g = sympify(g)
 243         if var is None:
 244             var = merge_var(f.atoms(type=Symbol), g.atoms(type=Symbol))
 245         f = Polynomial(f, var=var, order=order)
 246         g = Polynomial(g, var=var, order=order)
 247
 248     # Check if we need to keep an integer factor.
 249     if coeff == 'int':
 250         cf, f = f.as_primitive()
 251         cg, g = g.as_primitive()
 252         cf, cg = int(cf), int(cg)
 253         c = Integer(cf*cg/numbers.gcd(cf, cg))
 254     else:
 255         c = S.One
 256
 257     if len(f.var) == 0: # Constant result.
 258         return Polynomial(c, var=f.var, order=f.order)
 259     elif len(f.var) == 1: # Use gcd to get univariate lcm.
 260         gcd_temp = gcd(f, g)
 261         q, r = div(f*g, gcd_temp)
 262         assert r.sympy_expr is S.Zero
 263         lc, f = q.as_monic()
 264     else:
 265         # Compute a lexicographic Groebner base of the ideal generated
 266         # by t*f and (t-1)*g, with unrelated t.
 267         from sympy.polynomials import groebner_
 268
 269         t = Symbol('t', dummy=True)
 270         var = [t] + f.var
 271         G = groebner_.groebner([Polynomial(t*f.sympy_expr,
 272                                            var=var, order='1-el'),
 273                                 Polynomial((t-1)*g.sympy_expr,
 274                                            var=var, order='1-el')],
 275                                reduced=True)
 276         # Now intersect this result with the polynomial ring in the
 277         # var in `var', that is, eliminate t.
 278         I = filter(lambda p: t not in p.sympy_expr.atoms(type=Symbol), G)
 279         # The intersection should be a principal ideal, that is generated
 280         # by a single polynomial.
 281         if not len(I) == 1:
 282             raise PolynomialException("No single generator for intersection.")
 283         f = Polynomial(I[0].sympy_expr, var=f.var, order=f.order)
 284
 285     return Polynomial(coeffs=tuple([(c*t[0],) + t[1:] for t in f.coeffs]),
 286                       var=f.var, order=f.order)