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;
typedef long long ll;
typedef pair<int,int> Pii;
typedef pair<ll,ll> Pll;
#define rep(i,n) for (ll i=0;i<n;++i)
#define rep2(i,a,b) for (ll i=a;i<b;++i)
const ll MOD=1000000007;

template<class T> inline bool chmin(T &a,T b){
    if (a>b){
        a=b;
        return true;
    }
    return false;
}
template<class T> inline bool chmax(T &a,T b){
    if (a<b){
        a=b;
        return true;
    }
    return false;
}

typedef pair<double,double> Pdd;

int main(){
    int N;
    while(cin >> N,N>0){
        vector<Pdd> pos(N);
        rep(i,N){
            double x,y;
            cin >> x >> y;
            pos[i]=Pdd(x,y);
        }
        int ans=1;
        rep(i,N){
            rep2(j,i+1,N){
                double x1=pos[i].first,y1=pos[i].second;
                double x2=pos[j].first,y2=pos[j].second;
                double dist=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
                if (dist>2.0){
                    continue;
                }
                else {
                    int count=2;
                    double halfd=dist/2.0;
                    double len=sqrt(1.0-halfd*halfd);
                    double cx=(x1+x2)/2.0,cy=(y1+y2)/2.0;
                    double dx=cx+(y2-y1)*len/dist;
                    double dy=cy-(x2-x1)*len/dist;
                    rep(k,N){
                        if (k==i||k==j){
                            continue;
                        }
                        double x=pos[k].first;
                        double y=pos[k].second;
                        if ((x-dx)*(x-dx)+(y-dy)*(y-dy)<=1.0){
                            count++;
                        }
                    }
                    chmax(ans,count);
                    count=2;
                    dx=cx-(y2-y1)*len/dist,dy=cy+(x2-x1)*len/dist;
                    rep(k,N){
                        if (k==i||k==j){
                            continue;
                        }
                        double x=pos[k].first;
                        double y=pos[k].second;
                        if ((x-dx)*(x-dx)+(y-dy)*(y-dy)<=1.0){
                            count++;
                        }
                    }
                    chmax(ans,count);
                }
            }
        }
        cout << ans << endl;
    }
}

C++: #include <iostream>
#include <vector>
#include <string>
#include <complex>
#include <algorithm>
using namespace std;

typedef complex<double> p;
const double eps = 1e-9;

int n;
p ps[300];

int main() {
	while (cin >> n , n != 0) {
		for (int i = 0; i < n; i++) {
			double x, y;
			cin >> x >> y;
			ps[i] = p(x, y);
		}

		if (n == 1) {
			cout << 1 << endl;
			continue;
		}

		int ans = 0;
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				if (!(norm(ps[i] - ps[j]) < 1.0 + eps)) continue;
				p v = ps[i] - ps[j], ve = v / abs(v);
				double l = sqrt(1 - norm(v) * 0.25);
				p vn1 = ve * p(0, 1) * l, vn2 = ve * p(0, -1) * l;
				p c1 = ps[j] + v * 0.5 + vn1, c2 = ps[j] + v * 0.5 + vn2;
				int count1 = 2, count2 = 2;
				for (int k = 0; k < n; k++) {
					if (k == i || k == j) continue;
					if (norm(c1 - ps[k]) < 1.0 + eps) count1++;
					if (norm(c2 - ps[k]) < 1.0 + eps) count2++;
				}
				ans = max(ans, count1);
				ans = max(ans, count2);
			}
		}

		cout << ans << endl;
	}
}

C++: #include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <cstring>
#include <map>
#include <utility>
#include <set>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>
#include <sstream>
#include <complex>
#include <stack>
#include <queue>

#define FOR(i,b,n) for(int i=b;i<n;i++)
#define RFOR(i,b,n) for(int i=n-1;i>=b;i--)
#define CLR(mat) memset(mat, 0, sizeof(mat))
#define NCLR(mat) memset(mat, -1, sizeof(mat))

#define EPS 0.00001
#define EQ(a,b) (abs((a)-(b)) < EPS)
#define EQV(a,b) ( EQ((a).real(), (b).real()) && EQ((a).imag(), (b).imag()) )
#define PI (acos(-1.0))
#define POSARG(a) (arg((a)) > 0.0? arg((a)) : 2.0*PI + arg((a)))

using namespace std;

typedef complex<double> P;

int n;

int solve(vector<P> vc)
{	
	int res = 1;
	
	FOR(j, 0, n)
		FOR(i, j+1, n)
		{
			int tmpp = 2;
			int tmpn = 2;
			
			if(abs(vc[j]-vc[i]) > 2.0)
				continue;
			
			double cs = abs(vc[j]-vc[i])/2.0, sn = sqrt(1.0-cs*cs);	//線分IJと円の中心の距離
			P mid = (vc[j]+vc[i])/2.0;//線分IJの中点
			P nml = (vc[j]-vc[i])*P(0, 1);
			nml = nml/abs(nml);//線分IJの単位法線ベクトル
			P c1 = mid + nml*sn, c2 = mid - nml*sn;//円の中心2つ
						
			FOR(k, 0, n)
			{	
				if( k != j && k != i )
				{
					if( abs(vc[k]-c1) < 1.0+EPS)
					{
						tmpp++;
					}
					
					if( abs(vc[k]-c2) < 1.0+EPS)
					{
						tmpn++;
					}
				}
			}
			
			res = max(res, max(tmpp, tmpn));
		}
		
	return res;
}

int main()
{
	while(cin >> n, n)
	{
		vector<P> vc = vector<P>();
		
		FOR(i, 0, n)
		{
			double x, y;
			cin >> x >> y;
			vc.push_back(P(x,y));
		}

		cout << solve(vc) << endl;
	}
	
	return 0;
}

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

using namespace std;

// ??????????????§?????????????????§???3????????????????????????????±?????????????????????¨?????????

double pow2(double d){
	return d * d;
}

struct Vector2{
	double x;
	double y;
};

Vector2 P[500];

int main(){
	int n;
	while(true){
		cin >> n;
		if(n == 0){
			break;
		}
		for(int i = 0; i < n; i++){
			cin >> P[i].x >> P[i].y;
		}
		double k, L_2, a, b;
		
		int max = 0;
		
		for(int num1 = 0; num1 < n; num1++){
			for(int num2 = num1 + 1; num2 < n; num2++){
				L_2 = pow2(P[num1].x - P[num2].x) + pow2(P[num1].y - P[num2].y);
				if(L_2 >= 4){
					continue;
				}
				k = sqrt((4 - L_2) / L_2);
				
				a = (k * (P[num1].y - P[num2].y) + P[num1].x + P[num2].x) / 2;
				b = (k * (P[num2].x - P[num1].x) + P[num1].y + P[num2].y) / 2;
				
				int count;
				
				count = 0;
				
				for(int i = 0; i < n; i++){
					if(i == num1 || i == num2){
						count++;
						continue;
					}
					if(pow2(P[i].x - a) + pow2(P[i].y - b) <= 1){
						count++;
					}
				}
				
				if(max < count){
					max = count;
				}
				
				a = (k * (P[num2].y - P[num1].y) + P[num1].x + P[num2].x) / 2;
				b = (k * (P[num1].x - P[num2].x) + P[num1].y + P[num2].y) / 2;
				
				count = 0;
				
				for(int i = 0; i < n; i++){
					if(i == num1 || i == num2){
						count++;
						continue;
					}
					if(pow2(P[i].x - a) + pow2(P[i].y - b) <= 1){
						count++;
					}
				}
				
				if(max < count){
					max = count;
				}
			}
		}
		printf("%d\n", max);
	}
	return 0;
}