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 <utility>
#include <cmath>
#include <algorithm>

using namespace std;
typedef pair<double, double> pt;
const double EPS = 1e-4;
#define x first
#define y second

double sq(double a) { return a * a; }

double dist(pt p, pt q) {
  return sqrt(sq(p.x - q.x) + sq(p.y - q.y));
}

bool lessthan(double a, double b) {
    if (a - b <= EPS)
        return true;
    return false;
}

pt midpoint(pt a, pt b) {
    return pt((a.x + b.x)/2, (a.y + b.y)/2);
}

pt operator+(pt a, pt b) {
    return pt(a.x + b.x, a.y + b.y);
}

pt operator-(pt a, pt b) {
  return pt(a.x - b.x, a.y - b.y);
}

pt operator*(pt a, double b) {
    return pt(a.x * b, a.y * b);
}

pt rotate(pt p, double a) {
  pt res;
  res.x = p.x * cos(a) - p.y * sin(a);
  res.y = p.x * sin(a) + p.y * cos(a);
  return res;
}



const int maxn = 300;
pt pts[maxn + 1];

int main() {
    int N;
    
    cin >> N;

    int maxcount = -1;
    while (N != 0) {
        if (maxcount != -1)
            cout << maxcount << endl;
        for (int i = 0; i < N; ++i) {
            double x, y;
            cin >> x >> y;
            pts[i] = pt(x, y);
        }
        maxcount = 0;
        for (int i = 0; i < N ; i++){
            for (int j = 0; j < N; j++) {
                if (j == i) continue;
                // for every pair of points
                if (dist(pts[i], pts[j]) > 2) continue;
                pt A = pts[i];
                pt B = pts[j];
                cout << "for points " << A.x << " " << A.y << " b: " << B.x << " "<< B.y << endl;
                pt M = midpoint(A, B);
                //cout << " M is " << M.x << " " << M.y << endl;
                pt BM = M - B;
                double distBM = dist(B,M);
                //cout << "dist BM " << distBM << endl; 
                double alpha = acos(distBM);
                pt BK = rotate(BM, alpha);
                pt BC = BK * (1 / distBM);
                pt C = B + BC;
                cout << "C is  "<< C.x << " " << C.y << endl;
                int count = 2;
                for (int k = 0; k < N; k++) {
                    if (k == i || k == j) continue;
                    if (lessthan(dist(C, pts[k]), 1)) count++; 
                }
                maxcount = max(count, maxcount);
            }
        }
        
        cin >> N;
    }
    cout << maxcount;

    return 0;
}

C++: #include <iostream>
#include <complex>
#include <cmath>
#include <algorithm>
using namespace std;

#define EPS 1.0e-8
typedef complex<double> P;

int main(){
  int n;
  P p[301];

  while( cin>>n && n ){
    for( int i=0;i<n;i++ )
      cin >> p[i].real() >> p[i].imag();

    int ans=0;
    for( int i=0;i<n;i++ ){
      for( int j=0;j<n;j++ ){
	if( i==j ) 	  continue;
	if( abs( p[i]-p[j] )+EPS > 2.0 ) continue;

	P m  = (p[i]+p[j]) / 2.0;         // p[i]Æp[j]Ì_
	P nn = (p[i]-p[j]) * P(0.0,1.0);  // @üxNg
	nn = nn / abs(nn);

	double x = abs( p[i]-p[j] )/2.0;
	P c = m + nn*sqrt( 1.0-x*x );     // ~ÌS

	int num=0;
	for( int k=0;k<n;k++ )
	  if( abs( c-p[k] ) < 1.0 + EPS  )
	    num++;
	ans = max( ans,num );
      }
    }
    cout << ans << endl;
  }

  return 0;
}

C++: #define _USE_MATH_DEFINES
#define INF 100000000

#include <iostream>
#include <sstream>
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <queue>
#include <stack>
#include <limits>
#include <map>
#include <string>
#include <cstring>
#include <set>
#include <deque>
#include <bitset>
#include <list>
#include <cctype>
#include <utility>

using namespace std;

typedef long long ll;
typedef pair <int,int> P;
typedef pair <int,P> PP;

static const double EPS = 1e-8;

int tx[] = {0,1,0,-1};
int ty[] = {-1,0,1,0};

class Point{
public:
	double x;
	double y;
	Point() : x(0),y(0){}
	Point(double _x,double _y){
		x = _x;
		y = _y;
	}
};

class Region{
public:
	int larger_num;
	int smaller_num;
	double size;
	double x;
	double y;
	Region(double _x,double _y,double _size,int _larger_num,int _smaller_num){
		x = _x;
		y = _y;
		size = _size;
		larger_num = _larger_num;
		smaller_num = _smaller_num;
	}

	bool operator<(const Region& r) const{
		return larger_num < r.larger_num;
	}
};

void save(double x,double y,double size,const vector<Point>& points,priority_queue<Region>& que,int maxv){
	int larger_num=0;
	int smaller_num=0;
	double cx = x + (size / 2.0);
	double cy = y + (size / 2.0);

	for(int i=0;i<points.size();i++){
		double d = sqrt((points[i].x - cx) * (points[i].x - cx)
			+ (points[i].y - cy) * (points[i].y - cy));
		if(d <= 1.0) smaller_num++;
		if(d <= 1.0 + (size/2.0)*sqrt(2.0)) larger_num++;
	}

	if(maxv >= larger_num) return;

	que.push(Region(x,y,size,larger_num,smaller_num));
}

int main(){
	int n;
	while(~scanf("%d",&n)){
		if(n==0) break;
		vector<Point> points;
		for(int i=0;i<n;i++){
			double x,y;
			scanf("%lf %lf",&x,&y);
			points.push_back(Point(x,y));
		}

		priority_queue<Region> que;

		int maxv = 1;
		save(0,0,10,points,que,maxv);

		while(!que.empty()){
			Region re = que.top();
			int larger = re.larger_num;
			int smaller = re.smaller_num;
			double size = re.size;
			double x = re.x;
			double y = re.y;

			que.pop();

			if(maxv > larger) continue;
			maxv = max(smaller,maxv);
		
			save(x,y,size/2.0,points,que,maxv);
			save(x+size/2.0,y,size/2.0,points,que,maxv);
			save(x,y+size/2.0,size/2.0,points,que,maxv);
			save(x+size/2.0,y+size/2.0,size/2.0,points,que,maxv);
		}

		printf("%d\n",maxv);
	}
}

C++: #include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <string>
#include <deque>
#include <utility>

#define REP(i, n) FOR(i, 0, n)
#define FOR(i, a, b) for(int i = int(a); i < int(b); ++i)
#define RREP(i, n) RFOR(i, 0, n)
#define RFOR(i, a, b) for(int i = int(b) - 1; i >= int(a); --i)

constexpr double EPS = 1e-7;

signed main() {

	while (true) {
		int n;
		std::cin >> n;
		if (n == 0) break;

		std::vector<double> x(n), y(n);
		REP(i, n) std::cin >> x[i] >> y[i];

		int ans = 0;

		REP(i, n) FOR(j, i + 1, n) {
			double dist = std::hypot(x[i] - x[j], y[i] - y[j]);
			if (dist > 2 + EPS) continue;

			// iとjの中点m
			double mx = (x[i] + x[j]) / 2;
			double my = (y[i] + y[j]) / 2;

			// 線分ijから中心への距離
			double r = std::sqrt(1 - dist / 2);

			// 中心c
			for (int sign = -1; sign <= 1; sign += 2) {
				// 移動のベクトルv
				double vx = sign * r * (y[j] - y[i]) / dist;
				double vy = sign * r * -(x[j] - x[i]) / dist;

				double cx = mx + vx;
				double cy = my + vy;

				int cnt = 0;
				REP(k, n) {
					if (std::hypot(cx - x[k], cy - y[k]) < 1 + EPS) ++cnt;
				}
				ans = std::max(ans, cnt);
			}
		}
		std::cout << ans << "\n";
	}

	return 0;
}