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

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 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:

 1 2 3 4 5 6 7 8 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.

 1 2 3 4 5 6 7 8 9 @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 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.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 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$ = True, i.e., an empty centrifuge is balanced.

The whole function is extremely short:

 1 2 3 4 5 6 7 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.

 1 plot_centrifuge(6, "6-hole-centrifuge.svg") 1 plot_centrifuge(10, "10-hole-centrifuge.svg") 1 plot_centrifuge(12, "12-hole-centrifuge.svg") 1 plot_centrifuge(18, "18-hole-centrifuge.svg") 1 plot_centrifuge(20, "20-hole-centrifuge.svg") 1 plot_centrifuge(24, "24-hole-centrifuge.svg") 1 plot_centrifuge(33, "33-hole-centrifuge.svg") ### Python code¶

The code to generate the plots above:

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 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") 

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.