Hey everyone, Vitaliy suggested I post this. Explanation and code included!
The basic idea is twofold:
For most polynomials $P(x)$ and the series of polynomials $Q_ n(x)=x^n+P(x)$, if $n\to\infty$ then $Q_n(x)=0$ will have some solutions approaching the unit circle.
The structure of the roots of $Q_n(x)=x^n-1$ has the same structure as the set of rational numbers in $[0,1]$, since its solutions are $x=e^{2\pi i \frac{a}{n}}$.
So, you make some polynomial in x called $P(x,\theta)$ and animate over theta. Better yet, add a random polynomial of finite degree to start! Circle's sizes are decided by the leading coefficient of the polynomial, $n$.
This is another one I made without any randomness. I don't have the source code/specific polynomial that I used, but it looks like a cubic to me. Also it looks like I didn't set "AnimationRepetitions"->Infinity on this one, so you'll have to refresh if you don't see it animating. Keep an eye out for Euclid's orchard!
CODE
coef[n_, k_] := coef[n, k] = (RandomReal[{-1, 1}] + RandomReal[{-1, 1}] I);
nframes = 240;
solve[theta_, n_] := NSolve[
Sum[coef[n, k] z^k, {k, 1, n}] + 0.2 E^(I theta) +
0.9 E^(I 2 theta) z^3 If[n >= 3, 1, 0] +
0.9 E^(I 3 theta) z^6 If[n >= 6, 1, 0] == 0,
z];
disks[theta_, n_] := Disk @@@ Flatten[Table[
{{Re[z], Im[z]}, Max[0.3/n, 0.004]} /. solve[theta, n],
{n, 4, 120}], 1];
Export["algebraicAnimation.gif", Table[
Graphics[{Black, disks[theta, n]}, PlotRange -> 1.3],
{theta, 0, 2 Pi (1 - 1/nframes), 2 Pi/nframes}],
"DisplayDurations" -> 1/60, "AnimationRepetitions" -> Infinity]