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 <cmath>

using namespace std;
const double INF = 1e10;
const double EPS = 1e-8;

struct Point {
	double x, y;
	Point() : x(0), y(0) { }
	Point(double x, double y) : x(x), y(y) { }
	Point operator+(const Point &p) const { return Point(x + p.x, y + p.y); }
	Point operator-(const Point &p) const { return Point(x - p.x, y - p.y); }
	Point operator*(double s) const { return Point(x * s, y * s); }
	Point operator/(double s) const { return Point(x / s, y / s); }
	Point &operator+=(const Point &p){ x += p.x; y += p.y; return *this; }
	Point &operator-=(const Point &p){ x -= p.x; y -= p.y; return *this; }
	Point &operator*=(double s){ x *= s; y *= s; return *this; }
	Point &operator/=(double s){ x /= s; y /= s; return *this; }
	bool operator==(const Point &p) const { return x == p.x && y == p.y; }
	bool operator!=(const Point &p) const { return x != p.x || y != p.y; }
	bool operator<(const Point &p) const {
		return (x == p.x) ? (y < p.y) : (x < p.x);
	}
};
inline Point operator*(double s, const Point &p){ return p * s; }

inline double norm(const Point &p){ return p.x * p.x + p.y * p.y; }
inline double abs(const Point &p){ return sqrt(norm(p)); }
inline double arg(const Point &p){ return atan2(p.y, p.x); }
inline Point unit(const Point &p){ return p / abs(p); }
inline Point ortho(const Point &p){ return Point(-p.y, p.x); }

int main(){
	ios_base::sync_with_stdio(false);
	while(true){
		int n;
		cin >> n;
		if(n == 0){ break; }
		vector<Point> p(n);
		for(int i = 0; i < n; ++i){ cin >> p[i].x >> p[i].y; }
		const double limit = (1.0 + EPS) * (1.0 + EPS);
		int answer = 1;
		for(int i = 0; i < n; ++i){
			const Point a = p[i];
			for(int j = i + 1; j < n; ++j){
				const Point b = p[j];
				if(abs(b - a) > 2.0){ continue; }
				const Point c = (a + b) * 0.5;
				const Point v = unit(ortho(b - a));
				const double t = sqrt(1.0 - norm(b - c));
				const Point q[2] = { c + v * t, c + v * -t };
				for(int k = 0; k < 2; ++k){
					int count = 0;
					for(int l = 0; l < n; ++l){
						if(norm(p[l] - q[k]) <= limit){ ++count; }
					}
					answer = max(answer, count);
				}
			}
		}
		cout << answer << endl;
	}
	return 0;
}

C++: #include <iostream>
#include <algorithm>
#include <cmath>
#define COS45 0.707107
using namespace std;

class Point {
public:
  double x, y;
  Point(double x = 0, double y = 0) : x(x), y(y) {}
  Point operator +(Point p) {
    return Point(x+p.x, y+p.y);
  }
  Point operator -(Point p) {
    return Point(x-p.x, y-p.y);
  }
  Point operator *(double a) {
    return Point(x*a, y*a);
  }
  Point operator /(double a) {
    return Point(x/a, y/a);
  }
  bool operator ==(const Point &p) const {
    return x == p.x && y == p.y;
  }
  bool operator <(const Point &p) const {
    double l1 = x+y;
    double l2 = p.x+p.y;
    if(l1 != l2) return l1 < l2;
    if(y != p.y) return y < p.y;
    return x < p.x;
  }
};

typedef Point Vector;

double norm(Point p1, Point p2) {
  return (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y);
}

int n;
Point p[300];

int calc(Point p1, Point p2) {
  Point m = Point((p1.x+p2.x)/2.0, (p1.y+p2.y)/2.0);
  Point o;
  Vector moto = Vector(p2.y-p1.y, p1.x-p2.x);
  int cnt = 0;
  int newCnt;
  double p1p2 = norm(p1,p2);
  if(p1p2 > 4) return 1;
  double l = sqrt(4.0/p1p2-1)/2.0;

  for(int t = 0; t < 2; t++) {
    newCnt = 0;
    if(t == 1) {
      moto = moto * (-1);
    }
    o = m + moto * l;

    for(int i = 0; i < n; i++) {
      if(p[i].x > o.x + COS45 && p[i].y > o.y + COS45) break;
      if(newCnt + n-i <= cnt) break;
      if(p[i] == p1 || p[i] == p2) continue;
      if((p[i].x-o.x)*(p[i].x-o.x) + (p[i].y-o.y)*(p[i].y-o.y) <= 1.0) {
	newCnt++;
      }
    }
    if(newCnt > cnt) cnt = newCnt;
  }
  return cnt + 2;
}

main() {
  while(1) {
    int ans;
    cin >> n;
    if(n == 0) break;
    for(int i = 0; i < n; i++) {
      cin >> p[i].x >> p[i].y;
    }
    sort(p, &p[n]);
    ans = -1;
    for(int i = 0; i < n; i++) {
      for(int j = i+1; j < n; j++) {
	ans = max(ans, calc(p[i],p[j]));
      }
    }
    cout << ans << endl;
  }
}

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

using namespace std;

const double EPS = 1e-8;

bool equals(double x, double y) {
	return fabs(x-y) < EPS;
}

#define x first
#define y second
typedef pair<double, double> pt;


pt operator-(pt u, pt v) { return pt(u.x - v.x, u.y - v.y); }
pt operator+(pt u, pt v) { return pt(u.x + v.x, u.y + v.y); }
pt operator*(pt u, double d) { return pt(u.x * d, u.y * d); }

// dot product
double operator*(pt u, pt v) { return u.x * v.x + u.y * v.y; }
// norm
double operator!(pt p) { return sqrt(p * p); }
// "cross product"
double operator^(pt u, pt v) { return u.x * v.y - u.y * v.x; }

double operator==(pt a, pt b) { return equals(a.x, b.x) && equals(a.y, b.y);}

// angle with positive x axis
double operator~(pt v)
{
	return atan2(v.y, v.x);
}

int pointsInCircle(pt p, vector<pt>& points)
{
	vector<pair<double, bool>> angles;
	
	for(auto point : points)
	{
		pt x = point - p;
		double dSquared = x * x;
		if(p == point || dSquared > 4) continue;
		double a = acos(sqrt(dSquared) / 2.0); 
		double b = ~x;
		angles.emplace_back(b - a, true);
		angles.emplace_back(a + b, false);
	}
	
	int count = 1;
	int result = 1;
	sort(angles.begin(), angles.end());
	for(auto a : angles)
	{
		if(a.second) 
		{
			count++;
			result = max(result, count);
		}
		else count--;
	}
	return result;
}

int solve(vector<pt>& points)
{
	int sln = 0;
	for(auto p : points)
		sln = max(sln, pointsInCircle(p, points));
	return sln;
}

int main(void) {
	
	ios::sync_with_stdio(false);
	int N;
	cin >> N;

	while(N)
	{
		vector<pt> points;
		for(int i  = 0; i < N; ++i)
		{
			double x, y;
			cin >> x >> y;
			points.emplace_back(x, y);
		}
		cout << solve(points) << '\n';
		cin >> N;
	}
	return 0;
}