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>

#define rep(n) for(int i=0;i<n;i++)
#define repp(j, n) for(int j=0;j<n;j++)
#define reppp(i, m, n) for(int i=m;i<n;i++)
#define all(c) c.begin(), c.end()
#define rall(c) c.rbegin(), c.rend()
#define debug(x) cerr << #x << ": " << x << endl

using namespace std;

typedef long long ll;
typedef pair<ll, ll> Pll;
typedef pair<int, int> Pii;

const long double EPS = 10e-10;

class Point{
    public: 
    long double x,y;
    
    Point(long double x = 0.0, long double y = 0.0): x(x), y(y) {}
    
    Point operator + (Point p) {
        return Point(x+p.x, y+p.y);
    }
    Point operator - (Point p) {
        return Point(x-p.x, y-p.y);
    }
    Point operator * (long double a) {
        return Point(x*a, y*a);
    }
    Point operator / (long double a) {
        return Point(x/a, y/a);
    }
    bool operator < (const Point &p) const {
        return (x != p.x ? x < p.x : y < p.y);
    }
    bool operator == (const Point &p) const {
        return fabsl(x - p.x) < EPS && fabsl(y - p.y) < EPS;
    }
    
    long double abs() {
        return sqrt(norm());
    }
    
    long double norm() {
        return x*x + y*y;
    }
    
    long double dot(Point p) {
        return x*p.x + y*p.y;
    }
    
    long double cross(Point p) {
        return x*p.y - y*p.x;
    }
};

int main() {
    int n;
    while(scanf("%d", &n) && n){
        vector<Point> points(n);
        rep(n) cin >> points[i].x >> points[i].y;

        int ans = 1;
        rep(n){ /* 円周上の点 */
            repp(j, n){ /* 直径上の点 */
                if(j == i) continue;
                Point p = points[j] - points[i];
                if(p.abs() > 2.0) continue;
                int tmp = 2;
                repp(k, n){ /* 円の中に入る？ */
                    if(k == i || k == j) continue;
                    Point q = points[k] - points[i];
                    long double cos_ = p.dot(q) / (p.abs() * q.abs());
                    // printf("%d %d %d\n", i, j, k);
                    // printf("%Lf %Lf\n", q.abs(), 2.0*cos_);
                    if(q.abs() <= 2.0*abs(cos_)) tmp++;
                }
                ans = max(ans, tmp);
            }
        }
        cout << ans << endl;
    }
}

C++: #include <bits/stdc++.h>
#include<stdio.h>
using namespace std;
#define FOR(i,k,n) for(int i = (k); i < (n); i++)
#define REP(i,n) FOR(i,0,n)
#define ALL(a) begin(a),end(a)
#define MS(m,v) memset(m,v,sizeof(m))
#define D10  fixed<<setprecision(10)
typedef vector<int> vi;
typedef vector<string> vs;
typedef pair<int, int> P;
typedef long long ll;
const ll INF = 114514810;
const ll MOD = 1000000007;
const double EPS = 1e-8;
const double PI = acos(-1.0);
struct edge
{
	int to, cost;
	bool operator < (const edge& e) const { return cost < e.cost; }
	bool operator >(const edge& e) const { return cost > e.cost; }
};
int dx[] = { -1, 0, 0, 1 }; int dy[] = { 0, -1, 1, 0 };
template<class T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template<class T> T &chmax(T &a, const T &b) { return a = max(a, b); }
bool valid(int x, int y, int h, int w) { return (x >= 0 && y >= 0 && x < h&&y < w); }
///*************************************************************************************///
///*************************************************************************************///
///*************************************************************************************///

bool eq(double a, double b) { return fabs(a - b) < EPS; }
typedef complex<double> Point;
typedef vector<Point> Polygon;

namespace std
{
	bool operator < (const Point& a, const Point& b)
	{
		return real(a) != real(b) ? real(a) < real(b) : imag(a) < imag(b);
	}
}

struct Line
{
	Point a, b;
	Line(Point p, Point q) :a(p), b(q) {};
};

struct Circle
{
	Point p; double r;
	Circle(Point a, double b) :p(a), r(b) {};
};

double dot(Point a, Point b)
{
	return real(conj(a)*b);
}

double cross(Point a, Point b)
{
	return imag(conj(a)*b);
}

int ccw(Point a, Point b, Point c)
{
	b -= a; c -= a;
	if (cross(b, c) > EPS) return 1;	//counter cloclwise
	if (cross(b, c) < -EPS) return -1;  //cloclwise
	if (dot(b, c) < 0) return 2;		//c--a--b on line 
	if (norm(b) < norm(c)) return -2;   //a--b--c on line
	return 0;							//a--c--b on line
}

vector<Point> is_cc(Circle c1, Circle c2)
{
	vector<Point> res;
	double d = abs(c1.p - c2.p);
	double rc = (d*d + pow(c1.r, 2) - pow(c2.r, 2)) / (2 * d);
	double dfr = pow(c1.r, 2) - rc*rc;
	if (abs(dfr) < EPS) dfr = 0;
	if (dfr < 0.0) return res;
	double rs = sqrt(dfr);
	Point diff = (c2.p - c1.p) / d;
	res.push_back(c1.p + diff*Point(rc, rs));
	if (dfr != 0.0) res.push_back(c1.p + diff*Point(rc, -rs));
	return res;
}




int main()
{
	int n;
	while (cin >> n, n)
	{
		int ans = 1;
		vector<Point> v;
		REP(i, n)
		{
			double x, y;
			cin >> x >> y;
			v.push_back(Point(x, y));
		}
		REP(i, n - 1)
		{
			FOR(j, 1 + 1, n)
			{
				Circle c1(v[i], 1), c2(v[j], 1);
				auto vp = is_cc(c1, c2);

				REP(k, vp.size())
				{
					int cnt = 0;
					Point tmp = vp[k];
					REP(l, n)
					{
						if (abs(tmp - v[l])-1<= EPS) cnt++;
					}
					chmax(ans, cnt);
				}
			}
		}
		cout << ans << endl;
	}
	return 0;
}