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<string>
#include<algorithm>
#include<vector>
#include<iomanip>
#include<math.h>
#include<complex>
#include<queue>
#include<deque>
#include<map>
#include<set>
#include<bitset>
using namespace std;
#define REP(i,m,n) for(int i=(int)m ; i < (int) n ; i++ )
#define rep(i,n) REP(i,0,n)
typedef long long ll;
typedef pair<int,int> pint;
const int inf=1e9+7;
const ll longinf=1LL<<60 ;
const ll mod=1e9+7 ;
int dx[4]={1,0,-1,0} , dy[4]={0,1,0,-1} ;

const double eps=1e-10;
double add(double a,double b){
    if(abs(a+b)<eps*(abs(a)+abs(b)))return 0;
    return a+b;
}

typedef pair<double,double> P;
P operator+(const P& x,const P& y){
    return {add(x.first,y.first),add(x.second,y.second)};
}

P operator-(const P& x,const P& y){
    return {add(x.first,-y.first),add(x.second,-y.second)};

}

P operator*(const P& x,const double& y){
    return {x.first*y,x.second*y};
}
double dot(P x,P y){
    return add(x.first*y.first,x.second*y.second);


}

double det(P x,P y){
    return add(x.first*y.second,-x.second*y.first);
}

bool onseg(P p1,P p2,P q){
    return det(p1-q,p2-q)==0 &&dot(p1-q,p2-q)<=0; 
}

P intersection(P p1,P p2,P q1,P q2){
    return p1+(p2-p1)*(det(q2-q1,q1-p1)/det(q2-q1,p2-p1));
}

bool cross(P p1,P p2,P q1,P q2){
    if(det(p1-p2,q1-q2)==0)return onseg(p1,p2,q1)||onseg(p1,p2,q2)||onseg(q1,q2,p1)||onseg(q1,q2,p2);
    P ret=intersection(p1,p2,q1,q2);
    return onseg(p1,p2,ret)&&onseg(q1,q2,ret);
}

double dist(P p1,P p2,P q){
    if(dot(q-p1,p2-p1)<=0)return sqrt(dot(q-p1,q-p1));
    if(dot(q-p2,p1-p2)<=0)return sqrt(dot(q-p2,q-p2));
    P d=p1-p2;
    P ret=q+P{d.second,-d.first};
    P h=intersection(p1, p2, q, ret);
    return sqrt(dot(q-h,q-h));
}
double d(P p,P q){
    return sqrt(dot(p-q,p-q));
}

int main(){
    int n;
    while(cin>>n,n!=0){
        P a[n];
        rep(i,n)cin>>a[i].first>>a[i].second;
        int ans=1;
        rep(i,n)rep(j,n){
            if(i==j)continue;
            if(d(a[i],a[j])>2)continue;
            P v=(a[i]-a[j])*(1.0/d(a[i],a[j]));
            P h=P{v.second,-v.first};
            P c=(a[i]+a[j])*(1.0/2)+h*sqrt(1.0-dot(a[i]-a[j],a[i]-a[j])/4);
            int ret=2;
            rep(k,n){
                if(k==i||k==j)continue;
                if(d(a[k],c)<=1)ret++;
            }
            ans=max(ans,ret);
            c=(a[i]+a[j])*(1.0/2)-h*sqrt(1.0-dot(a[i]-a[j],a[i]-a[j])/4);
            ret=2;
            rep(k,n){
                if(k==i||k==j)continue;
                if(d(a[k],c)<=1)ret++;
            }
            ans=max(ans,ret);
        }
        cout<<ans<<endl;
    }
    return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;

typedef ostringstream OSS;
typedef istringstream ISS;

typedef long long LL;
typedef pair<int, int> PII;
typedef pair<double, double> PDD;

typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<LL> VLL;
typedef vector<VLL> VVLL;
typedef vector<double> VD;
typedef vector<VD> VVD;
typedef vector<string> VS;
typedef vector<VS> VVS;
typedef vector<bool> VB;
typedef vector<VB> VVB;
typedef vector<PII> VPII;
typedef vector<PDD> VPDD;

#define fst first
#define snd second
#define Y first
#define X second
#define MP make_pair
#define PB push_back
#define EB emplace_back 
#define ALL(x) (x).begin(),(x).end()
#define RANGE(x,y,maxX,maxY) (0 <= (x) && 0 <= (y) && (x) < (maxX) && (y) < (maxY))
#define DUMP( x ) cerr << #x << " = " << ( x ) << endl

template < typename T > inline T fromString(const string &s) { T res; ISS iss(s); iss >> res; return res; };
template < typename T > inline string toString(const T &a) { OSS oss; oss << a; return oss.str(); };

const int INF = 0x3f3f3f3f;
const LL INFL = 0x3f3f3f3f3f3f3f3fLL;
const double DINF = 0x3f3f3f3f;
const int DX[]={1,0,-1,0},DY[]={0,-1,0,1};

const double EPS = 1e-6;

VD xs;
VD ys;
int N;

double distance(int i, int j) {
    double x = xs[i] - xs[j];
    double y = ys[i] - ys[j];
    return sqrt(x * x + y * y);
}

bool is_over_circle(int i, int j) {
    return distance(i, j) < 2.0;
}

VPDD calc_circle(int i, int j) {
    double d = distance(i, j);
    double a = acos(d / 2.0);
    double b = acos(abs(xs[j] - xs[i]) / d);
    double sita = a + b;
    PDD p, q;
    p.X = xs[i] + cos(sita);
    p.Y = ys[i] + sin(sita);
    q.X = (xs[i] + xs[j]) - p.X;
    q.Y = (ys[i] + ys[j]) - p.Y;
    return VPDD({p, q});
}

int count_inner_points(PDD circle) {
    int cnt = 0;

    for (int i = 0; i < N; i++) {
        double x = circle.X - xs[i];
        double y = circle.Y - ys[i];
//        printf("--- %.3f %.3f %.3f %.3f\n", x, y, x * x + y * y, sqrt(x * x + y * y));
        cnt += x * x + y * y < 1.0 + EPS;
    }

    return cnt;
}

int solve() {
    int res = 1;

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            if (is_over_circle(i, j)) {
                VPDD ps = calc_circle(i, j);

                for (auto p : ps) {
//                    printf("%.3f %3f\n", p.X, p.Y);
                    int cnt = count_inner_points(p);
                    res = max(res, cnt);
                }
            }
        }
    }

    return res;
}

int main(void) {
    while (cin >> N, N) {
        xs.clear();
        ys.clear();
        xs.resize(N);
        ys.resize(N);
        for (int i = 0; i < N; i++) {
            cin >> xs[i] >> ys[i];
        }
        
        cout << solve() << endl;
    }

	return 0;
}