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;
#define int long long   // <-----!!!!!!!!!!!!!!!!!!!

#define rep(i,n) for (int i=0;i<(n);i++)
#define rep2(i,a,b) for (int i=(a);i<(b);i++)
#define rrep(i,n) for (int i=(n)-1;i>=0;i--)
#define rrep2(i,a,b) for (int i=(a)-1;i>=b;i--)
#define all(a) (a).begin(),(a).end()
#define rall(a) (a).rbegin(),(a).rend()
#define printV(_v) for(auto _x:_v){cout<<_x<<" ";}cout<<endl
#define printVS(vs) for(auto x : vs){cout << x << endl;}
#define printVV(_vv) for(auto _v:_vv){for(auto _x:_v){cout<<_x<<" ";}cout<<endl;}
#define printP(p) cout << p.first << " " << p.second << endl
#define printVP(vp) for(auto p : vp) printP(p);

typedef long long ll;
typedef pair<int, int> Pii;
typedef tuple<int, int, int> TUPLE;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef vector<vvi> vvvi;
typedef vector<Pii> vp;
const int inf = 1e9;
const int mod = 1e9 + 7;

typedef complex<double> P;
typedef vector<P> G;
#define here(g, i) g[i]
#define next(g, i) g[(i + 1) % g.size()]
#define prev(g, i) g[(i - 1 + g.size()) % g.size()]
const double EPS = 1e-10;
const double INF = 1e12;
const double PI = acos(-1);

P readP() {
    double x, y;
    cin >> x >> y;
    return P(x, y);
}

// rot p around q by theta (counter-clockwise)
P rotP(P p, P q, double theta) {
    p -= q;
    double x = p.real(), y = p.imag();
    p = P(x * cos(theta) - y * sin(theta), x * sin(theta) + y * cos(theta));
    p += q;
    return p;
}

P center(P p, P q) {
    P v = p - q;
    double theta = asin(abs(v) / 2);
    v *= (1 / abs(v));
    P u = rotP(v, p, theta);
    return p + u;
}

signed main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(0);

    int n;
    while (cin >> n, n) {
        vector<P> p(n);
        rep(i, n) p[i] = readP();

        int ans = 1;
        rep(i, n) {
            rep(j, n) {
                if (i == j) continue;
                if (abs(p[i] - p[j]) < 2 + EPS) {
                    P c = center(p[i], p[j]);
                    int cnt = 2;
                    rep(k, n) {
                        if (k == i || k == j) continue;
                        if (abs(p[k] - c) < 1 + EPS) {
                            cnt++;
                        }
                    }
                    ans = max(ans, cnt);
                }
            }
        }

        cout << ans << endl;
    }
}

C++: #include <cstdio>
#include <complex>
#define rep(i,n) for(int i=0;i<n;i++)
#define re real()
#define im imag()
using namespace std;
typedef complex<double> P;
int main(){
	while(true){
		int n,mx=1;
		scanf("%d",&n);
		if(n==0) break;
		complex<double> p[300];
		rep(i,n){
			double x,y;
			scanf("%lf%lf",&x,&y);
			p[i]=P(x,y);
		}
		rep(i,n) rep(j,i){
			double d=abs(p[i]-p[j]);
			if(d>2) continue;
			P c=(p[i]+p[j])/2.0,dis=(p[i]-p[j])/d*sqrt(1-d*d/4);
			dis=P(dis.im,-dis.re);
			P cen=c+dis;
			int cnt=0;
			rep(k,n) if(abs(cen-p[k])<=1.0000000001) cnt++;
			mx=max(mx,cnt);
			cen=c-dis;
			cnt=0;
			rep(k,n) if(abs(cen-p[k])<=1.0000000001) cnt++;
			mx=max(mx,cnt);
		}
		printf("%d\n",mx);
	}
	return 0;
}

C++: #include <iostream>
#include <cmath>
#include <algorithm>
#include <queue>
#include <map>
#define LP(i,n) for(int i = 0; i < n; i++)

using namespace std;
typedef pair<double,double> P;

double Dist(double x1[], double x2[]){
	return pow((x1[0]-x2[0])*(x1[0]-x2[0])+(x1[1]-x2[1])*(x1[1]-x2[1]),0.5);

}

int main(){
	int n;
	cin >> n;
	while(n){
		int mx = 0;
		if(n >= 2){
			P p[300];
			LP(i,n){
				cin >> p[i].first >> p[i].second;
			}
			sort(p,p+n);
			double x[300][2];
			LP(i,n){
				x[i][0] = p[i].first;
				x[i][1] = p[i].second;
			}
			double dist[300][300];
			LP(i,n){
				LP(j,n){
					dist[i][j] = Dist(x[i],x[j]);
				}
			}
			LP(i,n){
				for(int j = i+1; j < n; j++){
					if(x[j][0]-x[i][0] > 2){
						break;
					}else{
						if(dist[i][j] <= 2){
							double m[] = {(x[j][0]+x[i][0])/2,(x[j][1]+x[i][1])/2};
							double d[] = {x[j][1]-x[i][1],-x[j][0]+x[i][0]};
							double O[] = {0,0};
							double _d = Dist(d,O);
							LP(k,2){
								d[k] *= (pow(_d,-2)-0.25);
							}
							double o[2][2];
							LP(k,2){
								LP(l,2){
									o[k][l] = m[l];
									if(k) o[k][l] += d[l];
									else o[k][l] -= d[l];
								}
							}
							LP(k,2){
								int s = 0;
								for(int l = i; l < j; l++){
									if(Dist(x[l],o[k]) <= 1) s++;
								}
								for(int l = i-1; 0<=l && x[i][0]-x[l][0] <= 1;l--){
									if(Dist(x[l],o[k]) <= 1) s++;
								}
								for(int l = j; l < n && x[l][0]-x[j][0] <= 1; l++){
									if(Dist(x[l],o[k]) <= 1) s++;
								}
								mx = max(mx,s);
							}
						}
					}
				}
			}
		}else mx = 1;
		cout << mx << endl;
		cin >> n;
	}
}