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 <iostream>
#include <complex>
#include <vector>
#include <algorithm>
#include <cmath>
using namespace std;

// xy平面上の点(ベクトル)を表現するには、complex型を利用するとよい
typedef complex<double> P;

// 辺の表現 (座標を2つ pair でもつ)
typedef pair<P, P> L;

// 円の表現 (座標 P と 半径 d で表現する)
typedef pair<P, double> C;

// 成分を取り出すのを簡単にする
#define X real()
#define Y imag()

// 誤差(epsilon)の定義
#define EPS (1e-10)

// 2つの要素が等しいかどうか
#define EQ(a,b) (abs((a) - (b)) < EPS)

// 2つのベクトルが等しいかどうか
#define EQV(a,b) ( EQ((a).X, (b).X) && EQ((a).Y, (b).Y) )

// m は n より大きい(以上)かどうか
#define LE(n, m) ((n) < (m) + EPS)
#define LEQ(n, m) ((n) <= (m) + EPS)

// m は n より小さい(以下)かどうか
#define GE(n, m) ((n) + EPS > (m))
#define GEQ(n, m) ((n) + EPS >= (m))

// 2つのベクトルの内積を求める
double dot(P a, P b) {
    return (a.X * b.X + a.Y * b.Y);
}

// 2つのベクトルの外積を求める
double cross(P a, P b) {
    return (a.X * b.Y - a.Y * b.X);
}

// 点 a と 点 b を通り、半径が r の円の中心を返す
vector<P> circlesPointsRadius(P a, P b, double r) {
    vector<P> cs;
    P abH = (b-a) * 0.5;
    double d = abs(abH);

    if(d == 0 || d > r) return cs; // 必要なら !LE(d,r) として円1つになる側へ丸める
    double dN = sqrt(r*r - d*d);   // 必要なら max(r*r - d*d, 0) とする
    P n = abH * P(0,1) * (dN / d);
    cs.push_back(a + abH + n);
    if(dN > 0) cs.push_back(a + abH - n);
    return cs;
}

int main() {
    int N;

    while(cin >> N, N) {
        vector<P> ps;
        for(int i=0; i<N; i++) {
            double x, y; cin >> x >> y;
            ps.push_back(P(x, y));
        }

        if(N == 1) cout << 1 << endl;
        else if(N == 2) {
            cout << 1 + (abs(ps[0] - ps[1]) < 2 + EPS) << endl;
        }
        else {
            int ans = 1;
            for(int i=0; i<N; i++) {
                for(int j=i+1; j<N; j++) {
                    // printf("i = %d, j = %d, abs = %f\n", i, j, abs(ps[i] - ps[j]));
                    if(abs(ps[i] - ps[j]) > 2 + EPS) continue;
                    // printf("valid\n");

                    vector<P> centers = circlesPointsRadius(ps[i], ps[j], 1);
                    for(auto c : centers) {
                        int sum = 2;
                        for(int k=0; k<N; k++) {
                            if(i == k || j == k) continue;
                            sum += (abs(c - ps[k]) < 1 + EPS);
                        }
                        ans = max(ans, sum);
                    }
                }
            }
            cout << ans << endl;
        }
    }
    return 0;
}

C++: #include <iostream>
#include <fstream>
#include <cassert>
#include <typeinfo>
#include <vector>
#include <stack>
#include <cmath>
#include <set>
#include <map>
#include <string>
#include <algorithm>
#include <cstdio>
#include <queue>
#include <iomanip>
#include <cctype>
#include <random>
#include <time.h>
#define syosu(x) fixed<<setprecision(x)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> P;
typedef pair<double,double> pdd;
typedef pair<ll,ll> pll;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef vector<double> vd;
typedef vector<vd> vvd;
typedef vector<ll> vl;
typedef vector<vl> vvl;
typedef vector<char> vc;
typedef vector<vc> vvc;
typedef vector<string> vs;
typedef vector<bool> vb;
typedef vector<vb> vvb;
typedef vector<P> vp;
typedef vector<vp> vvp;
typedef vector<pll> vpll;
typedef pair<P,int> pip;
typedef vector<pip> vip;
const int inf=1<<28;
const ll INF=1ll<<60;
const double pi=acos(-1);
const double eps=1e-8;
const ll mod=1e9+7;
const int dx[4]={0,1,0,-1},dy[4]={1,0,-1,0};

struct point{
	double x,y;
	point operator+(point p){
		return point{x+p.x,y+p.y};
	}
	point operator-(point p){
		return point{x-p.x,y-p.y};
	}
	point operator*(double p){
		return point{x*p,y*p};
	}
	point operator/(double p){
		if(!p) return point{0,0}; 
		return point{x/p,y/p};
	}
	bool operator==(point p){
		return fabs(x-p.x)<eps&&fabs(y-p.y)<eps;
	}
	bool operator<(point p){
		if(fabs(x-p.x)>eps) return x<p.x;
		return y<p.y;
	}
};

typedef pair<point,point> pp;
typedef vector<point> VP;
const point O{0,0};

double Length(point x,point y){
	point z=y-x;
	return sqrt(z.x*z.x+z.y*z.y);
}

point Normal(point p){
	return point{p.y,-p.x};
}

int n;
VP p;

int f(point q){
	int res=0;
	for(int i=0;i<n;i++) if(Length(p[i],q)<1+eps) res++;
	return res;
}

int main(){
	while(1){
		cin>>n;
		if(!n) break;
		p=VP(n);
		for(int i=0;i<n;i++) cin>>p[i].x>>p[i].y;
		int res=1;
		for(int i=0;i<n;i++) for(int j=i+1;j<n;j++){
			point c=(p[i]+p[j])/2,e=Normal(p[i]-p[j]);
			double l=Length(c,p[i]);
			e=e/Length(O,e)*sqrt(1.0-l*l);
			res=max(res,max(f(c+e),f(c-e)));
		}
		cout<<res<<endl;
	}
}

C++: #include<cmath>  
#include<iostream>  
#include<algorithm>  
#define EPS 1e-5

typedef std::pair<double,double> pt;
#define x first
#define y second


double dist_square(pt a, pt b){
    double x = a.x - b.x;
    double y = a.y - b.y;
    return x*x + y*y;
}

int solve(pt P,pt Q, std::vector<pt> pts){  
    pt tmp=pt((P.x+Q.x)/2,(P.y+Q.y)/2);  
    double angle=atan2(P.y-Q.y,P.x-Q.x)+acos(-1.0)/2;  
    double len=sqrt(1-dist_square(P,Q)/4);  
    tmp.x+=cos(angle)*len;  
    tmp.y+=sin(angle)*len;  
    int inside = 0;  
    for(auto it = pts.begin();it != pts.end();++it){
        if(dist_square(*it,tmp) <= 1+EPS) ++inside;
    }  

    return inside;  
}  
int main()  
{  
    int numPts;
    std::cin >> numPts;
    while(numPts != 0){  
        if(numPts == 1){
            std::cout << 1 << "\n";
            std::cin >> numPts;
            continue;
        }
        std::vector<pt> pts;

        for(int i=0;i<numPts;++i){
            double x,y;
            std::cin >> x >> y;
            pts.push_back(pt(x,y));
        }
        std::sort(pts.begin(),pts.end(),[](pt a, pt b){
            return a.x < b.x || (a.x == b.x && a.y < b.y);
            });  

        
        int ans = 0;
        for(auto it = pts.begin();it != pts.end();++it){
            for(auto it2 = it + 1;it2 != pts.end();++it2){
                if(dist_square(*it,*it2) <= 4){
                    ans = std::max(ans,solve(*it,*it2,pts));
                    ans = std::max(ans,solve(*it2,*it,pts));
                }
            }
        }
        std::cout << ans << "\n";
        std::cin >> numPts;
    }  
    return 0;  
}