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<algorithm>
#include<cmath>
using namespace std;

#define MAX_N 310

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;

				double X = (x[i] + x[j]) / 2.0;
				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);

				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.0001) {
						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.0001) {
						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>
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) {
                    ++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) {
                    ++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 <iostream>
#include <cstdio>
#include <vector>
#include <map>
#include <complex>
#include <cmath>
#include <algorithm>
using namespace std;

// * XY座標
#define X real()
#define Y imag()

// * 点の表現
typedef complex<double> P;

// * 許容する誤差ε
const double EPS = 1e-10;
// * Infinity
const double INF = 1e12;
// * 円周率
const double PI = acos(-1.0);
// * 
const P INF_P(INF,INF);

// * complex<double> の順序を定義する
namespace std {
	bool operator < (const P& a, const P& b) {
		return (real(a) != real(b))? real(a) < real(b) : imag(a) < imag(b);
	}
}

// * 点座標のデバッグ出力
void print(P p) {
	cout << "point" << p << ";" << endl;
}

// * 2乗する
double sq(double x) {
	return x * x;
}

// * 2点間の距離
double dist(P p1, P p2) {
	return abs(p1 - p2);
}

// * 2つのスカラーが等しいかどうか
bool equal(double a, double b) {
	return ( fabs( a - b ) < EPS );
}
// * 2つのベクトルが等しいかどうか
bool eq_v(P a, P b) {
	return ( equal( a.X, b.X ) && equal( a.Y, b.Y ) );
}

// * 内積 (dot product) : a・b = |a||b|cosΘ
double dot(P a, P b) {
	return real( conj(a) * b );
}

// * 外積 (cross product) : a×b = |a||b|sinΘ
double cross(P a, P b) {
	return imag( conj(a) * b ) ;
}

// ベクトル a の単位法線ベクトルを求める
vector<P> normal_unit(P a) {
	vector<P> vp;
	vp.push_back( (a * P(0,  1)) / abs(a) );
	vp.push_back( (a * P(0, -1)) / abs(a) ) ;
	return vp;
}

// 原点を軸に点 p を角度 a (ラジアン)だけ回転させた点を返す
P rot(P p, double a){
	double x = p.X * cos(a) - p.Y * sin(a);
	double y = p.X * sin(a) + p.Y * cos(a);
	return P(x,y);
}

int main(){
	int n;
	while( cin >> n , n ){
		vector<P> v;
		int ans=1;
		for(int i=0 ; i < n ; i++ ){
			double x, y;
			scanf("%lf %lf", &x, &y);
			v.push_back( P(x,y) );
		}
		for(int i=0 ; i < v.size() ; i++ ){
			for(int j=i+1 ; j < v.size() ; j++ ){
				if( i == j ) continue;
				P p1 = v[i];
				P p2 = v[j];
				double d = dist(p1,p2);
				if( d >= 2.0 ) continue;
				
				double angle = atan2( p2.Y-p1.Y , p2.X-p1.X );
				P mid = rot( p2-p1 , -angle ) / 2.0;
				double h = sqrt(1.0 - sq(d/2.0));
				P p3 = P( mid.X , h );
				P p4 = P( mid.X , -h );
				int cnt1=0, cnt2=0;
				for(int k=0 ; k < v.size() ; k++ ){
					P p = rot( v[k]-p1 , -angle );
					if( dist(p3,p) <= 1.0 ) cnt1++;
					if( dist(p4,p) <= 1.0 ) cnt2++;
				}
				ans = max( ans , max(cnt1,cnt2) );
			}
		}
		printf("%d\n", ans);
	}
}

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

#define loop(i,a,b) for(int i=a; i<int(b); i++)
#define rep(i,n) loop(i,0,n)

using namespace std;
typedef complex<double> vec;

double const pi = acos(-1);
double const eps = 1e-6;

int main(){
    int n;
    while (cin >> n && n){
        vector<vec> v(n);
        rep(i, n){
            double x, y; cin >> x >> y;
            vec p(x, y);
            v[i] = p;
        }

        int ans = 1;
        rep(i, n)rep(j, n)if (i != j){
            vec dir = v[i] - v[j];
            double d = abs(dir);
            if (d > 2) continue;
            ans = max(ans, 2);
            double h = sqrt(1 - d*d / 4);
            vec med = v[i] + v[j];
            med /= 2;
            vec c1 = polar(h, arg(dir) + pi / 2) + med;
            vec c2 = polar(h, arg(dir) - pi / 2) + med;
            int a1 = 0, a2 = 0;
            rep(k, n){
                if (abs(c1 - v[k]) <= 1) a1++;
                if (abs(c2 - v[k]) <= 1) a2++;
            }
            ans = max(ans, a1);
            ans = max(ans, a2);
        }
        cout << ans << endl;
    }
}