[jkafjanvnfaf]

After thinking about the problem a bit, I've come up with a completely different approach. I'm sure this is known, but I can't find anything so I'd appreciate if someone can point me to a citation.

Generate two independent doubles a and b the "usual way", meaning from [0, 1) with equal spacings of 2^-53. Then calculate 2^-53*b + a. Assuming for simplicity that the rounding mode is towards zero, this gives a new random double in [0, 1). I've only done the math in my head so I'm not 100% sure on all the numbers, but I believe it scores pretty well at the criteria in the blog post; tentative results are

Probability of zero: 2^-106

Smallest nonzero: 2^-106

Largest non-seeable: approx. 2^-54

Distinct values: 2^58.7

The usual rounding mode (round to nearest, ties to even) will generate numbers in [0, 1] with very slight biases that need closer analysis. One could also go further and do e.g. 2^-106*c + 2^-53*b + a, but something like that should only be needed for for single-precision floats.

The obvious downside is that you always need two random integers and not just in the slow path. However, SIMD becomes easier because you don't need AVX512 for lzcnt.