Skip to content

Instantly share code, notes, and snippets.

@pjt33

pjt33/A003122.py

Last active Feb 26, 2017
Embed
What would you like to do?
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

This comment has been minimized.

Copy link
Owner Author

@pjt33 pjt33 commented Feb 26, 2017

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.