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 <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <complex>
#include <utility>

using namespace std;

typedef complex<double> Point;

bool eqv(double a, double b) {
    return abs(a-b) < 1e-10;
}

bool inCircle(const Point &p, const Point &center, double radius) {
    double dist2 = norm(p-center);
    double rad2 = radius*radius;
    return dist2 < radius*radius || eqv(dist2, rad2);
}

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]);
    const double hlen = abs(v)/2;

    if(hlen >= 1) return 0;
    const Point unit = v / abs(v);
    const double angle = acos(hlen);
    return max(check(points[b] + unit*polar(1.0, angle), points), check(points[b] + unit*polar(1.0, -angle), 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<iostream>
#include<algorithm>
#include<cmath>
using namespace std;

#define MAX_N 310

long double x[MAX_N], y[MAX_N];
int n, maxn;

int main() {
	while (true) {
		cin >> n;
		if (n == 0) { break; }
		for (int i = 0; i < n; i++) {
			cin >> x[i] >> y[i];
		}
		maxn = 0;
		for (int i = 0; i < n; i++) {
			for (int j = i + 1; j < n; j++) {
				int sum = 0;
				long double X = (x[i] + x[j]) / 2.0;
				long double Y = (y[i] + y[j]) / 2.0;
				for (int k = 0; k < n; k++) {
					if (((x[k] - X)*(x[k] - X) + (y[k] - Y)*(y[k] - Y)) <= 1) {
						sum++;
					}
				}
				maxn = max(maxn, sum);
				long double X1, Y1, X2, Y2, X3, Y3, VX2, VY2, Dist, Dist2;
				Dist = sqrt((x[j] - x[i])*(x[j] - x[i]) + (y[j] - y[i])*(y[j] - y[i]));
				Dist2 = sqrt(4 - Dist*Dist);
				X2 = x[j] - x[i]; Y2 = y[j] - y[i];
				VX2 = (X2 / Dist)*Dist2 / 2; VY2 = (Y2 / Dist)*Dist2 / 2;
				X1 = X + VY2; Y1 = Y - VX2;
				X3 = X - VY2; Y3 = Y + VX2;

				sum = 0;
				for (int k = 0; k < n; k++) {
					if (((x[k] - X1)*(x[k] - X1) + (y[k] - Y1)*(y[k] - Y1)) <= 1) {
						sum++;
					}
				}
				maxn = max(maxn, sum);
				sum = 0;
				for (int k = 0; k < n; k++) {
					if (((x[k] - X3)*(x[k] - X3) + (y[k] - Y3)*(y[k] - Y3)) <= 1) {
						sum++;
					}
				}
				maxn = max(maxn, sum);
			}
		}
		for (int i = 0; i < n; i++) {
			int sum = 0;
			long double X = x[i];
			long double Y = y[i];
			for (int k = 0; k < n; k++) {
				if (((x[k] - X)*(x[k] - X) + (y[k] - Y)*(y[k] - Y)) <= 1) {
					sum++;
				}
			}
			maxn = max(maxn, sum);
		}
		cout << maxn << endl;
	}
	return 0;
}

C++: #include <bits/stdc++.h>
 
#define loop(n, i) for(int i=0;i<n;i++)
#define all(v) v.begin(),v.end()

using namespace std;
typedef pair<double, double> P;
const double EPS = 1e-4;

#include <bits/stdc++.h>
 
#define loop(n, i) for(int i=0;i<n;i++)
#define all(v) v.begin(),v.end()

using namespace std;

vector<P> center(P a, P b)
{
    auto t = atan2(b.second - a.second, b.first - a.first);
    auto x = hypot(a.first - b.first, a.second - b.second) / 2;
    auto y = sqrt(1.0 - x * x);
    auto cs = cos(t), sn = sin(t);
    return { 
        { cs * x - sn * y + a.first, sn * x + cs * y + a.second },
        { cs * x + sn * y + a.first, sn * x - cs * y + a.second }
    };
}

int main()
{
    while (1) {
        int n; cin >> n;
        if (!n) break;
        vector<P> pnt;
        loop (n, i) {
            double x, y; cin >> x >> y;
            pnt.emplace_back(x, y);
        }
        int ans = 1;
        loop (n, i) loop (n, j) {
            auto ps = center(pnt[i], pnt[j]);
            for (auto o : ps) {
                int cnt = 0;
                loop (n, l) {
                    auto x = o.first - pnt[l].first;
                    auto y = o.second - pnt[l].second;
                    if (x * x + y * y < 1 + EPS) cnt++;
                }
                ans = max(ans, cnt);
            }
        }
        cout << ans << endl;
    }
    return 0;
}

C++: #include <iostream>
#include <cmath>
#include <vector>
#include <utility>
#include <algorithm>

using namespace std;

#define rep(i, n) for (int i = 0; i < (int)n; i++)
#define vec vector
#define pdd pair<double, double>
#define mp make_pair

double dis(pdd a, pdd b)
{
	pdd x = mp(a.first - b.first, a.second - b.second);
	return sqrt(x.first * x.first + x.second * x.second);
}

pdd find(pdd a, pdd b)
{
	pdd v = mp(a.first - b.first, a.second - b.second);
	double d = sqrt(v.first * v.first + v.second * v.second);
	// double d = dis(a, b);
	if (d > 2.0) return mp(1000, 1000);
	// pdd o = mp(a.first + 0.5 * v.first, a.second + 0.5 * v.second);
	pdd o = mp((a.first + b.first) / 2.0, (a.second + b.second) / 2.0);
	pdd g = mp(v.second, - v.first);
	double s = sqrt(1.0 - d * d / 4.0);
	// double g_size = sqrt(g.first * g.first + g.second * g.second);
	double g_size = d;
	pdd ret = mp(o.first + g.first / g_size * s, o.second + g.second / g_size * s);
	return ret;
}

int main(void)
{
	while(1) {
		int n;
		cin >> n;
		if (n == 0) break;
		vec<pdd> a(n);
		rep(i, n) cin >> a[i].first >> a[i].second;
		// if (n == 1) {
		// 	cout << 1 << endl;
		// 	continue;
		// }

		int MAX = 0;
		rep(i, n) rep(j, n) {
			if (i == j) continue;
			pdd c = find(a[i], a[j]);
			// cout << i << " " << j << " " << c.first << " " << c.second << endl;

			int p = 0;
			if (c.first == 1000) p = 1;
			else {
					rep(k, n) {
					// cout << dis(a[k], c) << " ";
					if (dis(a[k], c) <= 1.0001) p++;
				}
			}
			// cout << p << endl;
			MAX = max(p, MAX);
		}

		cout << MAX << endl;
	}
	return 0;
}