Computes odds of Bob winning in Daniel Litt's diabolical puzzle

hh_vs_ht.py

# Flip a fair coin 100 times--it gives a sequence of heads (H) and tails (T).
# For each HH in the sequence of flips, Alice gets a point; for each HT,
# Bob does, so e.g. for the sequence THHHT Alice gets 2 points and Bob
# gets 1 point. Who is most likely to win?
#   -- https://twitter.com/littmath/status/1769044719034647001

# Based on https://gist.github.com/llllvvuu/02e08ee87b2191c0f3f504072be73a0f#file-hh_vs_ht-c,
# a _very_ clever algorithm by L (llllvvuu). See https://twitter.com/llllvvuu/status/1770557954372382765
# I implemented this to try to understand llllvvuu's algorithm and as a bonus to benefit from
# Python's bignum support.

def hh_vs_ht(maxflips, printIntermediates = False, printFinal = False):
    # Uses two pairs of lists H and T (along with their previous iteration).
    # These arrays represent the number of ways Bob can be ahead of Alice by point difference,
    # depending on whether the last flip was heads or tails. Because the difference can be negative
    # and python arrays can't have negative indexes, the arrays are indexed by
    # offset - (bob_points - alice_points)
    sz, origin = maxflips * 2, maxflips+1
    H, H_prev = [0]*sz, [0]*sz
    T, T_prev = [0]*sz, [0]*sz

    # This represents the initial condition which is as if a tail had already been flipped
    # resulting in a tie
    T[origin] = 1
    for n in range(1, maxflips+1):
        H, H_prev = H_prev, H
        T, T_prev = T_prev, T

        # Alice can be ahead (Bob can be behind, rather) by as much as n
        # but Bob can only be ahead by as much as floor(n/2)
        for d in range(-n, n//2+1):
            H[origin + d] = T_prev[origin + d] + H_prev[origin + d + 1];
            T[origin + d] = T_prev[origin + d] + H_prev[origin + d - 1];
        # print("State:", n, d, H, T)   # For those following along at home...

        begin, end = origin-maxflips, origin
        alice_wins = sum(H[begin:end]) + sum(T[begin:end])

        begin, end = origin+1, origin+(n//2)+1
        bob_wins = sum(H[begin:end]) + sum(T[begin:end])

        # Optionally print intermediate results as they occur
        if printIntermediates or (printFinal and n == maxflips):
            flips = 2**n
            if n <= 25:  # Numbers get ridiculously long
                print("{0:3d} a={1}({2:.5f}%) b={3}({4:.5f}%) b-a={5}({6:.5f}%)".format(n,
                    alice_wins, alice_wins/flips*100,
                    bob_wins, bob_wins/flips*100,
                    bob_wins-alice_wins, (bob_wins-alice_wins)/flips*100))
            else:
                print("{0:3d} a={1:.5f}% b={2:.5f}% b-a={3:.5f}%".format(n,
                    alice_wins/flips*100,
                    bob_wins/flips*100,
                    (bob_wins-alice_wins)/flips*100))
    return (alice_wins, bob_wins)

print(hh_vs_ht(100))
hh_vs_ht(100, printIntermediates=True)
hh_vs_ht(1000, printFinal=True)
hh_vs_ht(10000, printFinal=True)     # patience grasshopper
# hh_vs_ht(100000, printFinal=True)  # are you feeling lucky?