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 <bits/stdc++.h>
using namespace std;

typedef complex<double> C;

int solve(int N, vector<double>& x, vector<double>& y) {
    vector<C> point(N);

    for (int j = 0; j < N; ++j) {
        point[j] = C(x[j], y[j]);
    }

    int ret = 1;

    for (int j = 0; j < N; ++j) {
        for (int k = j+1; k < N; ++k) {
            C mid = (point[j] + point[k]) * 0.5;

            C normal = (point[j] - point[k]) * C(0, 1);
            normal /= abs(point[j] - point[k]);

            double a = sqrt(1 - norm(mid - point[j]));

            C centre1 = mid + a * normal;
            int count = 0;
            for (int l = 0; l < N; ++l) {
                if (abs(centre1 - point[l]) < 1.0 + 1e-9) {
                    ++count;
                }
            }
            ret = max(ret, count);

            C centre2 = mid - a * normal;
            count = 0;
            for (int l = 0; l < N; ++l) {
                if (abs(centre2 - point[l]) < 1.0 + 1e-9) {
                    ++count;
                }
            }
            ret = max(ret, count);

        }
    }

    return ret;
}

int main () {
    int N;
    vector<double> x, y;

    while (true) {
        cin >> N;
        if (N == 0) break;

        x.resize(N);
        y.resize(N);

        for (int j = 0; j < N; ++j) {
            cin >> x[j] >> y[j];
        }

        cout << solve(N, x, y) << endl;
    }

    return 0;
}

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

using namespace std;

#define X first
#define Y second

using R=long double;
using P=pair<R, R>;

const double EPS = 1e-10;

R dist2(P a, P b) {
	R dx = a.X - b.X;
	R dy = a.Y - b.Y;
	return dx*dx+dy*dy;
}

int main() {
	int n;
	while(cin >> n, n) {
		vector<P> p(n), pl;
		for(int i = 0; i < n; i++) cin >> p[i].X >> p[i].Y;
		for(int i = 0; i < p.size(); i++) {
			P p1 = p[i];
			for(int j = i+1; j < p.size(); j++) {
				P p2 = p[j];
				if(dist2(p1, p2)-EPS > 4.0) continue;

				R dxy2 = sqrt(dist2(p1, p2)/4.0);
				R l = sqrt(1-(dxy2*dxy2));
				R cx = (p1.X+p2.X)/2.0;
				R cy = (p1.Y+p2.Y)/2.0;
				R dx = (p2.X-p1.X);
				R dy = (p2.Y-p1.Y);

				R dlx = l/sqrt(dist2(p1, p2)) * dy;
				R dly = l/sqrt(dist2(p1, p2)) * -dx;

				pl.push_back(P(cx+dlx, cy+dly));
				pl.push_back(P(cx-dlx, cy-dly));
			}
		}

		int res = 1;
		for(P p1:pl) {
			int tmp = 2;
			for(P p2:p) {
				if(sqrt(dist2(p1, p2)) < 1.0-EPS) tmp++;
			}
			res = max(res, tmp);
		}

		cout << res << endl;
	}

	return 0;
}

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

using namespace std;

const double r = 1.0;
const double PI = acos(-1.0);

struct Point {
    double x, y;
    Point() {}
    Point(double x, double y) : x(x), y(y) {}
};
inline double Distance2(const Point &a, const Point &b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

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

    int n;
    std::random_device rnd;
    std::mt19937 mt(rnd());
    std::uniform_real_distribution<> rand100(0.0, 10.0);

    while (cin >> n, n) {
        // Input
        vector<Point> p(n);
        for (int i = 0; i < n; ++i)
            cin >> p[i].x >> p[i].y;

        // Solve
        int ans = 1;
        auto start = std::chrono::high_resolution_clock::now();
        auto end = std::chrono::high_resolution_clock::now();
        auto take_time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

        while (true) {
            end = std::chrono::high_resolution_clock::now();
            take_time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
            if (100 < take_time.count())
                break;

            for (int i = 0; i < 10; ++i) {
                Point c(rand100(mt), rand100(mt));
                int num = 0;
                for (int j = 0; j < n; ++j) {
                    double d = Distance2(p[j], c);
                    if (d <= 1.0)
                        ++num;
                }
                ans = max(ans, num);
            }
        }

        // Output
        cout << ans << endl;
    }

    return 0;
}

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

using namespace std;

#define rep(i,n) for(int i=0;i<(int)n;i++)
#define cd complex<double>

complex<double> ci = complex<double>(0.0,1.0);
complex<double> mci = complex<double>(0.0,-1.0);

int main(){

  while(1){
    int n;
    cin >> n;
    if (n==0)break;
    vector<complex<double> > p(n);
    rep(i,n){
      double real,imag;
      cin >> real >> imag;
      p[i] = complex<double>(real,imag);
    }
    int summax = 1;
    rep(i,n){
      rep(j,n){
        if (i==j)continue;
        if (abs(p[i]-p[j])>2)continue;
        complex<double> mid = (p[i]+p[j])/2.0;
        double l = abs(p[i]-mid);
        double h = sqrt(1-l*l);
        complex<double> c1 = mid + (mid-p[i])*ci*h/l;
        complex<double> c2 = mid + (mid-p[i])*mci*h/l;

        int sum1 = 2,sum2 = 2;
        rep(k,n){
          if(k==i||k==j)continue;
          if(abs(p[k]-c1)<1.0)sum1++;
          if(abs(p[k]-c2)<1.0)sum2++;
        }
        if (max(sum1,sum2)>summax)summax=max(sum1,sum2);
      }
    }
    cout << summax << endl;
  }

  return 0;
}