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 <iostream>
#include <cstdio>
#include <complex>
#include <vector>
#define X real()
#define Y imag()
#define EPS (1e-10)
using namespace std;
typedef complex<double> P;
typedef pair<P, double> C;

namespace std{
  bool operator<(const P& a, const P& b){
    return a.X != b.X ? a.X < b.X : a.Y < b.Y;
  }
}

P n_vector(P a){
  double l = abs(a);
  return P(-a.Y,a.X)/l;
}

vector<P> cpcc(C a, C b){
  double d = abs(a.first-b.first), ra = a.second, rb = b.second;
  double t = (ra*ra-rb*rb+d*d)/2/d, h = sqrt(ra*ra-t*t);
  P m = t/abs(b.first-a.first)*(b.first-a.first)+a.first;
  vector<P> cp(2,m);
  P n = n_vector(a.first-b.first);
  cp[0] -= h*n;
  cp[1] += h*n;
  return cp;
}


int main(){
  int N;
  double x, y;
  while(cin >> N, N){
    P ps[N];
    for(int i = 0; i < N; ++i){
      cin >> x >> y;
      ps[i] = P(x,y);
    }
    if(N == 1){
      cout << 1 << endl;
      continue;
    }
    vector<P> cps;
    for(int i = 0; i < N; ++i){
      for(int j = i+1; j < N; ++j){
	if(abs(ps[i]-ps[j]) <= 2 + EPS){
	  vector<P> cp = cpcc(C(ps[i],1),C(ps[j],1));
	  for(int k = 0; k < cp.size(); ++k)
	    cps.push_back(cp[k]);
	}
      }
    }
    int n = cps.size(), ans = 0;;
    for(int i = 0; i < n; ++i){
      int cnt = 0;
      for(int j = 0; j < N; ++j){
	if(norm(cps[i]-ps[j]) <= 1 + EPS)
	  ++cnt;
      }
      ans = max(ans, cnt);
    }
    cout << ans << endl;
  }
  return 0;
}

C++: #include <iostream>
#include <deque>
#include <string>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <string>
#include <cmath>
#include <complex>

#define loop(i,a,b) for(int i=a; i<int(b); i++)
#define rep(i,n) loop(i,0,n)

using namespace std;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef complex<double> vec;

double const pi = acos(-1);
double const eps = 1e-6;

int main(){
    int n;
    while (cin >> n && n){
        vector<vec> v(n);
        rep(i, n){
            double x, y; cin >> x >> y;
            vec p(x, y);
            v[i] = p;
        }

        if (n == 1) {
            puts("1");
            continue;
        }
        int ans = 1;
        rep(i, n)rep(j, i){
            double d = abs(v[i] - v[j]);
            if (d > 2) continue;
            ans = max(ans, 2);
            double h = 1 - d*d / 4;
            vec c1 = polar(h, arg(v[i] - v[j]) + pi / 2) + (v[i] + v[j]) / 2.;
            vec c2 = polar(h, arg(v[i] - v[j]) - pi / 2) + (v[i] + v[j]) / 2.;
            int a1 = 0, a2 = 0;
            rep(k, n){
                if (abs(c1 - v[k]) <= 1) a1++;
                if (abs(c2 - v[k]) <= 1) a2++;
            }
            ans = max(ans, a1);
            ans = max(ans, a2);
        }
        cout << ans << endl;
    }
}

C++: //#define __USE_MINGW_ANSI_STDIO 0
#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
#define int ll
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<ll> VL;
typedef vector<VL> VVL;
typedef pair<int, int> PII;

#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)
#define REP(i, n) FOR(i, 0, n)
#define ALL(x) x.begin(), x.end()
#define IN(a, b, x) (a<=x&&x<b)
#define MP make_pair
#define PB push_back
#ifdef int
const ll INF = (1LL<<60);
#else
const int INF = (1LL<<30);
#endif
const double PI = 3.14159265359;
const double EPS = 1e-12;
const int MOD = 1000000007;

template <typename T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template <typename T> T &chmax(T &a, const T &b) { return a = max(a, b); }

int dx[] = {0, 1, 0, -1}, dy[] = {1, 0, -1, 0};

// y???real(), x???imag()
//point
typedef complex<double> P;
namespace std {
	bool operator < (const P& a, const P& b) {
		return real(a) != real(b) ? real(a) < real(b) : imag(a) < imag(b);
	}
	bool cmp_y(const P &a, const P &b){
		return a.imag() != b.imag() ? a.imag() < b.imag() : a.real() < b.real();
	}
}
// circle
struct C {
	P p; double r;
	C(const P& p, double r) : p(p), r(r) {}
};

// 2???p1, p2?????????????????????r?????????2?????????
vector<C> calcCircle(P p1, P p2, double r) {
	if(abs(p1-p2) > 2*r) return {};
	P p3 = {(p1.real()+p2.real())/2, (p1.imag()+p2.imag())/2};
	// cout << p3.real() << " " << p3.imag() << endl;
	double l = abs(p1-p3);
	// cout << "l:" << l << endl;
	// ????????????p_1p_2
	P p1p2 = p2-p1;
	double a = p1p2.real(), b = p1p2.imag();
	double dx = b*sqrt((r*r-l*l)/(a*a+b*b)), dy = a*sqrt((r*r-l*l)/(a*a+b*b));
	// cout << "dx:" << dx << " dy:" << dy << endl;
	return {{{p3.real()+dx, p3.imag()-dy}, r}, {{p3.real()-dx, p3.imag()+dy}, r}};
}

// ???p??????c????????¨????????¨??????????????????
bool intersectCP(C c, P p) { return abs(p-c.p) <= c.r + EPS; }

P po[105];
signed main(void)
{
  while(true) {
    int n;
    cin >> n;
    if(!n) break;
    REP(i, n) {
      double x, y;
      cin >> x >> y;
      po[i] = {x, y};
    }

    int ans = 0;
    REP(i, n) FOR(j, i+1, n) {
      vector<C> ret = calcCircle(po[i], po[j], 1);
      if(ret.size() == 0) continue;
      // cout << "i:" << i << " j:" << j << endl;
      // cout << ret[0].p.real() << " " << ret[0].p.imag() << " " << ret[0].r << endl;
      // cout << ret[1].p.real() << " " << ret[1].p.imag() << " " << ret[1].r << endl;
      int num = 0;
      REP(k, n) {
        num += intersectCP(ret[0], po[k])?1:0;
        // if(intersectCP(ret[0], po[k])) cout << k << " ";
      }
      chmax(ans, num);
      // cout << "num:" << num << endl;
      num = 0;
      REP(k, n) {
        num += intersectCP(ret[1], po[k])?1:0;
        // if(intersectCP(ret[1], po[k])) cout << k << " ";
      }
      // cout << "num:" << num << endl;
      chmax(ans, num);
    }

    cout << ans << endl;
  }

  return 0;
}