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 <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;
	}
}

C++: #include<iostream>
#include<cmath>
#include<vector>
#include<algorithm>
#include<climits>
#include<cfloat>

using namespace std;

double EPS=1e-10;

double add(double a,double b){
  if(abs(a+b)<EPS * (abs(a)+abs(b)))return 0;
  return a+b;
}

struct point{
  double x,y;
  point(){}
  point(double x,double y) : x(x) , y(y) {}

  point operator + (point p) {
    return point(add(x,p.x),add(y,p.y));
  }

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

  point operator * (double d) {
    return point(x*d,y*d);
  }

  point operator / (double d) {
    return point(x/d,y/d);
  }

  bool operator == (const point &p)const{
    return abs(x-p.x)<EPS && abs(y-p.y)<EPS;
  }
};

struct circle{
  point p;
  double r;
  circle(){}
  circle(point p,double r):p(p),r(r){}
};

double dist(point a, point b){
  return sqrt(pow(a.x-b.x,2)+pow(a.y-b.y,2));
}

vector<point> circle_cross(circle a, circle b){
  vector<point>ret;

  double x=a.p.x-b.p.x,y=a.p.y-b.p.y;
  double L=x*x+y*y;

  if(L>pow(a.r+b.r,2))return ret;

  L=sqrt(L);
  double C=atan2(b.p.y-a.p.y,b.p.x-a.p.x);
  double alpha=acos((L*L+a.r*a.r-b.r*b.r)/(2*L*a.r));

  point r1,r2;
  r1.x=a.p.x+cos(C+alpha),r1.y=a.p.y+sin(C+alpha);
  r2.x=a.p.x+cos(C-alpha),r2.y=a.p.y+sin(C-alpha);

  ret.push_back(r1),ret.push_back(r2);
  return ret;
}

int main(void){

  int n;
  point a;
  vector<point>P;
  vector<circle>C;

  while(cin >> n,n){
    P.clear();
    C.clear();
    for(int i=0;i<n;i++){
      cin >> a.x >> a.y;
      P.push_back(a);
    }

    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
	if(dist(P[i],P[j])>2.0+EPS)continue;
	vector<point>res=circle_cross(circle(P[i],1.0),circle(P[j],1.0));
	for(int k=0;k<res.size();k++){
	  C.push_back(circle(res[k],1));
	}
      }
    }

    int ans=1;
    for(int i=0;i<C.size();i++){
      int cnt=0;
      for(int j=0;j<P.size();j++)
	if(dist(C[i].p,P[j])<1.0+EPS)cnt++;
	
      ans=max(ans,cnt);
    }

    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 = 1;
		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=i+1; 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*d1/4.0));

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

				int cnt = 0;
				for(int k=0; k<n; k++){
					if(k == i || k == j || (cx-p[k].first)*(cx-p[k].first) + (cy-p[k].second)*(cy-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));
}