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 <vector>
#include <limits>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <iomanip>

using namespace std;

const double EPS = 1e-9;

bool tolerant_eq(double a, double b){ return abs(a - b) < EPS; }

struct Point {
	double x, y;
	explicit Point(const double &x = 0.0, const double &y = 0.0) :
		x(x), y(y)
	{ }

	static Point invalid(){
		double qnan = numeric_limits<double>::quiet_NaN();
		return Point(qnan, qnan);
	}
	bool is_valid() const { return !(isnan(x) || isnan(y)); }

	Point operator+(const Point &p) const { return Point(x + p.x, y + p.y); }
	Point &operator+=(const Point &p){ return *this = *this + p; }
	Point operator-(const Point &p) const { return Point(x - p.x, y - p.y); }
	Point &operator-=(const Point &p){ return *this = *this - p; }
	Point operator*(double s) const { return Point(x * s, y * s); }
	Point &operator*=(double s){ return *this = *this * s; }
	Point operator/(double s) const { return Point(x / s, y / s); }
	Point &operator/=(double s){ return *this = *this / s; }

	bool operator==(const Point &p) const { return x == p.x && y == p.y; }
	bool operator!=(const Point &p) const { return x != p.x || y != p.y; }
	bool operator<(const Point &p) const {
		return (x == p.x) ? (y < p.y) : (x < p.x);
	}
};

Point operator*(double s, const Point &p){ return p * s; }

bool tolerant_eq(const Point &a, const Point &b){
	return tolerant_eq(a.x, b.x) && tolerant_eq(a.y, b.y);
}
double abs(const Point &p){ return sqrt(p.x * p.x + p.y * p.y); }
Point unit(const Point &p){ return p / abs(p); }
Point ortho(const Point &p){ return Point(-p.y, p.x); }
double cross(const Point &a, const Point &b){ return a.x * b.y - a.y * b.x; }
double dot(const Point &a, const Point &b){ return a.x * b.x + a.y * b.y; }
int ccw(const Point &a, const Point &b, const Point &c){
	Point d = b - a, e = c - a;
	if(cross(d, e) > 0.0){ return 1; }
	if(cross(d, e) < 0.0){ return -1; }
	if(dot(d, e) < 0.0){ return 2; }
	if(abs(d) < abs(e)){ return -2; }
	return 0;
}

struct Circle {
	Point c;
	double r;
	explicit Circle(const Point &c = Point(), double r = 0.0) :
		c(c), r(r)
	{ }

	static Circle invalid(){
		double qnan = numeric_limits<double>::quiet_NaN();
		return Circle(Point::invalid(), qnan);
	}
	bool is_valid() const { return c.is_valid() && !isnan(r); }

	bool operator<(const Circle &i) const {
		return (r == i.r) ? (c < i.c) : (r < i.r);
	}
};

bool intersect(const Circle &a, const Circle &b){
	return abs(a.c - b.c) <= a.r + b.r;
}

vector<Point> crosspoint(const Circle &a, const Circle &b){
	if(!intersect(a, b)){ return vector<Point>(); }
	double l = abs(b.c - a.c), s = (l + a.r + b.r) / 2.0;
	double S = sqrt(s * (s - l) * (s - a.r) * (s - b.r));
	double d = (S * 2) / l, t = sqrt(a.r * a.r - d * d);
	Point v = t * unit(b.c - a.c);
	Point w = d * unit(ortho(v));
	Point x = a.c + v + w, y = a.c + v - w;
	if(tolerant_eq(x, y)){ return vector<Point>(1, x); }
	vector<Point> ans(2); ans[0] = x; ans[1] = y;
	return ans;
}

int main(){
	ios_base::sync_with_stdio(false);
	cout << setiosflags(ios::fixed) << setprecision(10);
	while(true){
		int N;
		cin >> N;
		if(N == 0){ break; }
		vector<Point> points(N);
		for(int i = 0; i < N; ++i){ cin >> points[i].x >> points[i].y; }
		int answer = 1;
		for(int i = 0; i < N; ++i){
			Circle ci(points[i], 1.0);
			for(int j = i + 1; j < N; ++j){
				Circle cj(points[j], 1.0);
				if(!intersect(ci, cj)){ continue; }
				vector<Point> cp = crosspoint(ci, cj);
				for(int k = 0; k < cp.size(); ++k){
					Point c = cp[k];
					int count = 0;
					for(int l = 0; l < N; ++l){
						if(abs(c - points[l]) <= 1.0 + EPS){ ++count; }
					}
					answer = max(answer, count);
				}
			}
		}
		cout << answer << endl;
	}
	return 0;
}

C++: #include<iostream>
#include<vector>
#include<complex>
#define EPS 1.0e-5
using namespace std;
typedef double elem;
typedef complex<elem> point,vec;
inline elem dist(const point &a, const point &b){return abs(a-b);}
inline vec normr(const vec &v){return v*vec(0,-1);}
inline vec norml(const vec &v){return v*vec(0,1);}
inline vec unormr(const vec &v){return normr(v) / abs(v);}
inline vec unorml(const vec &v){return norml(v) / abs(v);}
int main(){
	while(true){
		int ans = 0;
		int n;
		scanf("%d", &n);
		if(n==0)break;
		vector< point > points;
		for(int i=0;i<n;++i){
			elem x, y;
			scanf("%lf%lf", &x, &y);
			points.push_back( point(x,y) );
		}
		for(unsigned int i = 0; i < points.size(); ++i){
			for(unsigned int j = 0; j < points.size(); ++j){
				if( i != j && dist( points[i], points[j] ) < 2 ){
					int count_p = 0;
					int count_q = 0;
					elem d = sqrt(1-pow(dist(points[i], points[j])/2,2 ) );
					point p,q;
					vec ab = points[j] - points[i];
					point t = 0.5 * ab + points[i];  
					p = d * unormr(ab) + t;
					q = d * unorml(ab) + t;
					for(unsigned int k = 0; k < points.size(); ++k){
						if( dist(p, points[k]) <= 1.0+EPS ){ count_p++; }
						if( dist(q, points[k]) <= 1.0+EPS ){ count_q++; }
					}
					ans = max(ans, max(count_q, count_p) );
				}
			}
		}
		printf("%d\n", ans);
	}
	return 0;
}

C++: #include<iostream>
#include<cmath>
#include<algorithm>
#include<utility>
using namespace std;

double X[310],Y[310];

bool check(double x1, double y1, double x2, double y2){
    return (x1-x2)*(x1-x2)+(y1-y2)*(y1-y2) <= 1.0;
}

int main(){
    int N;
    while(1){
        cin>>N;
        if(!N) return 0;
        for(int i=1;i<=N;++i){
            cin >>X[i]>>Y[i];
        }

        int ans = 1;
        for(int i=1;i<=N;++i){
            for(int j=i+1;j<=N;++j){
                double t = (X[i]-X[j])*(X[i]-X[j])+(Y[i]-Y[j])*(Y[i]-Y[j]);

                if(t>=4.) continue;

                double s = sqrt(1./t - 1./4.);
                double p = (X[i]+X[j])/2., q = (Y[i]+Y[j])/2.;

                double x = p + (Y[j]-Y[i])*s;
                double y = q - (X[j]-X[i])*s;
                int cnt = 2;
                for(int k=1;k<=N;++k){
                    if(k==i || k==j) continue;
                    if(check(X[k],Y[k],x,y)) ++cnt;
                }
                ans=max(ans,cnt);

                x = p - (Y[j]-Y[i])*s;
                y = q + (X[j]-X[i])*s;
                cnt = 2;
                for(int k=1;k<=N;++k){
                    if(k==i || k==j) continue;
                    if(check(X[k],Y[k],x,y)) ++cnt;
                }
                ans=max(ans,cnt);
            }
        }
        cout<<ans<<endl;
    }

}