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 _overload3(_1, _2, _3, name, ...) name
#define _rep(i, n) repi(i, 0, n)
#define repi(i, a, b) for (int i = (a); i < (b); ++i)
#define rep(...) _overload3(__VA_ARGS__, repi, _rep, )(__VA_ARGS__)
#define ALL(x) x.begin(), x.end()
#define chmax(x, y) x = max(x, y)
#define chmin(x, y) x = min(x, y)
using namespace std;
using ll = long long;
using lld = long double;
using VI = vector<int>;
using VVI = vector<VI>;
using VL = vector<ll>;
using VVL = vector<VL>;
using PII = pair<int, int>;
using Point = pair<lld, lld>;
const lld EPS = 1e-9;

random_device rnd;
mt19937 mt(rnd());

lld dist(Point &a, Point &b)
{
    return sqrtl((a.first - b.first) * (a.first - b.first) + (a.second - b.second) * (a.second - b.second));
}

void update(Point &a, Point &center)
{
    Point differ = {a.first - center.first,
                    a.second - center.second};
    lld d = dist(center, a);
    differ.first /= d;
    differ.second /= d;
    center.first += differ.first * (d - 1.0);
    center.second += differ.second * (d - 1.0);
}

pair<Point, Point> get_center_of_cirlce(Point &a, Point &b)
{
    Point center = {0, 0};
    while (abs((lld)1.0 - dist(center, a)) > EPS || abs((lld)1.0 - dist(center, b)) > EPS)
    {
        update(a, center);
        update(b, center);
    }
    //cerr << fixed << setprecision(18) << center.first << " " << center.second << endl;
    //cerr << fixed << setprecision(18) << dist(center, a) << endl;
    //cerr << fixed << setprecision(18) << dist(center, b) << endl;
    Point differ1 = {a.first - center.first,
                     a.second - center.second};
    Point differ2 = {b.first - center.first,
                     b.second - center.second};
    Point other = {
        -a.first - b.first + center.first,
        -a.second - b.second + center.second};
    return {center, other};
}

int calc(vector<Point> &v, Point &center)
{
    int ret = 0;
    for (auto &&p : v)
    {
        //cerr << fixed << setprecision(18) << dist(p, center) << endl;
        if (dist(p, center) <= 1.0 + 1e-6)
        {
            ret++;
        }
        else
        {
            //cerr << "*" << fixed << setprecision(18) << abs(1.0 - dist(p, center)) << endl;
        }
    }
    return ret;
}

void solve()
{
    int n;
    cin >> n;
    if (n == 0)
        exit(0);
    vector<Point> v(n);
    int ans = 1;
    rep(i, n) cin >> v[i].first >> v[i].second;
    rep(i, n)
    {
        rep(j, i + 1, n)
        {
            if (dist(v[i], v[j]) > 2.0 + EPS)
            {
                continue;
            }
            auto ret = get_center_of_cirlce(v[i], v[j]);
            //cerr << fixed << setprecision(18) << abs(1.0 - dist(v[i], ret.first)) << endl;
            chmax(ans, calc(v, ret.first));
            chmax(ans, calc(v, ret.second));
        }
    }
    cout << ans << endl;
    return;
}

int main()
{
    while (true)
    {
        solve();
    }
    return 0;
}

C++: #include <iostream>
#include <complex>
#include <vector>

using namespace std;

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

#define EPS 1e-8

typedef complex<double> P;

double dot(P a,P b){
	return real(a)*real(b)+imag(a)*imag(b);
}

double cross(P a,P b){
	return real(a)*imag(b)-imag(a)*real(b);
}

double sq(double x){
	return x*x;
}

int main(){
	int n;
	while(cin>>n,n){
		vector<P> v(n);
		rep(i,n){
			double a,b;
			cin>>a>>b;
			v[i]=P(a,b);
		}
		int ans = 1;
		rep(i,n){
			rep(j,i){
				// cout<<i<<" "<<j<<" "<<abs(v[j]-v[i])<<endl;
				if(abs(v[j]-v[i])-2>EPS) continue;
				P c = v[j]-v[i];
				P n1 = (c*P(0,+1))/abs(c);
				P n2 = (c*P(0,-1))/abs(c);
				double d = sqrt(1-sq(abs(v[j]-v[i])/2));
				// cout<<"d "<<d<<endl;
				n1*=d;
				n2*=d;
				n1 += (v[j]-v[i])/2.0;
				n2 += (v[j]-v[i])/2.0;
				// cout<<"n1 "<<n1.real()<<" "<<n1.imag()<<endl;
				int tmp1=0;
				int tmp2=0;
				rep(k,n){
					if(1-abs(n1-v[k]+v[i])>-EPS) tmp1++;
					if(1-abs(n2-v[k]+v[i])>-EPS) tmp2++;
				}
				ans=max(ans,tmp1);
				ans=max(ans,tmp2);
			}
		}
		cout<<ans<<endl;
		// break;
	}
}

C++: #include "bits/stdc++.h"

#define REP(i,n) for(ll i=0;i<n;++i)
#define RREP(i,n) for(ll i=n-1;i>=0;--i)
#define FOR(i,m,n) for(ll i=m;i<n;++i)
#define RFOR(i,m,n) for(ll i=n-1;i>=m;--i)
#define ALL(v) (v).begin(),(v).end()
#define UNIQUE(v) v.erase(unique(ALL(v)),v.end());
#define DUMP(v) REP(aa, (v).size()) { cout << v[aa]; if (aa != v.size() - 1)cout << " "; else cout << endl; }
#define INF 1000000001ll
#define MOD 1000000007ll
#define EPS 1e-9

const int dx[8] = { 1,1,0,-1,-1,-1,0,1 };
const int dy[8] = { 0,1,1,1,0,-1,-1,-1 };


using namespace std;
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef vector<vi> vvi;
typedef vector<vl> vvl;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
ll max(ll a, int b) { return max(a, ll(b)); }
ll max(int a, ll b) { return max(ll(a), b); }
ll min(ll a, int b) { return min(a, ll(b)); }
ll min(int a, ll b) { return min(ll(a), b); }
///(?´????????`)(?´????????`)(?´????????`)(?´????????`)(?´????????`)(?´????????`)///

int main() {
	cin.tie(0);
	ios::sync_with_stdio(false);
	int n;
	while (cin >> n, n) {
		vector<double> x(n), y(n);
		int ans = 1;
		REP(i, n)cin >> x[i] >> y[i];
		REP(i, n) {
			REP(j, n) {
				if (i == j)continue;
				double ax = x[j] - x[i], ay = y[j] - y[i];
				if (ax*ax + ay*ay+EPS > 4)continue;
				double px = x[i] + ax / 2 + sqrt((4 - (ax*ax + ay*ay)) / (4 * (ax*ax + ay*ay)))*ay,
				py = y[i] + ay / 2 - sqrt((4 - (ax*ax + ay*ay)) / (4 * (ax*ax + ay*ay)))*ax;
				int cnt = 0;
				REP(k, n) {
					if ((x[k] - px)*(x[k] - px) + (y[k] - py)*(y[k] - py) < 1+EPS)cnt++;
				}
				ans = max(ans, cnt);
			}
		}
		REP(i, n) {
			REP(j, n) {
				if (i == j)continue;
				double ax = x[j] - x[i], ay = y[j] - y[i];
				if (ax*ax + ay*ay+EPS > 4)continue;
				double px = x[i] + ax / 2 - sqrt((4 - (ax*ax + ay*ay)) / (4 * (ax*ax + ay*ay)))*ay,
					py = y[i] + ay / 2 + sqrt((4 - (ax*ax + ay*ay)) / (4 * (ax*ax + ay*ay)))*ax;
				int cnt = 0;
				REP(k, n) {
					if ((x[k] - px)*(x[k] - px) + (y[k] - py)*(y[k] - py) < 1+EPS)cnt++;
				}
			
				ans = max(ans, cnt);
			}
		}
		cout <<ans << endl;
	}
}