There is an [math]O(N^2 \log N)[/math] deterministic algorithm. This is usually (if not always) sufficient for competitive programming contests / online judges (e.g., POJ 1981 Circle and Points).
If there exists a circle C that encloses M points, then by slightly moving C it will touch at least two points. 
So, essentially we just need to check all circles which touch (at least) 2 points:
- For each pair of points in the given set, construct a circle with radius R that touches both points. For each point pair there are at most two such circles.
- For each constructed circle, check for each of the given points is inside it.
- Return the circle that encloses the maximum number of points.
The runtime of this algorithm is [math]O(N^3)[/math] because there are at most [math]O(N^2)[/math] such circles in step 1, and the linear scan at step 2 takes [math]O(N)[/math] time.
So, how can we achieve [math]O(N^2 \log N)[/math]?
The idea is as follows:
Pick an arbitrary point P from the given set. We rotate a circle C with radius R using P as the "axis of rotation" (by convention, in the counter-clockwise direction), i.e. we keep C to touch P any time during the rotation. While C is being rotated, we maintain a counter to count the number of points C encloses.
Note that this counter only changes when some point Q enters (or leaves) the region of the circle C. Our goal is to come up with an algorithm that will increment (or decrement) this counter, whenever some other point Q ≠ P enters (or leaves) the region of C.
The state of the (rotating) circle C can be described by a single parameter [math]\theta[/math], where [math](r, \theta)[/math] are the polar coordinates of the center of the circle C, if we pick P as the fixed point of the polar coordinate system. With this system, rotating C means increasing [math]\theta[/math].
(Figure from Wikipedia: Points in a polar coordinate system)
For each other point Q (≠ P), we can actually compute the range of [math]\theta[/math] for which C covers Q. Put more formally, C encloses Q whenever (iff) [math]\theta \in [\alpha, \beta][/math]. 
So, up to this point, the original problem has been reduced to:
What is an optimal value of [math]\theta[/math] that lies in the most number of [math][\alpha, \beta][/math] intervals?
The reduced problem can be solved with a pretty standard [math]O(N \log N)[/math] algorithm. This reduced problem has to be solved N times (one for each point P), hence the time complexity [math]O(N^2 \log N)[/math].
The above algorithms assume that a point is considered being enclosed in a circle even if it is on the circumference. If not, we could still tweak them a bit, for example, by changing the intervals of [math]\theta[/math] from being close-close to close-open. (Hmm, probably – didn't think too much to verify, but you get the idea.)
: Unless all pairwise distances are > 2R, where the answer to the maximum number of points will be 1. This can be checked in [math]O(N^2)[/math] time.
: We ignore Q's with distance > 2R from P.
: Probably need to be careful when handling these cyclic intervals.