Return to Answer

• added a proof of the above formula.

• fixed the time_limit.

  def solve(): # straightforward implementation

  from time import time from itertools import count def binom(n, r): return facts[n] // (facts[r] * facts[n - r]) def p(N): ans = 0 for i in range(1 + N // 2): t = binom(2 * (N - 2 * i), N - 2 * i) t *= binom(N, 2 * i) t *= binom(4 * i, 2 * i) ans += t e = (ans & -ans).bit_length() - 1 numer = ans >> e denom = 1 << (3 * N - 1 - e) return numer, denom facts = [1] time_limit = 10.0 + time() for i in count(1): facts.append(facts[-1] * (2 * i - 1)) facts.append(facts[-1] * (2 * i)) n, d = p(i) if time() > time_limit: break print("%d %d/%d" % (i, n, d)) 

  solve()

• added a proof of the above formula.
• fixed the time_limit.
• added a PARI/GP code.

Python

def solve(): # straightforward implementation from time import time from itertools import count def binom(n, r): return facts[n] // (facts[r] * facts[n - r]) def p(N): ans = 0 for i in range(1 + N // 2): t = binom(2 * (N - 2 * i), N - 2 * i) t *= binom(N, 2 * i) t *= binom(4 * i, 2 * i) ans += t e = (ans & -ans).bit_length() - 1 numer = ans >> e denom = 1 << (3 * N - 1 - e) return numer, denom facts = [1] time_limit = 10.0 + time() for i in count(1): facts.append(facts[-1] * (2 * i - 1)) facts.append(facts[-1] * (2 * i)) n, d = p(i) if time() > time_limit: break print("%d %d/%d" % (i, n, d)) solve() 

PARI/GP

p(n) = polcoeff( (exp(x/2) + 1) * besseli(0, x/4) ^ 2, n) * n!; 

Edit on 2015-06-11:

An exponential generating function of p(n) is

where I_0(x) is the modified Bessel function of the first kind.

Edit on 2015-06-11:

A closed formula of p(n) is

Then, c(n) = \sum_{i=0}^n a(n,i) * b(i) * b(n-i).

Since p(n) = c(n) / 8^n, we can get the above closed formula.

A closed-form formula of p(n) is

Then, c(n) = \sum_{i=0}^n a(n, i) * b(i) * b(n-i).

Since p(n) = c(n) / 8^n, we can get the above closed-form formula.