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;

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 (380 < 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 int long long
#define rep(i, a, b) for (int i = (a); i < (int)(b); ++i)
#define rrep(i, a, b) for (int i = (int)(b) - 1; i >= (int)(a); --i)
#define sz(c) ((int)c.size())
#define all(c) c.begin(), c.end()

using pii = pair<int, int>;
using vi = vector<int>;
using vvi = vector<vi>;

const int INF = 1e9;
const double EPS = 1e-6;
template <class T> bool chmin(T &a, const T &b) { if (a > b) { a = b; return true; } return false; }
template <class T> bool chmax(T &a, const T &b) { if (a < b) { a = b; return true; } return false; }
template <class T> using heap = priority_queue<T, vector<T>, greater<T>>;

using Point = complex<double>;

struct Circle {
    Point c;
    double r;
};

// 円と円の交点
pair<Point, Point> get_cross_points(Circle c1, Circle c2) {
    double d = abs(c1.c - c2.c);
    double a = acos((c1.r * c1.r - c2.r * c2.r + d * d) / (2 * c1.r * d));
    double t = arg(c2.c - c1.c);
    return {c1.c + polar(c1.r, t + a), c1.c + polar(c1.r, t - a)};
}

void solve(int N) {
    vector<Circle> C(N);
    rep(i, 0, N) {
        double x, y;
        cin >> x >> y;
        C[i] = {Point(x, y), 1};
    }

    auto check_intersect = [&](int i, int j) {
        return abs(C[i].c - C[j].c) < C[i].r + C[j].r;
    };

    auto check_inner = [&](Circle c, Point p) {
        return abs(c.c - p) <= c.r + EPS;
    };

    int ans = 0;
    rep(i, 0, N) {
        rep(j, i + 1, N) {
            if (check_intersect(i, j)) {
                auto cps = get_cross_points(C[i], C[j]);
                int cnt = 0;
                rep(k, 0, N) {
                    if (check_inner(C[k], cps.first)) {
                        cnt++;
                    }
                }
                chmax(ans, cnt);
                cnt = 0;
                rep(k, 0, N) {
                    if (check_inner(C[k], cps.second)) {
                        cnt++;
                    }
                }
                chmax(ans, cnt);
            }
        }
    }
    rep(i, 0, N) {
        int cnt = 0;
        rep(k, 0, N) {
            if (check_inner(C[k], C[i].c)) {
                cnt++;
            }
        }
        chmax(ans, cnt);
    }
    cout << ans << endl;
}

signed main() {
    while (1) {
        int N;
        cin >> N;
        if (N == 0) break;
        solve(N);
    }
}

C++: //???????????????????????????, ?????????????????????????????°????????????????????????????????????????????????2???????????¨???????????£?????\????????¨?????§????????????
//?????£???, 2??????????????????1?????????????????????????????????????????°??????????????????????£????, O(N^2)????????????????????????, ??¨?????§O(N^3)????????§?§£?????????
#include <iostream>
#include <complex>
#include <vector>
using namespace std;
typedef complex<double> Point;

int n;
double x[300], y[300];
vector<Point> points;	//????£???????

int main() {
	while (cin >> n) {
		if (n == 0) break;
		for (int i = 0; i < n; i++) cin >> x[i] >> y[i];
		if (n == 1) { cout << 1 << endl; continue; }
		
		for (int i = 0; i < n; i++) {
			for (int j = i + 1; j < n; j++) {
				Point p1 = Point(x[i], y[i]);
				Point p2 = Point(x[j], y[j]);
				if (norm(p1 - p2) < 4) {
					double h = sqrt(1.0 - norm(p1 - p2) / 4.0);
					Point t = p2 - p1; t /= abs(t);
					Point u = t * Point(0, 1);
					points.push_back(p1 + t * abs(p1 - p2) / 2.0 + u * h);
					points.push_back(p1 + t * abs(p1 - p2) / 2.0 - u * h);
				}
			}
		}
		
		double eps = 1e-10;
		int ans = 0;
		for (int i = 0; i < points.size(); i++) {
			double px = points[i].real();
			double py = points[i].imag();
			int cnt = 0;
			for (int j = 0; j < n; j++) {
				if ((x[j] - px) * (x[j] - px) + (y[j] - py) * (y[j] - py) < 1 + eps) {
					cnt++;
				}
			}
			if (ans < cnt) ans = cnt;
		}
		cout << ans << endl;
	}
	return 0;
}

C++: #include<iostream>
#include<cmath>

using namespace std;

int main(){

	int n;
	double x[301],y[301];
	int cnt,ans;

	while(1){

		cin>>n;
		if(n==0) break;

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

		ans=0;
		for(double i=0.5;i<=9.5;i+=0.01){
			for(double j=0.5;j<=9.5;j+=0.01){
				cnt=0;
				for(int k=0;k<n;k++){
					if(pow(i-x[k],2.0)+pow(j-y[k],2.0)<=1+1e-9) cnt++;
				}
				if(cnt>ans) ans=cnt;
				cnt=0;
			}
		}
		cout<<ans<<endl;
	}

	return 0;
}