Alice, Bob, Coin, and Dynamic Programming (or Daniel Litt)

abcd.py

"""This is a python code answers the question from Daniel Litt (@littmath) on coin tosses.
Link: https://x.com/littmath/status/1769044719034647001
One can count the exact number of cases where Alice wins / Bob wins / or draw for given n, using dynamic programming (in other words, recursive formula).
In particular, the keys of the dictionary `dp` have a form of `(n, c, l)`, where `dp[(n, c, l)]` counts the number of cases where
1) n is the number of tosses,
2) c is the (score of Alice - score of Bob),
3) l is the last flip ('H' or 'T').
Note that the number of draws is exactly the sequence A163493 of OEIS (https://oeis.org/A163493), which also appears as a Problem 11610 of AMM by Richard Stanley.
"""

from typing import List, Tuple
from collections import defaultdict
from timeit import default_timer

import numpy as np
import matplotlib.pyplot as plt

dp = defaultdict(int)


def c_bound(n: int, l: str) -> Tuple[int, int]:
    if n == 1:
        return 0, 0
    elif n == 2:
        if l == "H":
            return 0, 1
        else:
            return -1, 0
    if l == "H":
        max_c = n - 1
        if n % 2 == 1:
            min_c = - (n - 1) // 2
        else:
            min_c = - (n - 2) // 2
    else:  # l == "T"
        max_c = n - 3
        if n % 2 == 1:
            min_c = - (n - 1) // 2
        else:
            min_c = - n // 2
    return min_c, max_c


def check_c_range(c: int, n: int, l: str) -> bool:
    min_c, max_c = c_bound(n, l)
    return min_c <= c <= max_c


def run_dp(n: int):
    for i in range(1, n + 1):
        for l in ["H", "T"]:
            min_c, max_c = c_bound(i, l)
            # print(f"{min_c=}, {max_c=}")
            for c in range(min_c, max_c + 1):
                if i == 1:
                    dp[(i, c, l)] = 1
                elif i == 2:
                    if c == 1:
                        dp[(i, c, l)] = int(l == "H")
                    elif c == 0:
                        dp[(i, c, l)] = 1
                    elif c == -1:
                        dp[(i, c, l)] = int(l == "T")
                    else:
                        pass
                else:
                    if l == "H":
                        # a_n,c,H = a_n-1,c-1,H + a_n-1,c,T
                        if check_c_range(c-1, i-1, "H"):
                            dp[(i, c, l)] += dp[(i-1, c-1, "H")]
                        if check_c_range(c, i-1, "T"):
                            dp[(i, c, l)] += dp[(i-1, c, "T")]
                    else:  # l == "T"
                        # a_n,c,T = a_n-1,c+1,H + a_n-1,c,T
                        if check_c_range(c+1, i-1, "H"):
                            dp[(i, c, l)] += dp[(i-1, c+1, "H")]
                        if check_c_range(c, i-1, "T"):
                            dp[(i, c, l)] += dp[(i-1, c, "T")]


def result_cnt(n: int) -> Tuple[List[int], List[int], List[int]]:
    cnts_A = [0] * (n + 1)
    cnts_B = [0] * (n + 1)
    cnts_draw = [0] * (n + 1)
    for k, v in dp.items():
        if k[1] > 0:
            cnts_A[k[0]] += v
        elif k[1] == 0:
            cnts_draw[k[0]] += v
        else:
            cnts_B[k[0]] += v
    return cnts_A, cnts_B, cnts_draw


if __name__ == "__main__":
    N = 100
    st = default_timer()
    run_dp(N)
    et = default_timer()
    print(f"time elapsed: {et - st: .4f}s")  # less than a second when N = 1000 with M2 macbook.

    cnts_A, cnts_B, cnts_draw = result_cnt(N)
    # ratios
    cnts_A_ = [a / (2 ** i) for (i, a) in enumerate(cnts_A)]
    cnts_B_ = [b / (2 ** i) for (i, b) in enumerate(cnts_B)]
    cnts_draw_ = [d / (2 ** i) for (i, d) in enumerate(cnts_draw)]
    
    print(f"{N=}, Alice: {cnts_A[N]}, Bob: {cnts_B[N]}, Draw: {cnts_draw[N]}")
    # For N = 100, you get Alice: 580127949239420834381088427404, Bob: 615866238418960422359689555420, Draw: 71656412569848144755925222552

    # If you want to see plot the ratio as n grows, uncomment the following lines L112 - L118.
    # xs = list(range(1, N + 1))
    # plt.plot(xs, cnts_A_[1:], label='Alice')
    # plt.plot(xs, cnts_B_[1:], label='Bob')
    # plt.plot(xs, cnts_draw_[1:], label='Draw')
    # plt.xlabel('Number of tosses')
    # plt.legend()
    # plt.show()

    # I have two conjectures on the growth of number of wins and draws.
    # First, it seems that the ratio (Bob wins - Alice wins) / (Draws) converges to 1/2 when n -> infty.
    m = 1
    for n in range(3, N + 1, m):
        cnt_A = cnts_A[n]
        cnt_B = cnts_B[n]
        cnt_draw = cnts_draw[n]
        cnt_total = 2 ** n
        print(f"n: {n}, ratio: {(cnt_B - cnt_A) / cnt_draw}")
    
    # Second, the growth of the number of draws is asymptotically 2^(n - 5/6) / sqrt(n), as n -> infty.
    logd = [np.log(d) for d in cnts_draw_[2:]]
    xs = list(range(2, N + 1))
    logxs = [np.log(x) for x in xs]
    c = - (5 / 6) * np.log(2)
    lins = [-x / 2 + c for x in logxs]
    plt.plot(logxs, logd, label="log draw ratio")
    plt.plot(logxs, lins, label="linear, y = - x / 2 - (5 / 6) * log 2")
    plt.legend()
    plt.show()