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 REP(i,n) for(int i=0; i<(int)(n); ++i)

using namespace std;

const double EPS = 1e-8;

typedef complex<double> Point;

int calc(Point p, vector<Point>& ps){
    int ans = 0;
    for(int i = 0; i < ps.size(); i++){
        if(abs(ps[i] - p) < 1.0 + EPS) {
            ans++;
        }
    }
    return ans;
}
vector<Point> touching_circle(Point a, Point b, double r){
    vector<Point> res;

    double d = abs(b - a);
    if(d > 2 * r) return res;

    Point mid = 0.5 * (a + b);
    Point dir = polar(sqrt(r * r - d * d / 4), arg(b - a) + M_PI / 2);
    res.push_back(mid + dir);
    res.push_back(mid - dir);
    return res;
}
int main(){
    int N;
    while(cin >> N && N){
        vector<Point> ps(N);
        REP(i, N){
            double x, y;
            cin >> x >> y;
            ps[i] = Point(x, y);
        }
        int ans = 0;
        REP(i, N) REP(j, N) if(abs(ps[i] - ps[j]) > EPS) {
            vector<Point> tp = touching_circle(ps[i], ps[j], 1.0);
            for(Point c : tp){
            assert(abs(1.0 - abs(c - ps[i])) < EPS);
            assert(abs(1.0 - abs(c - ps[j])) < EPS);
                ans = max(ans, calc(c, ps));
            }
        }
        cout << ans << endl;
    }
    return 0;
}

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

double Dis(double x1, double y1, double x2, double y2);
int main(void){
	while(1){
		int n;
		int ans = 0;
		vector< pair<double, double> > p;

		cin >> n;
		if(n == 0) break;

		for(int i=0; i<n; i++){
			pair<double, double> t;
			cin >> t.first >> t.second;
			p.push_back(t);
		}

		for(int i=0; i<n; i++){
			for(int j=0; j<n; j++){
				//中点
				double mpx = (p[i].first + p[j].first) / 2.0;
				double mpy = (p[i].second + p[j].second) / 2.0;
				//2点間の距離
				double d1 = Dis(p[i].first, p[i].second, p[j].first, p[j].second);
				//p[i]からp[j]へ向かうベクトル
				double v1x = (p[i].first - p[j].first) / d1;
				double v1y = (p[i].second - p[j].second) / d1;
				//中点から円の中心へ向かう単位ベクトル
				double v2x = -v1y;
				double v2y = v1x;
				//中心から円の中心へ向かうベクトルの長さ
				double d2 = sqrt(1.0 - (d1/2.0)*(d1/2.0));

				//円の中心の座標
				double cx = mpx + v2x * d2;
				double cy = mpy + v2y * d2;

				int cnt = 0;
				for(int k=0; k<n; k++){
					if(Dis(cx, cy, p[k].first, p[k].second) <= 1.0)	cnt++;
				}
				if(cnt > ans){
					ans = cnt;
				}

				cx = mpx - v2x * d2;
				cy = mpy - v2y * d2;

				cnt = 0;
				for(int k=0; k<n; k++){
					if(Dis(cx, cy, p[k].first, p[k].second) <= 1.0)	cnt++;
				}
				if(cnt > ans){
					ans = cnt;
				}
			}
		}

		cout << ans << endl;


	}
	
	return 0;
}

double Dis(double x1, double y1, double x2, double y2){
	return sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2));
}

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

using namespace std;

complex<double> Z[300];

int main(void)
{
    int N, ans, count;
    double X, Y;
    while(true)
    {
        cin >> N;
        if (!N) break;
        for (int i = 0; i < N; i++)
        {
            cin >> X >> Y;
            Z[i].real(X);
            Z[i].imag(Y);
        }
        ans = 0;
        for (int i = 0; i < N; i++)
        {
            for (int j = i + 1; j < N; j++)
            {
                if (abs(Z[i] - Z[j]) >= 1.9999) continue;
                auto a = (Z[i] - Z[j]) / 2.;
                auto b = (Z[i] + Z[j]) / 2.;
                auto r = sqrt(1. - norm(a));
                a /= abs(a);
                a *= r * 1.i;
                complex<double> c[2] = {b + a, b - a};
                for (int k = 0; k < 2; k++)
                {
                    count = 2;
                    for (int l = 0; l < N; l++)
                    {
                        if (l == i || l == j) continue;
                        if (abs(c[k] - Z[l]) <= 1.) count++;
                    }
                    if (count > ans) ans = count;
                }
            }
        }
        cout << (ans ? ans : 1) << endl;
    }
    return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int,int> P;
typedef pair<int,P> P1;
typedef pair<P,P> P2;
#define pu push
#define pb push_back
#define mp make_pair
#define INF 1000000000
#define fi first
#define sc second
#define rep(i,x) for(int i=0;i<x;i++)
#define SORT(x) sort(x.begin(),x.end())
#define ERASE(x) x.erase(unique(x.begin(),x.end()),x.end())
#define POSL(x,v) (lower_bound(x.begin(),x.end(),v)-x.begin())
#define POSU(x,v) (upper_bound(x.begin(),x.end(),v)-x.begin())
typedef complex<double> pt;
typedef pair<pt,pt> L;
typedef vector<P> poly;
const double EPS = 1e-11;
#define x real()
#define y imag()
 
double dot(pt a,pt b){
	return (conj(a)*b).x;
}
double cross(pt a,pt b){
	return (conj(a)*b).y;
}
int ccw(pt a,pt b,pt c){
	b -= a; c -= a;
	if(cross(b,c) > EPS) return 1; // counter clockwise
	if(cross(b,c) < -EPS) return -1; // clockwise
	if(dot(b,c) < -EPS) return 2; //c-a-b
	if(norm(b) < norm(c)) return -2; //a-b-c
	return 0; //a-c-b
}
bool cmp(const pt& a,const pt& b){
	if(-EPS < a.x-b.x && a.x-b.x < EPS) return a.y < b.y;
	else return a.x < b.x;
}
vector<pt>convex_hull(vector<pt>ps)
{
	sort(ps.begin(),ps.end(),cmp);
	int k=0,n = ps.size();
	vector<pt>qs(n*2);
	for(int i=0;i<n;i++)
	{
		while(k>1 && ccw(qs[k-2],qs[k-1],ps[i]) == -1) k--;
		qs[k++] = ps[i];
	}
	for(int i=n-2,t=k;i>=0;i--)
	{
		while(k>t && ccw(qs[k-2],qs[k-1],ps[i]) == -1) k--;
		qs[k++] = ps[i];
	}
	qs.resize(k-1);
	return qs;
}
pt intersect(pt a,pt b,pt c,pt d){
	double dev = (b.y-a.y)*(d.x-c.x)-(b.x-a.x)*(d.y-c.y);
	
	double d1, d2;
 
	d1 = (c.y*d.x-c.x*d.y);
	d2 = (a.y*b.x-a.x*b.y);
 
	pt ap1;
	ap1.x = d1*(b.x-a.x) - d2*(d.x-c.x);
	ap1.x /= dev;
	ap1.y = d1*(b.y-a.y) - d2*(d.y-c.y);
	ap1.y /= dev;
 
	return ap1;
}
bool par(pt a,pt b,pt c,pt d){
	b.x -= a.x; b.y -= a.y;
	d.x -= c.x; d.y -= c.y;
	if(abs(b.y*d.x - b.x*d.y) < EPS){
		return 1;
	}
	return 0;
}
double PI = 3.1415926535897932;
int main(){
	while(1){
		int n; scanf("%d",&n); if(!n) return 0; if(n==1){ puts("1"); continue;}
		pt p[305];
		for(int i=1;i<=n;i++){
			double a,b; scanf("%lf%lf",&a,&b);
			p[i] = pt(a,b);
		}
		int ans = 0;
		for(int i=1;i<=n;i++){
			for(int j=i+1;j<=n;j++){
				pt q = (p[i]+p[j])/2.0;
				pt r = (p[i]-p[j]); if(abs(r) > 2.0+EPS) continue;
				pt p1 = q + r / abs(r) * polar(sqrt(1.0-abs(r)*abs(r)/4.0),PI/2);
				pt p2 = q - r / abs(r) * polar(sqrt(1.0-abs(r)*abs(r)/4.0),PI/2);
				int c=0,d=0;
				for(int k=1;k<=n;k++){
					if(k==i || k==j || abs(p1-p[k]) <= 1.0+EPS) c++;
					if(k==i || k==j || abs(p2-p[k]) <= 1.0+EPS) d++;
				}
				ans = max(ans,max(c,d));
			}
		}
		printf("%d\n",ans);
	}
}