1from sympy.core.basic import Basic 2from sympy.core.cache import cacheit 3from sympy.core.compatibility import is_sequence, iterable, ordered 4from sympy.core.containers import Tuple 5from sympy.core.decorators import call_highest_priority 6from sympy.core.parameters import global_parameters 7from sympy.core.function import AppliedUndef 8from sympy.core.mul import Mul 9from sympy.core.numbers import Integer 10from sympy.core.relational import Eq 11from sympy.core.singleton import S, Singleton 12from sympy.core.symbol import Dummy, Symbol, Wild 13from sympy.core.sympify import sympify 14from sympy.polys import lcm, factor 15from sympy.sets.sets import Interval, Intersection 16from sympy.simplify import simplify 17from sympy.tensor.indexed import Idx 18from sympy.utilities.iterables import flatten 19from sympy import expand 20 21 22############################################################################### 23# SEQUENCES # 24############################################################################### 25 26 27class SeqBase(Basic): 28 """Base class for sequences""" 29 30 is_commutative = True 31 _op_priority = 15 32 33 @staticmethod 34 def _start_key(expr): 35 """Return start (if possible) else S.Infinity. 36 37 adapted from Set._infimum_key 38 """ 39 try: 40 start = expr.start 41 except (NotImplementedError, 42 AttributeError, ValueError): 43 start = S.Infinity 44 return start 45 46 def _intersect_interval(self, other): 47 """Returns start and stop. 48 49 Takes intersection over the two intervals. 50 """ 51 interval = Intersection(self.interval, other.interval) 52 return interval.inf, interval.sup 53 54 @property 55 def gen(self): 56 """Returns the generator for the sequence""" 57 raise NotImplementedError("(%s).gen" % self) 58 59 @property 60 def interval(self): 61 """The interval on which the sequence is defined""" 62 raise NotImplementedError("(%s).interval" % self) 63 64 @property 65 def start(self): 66 """The starting point of the sequence. This point is included""" 67 raise NotImplementedError("(%s).start" % self) 68 69 @property 70 def stop(self): 71 """The ending point of the sequence. This point is included""" 72 raise NotImplementedError("(%s).stop" % self) 73 74 @property 75 def length(self): 76 """Length of the sequence""" 77 raise NotImplementedError("(%s).length" % self) 78 79 @property 80 def variables(self): 81 """Returns a tuple of variables that are bounded""" 82 return () 83 84 @property 85 def free_symbols(self): 86 """ 87 This method returns the symbols in the object, excluding those 88 that take on a specific value (i.e. the dummy symbols). 89 90 Examples 91 ======== 92 93 >>> from sympy import SeqFormula 94 >>> from sympy.abc import n, m 95 >>> SeqFormula(m*n**2, (n, 0, 5)).free_symbols 96 {m} 97 """ 98 return ({j for i in self.args for j in i.free_symbols 99 .difference(self.variables)}) 100 101 @cacheit 102 def coeff(self, pt): 103 """Returns the coefficient at point pt""" 104 if pt < self.start or pt > self.stop: 105 raise IndexError("Index %s out of bounds %s" % (pt, self.interval)) 106 return self._eval_coeff(pt) 107 108 def _eval_coeff(self, pt): 109 raise NotImplementedError("The _eval_coeff method should be added to" 110 "%s to return coefficient so it is available" 111 "when coeff calls it." 112 % self.func) 113 114 def _ith_point(self, i): 115 """Returns the i'th point of a sequence. 116 117 Explanation 118 =========== 119 120 If start point is negative infinity, point is returned from the end. 121 Assumes the first point to be indexed zero. 122 123 Examples 124 ========= 125 126 >>> from sympy import oo 127 >>> from sympy.series.sequences import SeqPer 128 129 bounded 130 131 >>> SeqPer((1, 2, 3), (-10, 10))._ith_point(0) 132 -10 133 >>> SeqPer((1, 2, 3), (-10, 10))._ith_point(5) 134 -5 135 136 End is at infinity 137 138 >>> SeqPer((1, 2, 3), (0, oo))._ith_point(5) 139 5 140 141 Starts at negative infinity 142 143 >>> SeqPer((1, 2, 3), (-oo, 0))._ith_point(5) 144 -5 145 """ 146 if self.start is S.NegativeInfinity: 147 initial = self.stop 148 else: 149 initial = self.start 150 151 if self.start is S.NegativeInfinity: 152 step = -1 153 else: 154 step = 1 155 156 return initial + i*step 157 158 def _add(self, other): 159 """ 160 Should only be used internally. 161 162 Explanation 163 =========== 164 165 self._add(other) returns a new, term-wise added sequence if self 166 knows how to add with other, otherwise it returns ``None``. 167 168 ``other`` should only be a sequence object. 169 170 Used within :class:`SeqAdd` class. 171 """ 172 return None 173 174 def _mul(self, other): 175 """ 176 Should only be used internally. 177 178 Explanation 179 =========== 180 181 self._mul(other) returns a new, term-wise multiplied sequence if self 182 knows how to multiply with other, otherwise it returns ``None``. 183 184 ``other`` should only be a sequence object. 185 186 Used within :class:`SeqMul` class. 187 """ 188 return None 189 190 def coeff_mul(self, other): 191 """ 192 Should be used when ``other`` is not a sequence. Should be 193 defined to define custom behaviour. 194 195 Examples 196 ======== 197 198 >>> from sympy import SeqFormula 199 >>> from sympy.abc import n 200 >>> SeqFormula(n**2).coeff_mul(2) 201 SeqFormula(2*n**2, (n, 0, oo)) 202 203 Notes 204 ===== 205 206 '*' defines multiplication of sequences with sequences only. 207 """ 208 return Mul(self, other) 209 210 def __add__(self, other): 211 """Returns the term-wise addition of 'self' and 'other'. 212 213 ``other`` should be a sequence. 214 215 Examples 216 ======== 217 218 >>> from sympy import SeqFormula 219 >>> from sympy.abc import n 220 >>> SeqFormula(n**2) + SeqFormula(n**3) 221 SeqFormula(n**3 + n**2, (n, 0, oo)) 222 """ 223 if not isinstance(other, SeqBase): 224 raise TypeError('cannot add sequence and %s' % type(other)) 225 return SeqAdd(self, other) 226 227 @call_highest_priority('__add__') 228 def __radd__(self, other): 229 return self + other 230 231 def __sub__(self, other): 232 """Returns the term-wise subtraction of ``self`` and ``other``. 233 234 ``other`` should be a sequence. 235 236 Examples 237 ======== 238 239 >>> from sympy import SeqFormula 240 >>> from sympy.abc import n 241 >>> SeqFormula(n**2) - (SeqFormula(n)) 242 SeqFormula(n**2 - n, (n, 0, oo)) 243 """ 244 if not isinstance(other, SeqBase): 245 raise TypeError('cannot subtract sequence and %s' % type(other)) 246 return SeqAdd(self, -other) 247 248 @call_highest_priority('__sub__') 249 def __rsub__(self, other): 250 return (-self) + other 251 252 def __neg__(self): 253 """Negates the sequence. 254 255 Examples 256 ======== 257 258 >>> from sympy import SeqFormula 259 >>> from sympy.abc import n 260 >>> -SeqFormula(n**2) 261 SeqFormula(-n**2, (n, 0, oo)) 262 """ 263 return self.coeff_mul(-1) 264 265 def __mul__(self, other): 266 """Returns the term-wise multiplication of 'self' and 'other'. 267 268 ``other`` should be a sequence. For ``other`` not being a 269 sequence see :func:`coeff_mul` method. 270 271 Examples 272 ======== 273 274 >>> from sympy import SeqFormula 275 >>> from sympy.abc import n 276 >>> SeqFormula(n**2) * (SeqFormula(n)) 277 SeqFormula(n**3, (n, 0, oo)) 278 """ 279 if not isinstance(other, SeqBase): 280 raise TypeError('cannot multiply sequence and %s' % type(other)) 281 return SeqMul(self, other) 282 283 @call_highest_priority('__mul__') 284 def __rmul__(self, other): 285 return self * other 286 287 def __iter__(self): 288 for i in range(self.length): 289 pt = self._ith_point(i) 290 yield self.coeff(pt) 291 292 def __getitem__(self, index): 293 if isinstance(index, int): 294 index = self._ith_point(index) 295 return self.coeff(index) 296 elif isinstance(index, slice): 297 start, stop = index.start, index.stop 298 if start is None: 299 start = 0 300 if stop is None: 301 stop = self.length 302 return [self.coeff(self._ith_point(i)) for i in 303 range(start, stop, index.step or 1)] 304 305 def find_linear_recurrence(self,n,d=None,gfvar=None): 306 r""" 307 Finds the shortest linear recurrence that satisfies the first n 308 terms of sequence of order `\leq` ``n/2`` if possible. 309 If ``d`` is specified, find shortest linear recurrence of order 310 `\leq` min(d, n/2) if possible. 311 Returns list of coefficients ``[b(1), b(2), ...]`` corresponding to the 312 recurrence relation ``x(n) = b(1)*x(n-1) + b(2)*x(n-2) + ...`` 313 Returns ``[]`` if no recurrence is found. 314 If gfvar is specified, also returns ordinary generating function as a 315 function of gfvar. 316 317 Examples 318 ======== 319 320 >>> from sympy import sequence, sqrt, oo, lucas 321 >>> from sympy.abc import n, x, y 322 >>> sequence(n**2).find_linear_recurrence(10, 2) 323 [] 324 >>> sequence(n**2).find_linear_recurrence(10) 325 [3, -3, 1] 326 >>> sequence(2**n).find_linear_recurrence(10) 327 [2] 328 >>> sequence(23*n**4+91*n**2).find_linear_recurrence(10) 329 [5, -10, 10, -5, 1] 330 >>> sequence(sqrt(5)*(((1 + sqrt(5))/2)**n - (-(1 + sqrt(5))/2)**(-n))/5).find_linear_recurrence(10) 331 [1, 1] 332 >>> sequence(x+y*(-2)**(-n), (n, 0, oo)).find_linear_recurrence(30) 333 [1/2, 1/2] 334 >>> sequence(3*5**n + 12).find_linear_recurrence(20,gfvar=x) 335 ([6, -5], 3*(5 - 21*x)/((x - 1)*(5*x - 1))) 336 >>> sequence(lucas(n)).find_linear_recurrence(15,gfvar=x) 337 ([1, 1], (x - 2)/(x**2 + x - 1)) 338 """ 339 from sympy.matrices import Matrix 340 x = [simplify(expand(t)) for t in self[:n]] 341 lx = len(x) 342 if d is None: 343 r = lx//2 344 else: 345 r = min(d,lx//2) 346 coeffs = [] 347 for l in range(1, r+1): 348 l2 = 2*l 349 mlist = [] 350 for k in range(l): 351 mlist.append(x[k:k+l]) 352 m = Matrix(mlist) 353 if m.det() != 0: 354 y = simplify(m.LUsolve(Matrix(x[l:l2]))) 355 if lx == l2: 356 coeffs = flatten(y[::-1]) 357 break 358 mlist = [] 359 for k in range(l,lx-l): 360 mlist.append(x[k:k+l]) 361 m = Matrix(mlist) 362 if m*y == Matrix(x[l2:]): 363 coeffs = flatten(y[::-1]) 364 break 365 if gfvar is None: 366 return coeffs 367 else: 368 l = len(coeffs) 369 if l == 0: 370 return [], None 371 else: 372 n, d = x[l-1]*gfvar**(l-1), 1 - coeffs[l-1]*gfvar**l 373 for i in range(l-1): 374 n += x[i]*gfvar**i 375 for j in range(l-i-1): 376 n -= coeffs[i]*x[j]*gfvar**(i+j+1) 377 d -= coeffs[i]*gfvar**(i+1) 378 return coeffs, simplify(factor(n)/factor(d)) 379 380class EmptySequence(SeqBase, metaclass=Singleton): 381 """Represents an empty sequence. 382 383 The empty sequence is also available as a singleton as 384 ``S.EmptySequence``. 385 386 Examples 387 ======== 388 389 >>> from sympy import EmptySequence, SeqPer 390 >>> from sympy.abc import x 391 >>> EmptySequence 392 EmptySequence 393 >>> SeqPer((1, 2), (x, 0, 10)) + EmptySequence 394 SeqPer((1, 2), (x, 0, 10)) 395 >>> SeqPer((1, 2)) * EmptySequence 396 EmptySequence 397 >>> EmptySequence.coeff_mul(-1) 398 EmptySequence 399 """ 400 401 @property 402 def interval(self): 403 return S.EmptySet 404 405 @property 406 def length(self): 407 return S.Zero 408 409 def coeff_mul(self, coeff): 410 """See docstring of SeqBase.coeff_mul""" 411 return self 412 413 def __iter__(self): 414 return iter([]) 415 416 417class SeqExpr(SeqBase): 418 """Sequence expression class. 419 420 Various sequences should inherit from this class. 421 422 Examples 423 ======== 424 425 >>> from sympy.series.sequences import SeqExpr 426 >>> from sympy.abc import x 427 >>> s = SeqExpr((1, 2, 3), (x, 0, 10)) 428 >>> s.gen 429 (1, 2, 3) 430 >>> s.interval 431 Interval(0, 10) 432 >>> s.length 433 11 434 435 See Also 436 ======== 437 438 sympy.series.sequences.SeqPer 439 sympy.series.sequences.SeqFormula 440 """ 441 442 @property 443 def gen(self): 444 return self.args[0] 445 446 @property 447 def interval(self): 448 return Interval(self.args[1][1], self.args[1][2]) 449 450 @property 451 def start(self): 452 return self.interval.inf 453 454 @property 455 def stop(self): 456 return self.interval.sup 457 458 @property 459 def length(self): 460 return self.stop - self.start + 1 461 462 @property 463 def variables(self): 464 return (self.args[1][0],) 465 466 467class SeqPer(SeqExpr): 468 """ 469 Represents a periodic sequence. 470 471 The elements are repeated after a given period. 472 473 Examples 474 ======== 475 476 >>> from sympy import SeqPer, oo 477 >>> from sympy.abc import k 478 479 >>> s = SeqPer((1, 2, 3), (0, 5)) 480 >>> s.periodical 481 (1, 2, 3) 482 >>> s.period 483 3 484 485 For value at a particular point 486 487 >>> s.coeff(3) 488 1 489 490 supports slicing 491 492 >>> s[:] 493 [1, 2, 3, 1, 2, 3] 494 495 iterable 496 497 >>> list(s) 498 [1, 2, 3, 1, 2, 3] 499 500 sequence starts from negative infinity 501 502 >>> SeqPer((1, 2, 3), (-oo, 0))[0:6] 503 [1, 2, 3, 1, 2, 3] 504 505 Periodic formulas 506 507 >>> SeqPer((k, k**2, k**3), (k, 0, oo))[0:6] 508 [0, 1, 8, 3, 16, 125] 509 510 See Also 511 ======== 512 513 sympy.series.sequences.SeqFormula 514 """ 515 516 def __new__(cls, periodical, limits=None): 517 periodical = sympify(periodical) 518 519 def _find_x(periodical): 520 free = periodical.free_symbols 521 if len(periodical.free_symbols) == 1: 522 return free.pop() 523 else: 524 return Dummy('k') 525 526 x, start, stop = None, None, None 527 if limits is None: 528 x, start, stop = _find_x(periodical), 0, S.Infinity 529 if is_sequence(limits, Tuple): 530 if len(limits) == 3: 531 x, start, stop = limits 532 elif len(limits) == 2: 533 x = _find_x(periodical) 534 start, stop = limits 535 536 if not isinstance(x, (Symbol, Idx)) or start is None or stop is None: 537 raise ValueError('Invalid limits given: %s' % str(limits)) 538 539 if start is S.NegativeInfinity and stop is S.Infinity: 540 raise ValueError("Both the start and end value" 541 "cannot be unbounded") 542 543 limits = sympify((x, start, stop)) 544 545 if is_sequence(periodical, Tuple): 546 periodical = sympify(tuple(flatten(periodical))) 547 else: 548 raise ValueError("invalid period %s should be something " 549 "like e.g (1, 2) " % periodical) 550 551 if Interval(limits[1], limits[2]) is S.EmptySet: 552 return S.EmptySequence 553 554 return Basic.__new__(cls, periodical, limits) 555 556 @property 557 def period(self): 558 return len(self.gen) 559 560 @property 561 def periodical(self): 562 return self.gen 563 564 def _eval_coeff(self, pt): 565 if self.start is S.NegativeInfinity: 566 idx = (self.stop - pt) % self.period 567 else: 568 idx = (pt - self.start) % self.period 569 return self.periodical[idx].subs(self.variables[0], pt) 570 571 def _add(self, other): 572 """See docstring of SeqBase._add""" 573 if isinstance(other, SeqPer): 574 per1, lper1 = self.periodical, self.period 575 per2, lper2 = other.periodical, other.period 576 577 per_length = lcm(lper1, lper2) 578 579 new_per = [] 580 for x in range(per_length): 581 ele1 = per1[x % lper1] 582 ele2 = per2[x % lper2] 583 new_per.append(ele1 + ele2) 584 585 start, stop = self._intersect_interval(other) 586 return SeqPer(new_per, (self.variables[0], start, stop)) 587 588 def _mul(self, other): 589 """See docstring of SeqBase._mul""" 590 if isinstance(other, SeqPer): 591 per1, lper1 = self.periodical, self.period 592 per2, lper2 = other.periodical, other.period 593 594 per_length = lcm(lper1, lper2) 595 596 new_per = [] 597 for x in range(per_length): 598 ele1 = per1[x % lper1] 599 ele2 = per2[x % lper2] 600 new_per.append(ele1 * ele2) 601 602 start, stop = self._intersect_interval(other) 603 return SeqPer(new_per, (self.variables[0], start, stop)) 604 605 def coeff_mul(self, coeff): 606 """See docstring of SeqBase.coeff_mul""" 607 coeff = sympify(coeff) 608 per = [x * coeff for x in self.periodical] 609 return SeqPer(per, self.args[1]) 610 611 612class SeqFormula(SeqExpr): 613 """ 614 Represents sequence based on a formula. 615 616 Elements are generated using a formula. 617 618 Examples 619 ======== 620 621 >>> from sympy import SeqFormula, oo, Symbol 622 >>> n = Symbol('n') 623 >>> s = SeqFormula(n**2, (n, 0, 5)) 624 >>> s.formula 625 n**2 626 627 For value at a particular point 628 629 >>> s.coeff(3) 630 9 631 632 supports slicing 633 634 >>> s[:] 635 [0, 1, 4, 9, 16, 25] 636 637 iterable 638 639 >>> list(s) 640 [0, 1, 4, 9, 16, 25] 641 642 sequence starts from negative infinity 643 644 >>> SeqFormula(n**2, (-oo, 0))[0:6] 645 [0, 1, 4, 9, 16, 25] 646 647 See Also 648 ======== 649 650 sympy.series.sequences.SeqPer 651 """ 652 653 def __new__(cls, formula, limits=None): 654 formula = sympify(formula) 655 656 def _find_x(formula): 657 free = formula.free_symbols 658 if len(free) == 1: 659 return free.pop() 660 elif not free: 661 return Dummy('k') 662 else: 663 raise ValueError( 664 " specify dummy variables for %s. If the formula contains" 665 " more than one free symbol, a dummy variable should be" 666 " supplied explicitly e.g., SeqFormula(m*n**2, (n, 0, 5))" 667 % formula) 668 669 x, start, stop = None, None, None 670 if limits is None: 671 x, start, stop = _find_x(formula), 0, S.Infinity 672 if is_sequence(limits, Tuple): 673 if len(limits) == 3: 674 x, start, stop = limits 675 elif len(limits) == 2: 676 x = _find_x(formula) 677 start, stop = limits 678 679 if not isinstance(x, (Symbol, Idx)) or start is None or stop is None: 680 raise ValueError('Invalid limits given: %s' % str(limits)) 681 682 if start is S.NegativeInfinity and stop is S.Infinity: 683 raise ValueError("Both the start and end value " 684 "cannot be unbounded") 685 limits = sympify((x, start, stop)) 686 687 if Interval(limits[1], limits[2]) is S.EmptySet: 688 return S.EmptySequence 689 690 return Basic.__new__(cls, formula, limits) 691 692 @property 693 def formula(self): 694 return self.gen 695 696 def _eval_coeff(self, pt): 697 d = self.variables[0] 698 return self.formula.subs(d, pt) 699 700 def _add(self, other): 701 """See docstring of SeqBase._add""" 702 if isinstance(other, SeqFormula): 703 form1, v1 = self.formula, self.variables[0] 704 form2, v2 = other.formula, other.variables[0] 705 formula = form1 + form2.subs(v2, v1) 706 start, stop = self._intersect_interval(other) 707 return SeqFormula(formula, (v1, start, stop)) 708 709 def _mul(self, other): 710 """See docstring of SeqBase._mul""" 711 if isinstance(other, SeqFormula): 712 form1, v1 = self.formula, self.variables[0] 713 form2, v2 = other.formula, other.variables[0] 714 formula = form1 * form2.subs(v2, v1) 715 start, stop = self._intersect_interval(other) 716 return SeqFormula(formula, (v1, start, stop)) 717 718 def coeff_mul(self, coeff): 719 """See docstring of SeqBase.coeff_mul""" 720 coeff = sympify(coeff) 721 formula = self.formula * coeff 722 return SeqFormula(formula, self.args[1]) 723 724 def expand(self, *args, **kwargs): 725 return SeqFormula(expand(self.formula, *args, **kwargs), self.args[1]) 726 727class RecursiveSeq(SeqBase): 728 """ 729 A finite degree recursive sequence. 730 731 Explanation 732 =========== 733 734 That is, a sequence a(n) that depends on a fixed, finite number of its 735 previous values. The general form is 736 737 a(n) = f(a(n - 1), a(n - 2), ..., a(n - d)) 738 739 for some fixed, positive integer d, where f is some function defined by a 740 SymPy expression. 741 742 Parameters 743 ========== 744 745 recurrence : SymPy expression defining recurrence 746 This is *not* an equality, only the expression that the nth term is 747 equal to. For example, if :code:`a(n) = f(a(n - 1), ..., a(n - d))`, 748 then the expression should be :code:`f(a(n - 1), ..., a(n - d))`. 749 750 yn : applied undefined function 751 Represents the nth term of the sequence as e.g. :code:`y(n)` where 752 :code:`y` is an undefined function and `n` is the sequence index. 753 754 n : symbolic argument 755 The name of the variable that the recurrence is in, e.g., :code:`n` if 756 the recurrence function is :code:`y(n)`. 757 758 initial : iterable with length equal to the degree of the recurrence 759 The initial values of the recurrence. 760 761 start : start value of sequence (inclusive) 762 763 Examples 764 ======== 765 766 >>> from sympy import Function, symbols 767 >>> from sympy.series.sequences import RecursiveSeq 768 >>> y = Function("y") 769 >>> n = symbols("n") 770 >>> fib = RecursiveSeq(y(n - 1) + y(n - 2), y(n), n, [0, 1]) 771 772 >>> fib.coeff(3) # Value at a particular point 773 2 774 775 >>> fib[:6] # supports slicing 776 [0, 1, 1, 2, 3, 5] 777 778 >>> fib.recurrence # inspect recurrence 779 Eq(y(n), y(n - 2) + y(n - 1)) 780 781 >>> fib.degree # automatically determine degree 782 2 783 784 >>> for x in zip(range(10), fib): # supports iteration 785 ... print(x) 786 (0, 0) 787 (1, 1) 788 (2, 1) 789 (3, 2) 790 (4, 3) 791 (5, 5) 792 (6, 8) 793 (7, 13) 794 (8, 21) 795 (9, 34) 796 797 See Also 798 ======== 799 800 sympy.series.sequences.SeqFormula 801 802 """ 803 804 def __new__(cls, recurrence, yn, n, initial=None, start=0): 805 if not isinstance(yn, AppliedUndef): 806 raise TypeError("recurrence sequence must be an applied undefined function" 807 ", found `{}`".format(yn)) 808 809 if not isinstance(n, Basic) or not n.is_symbol: 810 raise TypeError("recurrence variable must be a symbol" 811 ", found `{}`".format(n)) 812 813 if yn.args != (n,): 814 raise TypeError("recurrence sequence does not match symbol") 815 816 y = yn.func 817 818 k = Wild("k", exclude=(n,)) 819 degree = 0 820 821 # Find all applications of y in the recurrence and check that: 822 # 1. The function y is only being used with a single argument; and 823 # 2. All arguments are n + k for constant negative integers k. 824 825 prev_ys = recurrence.find(y) 826 for prev_y in prev_ys: 827 if len(prev_y.args) != 1: 828 raise TypeError("Recurrence should be in a single variable") 829 830 shift = prev_y.args[0].match(n + k)[k] 831 if not (shift.is_constant() and shift.is_integer and shift < 0): 832 raise TypeError("Recurrence should have constant," 833 " negative, integer shifts" 834 " (found {})".format(prev_y)) 835 836 if -shift > degree: 837 degree = -shift 838 839 if not initial: 840 initial = [Dummy("c_{}".format(k)) for k in range(degree)] 841 842 if len(initial) != degree: 843 raise ValueError("Number of initial terms must equal degree") 844 845 degree = Integer(degree) 846 start = sympify(start) 847 848 initial = Tuple(*(sympify(x) for x in initial)) 849 850 seq = Basic.__new__(cls, recurrence, yn, n, initial, start) 851 852 seq.cache = {y(start + k): init for k, init in enumerate(initial)} 853 seq.degree = degree 854 855 return seq 856 857 @property 858 def _recurrence(self): 859 """Equation defining recurrence.""" 860 return self.args[0] 861 862 @property 863 def recurrence(self): 864 """Equation defining recurrence.""" 865 return Eq(self.yn, self.args[0]) 866 867 @property 868 def yn(self): 869 """Applied function representing the nth term""" 870 return self.args[1] 871 872 @property 873 def y(self): 874 """Undefined function for the nth term of the sequence""" 875 return self.yn.func 876 877 @property 878 def n(self): 879 """Sequence index symbol""" 880 return self.args[2] 881 882 @property 883 def initial(self): 884 """The initial values of the sequence""" 885 return self.args[3] 886 887 @property 888 def start(self): 889 """The starting point of the sequence. This point is included""" 890 return self.args[4] 891 892 @property 893 def stop(self): 894 """The ending point of the sequence. (oo)""" 895 return S.Infinity 896 897 @property 898 def interval(self): 899 """Interval on which sequence is defined.""" 900 return (self.start, S.Infinity) 901 902 def _eval_coeff(self, index): 903 if index - self.start < len(self.cache): 904 return self.cache[self.y(index)] 905 906 for current in range(len(self.cache), index + 1): 907 # Use xreplace over subs for performance. 908 # See issue #10697. 909 seq_index = self.start + current 910 current_recurrence = self._recurrence.xreplace({self.n: seq_index}) 911 new_term = current_recurrence.xreplace(self.cache) 912 913 self.cache[self.y(seq_index)] = new_term 914 915 return self.cache[self.y(self.start + current)] 916 917 def __iter__(self): 918 index = self.start 919 while True: 920 yield self._eval_coeff(index) 921 index += 1 922 923 924def sequence(seq, limits=None): 925 """ 926 Returns appropriate sequence object. 927 928 Explanation 929 =========== 930 931 If ``seq`` is a sympy sequence, returns :class:`SeqPer` object 932 otherwise returns :class:`SeqFormula` object. 933 934 Examples 935 ======== 936 937 >>> from sympy import sequence 938 >>> from sympy.abc import n 939 >>> sequence(n**2, (n, 0, 5)) 940 SeqFormula(n**2, (n, 0, 5)) 941 >>> sequence((1, 2, 3), (n, 0, 5)) 942 SeqPer((1, 2, 3), (n, 0, 5)) 943 944 See Also 945 ======== 946 947 sympy.series.sequences.SeqPer 948 sympy.series.sequences.SeqFormula 949 """ 950 seq = sympify(seq) 951 952 if is_sequence(seq, Tuple): 953 return SeqPer(seq, limits) 954 else: 955 return SeqFormula(seq, limits) 956 957 958############################################################################### 959# OPERATIONS # 960############################################################################### 961 962 963class SeqExprOp(SeqBase): 964 """ 965 Base class for operations on sequences. 966 967 Examples 968 ======== 969 970 >>> from sympy.series.sequences import SeqExprOp, sequence 971 >>> from sympy.abc import n 972 >>> s1 = sequence(n**2, (n, 0, 10)) 973 >>> s2 = sequence((1, 2, 3), (n, 5, 10)) 974 >>> s = SeqExprOp(s1, s2) 975 >>> s.gen 976 (n**2, (1, 2, 3)) 977 >>> s.interval 978 Interval(5, 10) 979 >>> s.length 980 6 981 982 See Also 983 ======== 984 985 sympy.series.sequences.SeqAdd 986 sympy.series.sequences.SeqMul 987 """ 988 @property 989 def gen(self): 990 """Generator for the sequence. 991 992 returns a tuple of generators of all the argument sequences. 993 """ 994 return tuple(a.gen for a in self.args) 995 996 @property 997 def interval(self): 998 """Sequence is defined on the intersection 999 of all the intervals of respective sequences 1000 """ 1001 return Intersection(*(a.interval for a in self.args)) 1002 1003 @property 1004 def start(self): 1005 return self.interval.inf 1006 1007 @property 1008 def stop(self): 1009 return self.interval.sup 1010 1011 @property 1012 def variables(self): 1013 """Cumulative of all the bound variables""" 1014 return tuple(flatten([a.variables for a in self.args])) 1015 1016 @property 1017 def length(self): 1018 return self.stop - self.start + 1 1019 1020 1021class SeqAdd(SeqExprOp): 1022 """Represents term-wise addition of sequences. 1023 1024 Rules: 1025 * The interval on which sequence is defined is the intersection 1026 of respective intervals of sequences. 1027 * Anything + :class:`EmptySequence` remains unchanged. 1028 * Other rules are defined in ``_add`` methods of sequence classes. 1029 1030 Examples 1031 ======== 1032 1033 >>> from sympy import EmptySequence, oo, SeqAdd, SeqPer, SeqFormula 1034 >>> from sympy.abc import n 1035 >>> SeqAdd(SeqPer((1, 2), (n, 0, oo)), EmptySequence) 1036 SeqPer((1, 2), (n, 0, oo)) 1037 >>> SeqAdd(SeqPer((1, 2), (n, 0, 5)), SeqPer((1, 2), (n, 6, 10))) 1038 EmptySequence 1039 >>> SeqAdd(SeqPer((1, 2), (n, 0, oo)), SeqFormula(n**2, (n, 0, oo))) 1040 SeqAdd(SeqFormula(n**2, (n, 0, oo)), SeqPer((1, 2), (n, 0, oo))) 1041 >>> SeqAdd(SeqFormula(n**3), SeqFormula(n**2)) 1042 SeqFormula(n**3 + n**2, (n, 0, oo)) 1043 1044 See Also 1045 ======== 1046 1047 sympy.series.sequences.SeqMul 1048 """ 1049 1050 def __new__(cls, *args, **kwargs): 1051 evaluate = kwargs.get('evaluate', global_parameters.evaluate) 1052 1053 # flatten inputs 1054 args = list(args) 1055 1056 # adapted from sympy.sets.sets.Union 1057 def _flatten(arg): 1058 if isinstance(arg, SeqBase): 1059 if isinstance(arg, SeqAdd): 1060 return sum(map(_flatten, arg.args), []) 1061 else: 1062 return [arg] 1063 if iterable(arg): 1064 return sum(map(_flatten, arg), []) 1065 raise TypeError("Input must be Sequences or " 1066 " iterables of Sequences") 1067 args = _flatten(args) 1068 1069 args = [a for a in args if a is not S.EmptySequence] 1070 1071 # Addition of no sequences is EmptySequence 1072 if not args: 1073 return S.EmptySequence 1074 1075 if Intersection(*(a.interval for a in args)) is S.EmptySet: 1076 return S.EmptySequence 1077 1078 # reduce using known rules 1079 if evaluate: 1080 return SeqAdd.reduce(args) 1081 1082 args = list(ordered(args, SeqBase._start_key)) 1083 1084 return Basic.__new__(cls, *args) 1085 1086 @staticmethod 1087 def reduce(args): 1088 """Simplify :class:`SeqAdd` using known rules. 1089 1090 Iterates through all pairs and ask the constituent 1091 sequences if they can simplify themselves with any other constituent. 1092 1093 Notes 1094 ===== 1095 1096 adapted from ``Union.reduce`` 1097 1098 """ 1099 new_args = True 1100 while new_args: 1101 for id1, s in enumerate(args): 1102 new_args = False 1103 for id2, t in enumerate(args): 1104 if id1 == id2: 1105 continue 1106 new_seq = s._add(t) 1107 # This returns None if s does not know how to add 1108 # with t. Returns the newly added sequence otherwise 1109 if new_seq is not None: 1110 new_args = [a for a in args if a not in (s, t)] 1111 new_args.append(new_seq) 1112 break 1113 if new_args: 1114 args = new_args 1115 break 1116 1117 if len(args) == 1: 1118 return args.pop() 1119 else: 1120 return SeqAdd(args, evaluate=False) 1121 1122 def _eval_coeff(self, pt): 1123 """adds up the coefficients of all the sequences at point pt""" 1124 return sum(a.coeff(pt) for a in self.args) 1125 1126 1127class SeqMul(SeqExprOp): 1128 r"""Represents term-wise multiplication of sequences. 1129 1130 Explanation 1131 =========== 1132 1133 Handles multiplication of sequences only. For multiplication 1134 with other objects see :func:`SeqBase.coeff_mul`. 1135 1136 Rules: 1137 * The interval on which sequence is defined is the intersection 1138 of respective intervals of sequences. 1139 * Anything \* :class:`EmptySequence` returns :class:`EmptySequence`. 1140 * Other rules are defined in ``_mul`` methods of sequence classes. 1141 1142 Examples 1143 ======== 1144 1145 >>> from sympy import EmptySequence, oo, SeqMul, SeqPer, SeqFormula 1146 >>> from sympy.abc import n 1147 >>> SeqMul(SeqPer((1, 2), (n, 0, oo)), EmptySequence) 1148 EmptySequence 1149 >>> SeqMul(SeqPer((1, 2), (n, 0, 5)), SeqPer((1, 2), (n, 6, 10))) 1150 EmptySequence 1151 >>> SeqMul(SeqPer((1, 2), (n, 0, oo)), SeqFormula(n**2)) 1152 SeqMul(SeqFormula(n**2, (n, 0, oo)), SeqPer((1, 2), (n, 0, oo))) 1153 >>> SeqMul(SeqFormula(n**3), SeqFormula(n**2)) 1154 SeqFormula(n**5, (n, 0, oo)) 1155 1156 See Also 1157 ======== 1158 1159 sympy.series.sequences.SeqAdd 1160 """ 1161 1162 def __new__(cls, *args, **kwargs): 1163 evaluate = kwargs.get('evaluate', global_parameters.evaluate) 1164 1165 # flatten inputs 1166 args = list(args) 1167 1168 # adapted from sympy.sets.sets.Union 1169 def _flatten(arg): 1170 if isinstance(arg, SeqBase): 1171 if isinstance(arg, SeqMul): 1172 return sum(map(_flatten, arg.args), []) 1173 else: 1174 return [arg] 1175 elif iterable(arg): 1176 return sum(map(_flatten, arg), []) 1177 raise TypeError("Input must be Sequences or " 1178 " iterables of Sequences") 1179 args = _flatten(args) 1180 1181 # Multiplication of no sequences is EmptySequence 1182 if not args: 1183 return S.EmptySequence 1184 1185 if Intersection(*(a.interval for a in args)) is S.EmptySet: 1186 return S.EmptySequence 1187 1188 # reduce using known rules 1189 if evaluate: 1190 return SeqMul.reduce(args) 1191 1192 args = list(ordered(args, SeqBase._start_key)) 1193 1194 return Basic.__new__(cls, *args) 1195 1196 @staticmethod 1197 def reduce(args): 1198 """Simplify a :class:`SeqMul` using known rules. 1199 1200 Explanation 1201 =========== 1202 1203 Iterates through all pairs and ask the constituent 1204 sequences if they can simplify themselves with any other constituent. 1205 1206 Notes 1207 ===== 1208 1209 adapted from ``Union.reduce`` 1210 1211 """ 1212 new_args = True 1213 while new_args: 1214 for id1, s in enumerate(args): 1215 new_args = False 1216 for id2, t in enumerate(args): 1217 if id1 == id2: 1218 continue 1219 new_seq = s._mul(t) 1220 # This returns None if s does not know how to multiply 1221 # with t. Returns the newly multiplied sequence otherwise 1222 if new_seq is not None: 1223 new_args = [a for a in args if a not in (s, t)] 1224 new_args.append(new_seq) 1225 break 1226 if new_args: 1227 args = new_args 1228 break 1229 1230 if len(args) == 1: 1231 return args.pop() 1232 else: 1233 return SeqMul(args, evaluate=False) 1234 1235 def _eval_coeff(self, pt): 1236 """multiplies the coefficients of all the sequences at point pt""" 1237 val = 1 1238 for a in self.args: 1239 val *= a.coeff(pt) 1240 return val 1241