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 <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; }

vector<PD> solve(double ax, double ay, double bx, double by) {
  double sqdist = sq(ax - bx) + sq(ay - by);
  double dist = sqrt(sqdist);
  double ndist = sqrt(1 - sqdist / 4);
  double ex = (bx - ax) / dist, ey = (by - ay) / dist;
  vector<PD> ans(2);
  REP(i, 0, 2) {
    double coef = i == 0 ? 1 : -1;
    ans[i].first = (ax + bx) / 2 + ey * coef * ndist;
    ans[i].second = (ay + by) / 2 - ex * coef * ndist;
  }
  return ans;
}

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]);
    vector<PD> cand;
    REP(i, 0, n) {
      REP(j, 0, i) {
	for (PD e: solve(x[i], y[i], x[j], y[j])) {
	  cand.push_back(e);
	}
      }
    }
    int ma = 0;
    for (PD c: cand) {
      int cnt = 0;
      REP(i, 0, n) {
	double sqdist = sq(x[i] - c.first) + sq(y[i] - c.second);
	if (sqdist <= 1 + EPS) {
	  cnt += 1;
	}
      }
      ma = max(ma, cnt);
    }
    printf("%d\n", ma);
  }
}