Counting the socks¶
Introduction¶
Came across an amusing little probability problem the other day, from Fifty Challenging Problems in Probability with Solutions, by Frederick Mostell, with connections to some deep mathematics
The question is: given you have red and black socks in a drawer, how many of each do you need, to have a 50% chance of getting two red socks when picking at random from the drawer (sampling without replacement)r
import math
import numpy as np
import sympy
Analysis¶
Given we have $\require{cancel}$ $r$ red socks and $b$ black socks, with $n$ socks in total, the probability of getting two red socks is:
$$\frac{r * (r-1)} { n * (n-1)}$$
so we want integers $r$ and $n$ such that
$$\frac{r * (r-1)} { n * (n-1)} = \frac{1}{2}$$
Re-arranging terms, we get
$$ 2*r^2 -2*r = n*(n-1)$$
If $A*x^2+B*x+C = 0$, we have solutions $\frac{-B \pm \sqrt{B^2-4*A*C}}{2*A}$
Here, we have $A=2, B=2, C=-n*(n-1)$, so
$$ r = \frac{2 \pm \sqrt{4+4*2*n(n-1)}}{2*2}$$ or $$ r = \frac{2 \pm \sqrt{4+8*n(n-1)}}{2*2}$$
so for $r$ to be an integer, we must have $[1+2*n*(n-1)] $ to be a square
Mathematical connections¶
This is an example of a Diophantine equation:
a Diophantine equation is a polynomial equation with integer coefficients, for which only integer solutions are of interest.Diophantine equations have been studied for very long time: Diophantus
may have lived anywhere between 170 BCE, roughly contemporaneous with Hypsicles, the latest author he quotes from, and 350 CE, when Theon of Alexandria quotes from himPierre de Fermat famously stated his "Last Theorem" around 1637 in the margin of a copy of the ancient Greek text Arithmetica by Diophantus. This problem is also connected to Pell's equation: from OEIS-
... this is equivalent to the Pell equation A(n)^2-2*B(n)^2 = -1 Exploration¶
First I search out to find the first few candidate values, where $[1+2*n*(n-1)] $ is a a square
for i in range(2,5000):
x = math.sqrt(1+2*i*(i-1))
# test for integer solution
if( int(x) == x):
red = int((1+int(x))/2)
print(f'Sock count {i:4}, with {red:3} red socks')
#end if
#end for
Sock count 4, with 3 red socks Sock count 21, with 15 red socks Sock count 120, with 85 red socks Sock count 697, with 493 red socks Sock count 4060, with 2871 red socks
Checking, we get the probability of getting 2 red socks from 3 red socks out of 4 in the drawer $$ \frac{\cancel{3}}{4}*\frac{2}{\cancel{3}} = \frac{1}{2}$$
Checking, we get the probability of getting 2 red socks from $15 ( = 3*5)$ red socks out of $21 ( = 3*7)$ in the drawer
$$ \frac{\cancel{3}*\cancel{5}}{\cancel{3}*\cancel{7}}*\frac{2*\cancel{7}}{4*\cancel{5}} = \frac{1}{2}$$
Now it turns out that $4, 21, 120, 697, ...$ is a well known sequence, listed in The On-Line Encyclopedia of Integer Sequences [https://oeis.org/A046090], as a solution to this very problem.
Recurrance relations¶
The OEIS gives a recurrance relation for computing the next valid sock count (given the two previous sock count), but I thought I would find it for myself. Let us assume that any valid value depends only on the previous two values (a la Fibonacci) and a constant. If $a_{n}$ is the n-th value in the series:
$$ a_{n} = A*a_{n-1} + B*a_{n-2} + C$$
so we have
$$ \begin{align*}4060 &= A*697 + B*120 + C\\ 697 &= A*120 + B*21 + C\\ 120 &= A*21 + B*4 + C\\ \end{align*} $$
To solve for $A, B, C$, we use numpy linear algebra routines
b = [4060, 697, 120]
matrix = np.array([ [697,120,1], [120,21,1], [21,4,1] ])
sol = np.linalg.solve(matrix,b)
print(f'Recurrance relation: a(n) = {round(sol[0]):+}*a(n-1) {round(sol[1]):+}*a(n-2) {round(sol[2]):+}')
Recurrance relation: a(n) = +6*a(n-1) -1*a(n-2) -2
So to get the first few values of valid sock count values
# initial value in series
a_n_1 = 4
a_n = 21
values = [4,21]
# find next 30 values
for i in range(30):
# compute next value in series
x = 6*a_n - a_n_1 -2
# add it to out list
values.append(x)
# update recurrance values
a_n, a_n_1 = x,a_n
#end for
# display values nicely
for i,v in enumerate(values):
print(f' {i+1:2}: {v:40,}')
#end for
1: 4 2: 21 3: 120 4: 697 5: 4,060 6: 23,661 7: 137,904 8: 803,761 9: 4,684,660 10: 27,304,197 11: 159,140,520 12: 927,538,921 13: 5,406,093,004 14: 31,509,019,101 15: 183,648,021,600 16: 1,070,379,110,497 17: 6,238,626,641,380 18: 36,361,380,737,781 19: 211,929,657,785,304 20: 1,235,216,565,974,041 21: 7,199,369,738,058,940 22: 41,961,001,862,379,597 23: 244,566,641,436,218,640 24: 1,425,438,846,754,932,241 25: 8,308,066,439,093,374,804 26: 48,422,959,787,805,316,581 27: 282,229,692,287,738,524,680 28: 1,644,955,193,938,625,831,497 29: 9,587,501,471,344,016,464,300 30: 55,880,053,634,125,472,954,301 31: 325,692,820,333,408,821,261,504 32: 1,898,276,868,366,327,454,614,721
We can see that we get some very large integers very quickly! In manipulating these numbers (especially taking square roots), we move to a representation as powers of prime factors.
Checking a large value¶
We are dealing with very large integers, so it is more convenient to deal with numbers encoded as powers of their prime factors. We define some helper functions.
We use sympy to factor our numbers
def sqrt_factors(d):
'''
sqrt_factors: given a dictionary of prime factors of a number, return the prime factors of the sqrt
'''
d2 = {}
for k in d.keys():
d2[k] = round(d[k]/2)
#end for
return d2
#end sqrt_factors
def sqrt_value(d):
'''
sqrt_value: given a dictionary of prime factor of a number, return the sqrt value
'''
d2 = sqrt_factors(d)
value = 1
for k in d2.keys():
for p in range(d2[k]):
value = value * k
#end for
#end for
return value
#end sqrt_value
def div_by_4(n):
'''
div_by_4: given a number, return (number / 4 )
assumes number is divisible by 4, so reduce power of 2 in factor list by 2
'''
d = sympy.factorint(n)
d[2] = d[2]-2
# multiply out factors
value = 1
for k in d.keys():
value = value * k**d[k]
#end for
return value
#end div_by_4
def show_factors(factors):
'''
show_factors: give a dictionary of prime factors, return a string depiecting products of powers
'''
s = ""
first = True
for k in factors.keys():
# * only appears between factors
if(first):
s = s+ f'{k}^{factors[k]}'
first = False
else:
s = s+ f' * {k}^{factors[k]}'
#end if
#end for
return s
#end show_factors
def get_r_value(n):
'''
get_r_value: given the total sock count, find red sock count
'''
x = 4 + 8*(n*n-n)
y = sqrt_value(sympy.factorint(x))
r = div_by_4(2+y)
return r
#end get_r_value
def show_prob_calculation(v):
'''
show_prob_calculation: given a valid sock count, show the factors of (r/n)*((r-1)/(n-1))
'''
z = v
# get red sock count
r = get_r_value(z)
print ( f'Red sock count {r:,}, Total sock count {z:,}')
# show (r/n)*((r-1)/(n-1)) as product of powers of prime factores
print(
'(' + show_factors(sympy.factorint(r-1)) + ') * (' +
show_factors(sympy.factorint(r)) + ')'
)
print('-'*80)
print(
'(' + show_factors(sympy.factorint(z-1)) + ') * (' +
show_factors(sympy.factorint(z)) + ')'
)
print('\n\n')
#end for
#end show_prob_calculation
Now we can show how the factors of $(r/n)/((r-1)/(n-1))$ cancel out, leaving $\frac{1}{2}$, for our largest computed valid total sock count
show_prob_calculation(values[-1])
Red sock count 1,342,284,446,191,393,385,645,665, Total sock count 1,898,276,868,366,327,454,614,721 (2^5 * 3^1 * 7^1 * 17^1 * 23^1 * 353^1 * 577^1 * 665857^1 * 37667521^1) * (5^1 * 257^1 * 1409^1 * 5741^1 * 2448769^1 * 52734529^1) -------------------------------------------------------------------------------- (2^6 * 3^1 * 5^1 * 17^1 * 577^1 * 5741^1 * 665857^1 * 52734529^1) * (7^1 * 23^1 * 257^1 * 353^1 * 1409^1 * 2448769^1 * 37667521^1)
Conclusion¶
I wouldn't have guessed that such a simpe-sounding problem would have such deep mathematics roots
There is a saying (when you want to decline an invitation) "That's the day I'm re-organizing my sock drawer": based upon the above, this could indeed take some time!
Reproducability information¶
%load_ext watermark
%watermark
Last updated: 2026-08-20T19:55:35.266650+10:00 Python implementation: CPython Python version : 3.11.7 IPython version : 8.20.0 Compiler : MSC v.1916 64 bit (AMD64) OS : Windows Release : 10 Machine : AMD64 Processor : Intel64 Family 6 Model 170 Stepping 4, GenuineIntel CPU cores : 22 Architecture: 64bit
%watermark -h -iv -co
conda environment: base Hostname: INSPIRON16 numpy : 1.26.4 json : 2.0.9 pandas : 2.1.4 ipywidgets: 7.6.5 xarray : 2025.4.0 sympy : 1.12 sys : 3.11.7 | packaged by Anaconda, Inc. | (main, Dec 15 2023, 18:05:47) [MSC v.1916 64 bit (AMD64)]