Source code for sympy.polys.ring_series

"""Power series evaluation and manipulation using sparse Polynomials

Implementing a new function
---------------------------

There are a few things to be kept in mind when adding a new function here::

    - The implementation should work on all possible input domains/rings.
      Special cases include the ``EX`` ring and a constant term in the series
      to be expanded. There can be two types of constant terms in the series:

        + A constant value or symbol.
        + A term of a multivariate series not involving the generator, with
          respect to which the series is to expanded.

      Strictly speaking, a generator of a ring should not be considered a
      constant. However, for series expansion both the cases need similar
      treatment (as the user doesn't care about inner details), i.e, use an
      addition formula to separate the constant part and the variable part (see
      rs_sin for reference).

    - All the algorithms used here are primarily designed to work for Taylor
      series (number of iterations in the algo equals the required order).
      Hence, it becomes tricky to get the series of the right order if a
      Puiseux series is input. Use rs_puiseux? in your function if your
      algorithm is not designed to handle fractional powers.

Extending rs_series
-------------------

To make a function work with rs_series you need to do two things::

    - Many sure it works with a constant term (as explained above).
    - If the series contains constant terms, you might need to extend its ring.
      You do so by adding the new terms to the rings as generators.
      ``PolyRing.compose`` and ``PolyRing.add_gens`` are two functions that do
      so and need to be called every time you expand a series containing a
      constant term.

Look at rs_sin and rs_series for further reference.

"""

from sympy.polys.domains import QQ, EX
from sympy.polys.rings import PolyElement, ring, sring
from sympy.polys.polyerrors import DomainError
from sympy.polys.monomials import (monomial_min, monomial_mul, monomial_div,
                                   monomial_ldiv)
from mpmath.libmp.libintmath import ifac
from sympy.core import PoleError, Function, Expr
from sympy.core.numbers import Rational, igcd
from sympy.core.compatibility import as_int, range
from sympy.functions import sin, cos, tan, atan, exp, atanh, tanh, log, ceiling
from mpmath.libmp.libintmath import giant_steps
import math


def _invert_monoms(p1):
    """
    Compute ``x**n * p1(1/x)`` for a univariate polynomial ``p1`` in ``x``.

    Examples
    ========

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.rings import ring
    >>> from sympy.polys.ring_series import _invert_monoms
    >>> R, x = ring('x', ZZ)
    >>> p = x**2 + 2*x + 3
    >>> _invert_monoms(p)
    3*x**2 + 2*x + 1

    See Also
    ========

    sympy.polys.densebasic.dup_reverse
    """
    terms = list(p1.items())
    terms.sort()
    deg = p1.degree()
    R = p1.ring
    p = R.zero
    cv = p1.listcoeffs()
    mv = p1.listmonoms()
    for i in range(len(mv)):
        p[(deg - mv[i][0],)] = cv[i]
    return p

def _giant_steps(target):
    """Return a list of precision steps for the Newton's method"""
    res = giant_steps(2, target)
    if res[0] != 2:
        res = [2] + res
    return res

[docs]def rs_trunc(p1, x, prec): """ Truncate the series in the ``x`` variable with precision ``prec``, that is, modulo ``O(x**prec)`` Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_trunc >>> R, x = ring('x', QQ) >>> p = x**10 + x**5 + x + 1 >>> rs_trunc(p, x, 12) x**10 + x**5 + x + 1 >>> rs_trunc(p, x, 10) x**5 + x + 1 """ R = p1.ring p = R.zero i = R.gens.index(x) for exp1 in p1: if exp1[i] >= prec: continue p[exp1] = p1[exp1] return p
[docs]def rs_is_puiseux(p, x): """ Test if ``p`` is Puiseux series in ``x``. Raise an exception if it has a negative power in ``x``. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_is_puiseux >>> R, x = ring('x', QQ) >>> p = x**QQ(2,5) + x**QQ(2,3) + x >>> rs_is_puiseux(p, x) True """ index = p.ring.gens.index(x) for k in p: if k[index] != int(k[index]): return True if k[index] < 0: raise ValueError('The series is not regular in %s' % x) return False
[docs]def rs_puiseux(f, p, x, prec): """ Return the puiseux series for `f(p, x, prec)`. To be used when function ``f`` is implemented only for regular series. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_puiseux, rs_exp >>> R, x = ring('x', QQ) >>> p = x**QQ(2,5) + x**QQ(2,3) + x >>> rs_puiseux(rs_exp,p, x, 1) 1/2*x**(4/5) + x**(2/3) + x**(2/5) + 1 """ index = p.ring.gens.index(x) n = 1 for k in p: power = k[index] if isinstance(power, Rational): num, den = power.as_numer_denom() n = int(n*den // igcd(n, den)) elif power != int(power): num, den = power.numerator, power.denominator n = int(n*den // igcd(n, den)) if n != 1: p1 = pow_xin(p, index, n) r = f(p1, x, prec*n) n1 = QQ(1, n) if isinstance(r, tuple): r = tuple([pow_xin(rx, index, n1) for rx in r]) else: r = pow_xin(r, index, n1) else: r = f(p, x, prec) return r
[docs]def rs_puiseux2(f, p, q, x, prec): """ Return the puiseux series for `f(p, q, x, prec)`. To be used when function ``f`` is implemented only for regular series. """ index = p.ring.gens.index(x) n = 1 for k in p: power = k[index] if isinstance(power, Rational): num, den = power.as_numer_denom() n = n*den // igcd(n, den) elif power != int(power): num, den = power.numerator, power.denominator n = n*den // igcd(n, den) if n != 1: p1 = pow_xin(p, index, n) r = f(p1, q, x, prec*n) n1 = QQ(1, n) r = pow_xin(r, index, n1) else: r = f(p, q, x, prec) return r
[docs]def rs_mul(p1, p2, x, prec): """ Return the product of the given two series, modulo ``O(x**prec)``. ``x`` is the series variable or its position in the generators. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_mul >>> R, x = ring('x', QQ) >>> p1 = x**2 + 2*x + 1 >>> p2 = x + 1 >>> rs_mul(p1, p2, x, 3) 3*x**2 + 3*x + 1 """ R = p1.ring p = R.zero if R.__class__ != p2.ring.__class__ or R != p2.ring: raise ValueError('p1 and p2 must have the same ring') iv = R.gens.index(x) if not isinstance(p2, PolyElement): raise ValueError('p1 and p2 must have the same ring') if R == p2.ring: get = p.get items2 = list(p2.items()) items2.sort(key=lambda e: e[0][iv]) if R.ngens == 1: for exp1, v1 in p1.items(): for exp2, v2 in items2: exp = exp1[0] + exp2[0] if exp < prec: exp = (exp, ) p[exp] = get(exp, 0) + v1*v2 else: break else: monomial_mul = R.monomial_mul for exp1, v1 in p1.items(): for exp2, v2 in items2: if exp1[iv] + exp2[iv] < prec: exp = monomial_mul(exp1, exp2) p[exp] = get(exp, 0) + v1*v2 else: break p.strip_zero() return p
[docs]def rs_square(p1, x, prec): """ Square the series modulo ``O(x**prec)`` Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_square >>> R, x = ring('x', QQ) >>> p = x**2 + 2*x + 1 >>> rs_square(p, x, 3) 6*x**2 + 4*x + 1 """ R = p1.ring p = R.zero iv = R.gens.index(x) get = p.get items = list(p1.items()) items.sort(key=lambda e: e[0][iv]) monomial_mul = R.monomial_mul for i in range(len(items)): exp1, v1 = items[i] for j in range(i): exp2, v2 = items[j] if exp1[iv] + exp2[iv] < prec: exp = monomial_mul(exp1, exp2) p[exp] = get(exp, 0) + v1*v2 else: break p = p.imul_num(2) get = p.get for expv, v in p1.items(): if 2*expv[iv] < prec: e2 = monomial_mul(expv, expv) p[e2] = get(e2, 0) + v**2 p.strip_zero() return p
[docs]def rs_pow(p1, n, x, prec): """ Return ``p1**n`` modulo ``O(x**prec)`` Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_pow >>> R, x = ring('x', QQ) >>> p = x + 1 >>> rs_pow(p, 4, x, 3) 6*x**2 + 4*x + 1 """ R = p1.ring p = R.zero if isinstance(n, Rational): np = int(n.p) nq = int(n.q) if nq != 1: res = rs_nth_root(p1, nq, x, prec) if np != 1: res = rs_pow(res, np, x, prec) else: res = rs_pow(p1, np, x, prec) return res n = as_int(n) if n == 0: if p1: return R(1) else: raise ValueError('0**0 is undefined') if n < 0: p1 = rs_pow(p1, -n, x, prec) return rs_series_inversion(p1, x, prec) if n == 1: return rs_trunc(p1, x, prec) if n == 2: return rs_square(p1, x, prec) if n == 3: p2 = rs_square(p1, x, prec) return rs_mul(p1, p2, x, prec) p = R(1) while 1: if n & 1: p = rs_mul(p1, p, x, prec) n -= 1 if not n: break p1 = rs_square(p1, x, prec) n = n // 2 return p
[docs]def rs_subs(p, rules, x, prec): """ Substitution with truncation according to the mapping in ``rules``. Return a series with precision ``prec`` in the generator ``x`` Note that substitutions are not done one after the other >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_subs >>> R, x, y = ring('x, y', QQ) >>> p = x**2 + y**2 >>> rs_subs(p, {x: x+ y, y: x+ 2*y}, x, 3) 2*x**2 + 6*x*y + 5*y**2 >>> (x + y)**2 + (x + 2*y)**2 2*x**2 + 6*x*y + 5*y**2 which differs from >>> rs_subs(rs_subs(p, {x: x+ y}, x, 3), {y: x+ 2*y}, x, 3) 5*x**2 + 12*x*y + 8*y**2 Parameters ---------- p : :class:`PolyElement` Input series. rules : :class:`dict` with substitution mappings. x : :class:`PolyElement` in which the series truncation is to be done. prec : :class:`Integer` order of the series after truncation. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_subs >>> R, x, y = ring('x, y', QQ) >>> rs_subs(x**2+y**2, {y: (x+y)**2}, x, 3) 6*x**2*y**2 + x**2 + 4*x*y**3 + y**4 """ R = p.ring ngens = R.ngens d = R(0) for i in range(ngens): d[(i, 1)] = R.gens[i] for var in rules: d[(R.index(var), 1)] = rules[var] p1 = R(0) p_keys = sorted(p.keys()) for expv in p_keys: p2 = R(1) for i in range(ngens): power = expv[i] if power == 0: continue if (i, power) not in d: q, r = divmod(power, 2) if r == 0 and (i, q) in d: d[(i, power)] = rs_square(d[(i, q)], x, prec) elif (i, power - 1) in d: d[(i, power)] = rs_mul(d[(i, power - 1)], d[(i, 1)], x, prec) else: d[(i, power)] = rs_pow(d[(i, 1)], power, x, prec) p2 = rs_mul(p2, d[(i, power)], x, prec) p1 += p2*p[expv] return p1
def _has_constant_term(p, x): """ Check if ``p`` has a constant term in ``x`` Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import _has_constant_term >>> R, x = ring('x', QQ) >>> p = x**2 + x + 1 >>> _has_constant_term(p, x) True """ R = p.ring iv = R.gens.index(x) zm = R.zero_monom a = [0]*R.ngens a[iv] = 1 miv = tuple(a) for expv in p: if monomial_min(expv, miv) == zm: return True return False def _get_constant_term(p, x): """Return constant term in p with respect to x Note that it is not simply `p[R.zero_monom]` as there might be multiple generators in the ring R. We want the `x`-free term which can contain other generators. """ R = p.ring zm = R.zero_monom i = R.gens.index(x) zm = R.zero_monom a = [0]*R.ngens a[i] = 1 miv = tuple(a) c = 0 for expv in p: if monomial_min(expv, miv) == zm: c += R({expv: p[expv]}) return c def _check_series_var(p, x, name): index = p.ring.gens.index(x) m = min(p, key=lambda k: k[index])[index] if m < 0: raise PoleError("Asymptotic expansion of %s around [oo] not " "implemented." % name) return index, m def _series_inversion1(p, x, prec): """ Univariate series inversion ``1/p`` modulo ``O(x**prec)``. The Newton method is used. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import _series_inversion1 >>> R, x = ring('x', QQ) >>> p = x + 1 >>> _series_inversion1(p, x, 4) -x**3 + x**2 - x + 1 """ if rs_is_puiseux(p, x): return rs_puiseux(_series_inversion1, p, x, prec) R = p.ring zm = R.zero_monom c = p[zm] # giant_steps does not seem to work with PythonRational numbers with 1 as # denominator. This makes sure such a number is converted to integer. if prec == int(prec): prec = int(prec) if zm not in p: raise ValueError("No constant term in series") if _has_constant_term(p - c, x): raise ValueError("p cannot contain a constant term depending on " "parameters") one = R(1) if R.domain is EX: one = 1 if c != one: # TODO add check that it is a unit p1 = R(1)/c else: p1 = R(1) for precx in _giant_steps(prec): t = 1 - rs_mul(p1, p, x, precx) p1 = p1 + rs_mul(p1, t, x, precx) return p1
[docs]def rs_series_inversion(p, x, prec): """ Multivariate series inversion ``1/p`` modulo ``O(x**prec)``. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_series_inversion >>> R, x, y = ring('x, y', QQ) >>> rs_series_inversion(1 + x*y**2, x, 4) -x**3*y**6 + x**2*y**4 - x*y**2 + 1 >>> rs_series_inversion(1 + x*y**2, y, 4) -x*y**2 + 1 >>> rs_series_inversion(x + x**2, x, 4) x**3 - x**2 + x - 1 + x**(-1) """ R = p.ring if p == R.zero: raise ZeroDivisionError zm = R.zero_monom index = R.gens.index(x) m = min(p, key=lambda k: k[index])[index] if m: p = mul_xin(p, index, -m) prec = prec + m if zm not in p: raise NotImplementedError("No constant term in series") if _has_constant_term(p - p[zm], x): raise NotImplementedError("p - p[0] must not have a constant term in " "the series variables") r = _series_inversion1(p, x, prec) if m != 0: r = mul_xin(r, index, -m) return r
def _coefficient_t(p, t): r"""Coefficient of `x\_i**j` in p, where ``t`` = (i, j)""" i, j = t R = p.ring expv1 = [0]*R.ngens expv1[i] = j expv1 = tuple(expv1) p1 = R(0) for expv in p: if expv[i] == j: p1[monomial_div(expv, expv1)] = p[expv] return p1
[docs]def rs_series_reversion(p, x, n, y): r""" Reversion of a series. ``p`` is a series with ``O(x**n)`` of the form `p = a*x + f(x)` where `a` is a number different from 0. `f(x) = sum( a\_k*x\_k, k in range(2, n))` a_k : Can depend polynomially on other variables, not indicated. x : Variable with name x. y : Variable with name y. Solve `p = y`, that is, given `a*x + f(x) - y = 0`, find the solution x = r(y) up to O(y**n) Algorithm: If `r\_i` is the solution at order i, then: `a*r\_i + f(r\_i) - y = O(y**(i + 1))` and if r_(i + 1) is the solution at order i + 1, then: `a*r\_(i + 1) + f(r\_(i + 1)) - y = O(y**(i + 2))` We have, r_(i + 1) = r_i + e, such that, `a*e + f(r\_i) = O(y**(i + 2))` or `e = -f(r\_i)/a` So we use the recursion relation: `r\_(i + 1) = r\_i - f(r\_i)/a` with the boundary condition: `r\_1 = y` Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_series_reversion, rs_trunc >>> R, x, y, a, b = ring('x, y, a, b', QQ) >>> p = x - x**2 - 2*b*x**2 + 2*a*b*x**2 >>> p1 = rs_series_reversion(p, x, 3, y); p1 -2*y**2*a*b + 2*y**2*b + y**2 + y >>> rs_trunc(p.compose(x, p1), y, 3) y """ if rs_is_puiseux(p, x): raise NotImplementedError R = p.ring nx = R.gens.index(x) y = R(y) ny = R.gens.index(y) if _has_constant_term(p, x): raise ValueError("p must not contain a constant term in the series " "variable") a = _coefficient_t(p, (nx, 1)) zm = R.zero_monom assert zm in a and len(a) == 1 a = a[zm] r = y/a for i in range(2, n): sp = rs_subs(p, {x: r}, y, i + 1) sp = _coefficient_t(sp, (ny, i))*y**i r -= sp/a return r
[docs]def rs_series_from_list(p, c, x, prec, concur=1): """ Return a series `sum c[n]*p**n` modulo `O(x**prec)`. It reduces the number of multiplications by summing concurrently. `ax = [1, p, p**2, .., p**(J - 1)]` `s = sum(c[i]*ax[i]` for i in `range(r, (r + 1)*J))*p**((K - 1)*J)` with `K >= (n + 1)/J` Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_series_from_list, rs_trunc >>> R, x = ring('x', QQ) >>> p = x**2 + x + 1 >>> c = [1, 2, 3] >>> rs_series_from_list(p, c, x, 4) 6*x**3 + 11*x**2 + 8*x + 6 >>> rs_trunc(1 + 2*p + 3*p**2, x, 4) 6*x**3 + 11*x**2 + 8*x + 6 >>> pc = R.from_list(list(reversed(c))) >>> rs_trunc(pc.compose(x, p), x, 4) 6*x**3 + 11*x**2 + 8*x + 6 See Also ======== sympy.polys.ring.compose """ R = p.ring n = len(c) if not concur: q = R(1) s = c[0]*q for i in range(1, n): q = rs_mul(q, p, x, prec) s += c[i]*q return s J = int(math.sqrt(n) + 1) K, r = divmod(n, J) if r: K += 1 ax = [R(1)] b = 1 q = R(1) if len(p) < 20: for i in range(1, J): q = rs_mul(q, p, x, prec) ax.append(q) else: for i in range(1, J): if i % 2 == 0: q = rs_square(ax[i//2], x, prec) else: q = rs_mul(q, p, x, prec) ax.append(q) # optimize using rs_square pj = rs_mul(ax[-1], p, x, prec) b = R(1) s = R(0) for k in range(K - 1): r = J*k s1 = c[r] for j in range(1, J): s1 += c[r + j]*ax[j] s1 = rs_mul(s1, b, x, prec) s += s1 b = rs_mul(b, pj, x, prec) if not b: break k = K - 1 r = J*k if r < n: s1 = c[r]*R(1) for j in range(1, J): if r + j >= n: break s1 += c[r + j]*ax[j] s1 = rs_mul(s1, b, x, prec) s += s1 return s
[docs]def rs_diff(p, x): """ Return partial derivative of ``p`` with respect to ``x``. Parameters ---------- x : :class:`PolyElement` with respect to which ``p`` is differentiated. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_diff >>> R, x, y = ring('x, y', QQ) >>> p = x + x**2*y**3 >>> rs_diff(p, x) 2*x*y**3 + 1 """ R = p.ring n = R.gens.index(x) p1 = R.zero mn = [0]*R.ngens mn[n] = 1 mn = tuple(mn) for expv in p: if expv[n]: e = monomial_ldiv(expv, mn) p1[e] = R.domain_new(p[expv]*expv[n]) return p1
[docs]def rs_integrate(p, x): """ Integrate ``p`` with respect to ``x``. Parameters ---------- x : :class:`PolyElement` with respect to which ``p`` is integrated. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_integrate >>> R, x, y = ring('x, y', QQ) >>> p = x + x**2*y**3 >>> rs_integrate(p, x) 1/3*x**3*y**3 + 1/2*x**2 """ R = p.ring p1 = R.zero n = R.gens.index(x) mn = [0]*R.ngens mn[n] = 1 mn = tuple(mn) for expv in p: e = monomial_mul(expv, mn) p1[e] = R.domain_new(p[expv]/(expv[n] + 1)) return p1
[docs]def rs_fun(p, f, *args): r""" Function of a multivariate series computed by substitution. The case with f method name is used to compute `rs\_tan` and `rs\_nth\_root` of a multivariate series: `rs\_fun(p, tan, iv, prec)` tan series is first computed for a dummy variable _x, i.e, `rs\_tan(\_x, iv, prec)`. Then we substitute _x with p to get the desired series Parameters ---------- p : :class:`PolyElement` The multivariate series to be expanded. f : `ring\_series` function to be applied on `p`. args[-2] : :class:`PolyElement` with respect to which, the series is to be expanded. args[-1] : Required order of the expanded series. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_fun, _tan1 >>> R, x, y = ring('x, y', QQ) >>> p = x + x*y + x**2*y + x**3*y**2 >>> rs_fun(p, _tan1, x, 4) 1/3*x**3*y**3 + 2*x**3*y**2 + x**3*y + 1/3*x**3 + x**2*y + x*y + x """ _R = p.ring R1, _x = ring('_x', _R.domain) h = int(args[-1]) args1 = args[:-2] + (_x, h) zm = _R.zero_monom # separate the constant term of the series # compute the univariate series f(_x, .., 'x', sum(nv)) if zm in p: x1 = _x + p[zm] p1 = p - p[zm] else: x1 = _x p1 = p if isinstance(f, str): q = getattr(x1, f)(*args1) else: q = f(x1, *args1) a = sorted(q.items()) c = [0]*h for x in a: c[x[0][0]] = x[1] p1 = rs_series_from_list(p1, c, args[-2], args[-1]) return p1
[docs]def mul_xin(p, i, n): r""" Return `p*x_i**n`. `x\_i` is the ith variable in ``p``. """ R = p.ring q = R(0) for k, v in p.items(): k1 = list(k) k1[i] += n q[tuple(k1)] = v return q
[docs]def pow_xin(p, i, n): """ >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import pow_xin >>> R, x, y = ring('x, y', QQ) >>> p = x**QQ(2,5) + x + x**QQ(2,3) >>> index = p.ring.gens.index(x) >>> pow_xin(p, index, 15) x**15 + x**10 + x**6 """ R = p.ring q = R(0) for k, v in p.items(): k1 = list(k) k1[i] *= n q[tuple(k1)] = v return q
def _nth_root1(p, n, x, prec): """ Univariate series expansion of the nth root of ``p``. The Newton method is used. """ if rs_is_puiseux(p, x): return rs_puiseux2(_nth_root1, p, n, x, prec) R = p.ring zm = R.zero_monom if zm not in p: raise NotImplementedError('No constant term in series') n = as_int(n) assert p[zm] == 1 p1 = R(1) if p == 1: return p if n == 0: return R(1) if n == 1: return p if n < 0: n = -n sign = 1 else: sign = 0 for precx in _giant_steps(prec): tmp = rs_pow(p1, n + 1, x, precx) tmp = rs_mul(tmp, p, x, precx) p1 += p1/n - tmp/n if sign: return p1 else: return _series_inversion1(p1, x, prec)
[docs]def rs_nth_root(p, n, x, prec): """ Multivariate series expansion of the nth root of ``p``. Parameters ---------- n : `p**(1/n)` is returned. x : :class:`PolyElement` prec : Order of the expanded series. Notes ===== The result of this function is dependent on the ring over which the polynomial has been defined. If the answer involves a root of a constant, make sure that the polynomial is over a real field. It can not yet handle roots of symbols. Examples ======== >>> from sympy.polys.domains import QQ, RR >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_nth_root >>> R, x, y = ring('x, y', QQ) >>> rs_nth_root(1 + x + x*y, -3, x, 3) 2/9*x**2*y**2 + 4/9*x**2*y + 2/9*x**2 - 1/3*x*y - 1/3*x + 1 >>> R, x, y = ring('x, y', RR) >>> rs_nth_root(3 + x + x*y, 3, x, 2) 0.160249952256379*x*y + 0.160249952256379*x + 1.44224957030741 """ p0 = p n0 = n if n == 0: if p == 0: raise ValueError('0**0 expression') else: return p.ring(1) if n == 1: return rs_trunc(p, x, prec) R = p.ring zm = R.zero_monom index = R.gens.index(x) m = min(p, key=lambda k: k[index])[index] p = mul_xin(p, index, -m) prec -= m if _has_constant_term(p - 1, x): zm = R.zero_monom c = p[zm] if R.domain is EX: c_expr = c.as_expr() const = c_expr**QQ(1, n) elif isinstance(c, PolyElement): try: c_expr = c.as_expr() const = R(c_expr**(QQ(1, n))) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") else: try: # RealElement doesn't support const = R(c**Rational(1, n)) # exponentiation with mpq object except ValueError: # as exponent raise DomainError("The given series can't be expanded in " "this domain.") res = rs_nth_root(p/c, n, x, prec)*const else: res = _nth_root1(p, n, x, prec) if m: m = QQ(m, n) res = mul_xin(res, index, m) return res
[docs]def rs_log(p, x, prec): """ The Logarithm of ``p`` modulo ``O(x**prec)``. Notes ===== Truncation of ``integral dx p**-1*d p/dx`` is used. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_log >>> R, x = ring('x', QQ) >>> rs_log(1 + x, x, 8) 1/7*x**7 - 1/6*x**6 + 1/5*x**5 - 1/4*x**4 + 1/3*x**3 - 1/2*x**2 + x >>> rs_log(x**QQ(3, 2) + 1, x, 5) 1/3*x**(9/2) - 1/2*x**3 + x**(3/2) """ if rs_is_puiseux(p, x): return rs_puiseux(rs_log, p, x, prec) R = p.ring if p == 1: return R.zero c = _get_constant_term(p, x) if c: const = 0 if c == 1: pass else: c_expr = c.as_expr() if R.domain is EX: const = log(c_expr) elif isinstance(c, PolyElement): try: const = R(log(c_expr)) except ValueError: R = R.add_gens([log(c_expr)]) p = p.set_ring(R) x = x.set_ring(R) c = c.set_ring(R) const = R(log(c_expr)) else: try: const = R(log(c)) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") dlog = p.diff(x) dlog = rs_mul(dlog, _series_inversion1(p, x, prec), x, prec - 1) return rs_integrate(dlog, x) + const else: raise NotImplementedError
[docs]def rs_LambertW(p, x, prec): """ Calculate the series expansion of the principal branch of the Lambert W function. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_LambertW >>> R, x, y = ring('x, y', QQ) >>> rs_LambertW(x + x*y, x, 3) -x**2*y**2 - 2*x**2*y - x**2 + x*y + x See Also ======== LambertW """ if rs_is_puiseux(p, x): return rs_puiseux(rs_LambertW, p, x, prec) R = p.ring p1 = R(0) if _has_constant_term(p, x): raise NotImplementedError("Polynomial must not have constant term in " "the series variables") if x in R.gens: for precx in _giant_steps(prec): e = rs_exp(p1, x, precx) p2 = rs_mul(e, p1, x, precx) - p p3 = rs_mul(e, p1 + 1, x, precx) p3 = rs_series_inversion(p3, x, precx) tmp = rs_mul(p2, p3, x, precx) p1 -= tmp return p1 else: raise NotImplementedError
def _exp1(p, x, prec): r"""Helper function for `rs\_exp`. """ R = p.ring p1 = R(1) for precx in _giant_steps(prec): pt = p - rs_log(p1, x, precx) tmp = rs_mul(pt, p1, x, precx) p1 += tmp return p1
[docs]def rs_exp(p, x, prec): """ Exponentiation of a series modulo ``O(x**prec)`` Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_exp >>> R, x = ring('x', QQ) >>> rs_exp(x**2, x, 7) 1/6*x**6 + 1/2*x**4 + x**2 + 1 """ if rs_is_puiseux(p, x): return rs_puiseux(rs_exp, p, x, prec) R = p.ring c = _get_constant_term(p, x) if c: if R.domain is EX: c_expr = c.as_expr() const = exp(c_expr) elif isinstance(c, PolyElement): try: c_expr = c.as_expr() const = R(exp(c_expr)) except ValueError: R = R.add_gens([exp(c_expr)]) p = p.set_ring(R) x = x.set_ring(R) c = c.set_ring(R) const = R(exp(c_expr)) else: try: const = R(exp(c)) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") p1 = p - c # Makes use of sympy fuctions to evaluate the values of the cos/sin # of the constant term. return const*rs_exp(p1, x, prec) if len(p) > 20: return _exp1(p, x, prec) one = R(1) n = 1 k = 1 c = [] for k in range(prec): c.append(one/n) k += 1 n *= k r = rs_series_from_list(p, c, x, prec) return r
def _atan(p, iv, prec): """ Expansion using formula. Faster on very small and univariate series. """ R = p.ring mo = R(-1) c = [-mo] p2 = rs_square(p, iv, prec) for k in range(1, prec): c.append(mo**k/(2*k + 1)) s = rs_series_from_list(p2, c, iv, prec) s = rs_mul(s, p, iv, prec) return s
[docs]def rs_atan(p, x, prec): """ The arctangent of a series Return the series expansion of the atan of ``p``, about 0. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_atan >>> R, x, y = ring('x, y', QQ) >>> rs_atan(x + x*y, x, 4) -1/3*x**3*y**3 - x**3*y**2 - x**3*y - 1/3*x**3 + x*y + x See Also ======== atan """ if rs_is_puiseux(p, x): return rs_puiseux(rs_atan, p, x, prec) R = p.ring const = 0 if _has_constant_term(p, x): zm = R.zero_monom c = p[zm] if R.domain is EX: c_expr = c.as_expr() const = atan(c_expr) elif isinstance(c, PolyElement): try: c_expr = c.as_expr() const = R(atan(c_expr)) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") else: try: const = R(atan(c)) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") # Instead of using a closed form formula, we differentiate atan(p) to get # `1/(1+p**2) * dp`, whose series expansion is much easier to calculate. # Finally we integrate to get back atan dp = p.diff(x) p1 = rs_square(p, x, prec) + R(1) p1 = rs_series_inversion(p1, x, prec - 1) p1 = rs_mul(dp, p1, x, prec - 1) return rs_integrate(p1, x) + const
[docs]def rs_asin(p, x, prec): """ Arcsine of a series Return the series expansion of the asin of ``p``, about 0. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_asin >>> R, x, y = ring('x, y', QQ) >>> rs_asin(x, x, 8) 5/112*x**7 + 3/40*x**5 + 1/6*x**3 + x See Also ======== asin """ if rs_is_puiseux(p, x): return rs_puiseux(rs_asin, p, x, prec) if _has_constant_term(p, x): raise NotImplementedError("Polynomial must not have constant term in " "series variables") R = p.ring if x in R.gens: # get a good value if len(p) > 20: dp = rs_diff(p, x) p1 = 1 - rs_square(p, x, prec - 1) p1 = rs_nth_root(p1, -2, x, prec - 1) p1 = rs_mul(dp, p1, x, prec - 1) return rs_integrate(p1, x) one = R(1) c = [0, one, 0] for k in range(3, prec, 2): c.append((k - 2)**2*c[-2]/(k*(k - 1))) c.append(0) return rs_series_from_list(p, c, x, prec) else: raise NotImplementedError
def _tan1(p, x, prec): r""" Helper function of `rs\_tan`. Return the series expansion of tan of a univariate series using Newton's method. It takes advantage of the fact that series expansion of atan is easier than that of tan. Consider `f(x) = y - atan(x)` Let r be a root of f(x) found using Newton's method. Then `f(r) = 0` Or `y = atan(x)` where `x = tan(y)` as required. """ R = p.ring p1 = R(0) for precx in _giant_steps(prec): tmp = p - rs_atan(p1, x, precx) tmp = rs_mul(tmp, 1 + rs_square(p1, x, precx), x, precx) p1 += tmp return p1
[docs]def rs_tan(p, x, prec): """ Tangent of a series. Return the series expansion of the tan of ``p``, about 0. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_tan >>> R, x, y = ring('x, y', QQ) >>> rs_tan(x + x*y, x, 4) 1/3*x**3*y**3 + x**3*y**2 + x**3*y + 1/3*x**3 + x*y + x See Also ======== _tan1, tan """ if rs_is_puiseux(p, x): r = rs_puiseux(rs_tan, p, x, prec) return r R = p.ring const = 0 c = _get_constant_term(p, x) if c: if R.domain is EX: c_expr = c.as_expr() const = tan(c_expr) elif isinstance(c, PolyElement): try: c_expr = c.as_expr() const = R(tan(c_expr)) except ValueError: R = R.add_gens([tan(c_expr, )]) p = p.set_ring(R) x = x.set_ring(R) c = c.set_ring(R) const = R(tan(c_expr)) else: try: const = R(tan(c)) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") p1 = p - c # Makes use of sympy fuctions to evaluate the values of the cos/sin # of the constant term. t2 = rs_tan(p1, x, prec) t = rs_series_inversion(1 - const*t2, x, prec) return rs_mul(const + t2, t, x, prec) if R.ngens == 1: return _tan1(p, x, prec) else: return rs_fun(p, rs_tan, x, prec)
[docs]def rs_cot(p, x, prec): """ Cotangent of a series Return the series expansion of the cot of ``p``, about 0. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_cot >>> R, x, y = ring('x, y', QQ) >>> rs_cot(x, x, 6) -2/945*x**5 - 1/45*x**3 - 1/3*x + x**(-1) See Also ======== cot """ # It can not handle series like `p = x + x*y` where the coefficient of the # linear term in the series variable is symbolic. if rs_is_puiseux(p, x): r = rs_puiseux(rs_cot, p, x, prec) return r i, m = _check_series_var(p, x, 'cot') prec1 = prec + 2*m c, s = rs_cos_sin(p, x, prec1) s = mul_xin(s, i, -m) s = rs_series_inversion(s, x, prec1) res = rs_mul(c, s, x, prec1) res = mul_xin(res, i, -m) res = rs_trunc(res, x, prec) return res
[docs]def rs_sin(p, x, prec): """ Sine of a series Return the series expansion of the sin of ``p``, about 0. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_sin >>> R, x, y = ring('x, y', QQ) >>> rs_sin(x + x*y, x, 4) -1/6*x**3*y**3 - 1/2*x**3*y**2 - 1/2*x**3*y - 1/6*x**3 + x*y + x >>> rs_sin(x**QQ(3, 2) + x*y**QQ(7, 5), x, 4) -1/2*x**(7/2)*y**(14/5) - 1/6*x**3*y**(21/5) + x**(3/2) + x*y**(7/5) See Also ======== sin """ if rs_is_puiseux(p, x): return rs_puiseux(rs_sin, p, x, prec) R = x.ring if not p: return R(0) c = _get_constant_term(p, x) if c: if R.domain is EX: c_expr = c.as_expr() t1, t2 = sin(c_expr), cos(c_expr) elif isinstance(c, PolyElement): try: c_expr = c.as_expr() t1, t2 = R(sin(c_expr)), R(cos(c_expr)) except ValueError: R = R.add_gens([sin(c_expr), cos(c_expr)]) p = p.set_ring(R) x = x.set_ring(R) c = c.set_ring(R) t1, t2 = R(sin(c_expr)), R(cos(c_expr)) else: try: t1, t2 = R(sin(c)), R(cos(c)) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") p1 = p - c # Makes use of sympy cos, sin fuctions to evaluate the values of the # cos/sin of the constant term. return rs_sin(p1, x, prec)*t2 + rs_cos(p1, x, prec)*t1 # Series is calculated in terms of tan as its evaluation is fast. if len(p) > 20 and R.ngens == 1: t = rs_tan(p/2, x, prec) t2 = rs_square(t, x, prec) p1 = rs_series_inversion(1 + t2, x, prec) return rs_mul(p1, 2*t, x, prec) one = R(1) n = 1 c = [0] for k in range(2, prec + 2, 2): c.append(one/n) c.append(0) n *= -k*(k + 1) return rs_series_from_list(p, c, x, prec)
[docs]def rs_cos(p, x, prec): """ Cosine of a series Return the series expansion of the cos of ``p``, about 0. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_cos >>> R, x, y = ring('x, y', QQ) >>> rs_cos(x + x*y, x, 4) -1/2*x**2*y**2 - x**2*y - 1/2*x**2 + 1 >>> rs_cos(x + x*y, x, 4)/x**QQ(7, 5) -1/2*x**(3/5)*y**2 - x**(3/5)*y - 1/2*x**(3/5) + x**(-7/5) See Also ======== cos """ if rs_is_puiseux(p, x): return rs_puiseux(rs_cos, p, x, prec) R = p.ring c = _get_constant_term(p, x) if c: if R.domain is EX: c_expr = c.as_expr() t1, t2 = sin(c_expr), cos(c_expr) elif isinstance(c, PolyElement): try: c_expr = c.as_expr() t1, t2 = R(sin(c_expr)), R(cos(c_expr)) except ValueError: R = R.add_gens([sin(c_expr), cos(c_expr)]) p = p.set_ring(R) x = x.set_ring(R) c = c.set_ring(R) else: try: t1, t2 = R(sin(c)), R(cos(c)) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") p1 = p - c # Makes use of sympy cos, sin fuctions to evaluate the values of the # cos/sin of the constant term. p_cos = rs_cos(p1, x, prec) p_sin = rs_sin(p1, x, prec) R = R.compose(p_cos.ring).compose(p_sin.ring) p_cos.set_ring(R) p_sin.set_ring(R) t1, t2 = R(sin(c_expr)), R(cos(c_expr)) return p_cos*t2 - p_sin*t1 # Series is calculated in terms of tan as its evaluation is fast. if len(p) > 20 and R.ngens == 1: t = rs_tan(p/2, x, prec) t2 = rs_square(t, x, prec) p1 = rs_series_inversion(1+t2, x, prec) return rs_mul(p1, 1 - t2, x, prec) one = R(1) n = 1 c = [] for k in range(2, prec + 2, 2): c.append(one/n) c.append(0) n *= -k*(k - 1) return rs_series_from_list(p, c, x, prec)
[docs]def rs_cos_sin(p, x, prec): r""" Return the tuple `(rs\_cos(p, x, prec)`, `rs\_sin(p, x, prec))`. Is faster than calling rs_cos and rs_sin separately """ if rs_is_puiseux(p, x): return rs_puiseux(rs_cos_sin, p, x, prec) t = rs_tan(p/2, x, prec) t2 = rs_square(t, x, prec) p1 = rs_series_inversion(1 + t2, x, prec) return (rs_mul(p1, 1 - t2, x, prec), rs_mul(p1, 2*t, x, prec))
def _atanh(p, x, prec): """ Expansion using formula Faster for very small and univariate series """ R = p.ring one = R(1) c = [one] p2 = rs_square(p, x, prec) for k in range(1, prec): c.append(one/(2*k + 1)) s = rs_series_from_list(p2, c, x, prec) s = rs_mul(s, p, x, prec) return s
[docs]def rs_atanh(p, x, prec): """ Hyperbolic arctangent of a series Return the series expansion of the atanh of ``p``, about 0. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_atanh >>> R, x, y = ring('x, y', QQ) >>> rs_atanh(x + x*y, x, 4) 1/3*x**3*y**3 + x**3*y**2 + x**3*y + 1/3*x**3 + x*y + x See Also ======== atanh """ if rs_is_puiseux(p, x): return rs_puiseux(rs_atanh, p, x, prec) R = p.ring const = 0 if _has_constant_term(p, x): zm = R.zero_monom c = p[zm] if R.domain is EX: c_expr = c.as_expr() const = atanh(c_expr) elif isinstance(c, PolyElement): try: c_expr = c.as_expr() const = R(atanh(c_expr)) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") else: try: const = R(atanh(c)) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") # Instead of using a closed form formula, we differentiate atanh(p) to get # `1/(1-p**2) * dp`, whose series expansion is much easier to calculate. # Finally we integrate to get back atanh dp = rs_diff(p, x) p1 = - rs_square(p, x, prec) + 1 p1 = rs_series_inversion(p1, x, prec - 1) p1 = rs_mul(dp, p1, x, prec - 1) return rs_integrate(p1, x) + const
[docs]def rs_sinh(p, x, prec): """ Hyperbolic sine of a series Return the series expansion of the sinh of ``p``, about 0. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_sinh >>> R, x, y = ring('x, y', QQ) >>> rs_sinh(x + x*y, x, 4) 1/6*x**3*y**3 + 1/2*x**3*y**2 + 1/2*x**3*y + 1/6*x**3 + x*y + x See Also ======== sinh """ if rs_is_puiseux(p, x): return rs_puiseux(rs_sinh, p, x, prec) t = rs_exp(p, x, prec) t1 = rs_series_inversion(t, x, prec) return (t - t1)/2
[docs]def rs_cosh(p, x, prec): """ Hyperbolic cosine of a series Return the series expansion of the cosh of ``p``, about 0. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_cosh >>> R, x, y = ring('x, y', QQ) >>> rs_cosh(x + x*y, x, 4) 1/2*x**2*y**2 + x**2*y + 1/2*x**2 + 1 See Also ======== cosh """ if rs_is_puiseux(p, x): return rs_puiseux(rs_cosh, p, x, prec) t = rs_exp(p, x, prec) t1 = rs_series_inversion(t, x, prec) return (t + t1)/2
def _tanh(p, x, prec): r""" Helper function of `rs\_tanh` Return the series expansion of tanh of a univariate series using Newton's method. It takes advantage of the fact that series expansion of atanh is easier than that of tanh. See Also ======== _tanh """ R = p.ring p1 = R(0) for precx in _giant_steps(prec): tmp = p - rs_atanh(p1, x, precx) tmp = rs_mul(tmp, 1 - rs_square(p1, x, prec), x, precx) p1 += tmp return p1
[docs]def rs_tanh(p, x, prec): """ Hyperbolic tangent of a series Return the series expansion of the tanh of ``p``, about 0. Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_tanh >>> R, x, y = ring('x, y', QQ) >>> rs_tanh(x + x*y, x, 4) -1/3*x**3*y**3 - x**3*y**2 - x**3*y - 1/3*x**3 + x*y + x See Also ======== tanh """ if rs_is_puiseux(p, x): return rs_puiseux(rs_tanh, p, x, prec) R = p.ring const = 0 if _has_constant_term(p, x): zm = R.zero_monom c = p[zm] if R.domain is EX: c_expr = c.as_expr() const = tanh(c_expr) elif isinstance(c, PolyElement): try: c_expr = c.as_expr() const = R(tanh(c_expr)) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") else: try: const = R(tanh(c)) except ValueError: raise DomainError("The given series can't be expanded in " "this domain.") p1 = p - c t1 = rs_tanh(p1, x, prec) t = rs_series_inversion(1 + const*t1, x, prec) return rs_mul(const + t1, t, x, prec) if R.ngens == 1: return _tanh(p, x, prec) else: return rs_fun(p, _tanh, x, prec)
[docs]def rs_newton(p, x, prec): """ Compute the truncated Newton sum of the polynomial ``p`` Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_newton >>> R, x = ring('x', QQ) >>> p = x**2 - 2 >>> rs_newton(p, x, 5) 8*x**4 + 4*x**2 + 2 """ deg = p.degree() p1 = _invert_monoms(p) p2 = rs_series_inversion(p1, x, prec) p3 = rs_mul(p1.diff(x), p2, x, prec) res = deg - p3*x return res
[docs]def rs_hadamard_exp(p1, inverse=False): """ Return ``sum f_i/i!*x**i`` from ``sum f_i*x**i``, where ``x`` is the first variable. If ``invers=True`` return ``sum f_i*i!*x**i`` Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_hadamard_exp >>> R, x = ring('x', QQ) >>> p = 1 + x + x**2 + x**3 >>> rs_hadamard_exp(p) 1/6*x**3 + 1/2*x**2 + x + 1 """ R = p1.ring if R.domain != QQ: raise NotImplementedError p = R.zero if not inverse: for exp1, v1 in p1.items(): p[exp1] = v1/int(ifac(exp1[0])) else: for exp1, v1 in p1.items(): p[exp1] = v1*int(ifac(exp1[0])) return p
[docs]def rs_compose_add(p1, p2): """ compute the composed sum ``prod(p2(x - beta) for beta root of p1)`` Examples ======== >>> from sympy.polys.domains import QQ >>> from sympy.polys.rings import ring >>> from sympy.polys.ring_series import rs_compose_add >>> R, x = ring('x', QQ) >>> f = x**2 - 2 >>> g = x**2 - 3 >>> rs_compose_add(f, g) x**4 - 10*x**2 + 1 References ========== A. Bostan, P. Flajolet, B. Salvy and E. Schost "Fast Computation with Two Algebraic Numbers", (2002) Research Report 4579, Institut National de Recherche en Informatique et en Automatique """ R = p1.ring x = R.gens[0] prec = p1.degree() * p2.degree() + 1 np1 = rs_newton(p1, x, prec) np1e = rs_hadamard_exp(np1) np2 = rs_newton(p2, x, prec) np2e = rs_hadamard_exp(np2) np3e = rs_mul(np1e, np2e, x, prec) np3 = rs_hadamard_exp(np3e, True) np3a = (np3[(0,)] - np3)/x q = rs_integrate(np3a, x) q = rs_exp(q, x, prec) q = _invert_monoms(q) q = q.primitive()[1] dp = p1.degree() * p2.degree() - q.degree() # `dp` is the multiplicity of the zeroes of the resultant; # these zeroes are missed in this computation so they are put here. # if p1 and p2 are monic irreducible polynomials, # there are zeroes in the resultant # if and only if p1 = p2 ; in fact in that case p1 and p2 have a # root in common, so gcd(p1, p2) != 1; being p1 and p2 irreducible # this means p1 = p2 if dp: q = q*x**dp return q
_convert_func = { 'sin': 'rs_sin', 'cos': 'rs_cos', 'exp': 'rs_exp', 'tan': 'rs_tan', 'log': 'rs_log' } def rs_min_pow(expr, series_rs, a): """Find the minimum power of `a` in the series expansion of expr""" series = 0 n = 2 while series == 0: series = _rs_series(expr, series_rs, a, n) n *= 2 R = series.ring a = R(a) i = R.gens.index(a) return min(series, key=lambda t: t[i])[i] def _rs_series(expr, series_rs, a, prec): # TODO Use _parallel_dict_from_expr instead of sring as sring is # inefficient. For details, read the todo in sring. args = expr.args R = series_rs.ring # expr does not contain any function to be expanded if not any(arg.has(Function) for arg in args) and not expr.is_Function: return series_rs if not expr.has(a): return series_rs elif expr.is_Function: arg = args[0] if len(args) > 1: raise NotImplementedError R1, series = sring(arg, domain=QQ, expand=False, series=True) series_inner = _rs_series(arg, series, a, prec) # Why do we need to compose these three rings? # # We want to use a simple domain (like ``QQ`` or ``RR``) but they don't # support symbolic coefficients. We need a ring that for example lets # us have `sin(1)` and `cos(1)` as coefficients if we are expanding # `sin(x + 1)`. The ``EX`` domain allows all symbolic coefficients, but # that makes it very complex and hence slow. # # To solve this problem, we add only those symbolic elements as # generators to our ring, that we need. Here, series_inner might # involve terms like `sin(4)`, `exp(a)`, etc, which are not there in # R1 or R. Hence, we compose these three rings to create one that has # the generators of all three. R = R.compose(R1).compose(series_inner.ring) series_inner = series_inner.set_ring(R) series = eval(_convert_func[str(expr.func)])(series_inner, R(a), prec) return series elif expr.is_Mul: n = len(args) for arg in args: # XXX Looks redundant if not arg.is_Number: R1, _ = sring(arg, expand=False, series=True) R = R.compose(R1) min_pows = list(map(rs_min_pow, args, [R(arg) for arg in args], [a]*len(args))) sum_pows = sum(min_pows) series = R(1) for i in range(n): _series = _rs_series(args[i], R(args[i]), a, prec - sum_pows + min_pows[i]) R = R.compose(_series.ring) _series = _series.set_ring(R) series = series.set_ring(R) series *= _series series = rs_trunc(series, R(a), prec) return series elif expr.is_Add: n = len(args) series = R(0) for i in range(n): _series = _rs_series(args[i], R(args[i]), a, prec) R = R.compose(_series.ring) _series = _series.set_ring(R) series = series.set_ring(R) series += _series return series elif expr.is_Pow: R1, _ = sring(expr.base, domain=QQ, expand=False, series=True) R = R.compose(R1) series_inner = _rs_series(expr.base, R(expr.base), a, prec) return rs_pow(series_inner, expr.exp, series_inner.ring(a), prec) # The `is_constant` method is buggy hence we check it at the end. # See issue #9786 for details. elif isinstance(expr, Expr) and expr.is_constant(): return sring(expr, domain=QQ, expand=False, series=True)[1] else: raise NotImplementedError def rs_series(expr, a, prec): """Return the series expansion of an expression about 0. Parameters ---------- expr : :class:`Expr` a : :class:`Symbol` with respect to which expr is to be expanded prec : order of the series expansion Currently supports multivariate Taylor series expansion. This is much faster that Sympy's series method as it uses sparse polynomial operations. It automatically creates the simplest ring required to represent the series expansion through repeated calls to sring. Examples ======== >>> from sympy.polys.ring_series import rs_series >>> from sympy.functions import sin, cos, exp, tan >>> from sympy.core import symbols >>> from sympy.polys.domains import QQ >>> a, b, c = symbols('a, b, c') >>> rs_series(sin(a) + exp(a), a, 5) 1/24*a**4 + 1/2*a**2 + 2*a + 1 >>> series = rs_series(tan(a + b)*cos(a + c), a, 2) >>> series.as_expr() -a*sin(c)*tan(b) + a*cos(c)*tan(b)**2 + a*cos(c) + cos(c)*tan(b) >>> series = rs_series(exp(a**QQ(1,3) + a**QQ(2, 5)), a, 1) >>> series.as_expr() a**(11/15) + a**(4/5)/2 + a**(2/5) + a**(2/3)/2 + a**(1/3) + 1 """ R, series = sring(expr, domain=QQ, expand=False, series=True) if a not in R.symbols: R = R.add_gens([a, ]) series = series.set_ring(R) series = _rs_series(expr, series, a, prec) R = series.ring gen = R(a) prec_got = series.degree(gen) + 1 if prec_got >= prec: return rs_trunc(series, gen, prec) else: # increase the requested number of terms to get the desired # number keep increasing (up to 9) until the received order # is different than the original order and then predict how # many additional terms are needed for more in range(1, 9): p1 = _rs_series(expr, series, a, prec=prec + more) gen = gen.set_ring(p1.ring) new_prec = p1.degree(gen) + 1 if new_prec != prec_got: prec_do = ceiling(prec + (prec - prec_got)*more/(new_prec - prec_got)) p1 = _rs_series(expr, series, a, prec=prec_do) while p1.degree(gen) + 1 < prec: p1 = _rs_series(expr, series, a, prec=prec_do) gen = gen.set_ring(p1.ring) prec_do *= 2 break else: break else: raise ValueError('Could not calculate %s terms for %s' % (str(prec), expr)) return rs_trunc(p1, gen, prec)