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

// double
struct P{
    P(){}
    P(double _r, double _i):r(_r), i(_i){}
    void real(const double& v){r = v;}
    void imag(const double& v){i = v;}
    double r, i;
};

double real(const P& p){return p.r;}
double imag(const P& p){return p.i;}
double norm(const P& p){double r = real(p), i = imag(p); return r * r + i * i;}
double dist(const P& p){return std::sqrt(norm(p));}

P operator+(const P& lhs, const P& rhs){
    return P(real(lhs)+real(rhs), imag(lhs)+imag(rhs));
}
P operator-(const P& lhs, const P& rhs){
    return P(real(lhs)-real(rhs), imag(lhs)-imag(rhs));
}
P operator-(const P& p){
    return P(-real(p), -imag(p));
}

double cross(const P& lhs, const P& rhs){
    return real(lhs)*imag(rhs) - imag(lhs)*real(rhs);
}

double dot(const P& lhs, const P& rhs){
    return real(lhs)*real(rhs) + imag(lhs)*imag(rhs);
}

const double EPS = 1e-9, D_THETA = 0.06, D_R = 0.06;

int N;
P ps[300];
std::vector<P> is;

void addIntersection(int i, int j){
    if(dist(ps[i]-ps[j]) > 2.0){return;}
    
    P m{(real(ps[j])-real(ps[i])) / 2, (imag(ps[j])-imag(ps[i])) / 2}, n{-imag(m), real(m)};
    {
        double d = dist(n);
        n.real(real(n) / d);
        n.imag(imag(n) / d);
    }

    double t = std::sqrt(1 - norm(m));
    for(int k=-1;k<=1;k+=2){
        P p{t * k * real(n), t * k * imag(n)};
        P q = ps[i] + m + p;
        // printf("%.6f, %.6f\n", real(q), imag(q));
        is.push_back(q);
    }
}

int count(P p){
    int res = 0;
    for(int i=0;i<N;i++){
        if(dist(ps[i]-p) < 1.0 + EPS){res++;}
    }
    return res;
}

int main(){
    while(scanf("%d", &N), N){
        is.clear();
        
        for(int i=0;i<N;i++){
            double x, y;
            scanf("%lf %lf", &x, &y);

            ps[i].real(x);
            ps[i].imag(y);
        }

        for(int i=0;i<N;i++){
            for(int j=i+1;j<N;j++){
                addIntersection(i, j);
            }
        }
        
        int res = 1;
        for(auto i : is){
            res = std::max(res, count(i));
        }

        printf("%d\n", res);
    }
}

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

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

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 x = dist(a, b) / 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) {
                    if (dist(o, pnt[l]) < 1 + EPS) cnt++;
                }
                ans = max(ans, cnt);
            }
        }
        cout << ans << endl;
    }
    return 0;
}

C++: #include <bits/stdc++.h>
 
using namespace std;
 
#define EPS 1e-8
typedef complex<double> Point;
 
int main(){
    int N;
    while(cin >> N,N){
	Point p[N];
	for(int i = 0 ; i < N ; i++){
	    cin >> p[i].real() >> p[i].imag();
	}
	int res = 1;
	for(int i = 0 ; i < N-1 ; i++){
	    for(int j = i+1 ; j < N ; j++){
		Point np = p[j]-p[i];
		double d = abs(np);
		if(d + EPS > 2.0){ continue; }
		np /= d;
		np *= Point(0,1);
		Point cp = (p[i]+p[j])/2.0;
		Point v = sqrt(1.0-d*d/4.0)*np;
		int cnt = 2;
		Point center = cp + v;
		for(int k = 0 ; k < N ; k++){
		    if(i == k || j == k) continue;
		    Point diff = center-p[k];
		    if(norm(diff) < 1.0 + EPS){
			cnt++;
		    }
		}
		res = max(res,cnt);
	    }
	}
	cout << res << endl;
    }
    return 0;
}

C++: #include <string>
#include <vector>
#include <sstream>
#include <iostream>
#include <algorithm>
#include <map>
#include <list>
#include <set>
#include <numeric>
#include <queue>
#include <stack>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <cstring>
#include <climits>
#include <cfloat>
#include <ctime>
#include <complex>
#include <cassert>
#include <tuple>
using namespace std;
typedef long long LL;
typedef pair<int,int> P;

int N;
double x[300];
double y[300];
inline double dist2(double a, double b)
{
	return a*a+b*b;
}
void exec()
{
	cin >> N;
	if(N==0){
		exit(0);
	}
	for(int i=0;i<N;i++){
		cin >> x[i] >> y[i];
	}
	int ret=0;
	for(int i=0;i<N;i++){
		for(int j=i+1;j<N;j++){
			double dx=(x[i]-x[j])/2;
			double dy=(y[i]-y[j])/2;
			double r2=dx*dx+dy*dy;
			double r=sqrt(r2);
			double R=sqrt(1-r2);
			double px=x[j]+dx-R*dy/r;
			double py=y[j]+dy+R*dx/r;
			double qx=x[j]+dx+R*dy/r;
			double qy=y[j]+dy-R*dx/r;
			int sp=0,sq=0;
			for(int k=0;k<N;k++){
				if(k==i||k==j){
					sp++;sq++;
					continue;
				}
				if(dist2(x[k]-px,y[k]-py)<1){
					sp++;
				}
				if(dist2(x[k]-qx,y[k]-qy)<1){
					sq++;
				}
			}
			ret=max(ret,sp);
			ret=max(ret,sq);
		}
	}
	cout << ret << endl;
}
int main() {
	while(1)
		exec();
}

C++: #include <bits/stdc++.h>
using namespace std;
#define rep(i, n) for(int i = 0; i < (n); i++)
typedef complex<double> P;  // Point
static const double EPS = 1e-6;

int main() {
	int N;
	while (cin >> N, N) {
		vector<P> points;
		priority_queue<tuple<int, double, double, double>, vector<tuple<int, double, double, double>>> pq;
		rep(n, N) {
			double x, y;
			cin >> x >> y;
			points.push_back(P(x, y));
		}
		tuple<int, double, double, double> s = make_tuple(N, -5.0, 5.0, 5.0);
		pq.push(s);
		while(!pq.empty()) {
			s = pq.top();
			//cout << get<1>(s)<< " " <<get<0>(s) << endl;
			if (-get<1>(s) < EPS) break;
			pq.pop();
			double r = -get<1>(s) / 2;
			rep(i, 2) {
				rep(j, 2) {
					double x1 = get<2>(s) - get<1>(s) / 2 + i * get<1>(s);
					double y1 = get<3>(s) - get<1>(s) / 2 + j * get<1>(s);
					int cnt = 0;
					P p1 = P(x1, y1);
					for(auto p2 : points) {
						if (abs(p1 - p2) < r * sqrt(2) + 1 + EPS) cnt++;
					}
					pq.push(make_tuple(cnt, -r, x1, y1));
				}
			}
		}
		cout << get<0>(s) << endl;
	}
}