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 <sstream>
#include <string>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <map>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cassert>

using namespace std;

#define FOR(i,k,n) for(int i=(k); i<(int)(n); ++i)
#define REP(i,n) FOR(i,0,n)
#define FORIT(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();++i)

template<class T> void debug(T begin, T end){ for(T i = begin; i != end; ++i) cerr<<*i<<" "; cerr<<endl; }
inline bool valid(int x, int y, int W, int H){ return (x >= 0 && y >= 0 && x < W && y < H); }

typedef long long ll;
const int INF = 100000000;
const double EPS = 1e-8;
const int MOD = 1000000007;
int dx[8] = {1, 0, -1, 0, 1, -1, -1, 1};
int dy[8] = {0, 1, 0, -1, 1, 1, -1, -1};
#include <complex>
typedef complex<double> Point;
int main(){
  int N;
  while(cin >> N && N){
    vector<Point> ps(N);
    REP(i, N) cin >> ps[i].real() >> ps[i].imag();
    int ans = 1;
    REP(i, N)REP(j, N) if(i != j){
      double d = abs(ps[i] - ps[j]) / 2;
      if(d > 1) continue;
      double t = sqrt(1 - d * d);
      Point start = ps[i] + (ps[j] - ps[i]) * 0.5;
      Point vect = (ps[j] - ps[i]) * Point(0.0, 1.0);
      vect /= abs(vect);
      for(int sign = -1; sign <= 1; sign += 2){
        Point center = start + vect * t * (double)sign;
        int sub_ans = 0;
        REP(k, N) if(abs(center - ps[k]) < 1.0 + EPS) sub_ans ++;
        ans = max(ans, sub_ans);
      }
    }
    cout << ans << endl;
  }
  return 0;
}

C++: #include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace std;
class P{
public:
	double x,y;
	P(){}
	P(double xx,double yy){
		x=xx;
		y=yy;
	}
	P operator +(P p){
		return P(x+p.x,y+p.y);
	}
	P operator -(P p){
		return P(x-p.x,y-p.y);
	}
	P operator /(double d){
		return P(x/d,y/d);
	}
};

int n;
P p[301];

int main(void){
	while(1){
		scanf("%d",&n);
		if(n==0)break;
		for(int i=0;i<n;i++){
			scanf("%lf%lf",&p[i].x,&p[i].y);
		}
		int res=1;
		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				P np=p[j]-p[i];
				//printf("%lf %lf\n",np.x,np.y);
				double blen=np.x*np.x/4.0+np.y*np.y/4.0;
				if(blen>1.0)continue;
				double len=sqrt(1-blen);
				double dir=atan2(np.y,np.x);
				P cp=(p[i]+p[j])/2+P(len*cos(dir+M_PI/2),len*sin(dir+M_PI/2));
				int cnt=0;
				for(int k=0;k<n;k++){
					double len2=sqrt((cp.x-p[k].x)*(cp.x-p[k].x)+(cp.y-p[k].y)*(cp.y-p[k].y));
					if(len2<=1.0001)cnt++;
				}
				res=max(res,cnt);
				cp=(p[i]+p[j])/2+P(len*cos(dir-M_PI/2),len*sin(dir-M_PI/2));
				cnt=0;
				for(int k=0;k<n;k++){
					double len2=sqrt((cp.x-p[k].x)*(cp.x-p[k].x)+(cp.y-p[k].y)*(cp.y-p[k].y));
					if(len2<=1.0001)cnt++;
				}
				res=max(res,cnt);
			}
		}
		printf("%d\n",res);
	}
	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 = 0;
		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=0; 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/2.0)*(d1/2.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 || Dis(cx, cy, p[k].first, 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));
}