Calculate a probability exactly and quickly [This is a partner question to Calculate a probability exactly ]

This task is about writing code to compute a probability exactly and quickly. The output should be a precise probability written as a fraction in its most reduced form. That is it should never output 4/8 but rather 1/2.

For some positive integer n, consider a uniformly random string of 1s and -1s of length n and call it A. Now concatenate to A its first value. That is A[1] = A[n+1] if indexing from 1.  A now has length n+1. Now also consider a second random string of length n whose first n values are -1, 0, or 1 with probability 1/4,1/2, 1/4 each and call it B.

Now consider the inner product of A[1,...,n] and B and the inner product of A[2,...,n+1] and B.

For example, consider n=3. Possible values for A and B could be A = [-1,1,1,-1] and B=[0,1,-1]. In this case the two inner products are 0 and 2.

Your code must output the probability that both inner products are zero.

Copying the table produced by Martin Büttner we have the following sample results.

n   P(n)
1   1/2
2   3/8
3   7/32
4   89/512
5   269/2048
6   903/8192
7   3035/32768
8   169801/2097152


Languages and libraries

You can use any freely available language and libraries you like. I must be able to run your code so please include a full explanation for how to run/compile your code in linux if at all possible.

The task

Your code must start with n=1 and give the correct output for each increasing n on a separate line. It should stop after 10 seconds.

The score 

The score is simply the highest n reached before your code stops after 10 seconds when run on my computer.  If there is a tie, the winner be the one to get to the highest score quickest.

Table of entries


n = 64 in Python. Version 1 by Mitch Schwartz
n = 106 in Python.  Version June 11 2015 by  Mitch Schwartz
n = 151 in C++. Port of Mitch Schwartz's answer by kirbyfan64sos
n = 165 in Python. Version June 11 2015 the "pruning" version by Mitch Schwartz with  N_MAX = 165.
n = 945 in Python by Min_25 using an exact formula. Amazing!
n = 1228 in Python by Mitch Schwartz using another exact formula (based on Min_25's previous answer).
n = 2761 in Python by Mitch Schwartz using a faster implementation of the same exact formula.
n = 3250 in Python using Pypy by Mitch Schwartz using the same implementation. This score needs pypy MitchSchwartz-faster.py |tail to avoid console scrolling overhead.
Python

A closed-form formula of p(n) is 



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:
- updated the Python code.

Edit on 2015-06-13:
- 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 &amp; -ans).bit_length() - 1
    numer = ans &gt;&gt; e
    denom = 1 &lt;&lt; (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() &gt; 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!;


Proof:
This problem is similar to a 2-dimensional (restricted) random walk problem.

If A[i] = A[i+1], we can move from (x, y) to (x+1, y+1) [1 way], (x, y) [2 ways] or (x-1, y-1) [1 way]. 

If A[i] != A[i+1], we can move from (x, y) to (x-1, y+1) [1 way], (x, y) [2 ways] or (x+1, y-1) [1 way].

Let a(n, m) = [x^m]((x+1)^n + (x-1)^n), b(n) = [x^n](1+x)^{2n} and c(n) be the number of ways to move from (0, 0) to (0, 0) with n steps.

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.