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 <queue>
#include <vector>
#include <algorithm>
#include <complex>
using namespace std;
#define rep(i,n) for(int i=0;i<n;i++)
static const double EPS = 1e-8;
typedef complex<double> P;
struct C {P p; double r;C(const P &p, double r) : p(p), r(r) { }};

int xor128(){ 
    static unsigned long x=123456789,y=362436069,z=521288629,w=88675123; 
    unsigned long t; 
    t=(x^(x<<11));x=y;y=z;z=w; return( w=(w^(w>>19))^(t^(t>>8)) ); 
} 

vector<P> data;

int check(double x,double y){
	int ret = 0;
	rep(i,data.size()){
		if( (data[i].real()-x) * (data[i].real()-x) + (data[i].imag()-y)*(data[i].imag()-y) <= 1.0+EPS){
			ret++;
		}
	}
	return ret;
}
int main(){
	int n;
	while(cin >> n , n){
		data.resize(n);
		rep(i,n){
			cin >> data[i].real() >> data[i].imag();
		}
		int ret = 0;
		rep(i,n){
			rep(x,300000/(n) ){
				double x = (xor128()%200000-100000)/100000.0;
				double y = (xor128()%200000-100000)/100000.0;
				ret = max( check(data[i].real()+x,data[i].imag()+y) , ret);
			}
		}
		cout << ret << endl;
	}
}

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

// xy平面上の点(ベクトル)を表現するには、complex型を利用するとよい
typedef complex<double> P;

// 辺の表現 (座標を2つ pair でもつ)
typedef pair<P, P> L;

// 円の表現 (座標 P と 半径 d で表現する)
typedef pair<P, double> C;

// 成分を取り出すのを簡単にする
#define X real()
#define Y imag()

// 誤差(epsilon)の定義
#define EPS (1e-10)

// 2つの要素が等しいかどうか
#define EQ(a,b) (abs((a) - (b)) < EPS)

// 2つのベクトルが等しいかどうか
#define EQV(a,b) ( EQ((a).X, (b).X) && EQ((a).Y, (b).Y) )

// m は n より大きい(以上)かどうか
#define LE(n, m) ((n) < (m) + EPS)
#define LEQ(n, m) ((n) <= (m) + EPS)

// m は n より小さい(以下)かどうか
#define GE(n, m) ((n) + EPS > (m))
#define GEQ(n, m) ((n) + EPS >= (m))

// 2つのベクトルの内積を求める
double dot(P a, P b) {
    return (a.X * b.X + a.Y * b.Y);
}

// 2つのベクトルの外積を求める
double cross(P a, P b) {
    return (a.X * b.Y - a.Y * b.X);
}

// 点 a と 点 b を通り、半径が r の円の中心を返す
vector<P> circlesPointsRadius(P a, P b, double r) {
    vector<P> cs;
    P abH = (b-a) * 0.5;
    double d = abs(abH);

    if(d == 0 || d > r) return cs; // 必要なら !LE(d,r) として円1つになる側へ丸める
    double dN = sqrt(r*r - d*d);   // 必要なら max(r*r - d*d, 0) とする
    P n = abH * P(0,1) * (dN / d);
    cs.push_back(a + abH + n);
    if(dN > 0) cs.push_back(a + abH - n);
    return cs;
}

int main() {
    int N;

    while(cin >> N, N) {
        vector<P> ps;
        for(int i=0; i<N; i++) {
            double x, y; cin >> x >> y;
            ps.push_back(P(x, y));
        }

        if(N == 1) cout << 1 << endl;
        else if(N == 2) {
            cout << (abs(ps[0] - ps[1]) < 2 + EPS) << endl;
        }
        else {
            int ans = 0;
            for(int i=0; i<N; i++) {
                for(int j=i+1; j<N; j++) {
                    // printf("i = %d, j = %d, abs = %f\n", i, j, abs(ps[i] - ps[j]));
                    if(abs(ps[i] - ps[j]) > 2 - EPS) continue;

                    vector<P> centers = circlesPointsRadius(ps[i], ps[j], 1);
                    for(auto c : centers) {
                        int sum = 2;
                        for(int k=0; k<N; k++) {
                            if(i == k || j == k) continue;
                            sum += (abs(c - ps[k]) < 1 + EPS);
                        }
                        ans = max(ans, sum);
                    }
                }
            }
            cout << ans << endl;
        }
    }
    return 0;
}

C++: #include <algorithm>
#include <cmath>
#include <climits>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <fstream>
#include <iostream>
#include <list>
#include <map>
#include <queue>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <vector>
#include <cassert>
#include <functional>

using namespace std;

#define LOG(...) printf(__VA_ARGS__)
//#define LOG(...)
#define FOR(i,a,b) for(int i=(int)(a);i<(int)(b);++i)
#define REP(i,n) for(int i=0;i<(int)(n);++i)
#define ALL(a) (a).begin(),(a).end()
#define RALL(a) (a).rbegin(),(a).rend()
#define EXIST(s,e) ((s).find(e)!=(s).end())
#define SORT(c) sort((c).begin(),(c).end())
#define RSORT(c) sort((c).rbegin(),(c).rend())
#define CLR(a) memset((a), 0 ,sizeof(a))

typedef long long ll;
typedef unsigned long long ull;
typedef vector<bool> vb;
typedef vector<int> vi;
typedef vector<ll> vll;
typedef vector<vb> vvb;
typedef vector<vi> vvi;
typedef vector<vll> vvll;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;

const int dx[] = { -1, 0, 1, 0 }; const int dy[] = { 0, 1, 0, -1 };

struct UnionFind {
	vector<int> v;
	UnionFind(int n) : v(n) { for (int i = 0; i < n; i++) v[i] = i; }
	int  find(int x) { return v[x] == x ? x : v[x] = find(v[x]); }
	void unite(int x, int y) { v[find(x)] = find(y); }
};
pair<double, double> nor(pair<double, double> a){
	double dist = hypot(a.first,a.second);
	return {a.first/dist,a.second/dist};
}
double dis(pair<double, double> a, pair<double, double> b) {
	return sqrt((a.first - b.first)*(a.first - b.first) + (a.second - b.second)*(a.second - b.second));
}

pair<double, double> sum(pair<double, double> a, pair<double, double> b) {
	return{ a.first + b.first,a.second + b.second };
}

pair<double, double> sub(pair<double, double> a, pair<double, double> b) {
	return{ a.first - b.first,a.second - b.second };
}

pair<double, double> mal(pair<double, double> a, double b) {
	return{ a.first *b,a.second *b };
}
int main() {
	int n;
	while (cin >> n, n) {
		vector<pair<double, double>> point(n);
		REP(i, n)
			cin >> point[i].first >> point[i].second;
		int ans = 1;
		REP(i, n)
			FOR(j, i + 1, n) {
			double d = dis(point[i], point[j]);
			if (d <2.0) {
				int cnt = 2;
				pair<double, double> ba = sub(point[i], point[j]);
				pair<double, double> vba = { -ba.second,ba.first };
				pair<double, double> cen1 = sum(point[j], sum(mal(ba, 0.5), mal(nor(vba), sqrt(1-d*d/4))));
				REP(k, n) {
					if (i == k || j == k)
						continue;
					if (dis(cen1, point[k]) <= 1.0)
						cnt++;
				}
				ans = max(ans, cnt);
				cnt = 2;
				pair<double, double> cen2 = sum(point[j], sum(mal(ba, 0.5), mal(nor(vba), -sqrt(1 - d*d / 4))));
				REP(k, n) {
					if (i == k || j == k)
						continue;
					if (dis(cen2, point[k]) <= 1.0)
						cnt++;
				}

				ans = max(ans, cnt);
			}
		}
		cout << ans << endl;
	}
}