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++: //#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[305];
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 = 1;
    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;
}

C++: #include <bits/stdc++.h>
using namespace std;

#define rep(i, s, n) for (Int i = s; i < (Int)(n); i++)
#define dump(x) cout << (x) << '\n'
#define Int int64_t
#define fi first
#define se second
#define all(x) (x).begin(), (x).end()
 
double EPS = 1e-10;
Int INF = 1e18;
int inf = 1e9;
Int mod = 1e9+7;

int main() {
    vector<Int> res;
    Int n;
    while (cin >> n and n) {
        vector<double> x(n), y(n);
        rep(i, 0, n) cin >> x[i] >> y[i];
        Int ma = 1;
        rep(i, 0, n - 1) rep(j, i + 1, n) {
            double d = sqrt(pow(x[i] - x[j], 2.0) + pow(y[i] - y[j], 2.0));
            if (d > 2.0) continue;
            complex<double> m((x[i] + x[j]) / 2.0, (y[i] + y[j]) / 2.0);
            complex<double> mv(1.0, -1.0 / ((y[i] - y[j]) / (x[i] - x[j])));
            mv /= abs(mv);
            complex<double> ml(abs(x[i] - x[j]) / 2.0, abs(y[i] - y[j]) / 2.0);
            mv *= sqrt(1.0 - pow(d / 2.0, 2.0));
            complex<double> a = m + mv, b = m - mv;
            // printf("a : %.12lf %.12lf\n", a.real(), a.imag());
            // printf("b : %.12lf %.12lf\n", b.real(), b.imag());
            Int cnt_a = 0;
            rep(k, 0, n) {
                double xx = pow(x[k] - a.real(), 2.0);
                double yy = pow(y[k] - a.imag(), 2.0);
                // printf("xx + yy = %.12lf\n", xx + yy);
                if (xx + yy <= 1.0 + EPS) cnt_a++;
            }
            Int cnt_b = 0;
            rep(k, 0, n) {
                double xx = pow(x[k] - b.real(), 2.0);
                double yy = pow(y[k] - b.imag(), 2.0);
                // printf("xx + yy = %.12lf\n", xx + yy);
                if (xx + yy <= 1.0 + EPS) cnt_b++;
            }
            ma = max(ma, max(cnt_a, cnt_b));
        }
        res.push_back(ma);
    }
    rep(i, 0, res.size()) dump(res[i]);
    return 0;
}