# Centrifuge Problem

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:

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`.
"""
= get_prime_divisors(n) # simple cached function, skipped
prime_divisors 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:
= get_prime_divisors(n) # cached
prime_divisors 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):
= get_prime_divisors(n)
prime_divisors # ...
# special case when n is power of prime
if len(prime_divisors) == 1:
= prime_divisors[0]
p 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):
= x + p * i # p_ <= n
p_ if not result[p_]:
# x + p*i has not been tested, and is a linear combination of given numbers
= True
result[p_] # check whether we can add multiples of remaining numbers
for n2 in numbers if n2 != p], result)
bootstrap(p_, n, [n2
def centrifuge_bootstrap(n: int, k: int) -> bool:
= get_prime_divisors(n) # cached, `prime_divisors` is sorted
prime_divisors # result[k] represents whether k is valid, k=0...n
= [True] + [False] * (n-1) + [True]
result 0, n, prime_divisors, result) # TODO: bootstrap only once for a given `n`
bootstrap(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:
= get_prime_divisors(n) # cached, `prime_divisors` is sorted
prime_divisors = [True] + [False] * n
f for p in prime_divisors: # TODO: DP only once for a given `n`
for i in range(p, n+1):
= f[i] or f[i-p]
f[i] 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.

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.

`6, "6-hole-centrifuge.svg") plot_centrifuge(`

`10, "10-hole-centrifuge.svg") plot_centrifuge(`

`12, "12-hole-centrifuge.svg") plot_centrifuge(`

`18, "18-hole-centrifuge.svg") plot_centrifuge(`

`20, "20-hole-centrifuge.svg") plot_centrifuge(`

`24, "24-hole-centrifuge.svg") plot_centrifuge(`

`33, "33-hole-centrifuge.svg") plot_centrifuge(`

### 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 = 2
p while p**2 <= n:
if n % p == 0:
primes.append(p)//= p
n else:
+= 1 if p % 2 == 0 else 2
p 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"""
= [True] + [False] * n
F for p in prime_divisors(n):
for i in range(p, n + 1):
= F[i] or F[i - p]
F[i] 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:
* (k // p))
res.extend([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):
* i)
res.extend([p] 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
= factorize(k, prime_divisors(n))
factors = [False] * n
pos
def c(factors: list, pos: list) -> bool:
if sum(pos) == k:
return True
if not factors:
return False
= factors.pop(0)
p = [n // p * i for i in range(p)]
pos_wanted for offset in range(n):
= [(i + offset) % n for i in pos_wanted]
pos_rotated # 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:
= True
pos[i] if not c(factors, pos):
# unclaim the positions
for i in pos_rotated:
= False
pos[i] 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"):
= max(int(n**0.5), 1) # minimum 1 column
ncols = n // ncols if n % ncols == 0 else n // ncols + 1
nrows = 3 if nrows == ncols else 2
height = 2
width = plt.subplots(nrows, ncols, figsize=(height * nrows, width * ncols))
fig, axes = np.exp(2 * np.pi * 1j / n)
z
= np.linspace(0, 2 * np.pi, 20)
theta = 1 / (ncols + nrows)
radius = radius * np.cos(theta)
a = radius * np.sin(theta)
b
= centrifuge(n)
cent for nr in range(nrows):
for nc in range(ncols):
= nr * ncols + nc + 1
k = axes[nr, nc] if ncols > 1 else axes[nr]
axis if k > n:
"off")
axis.axis(continue
# draw the n-holes
for i in [z**i for i in range(n)]:
+ i.real, b + i.imag, color="b" if cent[k] else "gray")
axis.plot(a # draw the k tubes
if cent[k]:
if k > n // 2:
= [not b for b in centrifuge_k(n, n - k)]
pos else:
= centrifuge_k(n, k)
pos for i, ok in enumerate(pos):
= z**i
i if ok:
+ i.real, b + i.imag, color="r")
axis.fill(a
1)
axis.set_aspect(set(xticklabels=[], yticklabels=[])
axis.set(xlabel=None)
axis.f"k={k}", rotation=0, labelpad=10)
axis.set_ylabel(=False, left=False)
axis.tick_params(bottom
f"$k$ Test Tubes to Balance a {n}-Hole Centrifuge")
fig.suptitle(0.1, 0.05, "Red dot represents the position of test tubes.")
fig.text(
plt.savefig(figname)
plt.close(fig)
if __name__ == "__main__":
for n in range(6, 51):
print(f"Balancing {n}-hole centrifuge...")
f"{n}-hole-centrifuge.png") plot_centrifuge(n,
```

### Download plots of balanced centrifuges

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.