# 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.