Source code for sympy.polys.ring_series

"""Power series manipulating functions acting on polys.ring.PolyElement()"""

from sympy.polys.domains import QQ
from sympy.polys.rings import ring, PolyElement
from sympy.polys.monomials import monomial_min, monomial_mul
from sympy.mpmath.libmp.libintmath import ifac
from sympy.core.numbers import Rational
from sympy.core.compatibility import as_int
from sympy.mpmath.libmp.libintmath import giant_steps
import math

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

    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()
    ring = p1.ring
    p = ring.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):
    """
    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 """ ring = p1.ring p = ring.zero i = ring.gens.index(x) for exp1 in p1: if exp1[i] >= prec: continue p[exp1] = p1[exp1] return p
[docs]def rs_mul(p1, p2, x, prec): """ product of 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 """ ring = p1.ring p = ring.zero if ring.__class__ != p2.ring.__class__ or ring != p2.ring: raise ValueError('p1 and p2 must have the same ring') iv = ring.gens.index(x) if not isinstance(p2, PolyElement): raise ValueError('p1 and p2 must have the same ring') if ring == p2.ring: get = p.get items2 = list(p2.items()) items2.sort(key=lambda e: e[0][iv]) if ring.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 = ring.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 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 """ ring = p1.ring p = ring.zero iv = ring.gens.index(x) get = p.get items = list(p1.items()) items.sort(key=lambda e: e[0][iv]) monomial_mul = ring.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): raise NotImplementedError('to be implemented') 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 _has_constant_term(p, x): """ test 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 """ ring = p.ring iv = ring.gens.index(x) zm = ring.zero_monom a = [0]*ring.ngens a[iv] = 1 miv = tuple(a) for expv in p: if monomial_min(expv, miv) == zm: return True return False 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 """ ring = p.ring zm = ring.zero_monom if zm not in p: raise ValueError('no constant term in series') if _has_constant_term(p - p[zm], x): raise ValueError('p cannot contain a constant term depending on parameters') if p[zm] != ring(1): # TODO add check that it is a unit p1 = ring(1)/p[zm] else: p1 = ring(1) for precx in _giant_steps(prec): tmp = p1.square() tmp = rs_mul(tmp, p, x, precx) p1 = 2*p1 - tmp 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 """ ring = p.ring zm = ring.zero_monom ii = ring.gens.index(x) m = min(p, key=lambda k: k[ii])[ii] if m: raise NotImplementedError('no constant term in series') 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') return _series_inversion1(p, x, prec)
[docs]def rs_series_from_list(p, c, x, prec, concur=1): """ series ``sum c[n]*p**n`` modulo ``O(x**prec)`` reduce the number of multiplication 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 """ ring = p.ring n = len(c) if not concur: q = ring(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 = [ring(1)] b = 1 q = ring(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 = ring(1) s = ring(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]*ring(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_integrate(self, x): """ integrate ``p`` with respect to ``x`` 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 """ ring = self.ring p1 = ring.zero n = ring.gens.index(x) mn = [0]*ring.ngens mn[n] = 1 mn = tuple(mn) for expv in self: e = monomial_mul(expv, mn) p1[e] = self[expv]/(expv[n] + 1) return p1
[docs]def rs_log(p, x, prec): """ 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 """ ring = p.ring if p == 1: return 0 if _has_constant_term(p - 1, x): raise NotImplementedError('p - 1 must not have a constant term in the series variables') dlog = p.diff(x) dlog = rs_mul(dlog, _series_inversion1(p, x, prec), x, prec - 1) return rs_integrate(dlog, x)
def _exp1(p, x, prec): """ helper function for ``rs_exp`` """ ring = p.ring p1 = ring(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 """ ring = p.ring if _has_constant_term(p, x): raise NotImplementedError if len(p) > 20: return _exp1(p, x, prec) one = ring(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
[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