Go: package main

import (
	"fmt"
	"math"
	"math/cmplx"
)

func main() {
	for {
		var n int
		fmt.Scan(&n)

		if n == 0 {
			return
		}

		arr := make([]complex128, n)

		for i := 0; i < n; i++ {
			var x, y float64

			fmt.Scanf("%f%f", &x, &y)

			arr[i] = complex(x, y)
		}

		max := 1
		for i := 0; i < n; i++ {
			for j := i + 1; j < n; j++ {
				sub := arr[j] - arr[i]

				if cmplx.Abs(sub) >= 2 {
					continue
				}
				sum := arr[j] + arr[i]

				v := 1
				{
					halfd := cmplx.Abs(sub / 2)

					center := sum/2 + sub*cmplx.Rect(math.Sqrt(1-halfd*halfd)/cmplx.Abs(sub), float64(v)*math.Pi/2)

					counter := 0
					for k := 0; k < n; k++ {
						if i == k || j == k || cmplx.Abs(arr[k]-center) <= 1 {
							counter++
						}
					}
					if max < counter {
						max = counter
					}
				}

				v = -1
				{
					halfd := cmplx.Abs(sub / 2)

					center := sum/2 + sub*cmplx.Rect(math.Sqrt(1-halfd*halfd)/cmplx.Abs(sub), float64(v)*math.Pi/2)

					counter := 0
					for k := 0; k < n; k++ {
						if i == k || j == k || cmplx.Abs(arr[k]-center) <= 1 {
							counter++
						}
					}
					if max < counter {
						max = counter
					}
				}

			}
		}

		fmt.Println(max)
	}
}

C++: #include <bits/stdc++.h>
using namespace std;
typedef complex<double> Point;
#define X real()
#define Y imag()
#define LE(n,m) ((n)<(m)+EPS)
#define GE(n,m) ((n)+EPS>(m))
#define EQ(n,m) (abs((n)-(m))<EPS)
const double EPS=1e-8;

//center of circle going through two given points of which radius equals to r
vector<Point> centergivenradius(const Point &a,const Point &b,const double &r){
    vector<Point> cs;
    Point ls=(b-a)*0.5;
    double dist=abs(ls);
    if (dist==0.0||dist>r) return cs;
    double len=sqrt(r*r-dist*dist);
    Point n=ls*Point(0.0,1.0)*(len/dist);
    cs.emplace_back(a+ls+n);
    if (GE(len,0.0)) cs.emplace_back(a+ls-n);
    return cs;
}

void solve(int N){
    vector<Point> P;
    for (int i=0;i<N;++i){
        double x,y; cin >> x >> y;
        P.emplace_back(x,y);
    }
    int ans=1;
    for (int i=0;i<N;++i){
        for (int j=i+1;j<N;++j){
            double dist=norm(P[i]-P[j]);
            if (GE(dist,4.0)) continue;
            vector<Point> center=centergivenradius(P[i],P[j],1.0);
            for (int l=0;l<2;++l){
                int cnt=2;
                for (int k=0;k<N;++k){
                    if (k==i||k==j) continue;
                    cnt+=(LE(norm(P[k]-center[l]),1.0));
                }
                ans=max(ans,cnt);
            }
        }
    }
    cout << ans << '\n';
}

int main(){
    cin.tie(0);
    ios::sync_with_stdio(false);
    int N;
    while(cin >> N,N){
        solve(N);
    }
}

C++: #include <iostream>
#include <complex>
#include <vector>
#include <stack>

using namespace std;

template<class T> bool operator<(const complex<T> &a, const complex<T> &b){
    return a.real() == b.real() ? a.imag() < b.imag() : a.real() < b.real();
}

typedef complex<double> P;
#define X real()
#define Y imag()

double dot(P a, P b){
    return (a * conj(b)).X;
}
double cross(P a, P b){
    return (a * conj(b)).Y;
}

vector<P> points;
deque<P> selected;
int max_count;

inline bool is_out_of_circle_for_three(P p1, P p2, P p3) {
    double d12 = abs(p1 - p2);
    P center12 = P((p1.X + p2.X) / 2, (p1.Y + p2.Y) / 2);
    double allow_distance = 1 + sqrt(1 - ((d12 * d12) / 4));
    return abs(center12 - p3)  > allow_distance + 1e-10;
}

void search(int index) {
    if(index > points.size()) return;
    int selected_size = selected.size();

    bool is_all_in_circle = true;
    if(selected_size >= 2 && selected_size > max_count) {
        // calc is all overlaps
        vector<P> copied_selected;
        copy(selected.begin(), selected.end(), back_inserter(copied_selected));

        for (int i = 0; i < selected_size ; ++i) {
            for (int j = 0; j < selected_size; ++j) {
                if(i == j) continue;
                // two point check
                double d12 = abs(copied_selected[i] - copied_selected[j]);
                if(d12 > 2.0 + 1e-10) {
                    is_all_in_circle = false;
                    break;
                }

                for (int k = 0; k < selected_size; ++k) {
                    if(i == k || j == k) continue;
                    // three point check
                    P p1 = copied_selected[i];
                    P p2 = copied_selected[j];
                    P p3 = copied_selected[k];
                    if(is_out_of_circle_for_three(p1, p2, p3)) {
                        is_all_in_circle = false;
                        break;
                    }
                }

            }
            if(!is_all_in_circle) {
                break;
            }
        }

        copied_selected.clear();
        copied_selected.shrink_to_fit();
    }

    if (is_all_in_circle) {
        // update max count
        if(selected_size > max_count) {
            max_count = selected_size;
        }

        P p = points[index];
        // remove bad situation
        if(selected.empty()) {
            // take
            selected.push_front(p);
            search(index + 1);
            selected.pop_front();
        } else if (selected.size() == 1 && abs(p - selected.front()) <= 2.0 + 1e-10) {
            // take
            selected.push_front(p);
            search(index + 1);
            selected.pop_front();
        } else if(selected.size() >= 2 && abs(p - selected.front()) <= 2.0 + 1e-10) {
            P pm1 = selected.front();
            selected.pop_front();
            P pm2 = selected.front();
            selected.push_front(pm1); // restore
            if(!is_out_of_circle_for_three(pm2, pm1, p)) {
                // take
                selected.push_front(p);
                search(index + 1);
                selected.pop_front();
            }
        }

        // for next dispatch
        // not take
        search(index + 1);
    }
}

int main() {
    int n;
    while (cin >> n) {
        if(n == 0) break;
        points.clear();
        points.shrink_to_fit();
        selected.clear();
        selected.shrink_to_fit();
        max_count = 0;

        // load points
        for (char i = 0; i < n; ++i) {
            double x;
            double y;
            cin >> x >> y;
            P p = P(x, y);
            points.push_back(p);
        }

        search(0);
        cout << max_count << endl;
    }
    return 0;
}