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;
#define int long long   // <-----!!!!!!!!!!!!!!!!!!!

#define rep(i,n) for (int i=0;i<(n);i++)
#define rep2(i,a,b) for (int i=(a);i<(b);i++)
#define rrep(i,n) for (int i=(n)-1;i>=0;i--)
#define rrep2(i,a,b) for (int i=(a)-1;i>=b;i--)
#define all(a) (a).begin(),(a).end()
#define rall(a) (a).rbegin(),(a).rend()
#define printV(_v) for(auto _x:_v){cout<<_x<<" ";}cout<<endl
#define printVS(vs) for(auto x : vs){cout << x << endl;}
#define printVV(_vv) for(auto _v:_vv){for(auto _x:_v){cout<<_x<<" ";}cout<<endl;}
#define printP(p) cout << p.first << " " << p.second << endl
#define printVP(vp) for(auto p : vp) printP(p);

typedef long long ll;
typedef pair<int, int> Pii;
typedef tuple<int, int, int> TUPLE;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef vector<vvi> vvvi;
typedef vector<Pii> vp;
const int inf = 1e9;
const int mod = 1e9 + 7;

// #define double long double

typedef complex<double> P;
typedef vector<P> G;
#define here(g, i) g[i]
#define next(g, i) g[(i + 1) % g.size()]
#define prev(g, i) g[(i - 1 + g.size()) % g.size()]
const double EPS = 1e-6;
const double INF = 1e12;
const double PI = acos(-1);

P readP() {
    double x, y;
    cin >> x >> y;
    return P(x, y);
}

// rot p around q by theta (counter-clockwise)
P rotP(P p, P q, double theta) {
    p -= q;
    double x = p.real(), y = p.imag();
    p = P(x * cos(theta) - y * sin(theta), x * sin(theta) + y * cos(theta));
    p += q;
    return p;
}

P center(P p, P q) {
    // P v = q - p;
    // double theta = asin(abs(v) / 2);
    // v *= (1 / abs(v));
    // P u = rotP(v, p, theta);
    // return p + u;

    // P v = q - p;
    // P h = v * P(0., 1.) / abs(v);
    // return p + (v / 2.) + h;

    P v = q - p;
    P h = P(v.imag(), -v.real());
    h /= abs(h);
    h *= sqrt(1 - norm(v / 2.));
    P c = p + (v / 2.) + h;
    return c;
}

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(0);

    int n;
    while (cin >> n, n) {
        vector<P> p(n);
        rep(i, n) p[i] = readP();

        int ans = 1;
        rep(i, n) {
            rep2(j, i + 1, n) {
                if (i == j) continue;
                if (norm(p[i] - p[j]) > 4. + EPS) continue;
                rep(s, 2) {
                    P c = (s ? center(p[i], p[j]) : center(p[j], p[i]));
                    int cnt = 2;
                    rep(k, n) {
                        if (k == i || k == j) continue;
                        if (norm(p[k] - c) < 1. + EPS) {
                            cnt++;
                        }
                    }
                    ans = max(ans, cnt);
                }
            }
        }

        cout << ans << endl;
    }
}

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

#define REP(i, a, b) for (int i = (a); i < (b); i++)
#define rep(i, n) REP(i, 0, n)

typedef complex<double> P;

#define EPS (1e-10)
#define EQ(a, b) (abs((a)-(b)) < EPS)
#define EQV(a, b) (EQ((a).real(), (b).real) && EQ((a).imag(), (b).imag()))

P unit_vector(P a) { return a/abs(a); }
P normal_vector1(P a) { return a * P(0, 1); }
P normal_vector2(P a) { return a * P(0, -1); }
double distance_p_p(P a, P b) { return abs(a-b); }

int N;

int main()
{
	while (cin >> N, N) {
		vector<P> ps;

		rep(i, N) {
			double x, y;
			cin >> x >> y;
			ps.push_back(P(x, y));
		}

		int ans = 1;
		rep(i, N) {
			REP(j, i+1, N) {
				P p, q;
				p = ps[i]; q = ps[j];

				double d = distance_p_p(p, q);
				if (d > 2) continue;

				P pq_c = (q - p)/2.0;
				P a = unit_vector(normal_vector1(q-p));
				P b = unit_vector(normal_vector2(q-p));

				P c1 = (p + pq_c) + a*(sqrt(1.0-d*d/4.0));
				P c2 = (p + pq_c) + b*(sqrt(1.0-d*d/4.0));

				int buf1, buf2;
				buf1 = buf2 = 2;
				rep(k, N) {
					if (k == i || k == j) continue;
					if (distance_p_p(ps[k], c1) <= 1)
						buf1++;
					if (distance_p_p(ps[k], c2) <= 1)
						buf2++;
				}
				ans = max(ans, max(buf1, buf2));
			}
		}
		cout << ans << endl;
	}
	return 0;
}

C++: #define _USE_MATH_DEFINES
#include <algorithm>
#include <cstdio>
#include <functional>
#include <iostream>
#include <cfloat>
#include <climits>
#include <cstring>
#include <cmath>
#include <map>
#include <queue>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <time.h>
#include <vector>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> i_i;
typedef pair<ll, int> ll_i;
typedef pair<double, int> d_i;
typedef pair<ll, ll> ll_ll;
typedef pair<double, double> d_d;
struct edge { int u, v; ll w; };

ll MOD = 1000000007;
ll _MOD = 1000000009;
double EPS = 1e-4;

int main() {
	for (;;) {
		int N; cin >> N;
		if (N == 0) break;
		vector<double> x(N), y(N);
		for (int i = 0; i < N; i++)
			cin >> x[i] >> y[i];
		int maxi = 1;
		for (int i = 0; i < N; i++)
			for (int j = 0; j < N; j++) {
				if (i == j) continue;
				double dx = x[j] - x[i], dy = y[j] - y[i];
				double d = sqrt(dx * dx + dy * dy);
				if (d > 2) continue;
				double h = sqrt(1 - d * d / 4);
				double ox = x[i] + dx / 2 - dy / d * h;
				double oy = y[i] + dy / 2 + dx / d * h;
				int cnt = 0;
				for (int k = 0; k < N; k++) {
					double dx = x[k] - ox, dy = y[k] - oy;
					double d2 = dx * dx + dy * dy;
					if (d <= 1 + EPS) cnt++;
				}
				maxi = max(maxi, cnt);
			}
		cout << maxi << endl;
	}
}