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<string>
#include<cstdio>
#include<vector>
#include<cmath>
#include<algorithm>
#include<functional>
#include<iomanip>
#include<queue>
#include<ciso646>
#include<random>
#include<map>
#include<set>
#include<complex>
#include<bitset>
#include<stack>
#include<unordered_map>
#include<utility>
#include<tuple>
using namespace std;
typedef long long ll;
typedef unsigned int ui;
const ll mod = 1000000007;
const ll INF = (ll)1000000007 * 1000000007;
typedef pair<int, int> P;
#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 Per(i,sta,n) for(int i=n-1;i>=sta;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 double ld;
const ld eps = 1e-6;
const ld pi = acos(-1.0);
typedef pair<ll, ll> LP;



ld torad(int deg) {return (ld)(deg) * pi / 180;}
ld todeg(ld ang) {return ang * 180 / pi;}

struct Point {
    ld x, y;
    Point(ld x = 0.0, ld y = 0.0) : x(x), y(y) {}
    friend ostream& operator << (ostream &s, const Point &p) {return s << '(' << p.x << ", " << p.y << ')';}
};

Point operator + (const Point &p, const Point &q) {return Point(p.x + q.x, p.y + q.y);}
Point operator - (const Point &p, const Point &q) {return Point(p.x - q.x, p.y - q.y);}
Point operator * (const Point &p, ld a) {return Point(p.x * a, p.y * a);}
Point operator * (ld a, const Point &p) {return Point(a * p.x, a * p.y);}
Point operator * (const Point &p, const Point &q) {return Point(p.x * q.x - p.y * q.y, p.x * q.y + p.y * q.x);}
Point operator / (const Point &p, ld a) {return Point(p.x / a, p.y / a);}
Point conj(const Point &p) {return Point(p.x, -p.y);}
Point rot(const Point &p, ld ang) {return Point(cos(ang) * p.x - sin(ang) * p.y, sin(ang) * p.x + cos(ang) * p.y);}
Point rot90(const Point &p) {return Point(-p.y, p.x);}
ld cross(const Point &p, const Point &q) {return p.x * q.y - p.y * q.x;}
ld dot(const Point &p, const Point &q) {return p.x * q.x + p.y * q.y;}
ld norm(const Point &p) {return dot(p, p);}
ld abs(const Point &p) {return sqrt(dot(p, p));}
ld amp(const Point &p) {ld res = atan2(p.y, p.x); if (res < 0) res += pi*2; return res;}
bool eq(const Point &p, const Point &q) {return abs(p - q) < eps;}
bool operator < (const Point &p, const Point &q) {return (abs(p.x - q.x) > eps ? p.x < q.x : p.y < q.y);}
bool operator > (const Point &p, const Point &q) {return (abs(p.x - q.x) > eps ? p.x > q.x : p.y > q.y);}
Point operator / (const Point &p, const Point &q) {return p * conj(q) / norm(q);}

int ccw(const Point &a, const Point &b, const Point &c) {
    if (cross(b-a, c-a) > eps) return 1;
    if (cross(b-a, c-a) < -eps) return -1;
    if (dot(b-a, c-a) < 0) return 2;
    if (norm(b-a) < norm(c-a)) return -2;
    return 0;
}

struct Line : vector<Point> {
    Line(Point a = Point(0.0, 0.0), Point b = Point(0.0, 0.0)) {
        this->push_back(a);
        this->push_back(b);
    }
};


Point proj(Point p, Line l) {//射影
    ld t = dot(p - l[0], l[1] - l[0]) / norm(l[1] - l[0]);
    return l[0] + (l[1] - l[0]) * t;
}
Point refl(Point p, Line l) {//対象な点
    return p + (proj(p, l) - p) * 2;
}
bool isinterPL(Point p, Line l) {//点pが直線lに含まれるか
    return (abs(p - proj(p, l)) < eps);
}
bool isinterPS(Point p, Line s) {//点pが線分sに含まれるか
    return (ccw(s[0], s[1], p) == 0);
}
bool isinterLL(Line l, Line m) {//直線l,mが交わるか
    return (abs(cross(l[1] - l[0], m[1] - m[0])) > eps ||
            abs(cross(l[1] - l[0], m[0] - l[0])) < eps);
}
bool isinterSS(Line s, Line t) {//線分s,tが交わるか
    return (ccw(s[0], s[1], t[0]) * ccw(s[0], s[1], t[1]) <= 0 && 
            ccw(t[0], t[1], s[0]) * ccw(t[0], t[1], s[1]) <= 0);
}
ld distancePL(Point p, Line l) {//点pと直線lの距離
    return abs(p - proj(p, l));
}
ld distancePS(Point p, Line s) {//点pと線分sの距離
    Point h = proj(p, s);
    if (isinterPS(h, s)) { return abs(p - h); }
    return min(abs(p - s[0]), abs(p - s[1]));
}
ld distanceLL(Line l, Line m) {//直線l,mの距離
    if (isinterLL(l, m)) return 0;
    else return distancePL(m[0], l);
}
ld distanceSS(Line s, Line t) {//線分s,tの距離
    if (isinterSS(s, t)) return 0;
    else return min(min(distancePS(s[0], t), distancePS(s[1], t)), min(distancePS(t[0], s), distancePS(t[1], s)));
}
struct Circle : Point{
    ld r;
    Circle(Point p = Point(0.0, 0.0), ld r = 0.0) : Point(p), r(r) {}
};


vector<Point> crosspoint(Line l, Line m) {
    vector<Point> res;
    ld d = cross(m[1] - m[0], l[1] - l[0]);
    if (abs(d) < eps) return vector<Point>();
    res.push_back(l[0] + (l[1] - l[0]) * cross(m[1] - m[0], m[1] - l[0]) / d);
    return res;
}
vector<Point> crosspoint(Circle e, Circle f) {
    vector<Point> res;
    ld d = abs(e - f);
    if (d < eps) return vector<Point>();
    ld rcos = (d * d + e.r * e.r - f.r * f.r) / (2.0 * d);
    ld rsin = sqrt(e.r * e.r - rcos * rcos);
    Point dir = (f - e) / d;
    Point p1 = e + dir * Point(rcos, rsin);
    Point p2 = e + dir * Point(rcos, -rsin);
    res.push_back(p1);
    if (!eq(p1, p2)) res.push_back(p2);
    return res;
}
vector<Point> crosspoint(Circle e, Line l) {
    vector<Point> res;
    Point p = proj(e, l);
    ld rcos = abs(e - p), rsin;
    if (rcos > e.r + eps) return vector<Point>();
    else if (e.r - rcos < eps) rsin = 0;
    else rsin = sqrt(e.r * e.r - rcos * rcos);
    Point dir = (l[1] - l[0]) / abs(l[1] - l[0]);
    Point p1 = p + dir * rsin;
    Point p2 = p - dir * rsin;
    res.push_back(p1);
    if (!eq(p1, p2)) res.push_back(p2);
    return res;
}

void solve(){
  while(true){
    int n;cin >> n;
    if(n==0) return;
    vector<Circle> cs;
    rep(i,n){
      ld x,y;
      cin >> x >> y;
      cs.push_back(Circle(Point(x,y),1));
    }
    int ans=0;
    rep(i,n){
      int res=0;
      rep(k,n){
        if(abs(cs[i]-cs[k])<1+eps){
          res+=1;
        }
      }
      ans=max(ans,res);
      rep(j,n){
        vector<Point> v=crosspoint(cs[i],cs[j]);
        for(Point p:v){
          res=0;
          rep(k,n){
            if(abs(p-cs[k])<1+eps){
              res+=1;
            }
          }
          ans=max(ans,res);
        }
      }
    }
    cout << ans << endl;
  }
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout << fixed << setprecision(50);
    solve();
}