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 <cstdio>
#include <algorithm>
#include <cmath>
#include <queue>
#include <set>
#include <functional>
#include <ctime>

using namespace std;

#define fst first
#define snd second
#define all(c) ((c).begin()), ((c).end())
#define TEST(x,a) { auto y=(x); if (sign(y-a)) { cout << "line " << __LINE__  << #x << " = " << y << " != " << a << endl; exit(-1); } }
double urand() { return rand() / (1.0 + RAND_MAX); }

const double PI = acos(-1.0);

// implementation note: use EPS only this function
// usage note: check sign(x) < 0, sign(x) > 0, or sign(x) == 0
// notice: should be normalize to O(1)
const double EPS = 1e-8;
int sign(double x) {
  if (x < -EPS) return -1;
  if (x > +EPS) return +1;
  return 0;
}
struct point {
  typedef double T;
  T x, y; 
  point &operator+=(point p) { x += p.x; y += p.y; return *this; }
  point &operator-=(point p) { x -= p.x; y -= p.y; return *this; }
  point &operator*=(T a)     { x *= a;   y *= a;   return *this; }
  point &operator/=(T a)     { return *this *= (1.0/a); }
  point operator-() const    { return {-x, -y}; }
  bool operator==(point p) const { return !sign(x-p.x) && !sign(y-p.y); }
  bool operator!=(point p) const { return !operator==(p); }
  bool operator<(point p) const { return x!=p.x ? x<p.x : y<p.y; } // for sort
};
point operator+(point p, point q) { return p += q; }
point operator-(point p, point q) { return p -= q; }
point operator*(point::T a, point p) { return p *= a; }
point operator*(point p, point::T a) { return p *= a; }
point operator/(point p, point::T a) { return p /= a; }
point::T dot(point p, point q) { return p.x*q.x+p.y*q.y; }
point::T cross(point p, point q) { return p.x*q.y-p.y*q.x; } // left turn > 0
point::T norm2(point p) { return dot(p,p); }
point::T norm(point p) { return sqrt(dot(p,p)); }
point::T dist(point p, point q) { return norm(p-q); }
point orth(point p) { return {-p.y, p.x}; }

istream &operator>>(istream &is, point &p) { is>>p.x>>p.y;return is; }
ostream &operator<<(ostream &os, const point &p) { os<<"("<<p.x<<","<<p.y<<")"; return os; }

int maximum_circle_cover(vector<point> ps, double r) {
  struct range {
    point p; // center
    double w; // width
    int hi;
    bool operator<(range r) const {
      return hi < r.hi;
    }
  };
  double w = 0;
  for (point p: ps) w = max({w, abs(p.x), abs(p.y)});
  priority_queue<range> que;
  que.push({{0,0}, w, (int)ps.size()});
  point best_p;
  int best = 0;
  while (!que.empty()) {
    range R = que.top(); que.pop();
    if (R.hi <= best) continue;
    //cout << "processing " << R.p << " " << R.w << " " << R.hi << "/" << best << endl;
    double dx[] = {1,-1,-1,1}, dy[] = {1,1,-1,-1};
    for (int i = 0; i < 4; ++i) {
      range S = {R.p, R.w/2, 0};
      S.p += S.w*point({dx[i], dy[i]});
      int lo = 0;
      for (point q: ps) {
        auto d = dist(S.p, q);
        if (sign(d - r) <= 0) ++lo;
        if (sign(d - S.w*sqrt(2) - r) <= 0) ++S.hi;
      }
      if (lo > best) { best = lo; best_p = S.p; }
      best = max(lo, best);
      if (S.hi > best) que.push(S);
    }
  }
  return best; //best_p;
}
int maximum_circle_cover3(vector<point> ps, double r) {

  point best_p;
  int best = 0;
  function<void(point,double,vector<point>&)> rec = [&](point p, double w, vector<point> &ps) {
    w /= 2;
    const double dx[] = {1,-1,-1,1}, dy[] = {1,1,-1,-1};
    point qs[4];
    vector<point> pss[4];
    for (int i = 0; i < 4; ++i) {
      qs[i] = p + w * point({dx[i], dy[i]});
      int lo = 0;
      for (point q: ps) {
        auto d = dist(qs[i], q);
        if (sign(d - r) <= 0) ++lo;
        if (sign(d - w*sqrt(2) - r) <= 0) pss[i].push_back(q);
      }
      if (lo > best) { best = lo; best_p = qs[i]; }
    }
    int a = 0, b = 1, c = 2, d = 3;
    auto SW = [&](int &a, int &b) { if (pss[a].size() > pss[b].size()) swap(a, b); };
    SW(a,b);SW(c,d);SW(b,d);SW(a,c);SW(b,c);
    if (pss[d].size() > best) rec(qs[d], w, pss[d]);
    if (pss[b].size() > best) rec(qs[b], w, pss[b]);
    if (pss[c].size() > best) rec(qs[c], w, pss[c]);
    if (pss[a].size() > best) rec(qs[a], w, pss[a]);
  };
  double w = 0;
  for (point p: ps) w = max({w, abs(p.x), abs(p.y)});
  rec({0,0}, w, ps);
  return best; //best_p;
}

int maximum_circle_cover2(vector<point> ps, double r) {
  int best = 0;
  for (point p: ps) {
    int count = 0;
    vector<pair<double,int>> aux;
    for (point q: ps) {
      auto d = dist(p, q);
      if (sign(d) == 0) ++count;
      else if (sign(d - 2*r) <= 0) {
        double theta = atan2(q.y-p.y, q.x-p.x);
        double phi   = acos(min(1., d/(2*r)));
        aux.push_back({theta-phi, -1});
        aux.push_back({theta+phi, +1});
      }
    }
    sort(all(aux));
    /*
    cout << "for point " << p << endl;
    for (auto a: aux) 
      cout << "(" << a.fst << "," << a.snd << ") ";
    cout << endl;
    */

    for (auto a: aux) 
      best = max(best, count -= a.snd);
  }
  return best;
}


void verify_maximum_circle_cover2() {
  for (int n; scanf("%d", &n) && n; ) {
    vector<point> ps(n);
    for (int i = 0; i < n; ++i) 
      scanf("%lf %lf", &ps[i].x, &ps[i].y);
    printf("%d\n", maximum_circle_cover(ps, 1.0));
  }
}


int main() {
  //verify_maximum_circle_cover();
  verify_maximum_circle_cover2();
}