"""Recurrence Operators""" from sympy.core.singleton import S from sympy.core.symbol import (Symbol, symbols) from sympy.printing import sstr from sympy.core.sympify import sympify def RecurrenceOperators(base, generator): """ Returns an Algebra of Recurrence Operators and the operator for shifting i.e. the `Sn` operator. The first argument needs to be the base polynomial ring for the algebra and the second argument must be a generator which can be either a noncommutative Symbol or a string. Examples ======== >>> from sympy import ZZ >>> from sympy import symbols >>> from sympy.holonomic.recurrence import RecurrenceOperators >>> n = symbols('n', integer=True) >>> R, Sn = RecurrenceOperators(ZZ.old_poly_ring(n), 'Sn') """ ring = RecurrenceOperatorAlgebra(base, generator) return (ring, ring.shift_operator) class RecurrenceOperatorAlgebra: """ A Recurrence Operator Algebra is a set of noncommutative polynomials in intermediate `Sn` and coefficients in a base ring A. It follows the commutation rule: Sn * a(n) = a(n + 1) * Sn This class represents a Recurrence Operator Algebra and serves as the parent ring for Recurrence Operators. Examples ======== >>> from sympy import ZZ >>> from sympy import symbols >>> from sympy.holonomic.recurrence import RecurrenceOperators >>> n = symbols('n', integer=True) >>> R, Sn = RecurrenceOperators(ZZ.old_poly_ring(n), 'Sn') >>> R Univariate Recurrence Operator Algebra in intermediate Sn over the base ring ZZ[n] See Also ======== RecurrenceOperator """ def __init__(self, base, generator): # the base ring for the algebra self.base = base # the operator representing shift i.e. `Sn` self.shift_operator = RecurrenceOperator( [base.zero, base.one], self) if generator is None: self.gen_symbol = symbols('Sn', commutative=False) else: if isinstance(generator, str): self.gen_symbol = symbols(generator, commutative=False) elif isinstance(generator, Symbol): self.gen_symbol = generator def __str__(self): string = 'Univariate Recurrence Operator Algebra in intermediate '\ + sstr(self.gen_symbol) + ' over the base ring ' + \ (self.base).__str__() return string __repr__ = __str__ def __eq__(self, other): if self.base == other.base and self.gen_symbol == other.gen_symbol: return True else: return False def _add_lists(list1, list2): if len(list1) <= len(list2): sol = [a + b for a, b in zip(list1, list2)] + list2[len(list1):] else: sol = [a + b for a, b in zip(list1, list2)] + list1[len(list2):] return sol class RecurrenceOperator: """ The Recurrence Operators are defined by a list of polynomials in the base ring and the parent ring of the Operator. Explanation =========== Takes a list of polynomials for each power of Sn and the parent ring which must be an instance of RecurrenceOperatorAlgebra. A Recurrence Operator can be created easily using the operator `Sn`. See examples below. Examples ======== >>> from sympy.holonomic.recurrence import RecurrenceOperator, RecurrenceOperators >>> from sympy import ZZ >>> from sympy import symbols >>> n = symbols('n', integer=True) >>> R, Sn = RecurrenceOperators(ZZ.old_poly_ring(n),'Sn') >>> RecurrenceOperator([0, 1, n**2], R) (1)Sn + (n**2)Sn**2 >>> Sn*n (n + 1)Sn >>> n*Sn*n + 1 - Sn**2*n (1) + (n**2 + n)Sn + (-n - 2)Sn**2 See Also ======== DifferentialOperatorAlgebra """ _op_priority = 20 def __init__(self, list_of_poly, parent): # the parent ring for this operator # must be an RecurrenceOperatorAlgebra object self.parent = parent # sequence of polynomials in n for each power of Sn # represents the operator # convert the expressions into ring elements using from_sympy if isinstance(list_of_poly, list): for i, j in enumerate(list_of_poly): if isinstance(j, int): list_of_poly[i] = self.parent.base.from_sympy(S(j)) elif not isinstance(j, self.parent.base.dtype): list_of_poly[i] = self.parent.base.from_sympy(j) self.listofpoly = list_of_poly self.order = len(self.listofpoly) - 1 def __mul__(self, other): """ Multiplies two Operators and returns another RecurrenceOperator instance using the commutation rule Sn * a(n) = a(n + 1) * Sn """ listofself = self.listofpoly base = self.parent.base if not isinstance(other, RecurrenceOperator): if not isinstance(other, self.parent.base.dtype): listofother = [self.parent.base.from_sympy(sympify(other))] else: listofother = [other] else: listofother = other.listofpoly # multiply a polynomial `b` with a list of polynomials def _mul_dmp_diffop(b, listofother): if isinstance(listofother, list): sol = [] for i in listofother: sol.append(i * b) return sol else: return [b * listofother] sol = _mul_dmp_diffop(listofself[0], listofother) # compute Sn^i * b def _mul_Sni_b(b): sol = [base.zero] if isinstance(b, list): for i in b: j = base.to_sympy(i).subs(base.gens[0], base.gens[0] + S.One) sol.append(base.from_sympy(j)) else: j = b.subs(base.gens[0], base.gens[0] + S.One) sol.append(base.from_sympy(j)) return sol for i in range(1, len(listofself)): # find Sn^i * b in ith iteration listofother = _mul_Sni_b(listofother) # solution = solution + listofself[i] * (Sn^i * b) sol = _add_lists(sol, _mul_dmp_diffop(listofself[i], listofother)) return RecurrenceOperator(sol, self.parent) def __rmul__(self, other): if not isinstance(other, RecurrenceOperator): if isinstance(other, int): other = S(other) if not isinstance(other, self.parent.base.dtype): other = (self.parent.base).from_sympy(other) sol = [] for j in self.listofpoly: sol.append(other * j) return RecurrenceOperator(sol, self.parent) def __add__(self, other): if isinstance(other, RecurrenceOperator): sol = _add_lists(self.listofpoly, other.listofpoly) return RecurrenceOperator(sol, self.parent) else: if isinstance(other, int): other = S(other) list_self = self.listofpoly if not isinstance(other, self.parent.base.dtype): list_other = [((self.parent).base).from_sympy(other)] else: list_other = [other] sol = [] sol.append(list_self[0] + list_other[0]) sol += list_self[1:] return RecurrenceOperator(sol, self.parent) __radd__ = __add__ def __sub__(self, other): return self + (-1) * other def __rsub__(self, other): return (-1) * self + other def __pow__(self, n): if n == 1: return self result = RecurrenceOperator([self.parent.base.one], self.parent) if n == 0: return result # if self is `Sn` if self.listofpoly == self.parent.shift_operator.listofpoly: sol = [self.parent.base.zero] * n + [self.parent.base.one] return RecurrenceOperator(sol, self.parent) x = self while True: if n % 2: result *= x n >>= 1 if not n: break x *= x return result def __str__(self): listofpoly = self.listofpoly print_str = '' for i, j in enumerate(listofpoly): if j == self.parent.base.zero: continue j = self.parent.base.to_sympy(j) if i == 0: print_str += '(' + sstr(j) + ')' continue if print_str: print_str += ' + ' if i == 1: print_str += '(' + sstr(j) + ')Sn' continue print_str += '(' + sstr(j) + ')' + 'Sn**' + sstr(i) return print_str __repr__ = __str__ def __eq__(self, other): if isinstance(other, RecurrenceOperator): if self.listofpoly == other.listofpoly and self.parent == other.parent: return True else: return False else: if self.listofpoly[0] == other: for i in self.listofpoly[1:]: if i is not self.parent.base.zero: return False return True else: return False class HolonomicSequence: """ A Holonomic Sequence is a type of sequence satisfying a linear homogeneous recurrence relation with Polynomial coefficients. Alternatively, A sequence is Holonomic if and only if its generating function is a Holonomic Function. """ def __init__(self, recurrence, u0=[]): self.recurrence = recurrence if not isinstance(u0, list): self.u0 = [u0] else: self.u0 = u0 if len(self.u0) == 0: self._have_init_cond = False else: self._have_init_cond = True self.n = recurrence.parent.base.gens[0] def __repr__(self): str_sol = 'HolonomicSequence(%s, %s)' % ((self.recurrence).__repr__(), sstr(self.n)) if not self._have_init_cond: return str_sol else: cond_str = '' seq_str = 0 for i in self.u0: cond_str += ', u(%s) = %s' % (sstr(seq_str), sstr(i)) seq_str += 1 sol = str_sol + cond_str return sol __str__ = __repr__ def __eq__(self, other): if self.recurrence == other.recurrence: if self.n == other.n: if self._have_init_cond and other._have_init_cond: if self.u0 == other.u0: return True else: return False else: return True else: return False else: return False