Coverage for sympy/polys/ring_series.py : 21%

"""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 

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 

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 

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 

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 

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 

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 

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 

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 

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): 

"""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 

def rs_series_reversion(p, x, n, y): 

""" 

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