Continually replace numbers in hat by absolute value of difference; what's the distribution of the final number in the hat?

Question

I was in a class recently, and was asked the following question:

Imagine you have a hat containing pieces of paper numbered from $1$ to $N$. You remove two pieces of paper at random from the hat, and replace them with the absolute value of the difference between the two numbers. You repeat this process until there is one piece of paper remaining. What can you tell about the final piece of paper?

When we draw two pieces of paper, we either reduce the amount of odd numbers by two (if two odd numbers are drawn, resulting in an even number), or the amount of odd numbers stays the same. Thus we can deduce that the final piece of paper will be even if the starting amount of odd numbers is even, and odd if the starting amount of odd numbers is odd.

As an extension to this, I wondered what the probability distribution of the final piece of paper in the hat is, i.e. for $X$ being the random variable representing the final piece of paper in the hat, what is $P_N(X = k)$, for $k \in \{0, 1, ..., N\}$?

In order to find some sort of pattern, I wrote the following Python program to simulate the game for any N, and return an array of the amount of times each number was left in the hat.

# Import necessary dependencies
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm

def finalPiece(n):
    ''' 
    Pick pieces out of the hat randomly, replace with the absolute value
    of the difference and return the final number left.
    '''
    numberOfPieces = n
    piecesInHat = list(range(1, n+1))
    while numberOfPieces > 1:
        # Pick random piece of paper
        choice1Index = np.random.randint(0, numberOfPieces)
        choice2Index = np.random.randint(0, numberOfPieces-1)
        # Remove pieces of paper from hat
        choice1 = piecesInHat.pop(choice1Index)
        choice2 = piecesInHat.pop(choice2Index)
        # Replace with new number
        piecesInHat.append(abs(choice1-choice2))
        numberOfPieces = numberOfPieces - 1
    return piecesInHat[0]

def experiment(numbersInHat, numberOfTrials, plot=False, save=False):
    ''' 
    Repeat the finalPiece function and count how many times each number
    is left in the hat. Plot the result if plot == True. Save the results
    array if save == True.
    '''
    results = np.zeros(numbersInHat+1, dtype=int)
    # Count number of times each number is left in the hat, with progress bar
    for _ in tqdm(range(numberOfTrials)):
        results[finalPiece(numbersInHat)] += 1
    # Make a plot if it is desired
    if plot:
        x = np.linspace(0, numbersInHat, numbersInHat+1, dtype=int)
        plt.figure(figsize=(8, 6), dpi=800)
        plt.xlabel('Final Number in the Hat')
        plt.ylabel('Percentage of Experiments')
        plt.title('Hat Numbers Experiment: ' + str(numbersInHat) + ', ' + str(numberOfTrials))
        plt.bar(x, results*100/numberOfTrials)
        plt.savefig('bar graph ' + str(numbersInHat) + ' ' + str(numberOfTrials) + '.png')
        #plt.show()
    # Save results to file if it is desired
    if save:
        np.savetxt('counts ' + str(numbersInHat) + ' ' + str(numberOfTrials) +'.txt', results, fmt='%d')
    # Return results array (counts of experiments)
    return results

This shows the probability decreasing as $k$ increases (with $k$ of appropriate parity, and $k\neq 0$), but I still haven't been able to work out what the distribution actually is. Any assistance would be much appreciated.

Edit: To clarify, I am seeking an explicit formula for $P_N(X = k)$ if possible. Using the code above, I have already explored the distribution stochastically for large $N$.

Answer 1

The number of ways the game can go is $\prod_{i=2}^n \binom{i}{2}=\dfrac{n!(n-1)!}{2^{n-1}}$ (see https://oeis.org/A006472). I wrote my own code in Wolfram Language (Mathematica) to just go down every branch to find the exact answers for starting numbers of slips up to 9 (10 takes longer than I'd like to wait using my naive method:

del[list_, n_] := del[list, n] = DeleteCases[list, n, 1, 1]; 
delpair[list_, pair_] := 
 delpair[list, pair] = del[del[list, pair[[1]]], pair[[2]]]; 
play[x_] := 
 play[x] = 
  Flatten[If[Length[x] == 1, x, 
    Map[play[Append[delpair[x, Sort@#], Abs[#[[1]] - #[[2]]]]] &, 
     Subsets[x, {2}]]]]; denom[n_] := n!*(n - 1)!/2^(n - 1); 
numer[n_] := Counts[Sort[play[Range[n]]]]; 
dist[n_] := numer[n]/denom[n]; Do[Print[i, ": ", dist[i]], {i, 9}]

Try it online!

The results are as follows: $$\begin{matrix}0&1\\0&1\\\dfrac23&0&\dfrac13\\\dfrac49&0&\dfrac49&0&\dfrac19\\0&\dfrac{19}{30}&0&\dfrac{29}{90}&0&\dfrac{2}{45}\\0&\dfrac{269}{450}&0&\dfrac{212}{675}&0&\dfrac{119}{1350}\\\dfrac{1444}{4725}&0&\dfrac{5881}{14175}&0&\dfrac{88}{405}&0&\dfrac{14}{225}\\\dfrac{57073}{198450}&0&\dfrac{4232}{11025}&0&\dfrac{22111}{99225}&0&\dfrac{6131}{66150}&0&\dfrac{431}{33075}\\0&\dfrac{3323063}{7144200}&0&\dfrac{2134871}{7144200}&0&\dfrac{286901}{1786050}&0&\dfrac{156479}{2381400}&0&\dfrac{923}{95256}\end{matrix}$$

I checked the numerators of the unsimplified fractions in the OEIS and found they didn't all appear in the OEIS at all.

As requested, here are the means and variances (Mean[play[Range[i]]] and Variance[play[Range[i]]]):

Means: $1,1,\dfrac23\approx0.67,\dfrac43\approx1.33,\dfrac{82}{45}\approx1.82,\dfrac{1337}{675}\approx1.98,\dfrac{29374}{14175}\approx2.07,\dfrac{230143}{99225}\approx2.32,\dfrac{322913}{119070}\approx2.71$

Variances:$0,0,\dfrac43\approx1.33,\dfrac{32}{17}\approx1.88,\dfrac{10724}{8055}\approx1.33,\dfrac{3107024}{1821825}\approx1.71,\dfrac{2476997696}{803708325}\approx3.08,\dfrac{47158935632}{12117654675}\approx3.89,\dfrac{866608104176}{226842634431}\approx3.82$