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;
typedef long long ll;
typedef pair<int,int> P;
typedef pair<int,P> P1;
typedef pair<P,P> P2;
#define pu push
#define pb push_back
#define mp make_pair
#define INF 1000000000
#define fi first
#define sc second
#define rep(i,x) for(int i=0;i<x;i++)
#define SORT(x) sort(x.begin(),x.end())
#define ERASE(x) x.erase(unique(x.begin(),x.end()),x.end())
#define POSL(x,v) (lower_bound(x.begin(),x.end(),v)-x.begin())
#define POSU(x,v) (upper_bound(x.begin(),x.end(),v)-x.begin())
typedef complex<double> pt;
typedef pair<pt,pt> L;
typedef vector<P> poly;
const double EPS = 1e-11;
#define x real()
#define y imag()
 
double dot(pt a,pt b){
	return (conj(a)*b).x;
}
double cross(pt a,pt b){
	return (conj(a)*b).y;
}
int ccw(pt a,pt b,pt c){
	b -= a; c -= a;
	if(cross(b,c) > EPS) return 1; // counter clockwise
	if(cross(b,c) < -EPS) return -1; // clockwise
	if(dot(b,c) < -EPS) return 2; //c-a-b
	if(norm(b) < norm(c)) return -2; //a-b-c
	return 0; //a-c-b
}
bool cmp(const pt& a,const pt& b){
	if(-EPS < a.x-b.x && a.x-b.x < EPS) return a.y < b.y;
	else return a.x < b.x;
}
vector<pt>convex_hull(vector<pt>ps)
{
	sort(ps.begin(),ps.end(),cmp);
	int k=0,n = ps.size();
	vector<pt>qs(n*2);
	for(int i=0;i<n;i++)
	{
		while(k>1 && ccw(qs[k-2],qs[k-1],ps[i]) == -1) k--;
		qs[k++] = ps[i];
	}
	for(int i=n-2,t=k;i>=0;i--)
	{
		while(k>t && ccw(qs[k-2],qs[k-1],ps[i]) == -1) k--;
		qs[k++] = ps[i];
	}
	qs.resize(k-1);
	return qs;
}
pt intersect(pt a,pt b,pt c,pt d){
	double dev = (b.y-a.y)*(d.x-c.x)-(b.x-a.x)*(d.y-c.y);
	
	double d1, d2;
 
	d1 = (c.y*d.x-c.x*d.y);
	d2 = (a.y*b.x-a.x*b.y);
 
	pt ap1;
	ap1.x = d1*(b.x-a.x) - d2*(d.x-c.x);
	ap1.x /= dev;
	ap1.y = d1*(b.y-a.y) - d2*(d.y-c.y);
	ap1.y /= dev;
 
	return ap1;
}
bool par(pt a,pt b,pt c,pt d){
	b.x -= a.x; b.y -= a.y;
	d.x -= c.x; d.y -= c.y;
	if(abs(b.y*d.x - b.x*d.y) < EPS){
		return 1;
	}
	return 0;
}
int main(){
	while(1){
		int n; scanf("%d",&n); if(!n) return 0; if(n==1){ puts("1"); continue;}
		pt p[305];
		for(int i=1;i<=n;i++){
			double a,b; scanf("%lf%lf",&a,&b);
			p[i] = pt(a,b);
		}
		int ans = 0;
		for(int i=1;i<=n;i++){
			for(int j=i+1;j<=n;j++){
				pt q = (p[i]+p[j])/2.0;
				pt r = (p[i]-p[j]); if(abs(r) > 2.0+EPS) continue;
				pt p1 = q + r / abs(r) * pt(0,sqrt(1.0-abs(r)*abs(r)/4.0));
				pt p2 = q - r / abs(r) * pt(0,sqrt(1.0-abs(r)*abs(r)/4.0));
				int c=0,d=0;
				for(int k=1;k<=n;k++){
					if(abs(p1-p[k]) <= 1.0+EPS) c++;
					if(abs(p2-p[k]) <= 1.0+EPS) d++;
				}
				ans = max(ans,max(c,d));
			}
		}
		printf("%d\n",ans);
	}
}

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 <complex>
#include <cmath>
#define EPS 1e-9
using namespace std;
typedef complex<double> xy_t;
int n;
xy_t points[301];
int count_in(xy_t center)
{
  int ans=0;
  for(int i=0;i<n;i++)
  {
    if (abs(center-points[i])<1.0+EPS)
      ans++;
  }
  return ans;
}
xy_t ppdcl(xy_t a)
{
  //return xy_t(a.imag(),-a.real());
  return a*xy_t(0,1);
}
xy_t delta_v(xy_t C1,xy_t C2)
{
  xy_t half=(C1-C2)*0.5;
  double l=abs(half);
  return ppdcl(half)*sqrt(1-l*l)*(1.0/l);
}
int main()
{
  while(1)
  {
    cin>>n;
    if (n==0) break;
    double x,y;
    for(int i=0;i<n;i++)
    {
      cin>>x>>y;
      points[i]=xy_t(x,y);
    }
    int maxn=-1;
    for(int i=0;i<n;i++)
    {
      for(int j=i+1;j<n;j++)
      {
        if(abs(points[i]-points[j])<2.0)
        {
          xy_t mid_p=(points[i]+points[j])*0.5;
          maxn=max(maxn,count_in(mid_p+delta_v(points[i],points[j])));
          maxn=max(maxn,count_in(mid_p-delta_v(points[i],points[j])));
        }
      }
    }
    cout<<maxn<<endl;
  }

}