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 <array>
#include <string>
#include <stack>
#include <queue>
#include <deque>
#include <map>
#include <unordered_map>
#include <set>
#include <unordered_set>
#include <tuple>
#include <bitset>
#include <memory>
#include <cmath>
#include <algorithm>
#include <functional>
#include <iomanip>
#include <numeric>
#include <climits>
#include <cfloat>




struct Point;
struct Circle;
struct Angle;
struct Vector;


Point operator+(const Point& p, const Vector& v);
Point operator+(const Vector& v, const Point& p);
Vector operator-(const Point& to, const Point& from);
Vector operator*(const Vector& v, const double scala);
Vector operator*(const double scala, const Vector& v);
Angle operator+(const Angle& a, const Angle& b);
Angle operator-(const Angle& a, const Angle& b);

struct Point {
	double x, y;
	double distance_from(const Point& p) const;
};
struct Circle {
	Point center;
	double radius;
	bool contains(const Point& p) const;
};
struct Angle {
	double sin, cos;
	static Angle from_cos(const double cos);
	Angle operator-() const;
};
struct Vector {
	double dx, dy;
	Vector rotate(const Angle& angle) const;
	double length() const;
};

int main() {
	while (true) {
		int n; std::cin >> n; if (n == 0) break;
		std::vector<Point> points(n); for (auto& p : points) std::cin >> p.x >> p.y;
		int max_covered_point = 1;
		for (auto i = 1; i < n; ++i) {
			const auto to = points[i];
			for (auto j = 0; j < i; ++j) if (to.distance_from(points[j]) <= 2) {
				const auto from = points[j];
				const auto angle = Angle::from_cos(to.distance_from(from) / 2);
				const auto vec = (to - from) * (1 / from.distance_from(to));
				const auto circle_a = Circle{ from + vec.rotate(angle), 1 };
				const auto circle_b = Circle{ from + vec.rotate(-angle), 1 };
				int count_a = 2, count_b = 2;
				for (auto k = 0; k < n; ++k) if (k != i && k != j) {
					if (circle_a.contains(points[k])) ++count_a;
					if (circle_b.contains(points[k])) ++count_b;
				}
				max_covered_point = std::max({ max_covered_point, count_a, count_b });
			}
		}
		std::cout << max_covered_point << '\n';
	}
}
Vector Vector::rotate(const Angle& angle) const
{
	return Vector{ angle.cos * dx - angle.sin * dy, angle.sin * dx + angle.cos * dy };
}

double Vector::length() const
{
	return std::sqrt(dx * dx + dy * dy);
}

Point operator+(const Point& p, const Vector& v)
{
	return Point{ p.x + v.dx, p.y + v.dy };
}

Point operator+(const Vector& v, const Point& p)
{
	return Point{ p.x + v.dx, p.y + v.dy };
}

Vector operator-(const Point& to, const Point& from)
{
	return Vector{ to.x - from.x, to.y - from.y };
}

Vector operator*(const Vector& v, const double scala)
{
	return Vector{ v.dx * scala, v.dy * scala };
}

Vector operator*(const double scala, const Vector& v)
{
	return Vector{ v.dx * scala, v.dy * scala };
}

Angle operator+(const Angle& a, const Angle& b)
{
	return Angle{ a.sin * b.cos + a.cos * b.sin, a.cos * b.cos - a.sin * b.sin };
}

Angle operator-(const Angle& a, const Angle& b)
{
	return Angle{ a.sin * b.cos - a.cos * b.sin, a.cos * b.cos + a.sin * b.sin };
}

double Point::distance_from(const Point& p) const
{
	return (*this - p).length();
}

bool Circle::contains(const Point& p) const
{
	return center.distance_from(p) <= radius;
}

Angle Angle::from_cos(const double cos)
{
	return Angle{ std::sqrt(1 - cos * cos), cos };
}

Angle Angle::operator-() const
{
	return Angle{ -sin, cos };
}

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

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

P center(P a, P b){
  P m  = (a+b) / 2.0;         // p[i]Æp[j]Ì_
  P n = (a-b) * P(0.0,1.0);  // @üxNg
  n = n / abs(n);
  
  double x = abs( a-b )/2.0;
  return m + n * sqrt( 1.0-x*x );     // ~ÌS
}

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 c = center( p[i],p[j] );

	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++: #include <cstdio>
#include <iostream>
#include <algorithm>
#include <string>
#include <cstring>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <cmath>
#include <iomanip>
#include <complex>
using namespace std;

//#define double long double
typedef complex<double> P;
#define rep(i, n) for (int i=0; i<(n); i++)
#define all(c) (c).begin(), (c).end()
#define uniq(c) c.erase(unique(all(c)), (c).end())
#define _1 first
#define _2 second
#define pb push_back
#define INF 1145141919
#define MOD 1000000007
const int R = 100000;

int N;
P X[300];

int f(P x, int a, int b) {
  int ctr = 2;
  rep(i, N) {
    if (i == a || i == b) continue;
    P p = X[i];
    if (abs(x-p) <= 1.0) ctr++;
  }
  return ctr;
}

signed main() {
  ios::sync_with_stdio(false); cin.tie(0);
  while (cin >> N && N != 0) {
    rep(i, N) {
      double x, y;
      cin >> x >> y;
      X[i] = P(x, y);
    }
    int ans = 1;
    rep(i, N) {
      rep(j, i) {
        P a = X[i], b = X[j];
        P base = (a+b)/(double)2.0;
        P v = b-a;
        double d = abs(v) / 2.0;
        if (d > 1) continue;
        P h1 = ((v*P(0, +1)) / abs(v)) * sqrt(1.0-d*d);
        P h2 = ((v*P(0, -1)) / abs(v)) * sqrt(1.0-d*d);
        ans = max(ans, f(base+h1, i, j));
        ans = max(ans, f(base+h2, i, j));
      }
    }
    cout << ans << "\n";
  }
  return 0;
}