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 eps = 1e-6;

double  caclDistance2(double x1, double y1, double x2, double y2){
    double ans;
    ans = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2);

    return ans;
}

int main(){
    while(1){
        int n;
        cin >> n;
        if(n == 0){
            break;
        }

        vector<pair<double, double> > point(n);

        int ans = 1;

        for(int i=0; i<n; i++){
            cin >> point[i].first >> point[i].second;
        }

        for(int i=0; i<n; i++){
            for(int j=i+1; j < n; j++){

                int point_num = 0;

                double distance_of_two_points = sqrt(caclDistance2(point[i].first, point[i].second, point[j].first, point[j].second));
                //cout << "distance = "<< distance_of_two_points << endl;
                if(distance_of_two_points > 2.0 + eps) continue;

                pair<double, double> midpoint;
                midpoint.first = (point[i].first + point[j].first) / 2;
                midpoint.second = (point[i].second + point[j].second) / 2;

                pair<double, double> vector;                
                pair<double, double> centor;
                
                vector.first = (-1) * (point[i].second - point[j].second);
                vector.second = (point[i].first - point[j].first);

                //cout << vector.first << " " << vector.second << endl;

                double vector_size = sqrt(vector.first * vector.first + vector.second * vector.second);

                double distance_from_centor = sqrt(1.0 - (distance_of_two_points / 2.0) * (distance_of_two_points / 2.0));

                centor.first = midpoint.first + (distance_from_centor/vector_size) * vector.first;
                centor.second = midpoint.second + (distance_from_centor/vector_size) * vector.second;

                //cout << centor.first << " " << centor.second << endl;

                for(int k = 0; k<n; k++){
                    if(caclDistance2(centor.first, centor.second, point[k].first, point[k].second) <= 1.00005){
                        point_num++;
                    }
                }

                ans = max(ans, point_num);

                point_num = 0;

                vector.first = (point[i].second - point[j].second);
                vector.second = (-1) * (point[i].first - point[j].first);

                centor.first = midpoint.first + (distance_from_centor/vector_size) * vector.first;
                centor.second = midpoint.second + (distance_from_centor/vector_size) * vector.second;

                //cout << centor.first << " " << centor.second << endl;

                for(int k = 0; k<n; k++){
                    if(caclDistance2(centor.first, centor.second, point[k].first, point[k].second) <= 1.00005){
                        point_num++;
                    }
                }

                ans = max(ans, point_num);                
            }
        }

        cout << ans << 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-6;

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

template <typename T>
class Mat
{
public:
    vector<vector<T>> mat;
    int rows, columns;

    Mat<T>(int r, int c) : rows(r), columns(c) {
        mat = vector<vector<T>>(rows, vector<T>(columns));
    }

    Mat<T>(initializer_list<initializer_list<T>> _mat)
        : rows(_mat.size()), columns(_mat.begin()->size()) {
        mat = vector<vector<T>>(rows, vector<T>(columns));
        int i = 0;
        for (auto&& row : _mat) {
            copy(row.begin(), row.end(), mat[i++].begin());
        }
    }

    vector<T> &operator [](int m) {
        return mat[m];
    }

    vector<T> operator *(const vector<T> v) const {
        vector<T> ret(rows);
        for (int m = 0; m < rows; m++) {
            for (int n = 0; n < columns; n++) {
                ret[m] += mat[m][n] * v[n];
            }
        }
        return ret;
    }

    Mat<T> operator *(Mat<T>& o) const {
        Mat<T> ret(rows, o.columns);

        for (int r = 0; r < rows; r++) {
            for (int c = 0; c < o.columns; c++) {
                for (int i = 0; i < columns; i++) {
                    ret[r][c] += mat[r][i] * o[i][c];
                }
            }
        }
        return ret;
    }

    void operator *=(Mat<T>& o) {
        *this = *this * o;
    }

    Mat<T> pow(int n) const {
        Mat ans(rows, columns);
        for (int i = 0; i < rows; i++) {
            ans[i][i] = 1;
        }

        Mat base = *this;

        for (int i = 0; n>>i; i++) {
            if (n>>i&1) {
                ans *= base;
            }
            base *= base;
        }
        return ans;
    }

    Mat<T> &t() {
        vector<vector<T>> _mat(columns, vector<T>(rows));
        for (int m = 0; m < rows; m++) for (int n = 0; n < columns; n++) {
            _mat[n][m] = mat[m][n];
        }
        mat = _mat;
        swap(rows, columns);
        return *this;
    }
};

double dist(P a, P b)
{
    return hypot(a.first - b.first, a.second - b.second);
}

vector<P> center(P a, P b)
{
    auto t = atan2(b.first - a.first, b.second - a.second);
    auto d = dist(a, b) / 2;
    vector<double> x1 = { d, sqrt(1.0 - d * d) };
    vector<double> x2 = { d, -sqrt(1.0 - d * d) };
    Mat<double> r = { { cos(t), -sin(t) }, { sin(t), cos(t) } };
    vector<P> ans;
    for (vector<double> v : { r * x1, r * x2 }) {
        ans.emplace_back(v[0] + a.first, v[1] + a.second);
    }
    return ans;
}

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) {
            if (dist(pnt[i], pnt[j]) < 2) ans = 2;
        }
        loop (n, i) loop (n, j) {
            auto ps = center(pnt[i], pnt[j]);
            for (auto o : ps) {
                int cnt = 0;
                loop (n, l) {
                    if (dist(o, pnt[l]) < 1 + EPS) cnt++;
                }
                ans = max(ans, cnt);
            }
        }
        cout << ans << endl;
    }
    return 0;
}