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 <algorithm>
#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <queue>
#include <set>
#include <tuple>
#include <vector>

using namespace std;

#define rep(i, n) for (int64_t i = 0; i < (n); i++)
#define irep(i, n) for (int64_t i = 0; i <= (n); i++)
#define rrep(i, n) for (int64_t i = (n)-1; i >= 0; i--)
#define rirep(i, n) for (int64_t i = n; i >= 0; i--)

class Vec;
Vec operator/(const Vec& v, const double& divisor);
Vec operator*(const double& scale, const Vec& v);

class Circle;
class Line;
template <typename G, typename H>
std::pair<std::vector<Vec>, bool> intersection(const G& lh, const H& rh,
                                               double eps);
template <>
std::pair<std::vector<Vec>, bool> intersection<Circle, Circle>(
    const Circle& circle, const Circle& line, double eps);
template <>
std::pair<std::vector<Vec>, bool> intersection<Circle, Line>(
    const Circle& circle, const Line& line, double eps);

class Vec {
  using F = double;
  vector<F> container;

 public:
  explicit Vec(size_t size) : container(size) {}
  explicit Vec(initializer_list<F> elements)
      : container(elements.begin(), elements.end()) {}
  size_t dim() const { return container.size(); }
  F& operator[](const size_t size) { return container.at(size); }
  const F& operator[](const size_t size) const { return container.at(size); }
  Vec operator+(const Vec& other) const {
    Vec ret(*this);
    for (size_t i = 0; i < dim(); i++) ret[i] += other[i];
    return ret;
  }
  Vec operator-(const Vec& other) const { return (*this) + (-other); }
  Vec operator-() const { return -1.0 * (*this); }
  F inner(const Vec& other) const {
    F acc = 0;
    for (size_t i = 0; i < dim(); i++) {
      acc += (*this)[i] * other[i];
    }
    return acc;
  }
  Vec& operator+=(const Vec& other) {
    for (size_t i = 0; i < dim(); i++) {
      (*this)[i] += other[i];
    }
    return *this;
  }
  F norm() const { return sqrt(inner(*this)); }
  Vec unit() const { return (*this) / this->norm(); }
  F ccw(const Vec& other) const {
    return (*this)[0] * other[1] - other[0] * (*this)[1];
  }
  Vec normal() const { return Vec({-(*this)[1], (*this)[0]}); }
};
Vec operator/(const Vec& v, const double& divisor) { return 1.0 / divisor * v; }
Vec operator*(const double& scale, const Vec& v) {
  Vec ret(v);
  for (size_t i = 0; i < ret.dim(); i++) ret[i] *= scale;
  return ret;
}

class Circle {
  Vec m_center;
  double m_radius;

 public:
  Circle(const Vec& center, double radius)
      : m_center(center), m_radius(radius) {}
  const Vec& center() const { return m_center; }
  double radius() const { return m_radius; }
};

class Line {
  Vec m_grad;
  double m_bias;

 public:
  Line(const Vec& grad, double bias) : m_grad(grad), m_bias(bias) {}
  static Line fromPoints(const Vec& p0, const Vec& p1) {
    Vec d = p0 - p1;
    Vec g({-d[1], d[0]});
    return Line(g, g.inner(p0));
  }
  const Vec& grad() const { return m_grad; }
  double bias() const { return m_bias; }
  double dist(const Vec& p) const {
    return abs(m_grad.inner(p) - m_bias) / m_grad.norm();
  }
  Vec projection(const Vec& p) const {
    const Vec g = m_grad.unit();
    const Vec h({-g[1], g[0]});
    return m_bias / m_grad.norm() * g + h.inner(p) * h;
  }
};

template <>
pair<vector<Vec>, bool> intersection<Circle, Circle>(const Circle& c0,
                                                     const Circle& c1,
                                                     double eps) {
  const Vec& p0 = c0.center();
  const Vec& p1 = c1.center();
  double r0 = c0.radius();
  double r1 = c1.radius();

  const Vec v = p0 - p1;

  if (r0 + r1 < v.norm()) {
    return make_pair(vector<Vec>(), false);
  } else if (r0 > r1 + v.norm() || r1 > r0 + v.norm()) {
    return make_pair(vector<Vec>(), false);
  } else {
    Vec grad = 2.0 * (p1 - p0);
    const double bias = (r0 - r1) * (r0 + r1) -
                        (p0[0] - p1[0]) * (p0[0] + p1[0]) -
                        (p0[1] - p1[1]) * (p0[1] + p1[1]);
    return move(intersection(c0, Line(grad, bias), eps));
  }
}
template <>
pair<vector<Vec>, bool> intersection<Circle, Line>(const Circle& circle,
                                                   const Line& line,
                                                   double eps) {
  const double g = line.grad().norm();

  if (g < eps) {
    return make_pair(vector<Vec>(), abs(line.bias()) < eps);
  } else {
    vector<Vec> result;
    const double dist = (line.bias() - line.grad().inner(circle.center())) / g;

    const double det = pow(circle.radius(), 2) - pow(dist, 2);
    if (det < 0) {
    } else {
      const double s = sqrt(det);
      const Vec u = line.grad().unit();
      const Vec p = circle.center() + dist * u;
      Vec dir = u.normal();

      result.push_back(p + s * dir);
      result.push_back(p - s * dir);
    }

    return make_pair(move(result), false);
  }
}

void solve(int N) {
  vector<Vec> p;
  rep(i, N) {
    double x, y;
    cin >> x >> y;
    p.push_back(Vec({x, y}));
  }

  int result = 1;
  rep(i, N) for (int j = i + 1; j < N; j++) {
    const double EPS = 1e-9;

    vector<Vec> s;
    bool eq;
    tie(s, eq) = intersection(Circle(p[i], 1), Circle(p[j], 1), EPS);
    for (auto c : s) {
      int cnt = 0;
      rep(k, N) {
        if (k == i || k == j) {
          cnt++;
        } else {
          if ((p[k] - c).norm() < 1.0) {
            cnt++;
          }
        }
      }

      result = max(result, cnt);
    }
  }
  cout << result << "\n";
}

int main() {
  cin.tie(nullptr);
  ios::sync_with_stdio(false);
  int N;
  cin >> N;
  while (N != 0) {
    solve(N);
    cin >> N;
  }
  cout << flush;
  return 0;
}