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++: 
//Name: Circle and Points
//Level: 3
//Category: 幾何,円
//Note:

/*
 * 半径が固定なので，円周上に2点が乗っているときが境界状態．
 * すべての2点の組み合わせについて中心の来る場所（2点）を求め，内包判定を行う．
 * N=1のケースと，全ての点が互いに距離2以上離れているケースに注意．
 *
 * オーダーは O(N^3)．
 */
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <complex>
#include <utility>

using namespace std;

typedef complex<double> Point;

bool inCircle(const Point &p, const Point &center, double radius) {
    double dist2 = norm(p-center);
    double rad2 = radius*radius;
    return dist2 < rad2+1e-9;
}

int check(const Point &center, const vector<Point> &points) {
    const int N = points.size();
    int cnt = 0;
    for(int i = 0; i < N; ++i) {
        if(inCircle(points[i], center, 1)) ++cnt;
    }
    return cnt;
}

int test(int a, int b, const vector<Point> &points) {
    const Point v(points[a]-points[b]);
    if(abs(v) >= 2) return 1;
    const double hlen = abs(v)/2;

    const double arg = acos(hlen);
    const Point unit = v / abs(v);
    return max(check(points[b] + unit*polar(1.0, arg), points), check(points[b] + unit*polar(1.0, -arg), points));
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(0);

    while(true) {
        int N;
        cin >> N;
        if(!N) break;

        vector<Point> points(N);
        for(int i = 0; i < N; ++i) {
            cin >> points[i].real() >> points[i].imag();
        }
        int ans = 0;
        for(int i = 0; i < N; ++i) {
            for(int j = 0; j < i; ++j) {
                ans = max(ans, test(i, j, points));
            }
        }
        if(N == 1) ans = 1;
        cout << ans << endl;
    }
    return 0;
}

C++: #include <algorithm>
#include <cassert>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <deque>
#include <functional>
#include <list>
#include <map>
#include <queue>
#include <set>
#include <sstream>
#include <string>
#include <utility>
#include <vector>

#define REP(i,s,n) for(int i=(int)(s);i<(int)(n);i++)
#define DEBUGP(val) cerr << #val << "=" << val << "\n"

using namespace std;
typedef long long int ll;
typedef vector<int> VI;
typedef vector<ll> VL;
typedef pair<int, int> PI;
typedef pair<double, double> PD;

const double EPS = 1e-11;
double sq(double x) { return x * x; }

void solve(double ax, double ay, double bx, double by, int &ma,
		 const vector<double> &x, const vector<double> &y) {
  double sqdist = sq(ax - bx) + sq(ay - by);
  if (sqdist >= 4 + EPS) return;
  double dist = sqrt(sqdist);
  double ndist = sqrt(1 - sqdist / 4);
  double ex = (bx - ax) / dist, ey = (by - ay) / dist;
  int n = x.size();
  REP(i, 0, 2) {
    double coef = i == 0 ? 1 : -1;
    double cx = (ax + bx) / 2 + ey * coef * ndist;
    double cy = (ay + by) / 2 - ex * coef * ndist;
    int cnt = 0;
    REP(j, 0, n) {
      double sqdist = sq(x[j] - cx) + sq(y[j] - cy);
      if (sqdist <= 1 + EPS) {
	cnt += 1;
      }
    }
    ma = max(ma, cnt);
  }
}

int main(void) {
  int n;
  while(scanf("%d", &n) && n) {
    vector<double> x(n), y(n);
    REP(i, 0, n) scanf("%lf%lf", &x[i], &y[i]);
    int ma = 0;
    REP(i, 0, n) {
      REP(j, 0, i) {
        solve(x[i], y[i], x[j], y[j], ma, x, y);
      }
    }
    printf("%d\n", ma);
  }
}