Probabilities - how high can you go? I previously asked a question for how to compute a probability quickly and accurately.  However, evidently it was too easy as a closed form solution was given!  Here is a more difficult version.

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 a copy of itself. That is A[1] = A[n+1] if indexing from 1, A[2] = A[n+2] and so on.  A now has length 2n. 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 B with A[1+j,...,n+j] for different j =0,1,2,....

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 first two inner products are 0 and 2.

Task

For each j, starting with j=1, your code should output the probability that all the first j+1 inner products are zero for every n=j,...,50.

Copying the table produced by Martin Büttner for j=1 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


Score

Your score is the largest j your code completes  in 1 minute on my computer.  To clarify a little, each j gets one minute. Note that the dynamic programming code in the previous linked question will do this easily for j=1.

Tie breaker

If two entries get the same j score then the winning entry will be the one that gets to the highest n in one minute on my machine for that j. If the two best entries are equal on this criterion too then the winner will be the answer submitted first.

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.

My Machine The timings will be run on my machine. This is a standard ubuntu install on an AMD FX-8350 Eight-Core Processor.  This also means I need to be able to run your code.

Winning entries


j=2 in Python by Mitch Schwartz.
j=2 in Python by feersum. Currently the fastest entry.
Python, j=2

The dynamic programming approach for j = 1 in my answer to the previous question does not need many modifications to work for higher j, but gets slow quickly. Table for reference:

n   p(n)

2   3/8
3   11/64
4   71/512
5   323/4096
6   501/8192
7   2927/65536
8   76519/2097152
9   490655/16777216
10  207313/8388608


And the code:

from time import*
from fractions import*
from collections import*

def main():
    N_MAX=50

    T=time()

    n=2
    Y=defaultdict(lambda:0)
    numer=0

    for a1 in [1]:
        for b1 in (1,0):
            for a2 in (1,-1):
                for b2 in (1,0,0,-1):
                    if not a1*b1+a2*b2 and not a2*b1+a1*b2:
                        numer+=1
                    Y[(a1,a2,b1,b2,a1*b1+a2*b2,a2*b1,0)]+=1

    thresh=N_MAX-1

    while time() &lt;= T+60:
        print('%d %s'%(n,Fraction(numer,8**n/4)))

        if thresh&lt;2:
            print('reached N_MAX with %.2f seconds remaining'%(T+60-time()))
            return

        n+=1
        X=Y
        Y=defaultdict(lambda:0)
        numer=0

        for a1,a2,b1,b2,s,t,u in X:
            if not ( abs(s)&lt;thresh and abs(t)&lt;thresh+1 and abs(u)&lt;thresh+2 ):
                continue

            c=X[(a1,a2,b1,b2,s,t,u)]

            # 1,1

            if not s+1 and not t+b2+a1 and not u+b1+a1*b2+a2: numer+=c
            Y[(a1,a2,b2,1,s+1,t+b2,u+b1)]+=c

            # -1,1

            if not s-1 and not t-b2+a1 and not u-b1+a1*b2+a2: numer+=c
            Y[(a1,a2,b2,1,s-1,t-b2,u-b1)]+=c

            # 1,-1

            if not s-1 and not t+b2-a1 and not u+b1+a1*b2-a2: numer+=c
            Y[(a1,a2,b2,-1,s-1,t+b2,u+b1)]+=c

            # -1,-1

            if not s+1 and not t-b2-a1 and not u-b1+a1*b2-a2: numer+=c
            Y[(a1,a2,b2,-1,s+1,t-b2,u-b1)]+=c

            # 1,0

            c+=c

            if not s and not t+b2 and not u+b1+a1*b2: numer+=c
            Y[(a1,a2,b2,0,s,t+b2,u+b1)]+=c

            # -1,0

            if not s and not t-b2 and not u-b1+a1*b2: numer+=c
            Y[(a1,a2,b2,0,s,t-b2,u-b1)]+=c

        thresh-=1

main()


Here we are keeping track of the first two elements of A, the last two elements of B (where b2 is the last element), and the inner products of (A[:n], B), (A[1:n], B[:-1]), and (A[2:n], B[:-2]).