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;
using ll=long long;
using vi=vector<int>;
using vvi=vector<vi>;
using pii=pair<int,int>;
#define rep(i,n) for(int i=0;i<n;i++)
#define range(i,a,n) for(int i=a;i<n;i++)
#define all(a) a.begin(),a.end()
#define rall(a) a.rbegin(),a.rend()
#define INF 1e9
#define EPS 1e-9
#define MOD (1e9+7)
void put(string d){}template<class H,class...T>void put(string d,H&h,T&...t){cout<<h;if(sizeof...(t))cout<<d;put(d,t...);}
template<class T>void puti(T&a,string d=" "){bool f=1;for(auto&_:a)cout<<(exchange(f,0)?"":d)<<_;cout<<endl;}
template<class T>void putii(T&a,string d=" "){for(auto&_:a)puti(_,d);}

using P=complex<double>;

int main(){
    int n;
    while(cin>>n,n){
        int ans=1;
        vector<P> c(n);
        rep(i,n) { double p,q; cin>>p>>q; c[i]=P(p,q); }
        rep(i,n)range(j,i+1,n){
            P v=c[i]-c[j];
            double d=abs(v);
            if(d>2.0) continue;
            P med=(c[i]+c[j])/2.0;
            P nml=v/d*P(0,1)*sqrt(1-d*d/4);
            P r1=med+nml;
            P r2=med-nml;
            int s1=0,s2=0;
            rep(k,n){
                s1+=(abs(r1-c[k])<=1.0);
                s2+=(abs(r2-c[k])<=1.0);
            }
            ans=max(ans,max(s1,s2));
        }
        cout<<ans<<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);
				}
			}
		}
		cout << ans << endl;
	}

	return 0;
}