https://twitter.com/bitshiftmask/status/1323809212875608066?s=20

dp.py

"""
https://twitter.com/bitshiftmask/status/1323809212875608066?s=20

Let P(i, j) be the probability that if there are `i` devices, then after one
round `j` devices selected a unique number, so they are out. The expected
value of the number of rounds required so that each device selected a unique number is:

E(n) = p(n, n) * E(0) + p(n, n - 1) * E(1) + ... + p(n, 0) * E(n)

E(n) * (1 - p(n, 0)) = p(n, n) * E(0) + p(n, n - 1) * E(1) + ...

       p(n, n) * E(0) + p(n, n - 1) * E(1) + ...
E(n) = -----------------------------------------
                      (1 - p(n, 0))
"""

from pprint import pprint

graph = {}


def build_graph(remaining, zero, one, two_or_more, prev, weight):
    node = remaining, zero, one, two_or_more

    edges = graph.get(prev, [])
    edges.append((node, weight))
    graph[prev] = edges

    if node in graph or remaining == 0:
        return

    if zero:
        build_graph(remaining - 1, zero - 1, one + 1, two_or_more, node, zero)

    if one:
        build_graph(remaining - 1, zero, one - 1, two_or_more + 1, node, one)

    if two_or_more:
        build_graph(remaining - 1, zero, one, two_or_more, node, two_or_more)


def transition(devices, colors):
    """
    Return a 1D table P, such that P[i] is the probability
    that if there are `devices` devices and `colors` numbers
    after one round the number of devices that selected a unique
    number is `i`.
    """
    graph.clear()

    build_graph(devices, colors, 0, 0, None, 1)

    degree = {}

    for node, neighbors in graph.items():
        if node is None:
            continue

        for next_node, _ in neighbors:
            degree[next_node] = degree.get(next_node, 0) + 1

    start = (devices, colors, 0, 0)
    ways = {start: 1}
    queue = [start]

    result = [0] * (devices + 1)

    pointer = 0

    while pointer < len(queue):
        node = queue[pointer]

        remaining, _, one, _ = node

        if remaining == 0:
            result[one] += ways[node]

        pointer += 1

        for next_node, weight in graph.get(node, []):
            ways[next_node] = ways.get(next_node, 0) + weight * ways[node]

            degree[next_node] -= 1

            if degree[next_node] == 0:
                queue.append(next_node)

    chances = sum(result)
    assert chances == colors**devices
    return [x / chances for x in result]


DEVICES = 32
COLORS = 256

matrix = []
for start in range(DEVICES + 1):
    row = transition(start, COLORS)
    matrix.append(row)


expected_value = [0.] * (DEVICES + 1)

for i in range(1, DEVICES + 1):
    rhs = 1.
    for j in range(i):
        rhs += expected_value[j] * matrix[i][i - j]
    rhs /= 1. - matrix[i][0]
    expected_value[i] = rhs

    print(f"Devices: {i} Colors: {COLORS} | Expected value: {rhs}")

simulation.py

from random import randint
from itertools import count


def go(devices, colors):
    rounds = 0
    while devices > 0:
        rounds += 1

        freq = [0] * colors
        for col in range(devices):
            freq[randint(0, colors-1)] += 1

        for tot in freq:
            if tot == 1:
                devices -= 1

    return rounds


acc = 0
for i in count(1):
    cur = go(32, 256)
    acc += cur
    if i % 1000 == 0:
        print("simulations", i, "mean", acc / i)