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;
const double EPS = 1e-9;   // 許容誤差。問題によって変える
const double pi=M_PI;
const double INFD=1e+9;
typedef complex<double> P; // Point
#define rep(i,n) for(int i=0;i<n;i++)
//点の座標は複素数で持つ
#define X real()
#define Y imag()
namespace std{
    bool operator<(const P p1, const P p2){
        return p1.X != p2.X ? p1.X < p2.X : p1.Y < p2.Y;
    }
    P operator+(const P p1, const P p2) { return P(p1.X + p2.X, p1.Y + p2.Y); }
    P operator-(const P p1, const P p2) { return P(p1.X - p2.X, p1.Y - p2.Y); }
    P operator*(const P p, double k) { return P(p.X * k, p.Y * k); }
    P operator/(const P p, double k) { return P(p.X / k, p.Y / k); }
} 
//Segment/L
struct Segment{
    P p1, p2;
    Segment() {}
    Segment(P p1, P p2) : p1(p1), p2(p2) {}
};
typedef Segment L;

double norm(P &p1, P &p2){
    return (p1.X - p2.X) * (p1.X - p2.X) + (p1.Y - p2.Y) * (p1.Y - p2.Y);
}

double cross(const P &a, const P &b){
    return (conj(a) * b).Y;
}

double dot(const P &a, const P &b){
    return (conj(a) * b).X;
}

P Projection(L l, P p){
    P base = l.p2 - l.p1;                        //始点
    double r = dot(p - l.p1, base) / norm(base); //線分に対する倍率の計算
    return l.p1 + base * r;
}

double distanceLP(const L &l, const P &p){
    return abs(p - Projection(l, p));
}
//Circle
struct Circle{
    P p;
    double r;
    Circle(const P &p,double r) : p(p),r(r){}
};
typedef Circle C;

vector<P> crosspointCL(const C &c, const L &l) {
    vector<P> res;
    double d = distanceLP(l, c.p);      // 中心と直線の距離
    if (abs(d - c.r) < EPS) {  // 触れている
        res.push_back(Projection(l, c.p));
        return res;
    }
    if (d > c.r) return res;  // 離れている
    P h = Projection(l, c.p);
    P u = sqrt(c.r * c.r - d * d) * (l.p2- l.p1) / abs(l.p2 - l.p1);
    res.push_back(h + u);
    res.push_back(h - u);
    return res;
}

int intersectCC(const C &c1, const C &c2) {
    double d = abs(c1.p - c2.p), r1 = c1.r, r2 = c2.r;
    if (r1 + r2 < d) return 0;                   // 離れている
    if (abs(r1 + r2 - d) < EPS) return -2;  // 外接
    if (abs(d + r1 - r2) < EPS) return +1;  // c1 が c2 の中で内接
    if (abs(d + r2 - r1) < EPS) return -1;  // c2 が c1 の中で内接
    if (d + r1 < r2) return +3;                  // c1 が c2 の中
    if (d + r2 < r1) return -3;                  // c2 が c1 の中
    return 2;                                    // 2つの交点を持つ
}

vector<P> crosspointCC(const C &c1, const C &c2) {
    vector<P> res;
    double r1 = c1.r, r2 = c2.r, d = abs(c1.p - c2.p);
    if (d<EPS) return res;  // 中心が同じ
    int i = intersectCC(c1, c2);
    if (i == +1 || i == -1) {  // 内接
        if (r2 < r1)
            res.push_back(c1.p + r1 / d * (c2.p - c1.p));
        else
            res.push_back(c2.p + r2 / d * (c1.p - c2.p));
        return res;
    }
    if (i == -2) {  // 外接
        res.push_back((c1.p * c2.r + c2.p * c1.r) / (c1.r + c2.r));
        return res;
    }
    if (i == 0 || i == +3 || i == -3) {  // 共通部分なし || 内部
        return res;
    }
    // 2つの交点を持つ
    P p = c1.p - c2.p;
    double A = -2 * p.real(), B = -2 * p.imag();
    double C = norm(c1.p) - norm(c2.p) - r1 * r1 + r2 * r2;
    L l;
    if(abs(A)<EPS){
        l = L(P(0, -C / B), P(1, -C / B));
    }else if(abs(B)<EPS){
        l = L(P(-C / A, 0), P(-C / A, 1));
    }else if(abs(C)<EPS){
        l = L(P(0, 0), P(1, -A / B));
    }else{
        l = L(P(0, -C / B), P(-C / A, 0));
    }
    return crosspointCL(c1, l);  // Ax + By + C = 0
}
vector<P> points;
signed main(){
    int n;
    while(cin>>n and n){
        rep(i,n){
            double x,y;cin>>x>>y;
            P a=P(x,y);
            points.push_back(a);
        }
        int ans=1;
        vector<P> cent;
        rep(i,n)for(int j=i+1;j<n;j++){
            C c1=C(points[i],1),c2=C(points[j],1);
            cent=crosspointCC(c1,c2);
            rep(l,cent.size()){
                int tmp=0;
                rep(k,n){
                    if(abs(points[k]-cent[l])<=1+EPS) tmp++;
                }
                ans=max(ans,tmp);
            }
        }
        cout<<ans<<endl;
        points.clear();
    }
    return 0;
}

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

class P
{
public:
	double x,y;
	P(double x=0.0,double y=0.0):x(x),y(y){}
	P operator+(const P p)const{return P(x+p.x,y+p.y);}
	P operator-(const P p)const{return P(x-p.x,y-p.y);}
};

double dis(P p1,P p2)
{
	return hypot(p1.x-p2.x,p1.y-p2.y);
}

P getc(P p1,P p2)
{
	double theta=acos(dis(p1,p2)/2.0);
	double alpha=atan2(p2.y-p1.y,p2.x-p1.x);
	return P(p1.x+cos(alpha-theta),p1.y+sin(alpha-theta));
}

bool in(P c,P p)
{
	return dis(c,p)<1.0+1e-9;
}

int main()
{
	int n;
	while(cin>>n,n)
	{
		P p[300];
		for(int i=0;i<n;i++)
		{
			cin>>p[i].x>>p[i].y;
		}
		if(n==1)
		{
			cout<<1<<endl;
			continue;
		}
		int m=0;
		for(int i=0;i<n;i++)
		{
			for(int j=i+1;j<n;j++)
			{
				P a=p[i],b=p[j];
				P c=getc(a,b);
				int cnt=0;
				for(int k=0;k<n;k++)
				{
					if(in(c,p[k]))cnt++;
				}
				m>cnt?:m=cnt;
			}
		}
		for(int i=0;i<n;i++)
		{
			P a=p[i];
			int cnt=0;
			for(int j=0;j<n;j++)
			{
				if(in(a,p[j]))cnt++;
			}
			m>cnt?:m=cnt;
		}
		cout<<m<<endl;
	}
}