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<math.h>
#include<vector>
using namespace std;

// 点
class Point {
  public:
    float x, y;
    Point(double ix = 0.0, double iy = 0.0) {
      set(ix, iy);
    }
    void set(double ix = 0.0, double iy = 0.0) {
      x = ix;
      y = iy;
    }
    void copy(Point p) {
      x = p.x;
      y = p.y;
    }
    double length() {
      return sqrt(x * x + y * y);
    }
    double distanceTo(Point p) {
      double dx = p.x - x, dy = p.y - y;
      return sqrt(dx * dx + dy * dy);
    }
    void setCentroidOf(Point p1, Point p2) {
      x = (p1.x + p2.x) / 2;
      y = (p1.y + p2.y) / 2;
    }
};

int main() {
  vector<int> results;
  while (true) {
    // 入力
    int n;
    cin >> n;
    if (n == 0) break;
    vector<Point> p(n);
    for (int i=0; i<n; i++) {
      double x, y;
      cin >> x >> y;
      p[i].set(x, y);
    }
    
    const int nC2 = n * (n - 1) / 2 + n;
    vector<Point> px(n); // 現時点で有力な点
    for (int i=0; i<n; i++) px[i].copy(p[i]);
    vector<Point> cp(nC2); // 中点たち
    int previousMaxScore = 0; // 前回のループでの最高スコア
    int maxScoreConst = 0; // 同じ最高スコアが連続して記録された回数
    while (true) {
      // 与えられた点の中点たち（+もとの点）を候補点とする
      int cpi = 0;
      for (int i=0; i<n-1; i++) {
	for (int j=i+1; j<n; j++) {
	  cp[cpi++].setCentroidOf(px[i], px[j]);
	}
      }
      for (int i=0; i<n; i++) cp[cpi++].copy(px[i]);
      // 各点のスコアを算出
      vector<int> scores(nC2);
      for (int i=0; i<nC2; i++) {
	scores[i] = 0;
	for (int j=0; j<n; j++) {
	  if (cp[i].distanceTo(p[j]) <= 1) scores[i]++;
	}
      }
      // cpをスコア降順にバブルソート
      for (int i=0; i<nC2; i++) {
	for (int j=0; j<nC2-1; j++) {
	  if (scores[j] < scores[j + 1]) {
	    double jx = cp[j].x, jy = cp[j].y;
	    cp[j].copy(cp[j + 1]);
	    cp[j + 1].set(jx, jy);
	    int js = scores[j];
	    scores[j] = scores[j + 1];
	    scores[j + 1] = js;
	  }
	}
      }
      // このループの最大スコア
      int maxScore = scores[0];
      //cout << maxScore << " -> ";
      if (maxScore == previousMaxScore) {
	// 最大スコアが3回連続で更新されなかったら確定
	previousMaxScore = maxScore;
	if (++maxScoreConst >= 3) break;
      } else {
	previousMaxScore = maxScore;
	maxScoreConst = 0;
      }
      // スコアが上位nの点を次のループの中点の素として採用
      for (int i=0; i<n; i++) px[i].copy(cp[i]);
    }
    
    results.push_back(previousMaxScore);
    //cout << endl;
  }
  
  int size = results.size();
  for (int i=0; i<size; i++) {
    cout << results[i] << 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=0;
		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 <iostream>
#include <vector>
#include <utility>
#include <cmath>

using namespace std;

typedef pair<double,double> pos;
const double r = 1;
const double eps = 0.0000000001;
int count(vector<pos> ps,double x,double y,int ban1,int ban2){
	int c=0;//2;
	for(int i=0;i<ps.size();i++){
		if(//i!=ban1&&i!=ban2&&
				pow(ps[i].first-x,2)+pow(ps[i].second-y,2)<=pow(1.0+eps,2)){
			c++;
		}
	}
	return c;
}

bool range(pos p1,pos p2,double r){
	return pow(p1.first-p2.first,2)+pow(p1.second-p2.second,2) <= r*r;
}

double len(pos p){
	return sqrt(pow(p.first,2)+pow(p.second,2));
}

double len2(pos p){
	return (pow(p.first,2)+pow(p.second,2));
}

vector<pos> point2(pos p1,pos p2){
	pos m = make_pair((p1.first+p2.first)/2,(p1.second+p2.second)/2);
	pos vec = make_pair(p1.first-p2.first,p1.second-p2.second);
	pos q = m;//make_pair(p2.first+vec.first/2,p2.second+vec.second/2);
	double l = sqrt(1-pow(len(vec)/2.0,2));
	double ll = len(vec);
	pos vec2 = make_pair(-vec.second,vec.first);
	vec2.first=vec2.first/ll*l;
	vec2.second = vec2.second/ll*l;
	vector<pos> ret;
	ret.push_back(make_pair(q.first+vec2.first,q.second+vec2.second));
	ret.push_back(make_pair(q.first-vec2.first,q.second-vec2.second));
	return ret;
}

int solve(vector<pos> ps){
	int mx=1;
	for(int i=0;i<ps.size()-1;i++)
		for(int k=i+1;k<ps.size();k++){
			if(range(ps[i],ps[k],2)){
				//cerr << "*****" << i << " " << k << endl;
				vector<pos> vp = point2(ps[i],ps[k]);
				mx = max(mx,count(ps,vp[0].first,vp[0].second,i,k));
				mx = max(mx,count(ps,vp[1].first,vp[1].second,i,k));
			}
		}
	return mx;
}

int main(){
	//cerr << point2(make_pair(0,0),make_pair(1,1))[1].first << endl;
	while(true){
		int n;
		vector<pos> v;
		cin >> n;
		if(n==0)return 0;
		for(int i=0;i<n;i++){
			double x,y;
			cin >> x >> y;
			v.push_back(make_pair(x,y));
		}
		cout << solve(v) << endl;
	}
}