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 loop(n, i) for(int i=0;i<n;i++)
#define all(v) v.begin(),v.end()

using namespace std;
typedef pair<double, double> P;
const double EPS = 1e-6;

#include <bits/stdc++.h>
 
#define loop(n, i) for(int i=0;i<n;i++)
#define all(v) v.begin(),v.end()

using namespace std;

double dist(P a, P b)
{
    return hypot(a.first - b.first, a.second - b.second);
}

vector<P> center(P a, P b)
{
    auto t = atan2(b.first - a.first, b.second - a.second);
    auto x = dist(a, b) / 2;
    auto y = sqrt(1.0 - x * x);
    return { 
        { cos(t) * x - sin(t) * y + a.first, sin(t) * x + cos(t) * y + a.second },
        { cos(t) * x + sin(t) * y + a.first, sin(t) * x - cos(t) * y  + a.second },
    };
}

int main()
{
    while (1) {
        int n; cin >> n;
        if (!n) break;
        vector<P> pnt;
        loop (n, i) {
            double x, y; cin >> x >> y;
            pnt.emplace_back(x, y);
        }
        int ans = 1;
        loop (n, i) loop (n, j) {
            if (dist(pnt[i], pnt[j]) < 2) ans = 2;
        }
        loop (n, i) loop (n, j) {
            auto ps = center(pnt[i], pnt[j]);
            for (auto o : ps) {
                int cnt = 0;
                loop (n, l) {
                    if (dist(o, pnt[l]) < 1 + EPS) cnt++;
                }
                ans = max(ans, cnt);
            }
        }
        cout << ans << endl;
    }
    return 0;
}

C++: #include <bits/stdc++.h>

using namespace std;

typedef pair<double, double> pii;
#define rep(i, n) REP(i, 0, (n))
#define REP(i, a, n) for(int i=(a); i<(n); i++)

struct Vec2D {
    double x, y;
    Vec2D() : x(0.0), y(0.0) {}
    Vec2D(double x, double y) : x(x), y(y) {}
    Vec2D operator=(const Vec2D& v) { x=v.x; y=v.y; return *this; }
    Vec2D operator-() { return Vec2D(-x, -y); }
    Vec2D operator+(const Vec2D& v) { return Vec2D(x+v.x, y+v.y); }
    Vec2D operator-(const Vec2D& v) { return Vec2D(x-v.x, y-v.y); }
    double operator*(const Vec2D& v) { return x*v.x + y*v.y; }
    Vec2D operator*(double k) { return Vec2D(x*k, y*k); }
    Vec2D operator/(double k) { return Vec2D(x/k, y/k); }
    double norm() { return sqrt(x*x + y*y); }
    Vec2D normalize() { return Vec2D(x, y)/norm(); }
    Vec2D normal() { return Vec2D(y, -x).normalize(); }
};
ostream& operator<<(ostream& os, const Vec2D& v) {
    os << "Vec2D [ x: " << v.x << ", y: " << v.y << " ]";
    return os;
}
istream& operator>>(std::istream& is, Vec2D& v) {
    is >> v.x >> v.y;
    return is;
}

int N;
Vec2D ps[300];

int main() {
    while(cin>>N, N) {
        rep(i, N) cin >> ps[i];
        
        int ans = 1;
        rep(i, N-1) {
            REP(j, i+1, N) {
                Vec2D d = ps[i]-ps[j];
                if(d.norm() > 2.0001) continue;
                double x = sqrt(1 - pow(d.norm()/2.0, 2));
                Vec2D n = d.normal(), m = (ps[i]+ps[j])/2.0;
                Vec2D o;
                int cnt_p = 0;
                o = m + (n*x);
                rep(k, N) if((ps[k]-o).norm()<1.0001) cnt_p++;
                int cnt_n = 0;
                o = m - (n*x);
                rep(k, N) if((ps[k]-o).norm()<1.0001) cnt_n++;
                ans = max(ans, max(cnt_p, cnt_n));
            }
        }
        cout << ans << endl;
    }
}

C++: #include <bits/stdc++.h>
#define PB push_back
#define MP make_pair
#define REP(i,n) for (int i=0;i<(n);i++)
#define FOR(i,a,b) for(int i=(a);i<(b);i++)
#define ALL(a) (a).begin(),(a).end()
int INF=1e9;
int MOD=1000000007;
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
double EPS = 1e-10;
double add(double a,double b){
	if(abs(a+b)<EPS*(abs(a)+abs(b)))return 0;
	return a+b;
}
struct P{
	double x,y;
	P(){}
	P(double x,double y):x(x),y(y){}
	P operator + (P p){
		return P(add(x,p.x),add(y,p.y));
	}
	P operator - (P p){
		return P(add(x,p.x),add(y,-p.y));
	}
	P operator * (P p){
		return P(x*p.x,y*p.y);
	}
};
double dist(P a,P b){
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int cnt(P pos,vector<P> dat){
	int res=0;
	REP(i,dat.size()){
		if(dist(pos,dat[i])<=1.0)res++;
	}
	return res;
}
int main(){
	int n;
	double x,y;
	while(scanf("%d",&n)&&n){
		vector<P> dat;
		REP(i,n){
			cin>>x>>y;
			dat.PB(P(x,y));
		}
		int ans=0;
		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				P a=dat[i],b=dat[j];
				if(dist(a,b)>2.0){
					ans=max(ans,1);
					continue;
				}
				P mid=P((a.x+b.x)/2,(a.y+b.y)/2);
				double dx=a.y-b.y,dy=b.x-a.x;
				double dd=sqrt(dx*dx+dy*dy);
				double ddd=sqrt(1-dist(a,b)/2*dist(a,b)/2);
				dx=dx/dd*ddd;dy=dy/dd*ddd;
				ans=max(ans,cnt(P(mid.x+dx,mid.y+dy),dat));
				ans=max(ans,cnt(P(mid.x-dx,mid.y-dy),dat));
			}
		}
		printf("%d\n",ans);
	}
	return 0;
}