Python code to solve the Lonesome King problem, http://fivethirtyeight.com/features/the-puzzle-of-the-lonesome-king/

main.py

import statistics
import random
import numpy as np
import time



def repeated_cull_sim(initial_size, runs):
    sum = 0
    init_time = time.time()

    for j in range(runs):
        sum += repeated_cull(initial_size)

    print("N=", initial_size, "  runs =",  runs, "  num_lonely=", runs-sum, "  prob_lonely=", (runs-sum)/runs,
          '    run time (hr) = {0:.3f}'.format((time.time() - init_time)/3600.0))


def repeated_cull(initial_size):
    size = int(initial_size)
    prob = [1.0, 0.0, 1.0, 0.25, 0.4074074074074074, 0.53125, 0.5838577777777778, 0.5611357992938068, 0.5109616766040821]

    while size >= len(prob):
        size = cull(size)

    return np.random.random() > prob[size]



def cull(start_num):
    subjects = np.ones(start_num, dtype=np.int)
    rand = np.random.randint(0, start_num-1, start_num)

    for i in range(start_num):
        if rand[i] >= i:
            rand[i] += 1

    subjects[rand] = 0

    return np.sum(subjects)




# main program
# first parameter is initial population size and the second is the number of simulations to run
repeated_cull_sim(56000, 1000*1000)

Running a million simulations with N = 56000 took about 10 hours on my machine.