Hey everyone, Vitaliy suggested I post this. Explanation and code included!

The basic idea is twofold:

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

2. 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]