Python code to solve the Secret Santa puzzle (Riddler Classic) from http://fivethirtyeight.com/features/can-you-unmask-the-secret-santas/

main.py

import numpy as np
import math
import matplotlib.pyplot as plot
import time


def gen_rand_derang(n):
    comparison = np.arange(n)
    sigma = np.random.permutation(n)

    # check if random permutation has one or more fixed points. If yes, discard and generate a new one
    while (np.count_nonzero(comparison-sigma) < n):
        sigma = np.random.permutation(n)

    return sigma


def sample_longest_cycle(n):
    derang = gen_rand_derang(n)

    sum = 0
    max_cycle_length = 0

    while (sum < n/2):
        head = 0
        while (derang[head] == 0): head += 1

        cycle_length = 2
        pointee = head
        pointee = derang[pointee]
        derang[head] = 0

        tmp = pointee
        pointee = derang[pointee]
        derang[tmp] = 0

        while (pointee != head):
            cycle_length += 1
            tmp = pointee
            pointee = derang[pointee]
            derang[tmp] = 0

        derang[pointee] = 0

        sum += cycle_length
        if cycle_length > max_cycle_length:
            max_cycle_length = cycle_length

        #print(derang, "    cycle_length=", cycle_length)

    return max_cycle_length


  
def estimate_prob_win(n, runs):
    init_time = time.time()
    win_count = 0

    for j in range(runs):
        if sample_longest_cycle(n) < (n / 2 + 1):
            win_count += 1

    prob_win = win_count/runs

    print("N =", runs, '   Run Time (hr) = {0:.3f}'.format((time.time() - init_time)/3600.0))
    print('     win prob = {0:.8f}'.format(prob_win))
    print(' 95% conf int = {0:.8f}'.format(1.96*math.sqrt(prob_win * (1.0-prob_win)/runs) ))


# Main program
# Estimates the probability of winning using the pointer chasing strategy, where winning is defined by everyone finding
# out who their secret Santa is.
# Under this strategy, a win occurs iff the randomly generated derangment contains cycles of length 21 or less
estimate_prob_win(41, 400*1000*1000)

After running 400 million simulations, the probability of winning is 0.318323 ± 0.000045 (95% confidence interval).