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 <queue>
#include <set>
#include <map>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <complex>

using namespace std;

#define all(c) ((c).begin()),((c).end())
#define debug(c) cerr<<"> "<<#c<<" = "<<(c)<<endl;
#define iter(c) __typeof((c).begin())
#define tr(i,c) for(iter(c) i=(c).begin();i!=(c).end();i++)
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define REP(i,a,b) for(int i=(int)(a);i<=(int)(b);i++)

#define mp make_pair
#define fst first
#define snd second
#define pb push_back

typedef long long ll;
typedef vector<int> vi;
typedef pair<int,int> pii;

#define EPS 1.0e-6

//ツ点
typedef complex<double> P;

struct C{
	P p;
	double r;
};

vector<P> isCC(const C &c1, const C &c2) {
	double d = abs(c1.p - c2.p);
	if (c1.r + c2.r - d < -EPS || fabs(c1.r - c2.r) - d > EPS) return vector<P>();
	double x = (d*d + c1.r * c1.r - c2.r * c2.r) / 2 / d;
	double y = sqrt(c1.r * c1.r - x * x);
	P diff = (c2.p - c1.p) / d;
	vector<P> res;
	res.pb(c1.p + diff * P(x, y));
	if (y != 0) res.pb(c1.p + diff * P(x, -y));
	return res;
}


int main() {
	for (int N; cin >> N, N != 0; ) {
		vector<P> ps;
		rep (i, N) {
			double x, y; cin >> x >> y;
			ps.pb(P(x, y));
		}
		int res = 0;
		rep (i, N) REP (j, i + 1, N-1) {
			vector<P> isp = isCC((C){ps[i],1},(C){ps[j],1});
			rep (k, isp.size()) {
				int cnt = 0;
				rep (l, N) if (abs(isp[k] - ps[l]) <= 1 + EPS) cnt++;
				res = max(res, cnt);
			}
		}
		if (res == 0) res++;
		cout << res <<endl;
	}
	
	return 0;
}

C++: #include <vector>
#include <map>
#include <set>
#include <stack>
#include <queue>
#include <algorithm>
#include <utility>
#include <functional>
#include <sstream>
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <string>
#include <cstring>
#include <ctime>
#include <climits>
#include <fstream>
using namespace std;
inline int toInt(string s) { int v; istringstream sin(s); sin >> v; return v;}
template<class T> inline string toStr(T x) { ostringstream sout; sout << x; return sout.str();}
typedef vector<int> vi;
typedef vector<vi>  vvi;
typedef vector<string> vs;
typedef pair<int, int> pii;
typedef long long ll;
#define ALL(a) (a).begin(),(a).end()
#define RALL(a) (a).rbegin(),(a).rend()
#define FOR(i,a,b) for(int i=(a);i<=(b);++i)
#define REP(i,n) FOR(i,0,(n)-1)
const double EPS = 1e-12;
const double PI = acos(-1.0);
const int INF = INT_MAX/10;

#define EQ(x,y) (abs((x)-(y))<EPS)
#define LE(x,y) ((x)-(y) < +EPS)

#include <complex>
typedef complex<double> P;
typedef vector<P> vp;

P rotate(P t, P p, double r) {
	double ta = cos(r)*(t.real()-p.real())-sin(r)*(t.imag()-p.imag())+p.real();
	double tb = sin(r)*(t.real()-p.real())+cos(r)*(t.imag()-p.imag())+p.imag();
	return P(ta, tb);
}

double to_radian(double degree) {
	return degree*PI/180;
}

int main() {
	int N;
	while(cin >> N, N) {
		vp points(N);
		double X, Y;
		REP(i, N) {
			cin >> X >> Y;
			points[i] = P(X, Y);
		}

		int ans = 0;
		REP(i, N) {
			FOR(j, i+1, N-1) {
				P a = points[i], b = points[j];
				if(LE(abs(a-b), 2.0)) {
					P m = P((a.real()+b.real())/2, (a.imag()+b.imag())/2);
					double l = abs(a-m);
					double d = sqrt(1-l*l);
					double c = l-d;
					P t = P((d*a.real()+c*m.real())/(c+d), (d*a.imag()+c*m.imag())/(c+d));
					P p1 = rotate(t, m, to_radian(90));
					int tmp1 = 0;
					REP(k, N) {
						if(LE(abs(p1-points[k]), 1.0)) {
							tmp1++;
						}
					}
					ans = max(ans, tmp1);
					P p2 = rotate(t, m, to_radian(-90));
					int tmp2 = 0;
					REP(k, N) {
						if(LE(abs(p2-points[k]), 1.0)) {
							tmp2++;
						}
					}
					ans = max(ans, tmp2);
				}
			}
		}

		if(N <= 2 && ans == 0) {
			cout << 1 << endl;	
		} else {
			cout << ans << endl;
		}
	}

	return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;
using ld = long double;
using P = complex<ld>;
const ld eps = 1e-7;

int main()
{
	int N;
	while (cin >> N, N) {
		vector<ld> x(N), y(N);
		vector<P> p(N);
		for (int i = 0; i < N; i++) {
			cin >> x[i] >> y[i];
			p[i] = P(x[i], y[i]);
		}
		vector<P> ko;
		for (int i = 0; i < N; i++) {
			for (int j = i + 1; j < N; j++) {
				if (norm(p[i] - p[j]) < 4.0 + eps) {
					auto d = p[i] - p[j];
					ko.push_back(p[j] + (d * (ld)0.5) + d * P(0, 1) / abs(d) * sqrtl(1 - norm(d * (ld)0.5)));
					ko.push_back(p[j] + (d * (ld)0.5) - d * P(0, 1) / abs(d) * sqrtl(1 - norm(d * (ld)0.5)));
				}
			}
		}
		int res = 0;
		for (auto c : ko) {
			int cnt = 0;
			for (int i = 0; i < N; i++) {
				if (norm(c - p[i]) < 1.0 + eps) {
					cnt++;
				}
			}
			res = max(res, cnt);
		}
		cout << res << endl;
	}
	return 0;
}