{{ message }}

Instantly share code, notes, and snippets.

# pjt33/A003122.py

Last active Feb 26, 2017
Calculate OEIS A003122 with heavy memoisation
 #!/usr/bin/python3 from functools import lru_cache class Partition: def __init__(self, maxPart, tail=None): self.part = maxPart self.tail = tail if tail is None: self.partFreq = 1 self.numParts = 1 tailTerm = 1 else: self.partFreq = tail.partFreq + 1 if tail.part == self.part else 1 self.numParts = tail.numParts + 1 tailTerm = tail.term self.term = tailTerm * r(self.part) * self.numParts // self.partFreq @lru_cache(maxsize=100) def partition(n): accum = [] for maxPart in range(1, n): for tail in partition(n - maxPart): if tail.part > maxPart: break accum.append(Partition(maxPart, tail)) accum.append(Partition(n)) return accum @lru_cache(maxsize=100) def r(n): catalan = binom(2 * n, n) // (n + 1) return catalan * catalan @lru_cache(maxsize=10000) def binom(n, k): if k < 0 or k > n: return 0 if k == 0 or k == n: return 1 return binom(n - 1, k - 1) + binom(n - 1, k) @lru_cache(maxsize=100) def p(n): # 3 (2n + 2)! / ((n + 3)! (n + 1)!) \binom{2n + 2}{n} val = 3 * binom(2 * n + 2, n) # Note that it's important to do this as (2n + 2)! / (n + 1)! and then # divide by (n + 3)! rather than (2n + 2)! / (n + 3)! and then divide by # (n + 1)! because the latter gives the wrong result for n == 0 for x in range(n + 2, 2 * n + 3): val *= x for y in range(2, n + 4): val //= y return val @lru_cache(maxsize=10000) def A(n, s): sum = 0 for p in partition(s): sum += binom(n, p.numParts) * p.term return sum @lru_cache(maxsize=100) def a(n): val = p(n) for s in range(0, n): val -= A(s + 3, n - s) * a(s) return val if __name__ == "__main__": for n in range(0, 40): print(a(n))

### pjt33 commented Feb 26, 2017

 NB version 2 following corrections for PEP8 as per http://pep8online.com/