I thought this was kinda fun, and I was trying to figure out why you couldn't just take a random x and use it to find a y that preserves x²+y²=r² (ie by setting y to a random number between 0 and ±√(r²-x²)). In the end I figured out that you need to generate x the same way, but if you have, say: n=rand(r), x=rand(±√(r²-n²)), y=rand(±√(r²-x²)) it seems to work.

Here's a version I extended to three dimensions, designed to be run in the JS console of the original article (it'll show up underneath the first sphere):

    new SphereSimulation(document.getElementById("spheres1"), () => {
      const rand = Math.random
      const root = (x) => rand() < 0.5 ? Math.sqrt(x) : -Math.sqrt(x)
      const invdist = (a=0, b=0) => root(RADIUS*RADIUS - a*a - b*b)

      const x1 = rand() * invdist()

      const x2 = rand() * invdist(x1)
      const y2 = rand() * invdist(x2)

      const x3 = rand() * invdist(x2, y2)
      const y3 = rand() * invdist(x3, y2)
      const z3 = rand() * invdist(x3, y3)

      return new THREE.Vector3(x3, y3, z3)
    })

It looks okay to me - am I missing something? Is this a reasonable approach?