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<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<algorithm>
#include<functional>
#include<iomanip>
#include<queue>
#include<ciso646>
#include<random>
#include<map>
#include<set>
#include<complex>
using namespace std;
typedef long long ll;
const ll MOD = 1000000007;
const ll INF = (ll)1000000007 * 1000000007;
const double EPS = 1e-9;
typedef pair<int, int> P;
typedef unsigned int ui;
#define stop char nyaa;cin>>nyaa;
#define rep(i,n) for(int i=0;i<n;i++)
#define per(i,n) for(int i=n-1;i>=0;i--)
#define Rep(i,sta,n) for(int i=sta;i<n;i++)
#define rep1(i,n) for(int i=1;i<=n;i++)
#define per1(i,n) for(int i=n;i>=1;i--)
#define Rep1(i,sta,n) for(int i=sta;i<=n;i++)
typedef long double ld;
typedef complex<ld> Point;
const ld eps = 1e-8;
const ld pi = acos(-1.0);
bool eq(ld a, ld b) {
	return (abs(a - b) < eps);
}
ld dot(Point a, Point b) {
	return real(conj(a)*b);
}
ld cross(Point a, Point b) {
	return imag(conj(a)*b);
}
class Line {
public:
	Point a, b;
};
class Circle {
public:
	Point p;
	ld r;
};
vector<Point> is_cc(Circle c1, Circle c2) {
	vector<Point> res;
	ld d = abs(c1.p - c2.p);
	ld rc = (d*d + c1.r*c1.r - c2.r*c2.r) / (2 * d);
	ld dfr = c1.r*c1.r - rc * rc;
	if (abs(dfr) < eps)dfr = 0.0;
	if (dfr < 0.0)return res;
	ld rs = sqrt(dfr);
	Point diff = (c2.p - c1.p) / d;
	res.push_back(c1.p + diff * Point(rc, rs));
	if (dfr != 0.0)res.push_back(c1.p + diff * Point(rc, -rs));
	return res;
}
ld dis(Point x, Point y) {
	return sqrt(dot(y - x, y - x));
}
int main() {
	int n;
	while (cin >> n, n) {
		Point a[300];
		rep(i, n) {
			ld x, y; cin >> x >> y; a[i] = { x,y };
		}
		int ma = 1;
		rep(i, n) {
			Rep(j, i + 1, n) {
				Circle s = { a[i],(ld)1.0 };
				Circle t = { a[j],(ld)1.0 };
				vector<Point> v = is_cc(s, t);
				rep(l, (int)v.size()) {
					int cnt = 0;
					rep(k, n) {
						if (dis(a[k], v[l]) <= (ld)1.0 + eps)cnt++;
					}
					ma = max(ma, cnt);
				}
			}
		}
		cout << ma << 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;
	*maxv = max(*maxv,smaller_num);

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