Epsilon to compare squared distances must be exactly 1.0 / (1 << 24) (at least for M*N algo). Even 23 and 25 leads to WA.

(std::numeric_limits< float >::epsilon() / 2) is absolute accuracy needed overall. You can invent twofold float-then-double algorithm to sieve bad points beforehand. For single precision algorithm part you have to use 20.1f (20 + small constant) "epsilon" to compare squared distances in case if you use <= or >= operator (say, vcmpge_oqps or vcmple_oqps instructions) and 28.2f (28 + small constant) in case of strict inequality. Surely you can invent adaptive algorithm to infer relative accuracy needed in particular test case, which takes into account max abs differences of input point coordinates.

epsilon should be explicitly specified in the problem statement, it is not hard to make a test to break even simplest O(n*m) double precision solution.

I totally agree.
The problem should be reformulated in integers or other way to make use of arbitrary precision numbers.
Also 10000*100000 is too small dimensionality to encourage participants to make submissions of O((N + M) * log(N + M)) solutions due to high constant factor of latter. Stupid algorithm with randomization and trivial vectorization is faster (and extremely easier to implement) then clever algorithm with Voronoi and point location.