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 <cmath>
#include <vector>

const double EPS = 1e-4; // epsilon
const double R = 1;      // radius of circle

/* Distance from (x1,y1) to (x2,y2) */
double dist(std::pair<double, double> a, std::pair<double, double> b) {
    
    double ax, ay, bx, by;
    ax = a.first;
    ay = a.second;
    bx = b.first;
    by = b.second;
    return std::sqrt(
                        (by - ay) * (by - ay) +
                        (bx - ax) * (bx - ax)
                    );
}

/* Midpoint of two points */
std::pair<double, double>
midpoint(std::pair<double, double> a, std::pair<double, double> b) {
    std::pair<double, double> mid;
    mid.first = (a.first + b.first) / 2;
    mid.second = (a.second + b.second) / 2;
    return mid;
}

/*
Find (By)
((By) - Ay)^2 + (Bx - Ax)^2 = 1
((By) - Ay)^2 = 1 - (Bx - Ax)^2
(By) - Ay = sqrt(1 - (Bx - Ax)^2)
(By) = sqrt(1 - (Bx - Ax)^2) + Ay
*/
int main() {
    int n;
    while (std::cin >> n) {
        if (n <= 0) break;
        int maxPts = 1;
        if (n > 1) {
            std::vector<std::pair<double, double>> pts;
        
            // insert points into vector
            for(int i = 0; i < n; ++i) {
                double x, y;
                std::cin >> x >> y;
                pts.emplace_back(x, y);
            }
            
            for (unsigned int i = 0; i < pts.size() - 1; ++i) {
                for (unsigned int j = i + 1; j < pts.size(); ++j) {
                    // std::pair<double, double> a, b;
                    auto a = pts.at(i);
                    auto b = pts.at(j);
                    if( dist(a, b) < (2*R-EPS)) {
                        int numPts = 0;
                        auto c = midpoint(a, b);
                        auto cy1 = std::sqrt(1 - (c.first-a.first)*(c.first-a.first)) +
                                   a.second;
                        auto cy2 = -std::sqrt(1 - (c.first-a.first)*(c.first-a.first)) +
                                   a.second;
                        c.second = cy1;
                        for (unsigned int k = 0; k < pts.size(); ++k) {
                            // std::pair<double, double> p;
                            auto p = pts.at(k);
                            if (dist(c, p) < (R + EPS)) {
                                numPts++;
                            }
                        }
                        maxPts = std::max(numPts, maxPts);
                        
                        numPts = 0;
                        c.second = cy2;
                        for (unsigned int k = 0; k < pts.size(); ++k) {
                            // std::pair<double, double> p;
                            auto p = pts.at(k);
                            if (dist(c, p) < (R + EPS)) {
                                numPts++;
                            }
                        }
                        maxPts = std::max(numPts, maxPts);
                    }
                }
            }
        
        }
        std::cout << maxPts << "\n";
    }
    return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;

double Dis(double x1, double y1, double x2, double y2);
int main(void){
	while(1){
		int n;
		int ans = 1;
		vector< pair<double, double> > p;

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

		for(int i=0; i<n; i++){
			pair<double, double> t;
			cin >> t.first >> t.second;
			p.push_back(t);
		}

		for(int i=0; i<n; i++){
			for(int j=i+1; j<n; j++){
				//中点
				double mpx = (p[i].first + p[j].first) / 2.0;
				double mpy = (p[i].second + p[j].second) / 2.0;
				//2点間の距離
				double d1 = Dis(p[i].first, p[i].second, p[j].first, p[j].second);
				//p[i]からp[j]へ向かうベクトル
				double v1x = (p[j].first - p[i].first) / d1;
				double v1y = (p[j].second - p[i].second) / d1;
				//中点から円の中心へ向かう単位ベクトル
				double v2x = -v1y;
				double v2y = v1x;
				//中心から円の中心へ向かうベクトルの長さ
				double d2 = sqrt(1.0 - (d1*d1/4.0));

				//円の中心の座標
				double cx = mpx + v2x * d2;
				double cy = mpy + v2y * d2;

				int cnt = 0;
				for(int k=0; k<n; k++){
					if(k == i || k == j || (cx-p[k].first)*(cx-p[k].first) + (cy-p[k].second)*(cy-p[k].second) <= 1.0)	cnt++;
				}
				if(cnt > ans){
					ans = cnt;
				}

				cx = mpx - v2x * d2;
				cy = mpy - v2y * d2;

				cnt = 0;
				for(int k=0; k<n; k++){
					if(k == i || k == j || (cx-p[k].first)*(cx-p[k].first) + (cy-p[k].second)*(cy-p[k].second) <= 1.0)	cnt++;
				}
				if(cnt > ans){
					ans = cnt;
				}
			}
		}

		cout << ans << endl;


	}
	
	return 0;
}

double Dis(double x1, double y1, double x2, double y2){
	return sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2));
}

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

using namespace std;

#define int long long
#define MKP make_pair
#define fi first
#define se second
typedef pair<int, int> pii;

constexpr int MOD = 1e9 + 7;
constexpr int MAX_N = 205;
constexpr int MAX_W = 10005;
constexpr int INF = 1LL << 60;

/* 基本要素 */

typedef double D;      // 座標値の型。doubleかlong doubleを想定
typedef complex<D> P;  // Point
typedef pair<P, P> L;  // Line
typedef vector<P> VP;
const D EPS = 1e-9;    // 許容誤差。問題によって変える
#define X real()
#define Y imag()
#define LE(n,m) ((n) < (m) + EPS)
#define GE(n,m) ((n) + EPS > (m))
#define EQ(n,m) (abs((n)-(m)) < EPS)

// 内積　dot(a,b) = |a||b|cosθ
D dot(P a, P b) {
  return (conj(a)*b).X;
}
// 外積　cross(a,b) = |a||b|sinθ
D cross(P a, P b) {
  return (conj(a)*b).Y;
}

VP crosspointCC(P a, D ar, P b, D br) {
  VP ps;
  P ab = b-a;
  D d = abs(ab);
  D crL = (norm(ab) + ar*ar - br*br) / (2*d);
  if (EQ(d, 0) || ar < abs(crL)) return ps;

  P abN = ab * P(0, sqrt(ar*ar - crL*crL) / d);
  P cp = a + crL/d * ab;
  ps.push_back(cp + abN);
  if (!EQ(norm(abN), 0)) ps.push_back(cp - abN);
  return ps;
}

int n;
vector<P> v;

bool solve() {
    v.clear();
    cin >> n;
    if (!n) return false;

    for (int i = 0; i < n; ++i) {
        double x, y;
        cin >> x >> y;
        v.emplace_back(P(x, y));
    }

    int res = 0;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            auto crs = crosspointCC(v[i], 1.0, v[j], 1.0);
            for (int k = 0; k < crs.size(); ++k) {
                int cnt = 0;
                for (int l = 0; l < n; ++l) {
                    double dist = abs(crs[k] - v[l]);
                    if (dist <= 1.0 + EPS) cnt++;
                }

                res = max(res, cnt);
            }
        }
    }

    cout << res << endl;

    return true;
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    while (solve());
}