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 <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <deque>
#include <queue>
#include <set>
#include <map>
#include <algorithm>
#include <functional>
#include <utility>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <cstdio>

using namespace std;

#define REP(i,n) for((i)=0;(i)<(int)(n);(i)++)
#define foreach(c,itr) for(__typeof((c).begin()) itr=(c).begin();itr!=(c).end();itr++)

struct point {double x,y;};

int N;
point P[310];

double dist(point A, point B){
	double dx = A.x - B.x, dy = A.y - B.y;
	return sqrt(dx*dx + dy*dy);
}

int func(void){
	int i,j,k,ans=1;
	
	REP(i,N) REP(j,N) if(i != j && dist(P[i],P[j]) < 2.0){
		double dx = (P[j].x - P[i].x) / 2.0, dy = (P[j].y - P[i].y) / 2.0;
		double d = sqrt(dx*dx + dy*dy);
		double t = sqrt(1.0 - d*d) / d;
		point O = {P[i].x + dx - t * dy, P[i].y + dy + t * dx};
		
		int tmp = 0;
		REP(k,N) if(dist(O,P[k]) < 1.00001) tmp++;
		ans = max(ans,tmp);
	}
	
	return ans;
}

int main(void){
	int i;
	
	while(1){
		cin >> N;
		if(N == 0) break;
		REP(i,N) cin >> P[i].x >> P[i].y;
		int ans = func();
		cout << ans << endl;
	}
	
	return 0;
}

C++: #include <bits/stdc++.h>
#define ll long long
#define INF 1000000005
#define MOD 1000000007
#define EPS 1e-10
#define rep(i,n) for(int i=0;i<n;++i)

using namespace std;

typedef pair<double,double>P;

const int MAX_N = 100005;

vector<P> a;
int n;

double dir(P p1,P p2)
{
    return (p1.first - p2.first) * (p1.first - p2.first) + (p1.second - p2.second) * (p1.second - p2.second);
}

vector<P> circle_center(P p1,P p2)
{
    vector<P> cs;
    P mid = P((p1.first + p2.first) / 2.0,(p1.second + p2.second) / 2.0);
    P vec = P(p1.first - p2.first,p1.second - p2.second);
    double len = dir(vec,P(0.0,0.0));
    double real_len = sqrt(len);
    double len2 = sqrt(1.0 - len / 4.0);
    cs.push_back(P(-vec.second * len2 / real_len + mid.first,vec.first * len2 / real_len + mid.second));
    cs.push_back(P(vec.second * len2 / real_len + mid.first,-vec.first * len2 / real_len + mid.second));
    return cs;
}

int main()
{
    while(1){
        scanf("%d",&n);
        if(n == 0){
            break;
        }
        a.clear();
        rep(i,n){
            double x,y;
            scanf("%lf%lf",&x,&y);
            a.push_back(P(x,y));
        }
        int res = -1;
        rep(i,n){
            for(int j = i+1; j < n; j++){
                if(dir(a[i],a[j]) < 4.0){
                    int cnt = 0;
                    vector<P> cs = circle_center(a[i],a[j]);
                    rep(k,n){
                        if(k != i && k != j && dir(a[k],cs[0]) <= 1.0){
                            cnt++;
                        }
                    }
                    res = max(res,cnt);
                    cnt = 0;
                    rep(k,n){
                        if(k != i && k != j && dir(a[k],cs[1]) <= 1.0){
                            cnt++;
                        }
                    }
                    res = max(res,cnt);
                }
            }
        }
        printf("%d\n",res+2);
    }
	return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;
#define rep(i,n) for(int i=0;i<n;++i)
struct Point{double x,y;};
double dist(Point a,Point b) {return sqrt(pow(a.x-b.x,2)+pow(a.y-b.y,2));}
#define eps 0.0001

int main(void){
    int N;
    while(cin>>N,N){
        vector<Point> p(N);
        int maxcnt=0;
        rep(i,N)cin>>p[i].x>>p[i].y;
        rep(i,N){
            for(int j=i+1;j<N;j++){
                double d = dist(p[i],p[j]);
                if(d>2.0)continue;
                Point C1,C2;
                C1.x = p[i].x + cos(atan2(p[j].y-p[i].y,p[j].x-p[i].x) + acos(d/2.0));
                C1.y = p[i].y + sin(atan2(p[j].y-p[i].y,p[j].x-p[i].x) + acos(d/2.0));
                C2.x = p[i].x + cos(atan2(p[j].y-p[i].y,p[j].x-p[i].x) - acos(d/2.0));
                C2.y = p[i].y + sin(atan2(p[j].y-p[i].y,p[j].x-p[i].x) - acos(d/2.0));
                int cnt=0;
                rep(k,N){
                    if(pow(C1.x-p[k].x,2) + pow(C1.y-p[k].y,2) <= 1.0+eps)cnt++;
                }
                if(maxcnt<cnt)maxcnt=cnt;
                cnt=0;
                rep(k,N){
                    if(pow(C2.x-p[k].x,2) + pow(C2.y-p[k].y,2) <= 1.0+eps)cnt++;
                }
                if(maxcnt<cnt)maxcnt=cnt;
            }
        }
        cout<<maxcnt<<endl;
    }

    return 0;
}

C++: #include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;

constexpr double EPS = 1e-6;

inline double sq(double x) { return x*x; }

struct point {
    double x, y;
    point() {}
    point(double x, double y): x(x), y(y) {}

    point operator + (const point& o) const { return point(x+o.x, y+o.y); }
    point operator - (const point& o) const { return point(x-o.x, y-o.y); }
    point operator * (const double a) const { return point(x*a, y*a); }
    point operator / (const double a) const { return point(x/a, y/a); }
};

double fabs(point a) { return sqrt(a.x*a.x + a.y*a.y); }
double dist(point a, point b) { return fabs(a - b); }

int N;
point P[300];

int solve() {
    int ans = 1;
    for (int i = 0; i < N; ++i) {
        for (int j = i+1; j < N; ++j) {
            if (dist(P[i], P[j]) > 2) continue;

            point d = (P[j] - P[i]) / dist(P[i], P[j]);
            point t = (P[i] + P[j]) / 2;
            double k = sqrt(1 - sq(dist(P[i], P[j]) / 2));
            point cs[2] = {t + point(d.y, -d.x) * k, t + point(-d.y, d.x) * k};

            for (point c : cs) {
                int count = 0;
                for (int i = 0; i < N; ++i) {
                    if (dist(c, P[i]) <= 1 + EPS) ++count;
                }
                ans = max(ans, count);
            }
        }
    }

    return ans;
}

int main(void) {
    while (true) {
        cin >> N; if (N == 0) break;
        for (int i = 0; i < N; ++i) cin >> P[i].x >> P[i].y;
        cout << solve() << endl;
    }
    return 0;
}