Why Using One Fixed Projection Will Not Work
For a projection to be of any use at all for distance-related queries, it and its inverse should be continuous almost everywhere. Consider, then, what happens when you pick one point--any point--and start to draw on the map a collection of routes that emanate from that point and move straight away from it, in all directions, on the earth. Initially these routes will fill out a local neighborhood of the point, because every location in that local neighborhood can (obviously) be reached by moving straight towards it. At any given distance, the set of locations that can be reached will have to fill out some kind of a blob on the map centered at the image of the starting point. As the distance grows, this blob has to continue expanding all along its perimeter. After a very long distance, though, these routes will all converge on a small collection of points that are diametrically opposite the starting point on the earth. But now there's no way they can converge on the map, because they have to form the perimeter of a (now very large) blob. Therefore there has to be some location where the projection creates an enormous distortion of distances. Any Euclidean map-based algorithm is guaranteed to fail when it is applied near any such location.
What to do
There are two problems you have to manage:
Changes in the scale of the projection from point to point on the earth.
Changes in the scale of the projection according to the direction relative to a fixed point on the earth.
Any projection where the scale changes with direction is called "non-conformal." When you use such a projection, a circle on the earth--even a very small one--will be non-circular on the map. That creates subtly bad results for your algorithm, because it will then have a tendency to select more points in some directions and fewer points in other directions. Therefore, although departures from conformality may be acceptable, they must be kept within narrow bounds.
Considerations like this lead to several possible strategies:
Use a collection of well-tuned projections blanketing the Earth. Each projection is allowed to extend to a distance--which in some cases can be around a thousand kilometers or so (I estimate)--over which its scale varies by an acceptably small amount. A standard such set of projections is UTM, the Universal Transverse Mercator. (Despite the name, it includes some Stereographic projections near the Poles.)
The price you pay in using this strategy is an initial step in which you identify an appropriate projection for the probe point and its surroundings and then reproject all target points. Then you apply your 2D Euclidean algorithm.
Use a very small number of conformal projections--two will do (such as two Polar Stereographic projections) or even just one in some cases--and modify the search algorithm to allow for the search radius to vary according to the location of the probe point.
This is easy to carry out--scale formulas for Mercator and Stereographic projections are simple (see the resources below)--but it works only for short distances, typically less than a few thousand kilometers, depending on the accuracy needed.
Perform the search in 3D. Use earth-centered X,Y,Z coordinates. Convert the distance on the spheroid to an (approximate) 3D distance (which will be a little shorter).
This might be simplest among the highly accurate solutions. There's a one-time precomputation cost associated with converting all point coordinates to 3D coordinates. The real price is that you need to implement a 3D Euclidean algorithm instead of a 2D one.
Use a crude approximation. For 50 m accuracy within 1000 km, simply rescaling the longitudes locally will work fine, as described at https://gis.stackexchange.com/a/108016. This is really an approximation to (1), but carries a smaller computational burden.
Resources
Scale factors for UTM distance distortions are given at Calculating areal distortion outside UTM zone?.
A discussion of these issues for the purpose of distance-weighted interpolation, along with some recommendations, appears at Should point data be equidistant projected when using ArcGIS IDW spatial interpolation?.
Mercator scale factors are given at Using SRTM Global DEM for Slope calculation?. If you don't need to search close to either Pole, this single projection will work with strategy (2) (varying the search radius according to the probe locations).
±delta
of the dividing line - wheredelta
is a distance that you'll need to determine by analysis of the error in the fast calculation.