Centrifuge Problem¶
Question
Given a centrifuge with \(n\) holes, can we balance it with \(k\) (\(1\le k \le n\)) identical test tubes?
This is a simple yet interesting problem, very well illustrated by Numberphile and discussed by Matt Baker's blog.
The now proved solution is that:
Note
You can balance \(k\) identical test tubes, \(1\le k\le n\), in an \(n\)-hole centrifuge if and only if both \(k\) and \(n-k\) can be expressed as a sum of prime divisors of \(n\).
Below is my attempt to programmatically answer the centrifuge problem.
Method 1: Naïve DFS¶
The very first method literally follows the solution. For a given \((n,k)\) pair, check if \(k\) and \(n-k\) can be written as a linear combination of the prime divisors of \(n\) (with non-negative coefficients).
def is_linear_combination(x: int, prime_numbers: list) -> bool:
"""Check if `x` can be written as a linear combination of prime numbers, i.e.,
x = b1*p1 + b2*p2 + b3*p3 + ... + bn*pn
where pi represents a prime number in `prime_numbers`, bi is a non-negative integer.
"""
# very naive and not optimized
for n in prime_numbers:
# n divides x
if x % n:
return True
# n does not divides x, check if the difference between x and multiples of n can be
# a linear combination of other remaining prime numbers
for i in range(x//n):
if is_linear_combination(x - i*n, [p for p in prime_numbers if p!=n]):
return True
return False
def centrifuge_naive(n: int, k: int) -> bool:
"""Check if a `n`-hole centrifuge can be balanced with `k` identical test tubes.
True if both `k` and `n-k` can be written as a linear combination of the prime divisors of `n`.
"""
prime_divisors = get_prime_divisors(n) # simple cached function, skipped
return is_linear_combination(k, prime_divisors) and is_linear_combination(n-k, prime_divisors)
Some Optimizations¶
The above method works just fine, but very slow if we want to compute the total number of solutions, instead of just checking whether a particular \(k\) works.
There can be a few optimizations, for example, we can compute only the lower half of \(k\)s:
from functools import lru_cache
@lru_cache(maxsize=None)
def centrifuge_naive(n: int, k: int) -> bool:
prime_divisors = get_prime_divisors(n) # cached
if k > n//2:
return centrifuge(n, n-k)
return is_linear_combination(k, prime_divisors) and is_linear_combination(n-k, prime_divisors)
Further, if \(n\) is a (large) prime number itself, we understand that all \(1\le k\lt n\) will not work. Similarly, if \(n\) is a power of prime number, we can bypass many values of \(k\) too.
@lru_cache(maxsize=None)
def centrifuge_naive(n, k):
prime_divisors = get_prime_divisors(n)
# ...
# special case when n is power of prime
if len(prime_divisors) == 1:
p = prime_divisors[0]
return (k % p == 0) and ((n - k) % p == 0)
# ...
At certain point, we will realize that it would be faster to simply compute all possible \(k\)s instead of checking one by one whether a certain \(k\) can balance the centrifuge. This leads us to the second approach, which I call "bootstrap".
Method 2: Bootstrap¶
The bootstrap method is a variant of DFS, which essentially generates all possible \(k\) for a given \(n\) by exhausting the values from linear combinations of \(n\)'s prime divisors. The generated values should be between 2 and \(n\). Then we can tell if \(k'\) can balance the \(n\)-hole centrifuge by checking whether \(k'\) and \(n-k'\) are in the generated values.
def bootstrap(x, n, numbers, result):
"""Compute all linear combinations of the given numbers smaller than n"""
for p in numbers:
if p+x > n:
break
for i in range((n-x) // p):
p_ = x + p * i # p_ <= n
if not result[p_]:
# x + p*i has not been tested, and is a linear combination of given numbers
result[p_] = True
# check whether we can add multiples of remaining numbers
bootstrap(p_, n, [n2 for n2 in numbers if n2 != p], result)
def centrifuge_bootstrap(n: int, k: int) -> bool:
prime_divisors = get_prime_divisors(n) # cached, `prime_divisors` is sorted
# result[k] represents whether k is valid, k=0...n
result = [True] + [False] * (n-1) + [True]
bootstrap(0, n, prime_divisors, result) # TODO: bootstrap only once for a given `n`
return result[k] and result[n-k]
This method invests some time in pre-computing all possible linear combinations of the prime divisors of \(n\). If we are only interested to see a particular \((n,k)\) pair, we can break out when we have done result[k] and result[n-k] in bootstrap().
Method 3: Dynamic Programming¶
The last method uses dynamic programming. We can use \(f[k]\)=True to represent that \(k\) is a linear combination of \(n\)'s prime divisors. A value \(i\) is either itself a prime divisor of \(n\) (and thus a linear combination of the prime divisors), or the sum of a \(n\)'s prime divisor \(p\) and \((i-p)\). In the latter case, if \((i-p)\) is a linear combination of \(n\)'s prime divisors, so is \(p+(i-p)=i\).
Hint
If \((i-p)\) is a linear combination of \(n\)'s prime divisors, i.e., \(i-p=a_1p_1+a_2p_2+...+a_np_n\), where \(\{p_i\}\) are the prime divisors of \(n\) and \(\{a_i\}\) are non-negative integers, then \(i-p+p\) is definitely a linear combination too: \(p\)'s coefficient becomes \(a+1\ge0\).
Hence,
- \(f[i] = f[i] \text{ or } f[i-p]\)
The boundary condition is \(f[0]\) = True, i.e., an empty centrifuge is balanced.
The whole function is extremely short:
def centrifuge_dp(n: int, k: int) -> bool:
prime_divisors = get_prime_divisors(n) # cached, `prime_divisors` is sorted
f = [True] + [False] * n
for p in prime_divisors: # TODO: DP only once for a given `n`
for i in range(p, n+1):
f[i] = f[i] or f[i-p]
return f[k] and f[n-k]
Performance Comparison¶
Obviously, the Method 2 and 3 are much faster than the naïve Method 1. Method 3 does not even use recursion and is the fastest.
Note
A note there is that if we are to check all \(1\le k\le n\), e.g., [i for i in range(1, n+1) if centrifuge(n,i)], we need to make some adjustment to the functions above so as to bootstrap or perform DP only once for each \(n\). This is trivial.
Visualization¶
Below are some plots of balanced centrifuges. Note that for a particular value of \(k\), there can be more than one way to balance the centrifuge. Here, I illustrate only one.
Python code¶
The code to generate the plots above:
from functools import lru_cache
import numpy as np
import matplotlib.pyplot as plt
@lru_cache(maxsize=None)
def prime_divisors(n):
"""Return list of n's prime divisors"""
primes = []
p = 2
while p**2 <= n:
if n % p == 0:
primes.append(p)
n //= p
else:
p += 1 if p % 2 == 0 else 2
if n > 1:
primes.append(n)
return primes
def centrifuge(n):
"""Return a list of which the k-th element represents if k tubes can balance the n-hole centrifuge"""
F = [True] + [False] * n
for p in prime_divisors(n):
for i in range(p, n + 1):
F[i] = F[i] or F[i - p]
return [F[k] and F[n - k] for k in range(n + 1)]
def factorize(k: int, nums: list) -> list:
"""Given k, return the list of numbers from the given numbers which add up to k.
The given numbers are guaranteed to be able to generate k via a linear combination.
Examples:
>>> factorize(5, [2, 3])
[2, 3]
>>> factorize(6, [2, 3])
[2, 2, 2]
>>> factorize(7, [2, 3])
[2, 2, 3]
"""
def _factorize(k, nums, res: list):
for p in nums:
if k % p == 0:
res.extend([p] * (k // p))
return True
else:
for i in range(1, k // p):
if _factorize(k - p * i, [n for n in nums if n != p], res):
res.extend([p] * i)
return True
return False
res = []
_factorize(k, nums, res)
return res
@lru_cache(maxsize=None)
def centrifuge_k(n, k):
"""Given (n, k) and that k balances a n-hole centrifuge, find the positions of k tubes"""
if n == k:
return [True] * n
factors = factorize(k, prime_divisors(n))
pos = [False] * n
def c(factors: list, pos: list) -> bool:
if sum(pos) == k:
return True
if not factors:
return False
p = factors.pop(0)
pos_wanted = [n // p * i for i in range(p)]
for offset in range(n):
pos_rotated = [(i + offset) % n for i in pos_wanted]
# the intended positions of the p tubes are all available
if not any(pos[i] for i in pos_rotated):
# claim the positions
for i in pos_rotated:
pos[i] = True
if not c(factors, pos):
# unclaim the positions
for i in pos_rotated:
pos[i] = False
else:
return True
# all rotated positions failed, add p back to factors to place later
factors.append(p)
c(factors, pos)
return pos
def plot_centrifuge(n, figname="centrifuge.svg"):
ncols = max(int(n**0.5), 1) # minimum 1 column
nrows = n // ncols if n % ncols == 0 else n // ncols + 1
height = 3 if nrows == ncols else 2
width = 2
fig, axes = plt.subplots(nrows, ncols, figsize=(height * nrows, width * ncols))
z = np.exp(2 * np.pi * 1j / n)
theta = np.linspace(0, 2 * np.pi, 20)
radius = 1 / (ncols + nrows)
a = radius * np.cos(theta)
b = radius * np.sin(theta)
cent = centrifuge(n)
for nr in range(nrows):
for nc in range(ncols):
k = nr * ncols + nc + 1
axis = axes[nr, nc] if ncols > 1 else axes[nr]
if k > n:
axis.axis("off")
continue
# draw the n-holes
for i in [z**i for i in range(n)]:
axis.plot(a + i.real, b + i.imag, color="b" if cent[k] else "gray")
# draw the k tubes
if cent[k]:
if k > n // 2:
pos = [not b for b in centrifuge_k(n, n - k)]
else:
pos = centrifuge_k(n, k)
for i, ok in enumerate(pos):
i = z**i
if ok:
axis.fill(a + i.real, b + i.imag, color="r")
axis.set_aspect(1)
axis.set(xticklabels=[], yticklabels=[])
axis.set(xlabel=None)
axis.set_ylabel(f"k={k}", rotation=0, labelpad=10)
axis.tick_params(bottom=False, left=False)
fig.suptitle(f"$k$ Test Tubes to Balance a {n}-Hole Centrifuge")
fig.text(0.1, 0.05, "Red dot represents the position of test tubes.")
plt.savefig(figname)
plt.close(fig)
if __name__ == "__main__":
for n in range(6, 51):
print(f"Balancing {n}-hole centrifuge...")
plot_centrifuge(n, f"{n}-hole-centrifuge.png")
Download plots of balanced centrifuges¶
Success
You can download the Python code and all plots of balanced \(n\)-hole centrifuge, \(6\le n\le50\), which I calculated using the code above.