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<string>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#include<stack>
#include<queue>
#include<cmath>
#include<numeric>
#include<cassert>
#include<limits>
#include<iomanip>
using namespace std;
typedef unsigned long long int ull;
typedef vector<char> vc;
typedef vector<vc> vvc;



typedef pair<double,double> Point;
#define F first
#define S second
#define MK make_pair
#define equal abs(a-b)<=1e-8

Point pls(Point a,Point b){
    return MK(a.F+b.F,a.S+b.S);
}

Point mns(Point a,Point b){
    return MK(a.F-b.F,a.S-b.S);
}

Point T(Point t){
    return MK(t.S,-1.0*t.F);
}

double abs(Point a){
    return sqrt(a.F*a.F+a.S*a.S);
}

double norm(Point a){
    return a.F*a.F+a.S*a.S;
}

Point mid(Point a,Point b){
    return MK((a.F+b.F)/2,(a.S+b.S)/2);
}

bool in(Point c,Point t){
    if(norm(mns(t,c))<=1){return true;}
    return false;
}


int main(){
    while(1){
        int n;
        cin>>n;
        if(n==0){break;}
        vector<Point> vp;
        int mx=0;
        for(int i=0;i<n;i++){
            double a,b;
            cin>>a>>b;
            vp.push_back(MK(a,b));
        }
        sort(vp.begin(),vp.end());
        for(int i=0;i<vp.size();i++){
            for(int t=i+1;t<vp.size() && vp[t].F-vp[i].F<=2;t++){
                if(abs(mns(vp[t],vp[i]))>2){continue;}
                Point m=mid(vp[i],vp[t]);
                Point l=mns(vp[t],vp[i]);
                Point lt=T(l);
                double h=sqrt(1-abs(l)*abs(l)/4.0);
                lt.F=h*lt.F/abs(l);
                lt.S=h*lt.S/abs(l);
                Point c=pls(m,lt);
                int count=2;
                for(int k=0;k<vp.size();k++){
                    if(k!=i && k!=t && in(c,vp[k])){count++;}
                }
                mx=max(mx,count);
                c=mns(m,lt);
                count=2;
                for(int k=0;k<vp.size();k++){
                    if(k!=i && k!=t && in(c,vp[k])){count++;}
                }
                mx=max(mx,count);
            }
        }
        cout<<mx<<endl;
    }
    return 0;
}

C++: #include<iostream>
#include<vector>
#include<complex>
#define EPS 0.0001
using namespace std;
typedef double elem;
typedef complex<elem> point,vec;
inline elem dist(const point &a, const point &b){return abs(a-b);}
inline vec normr(const vec &v){return v*vec(0,-1);}
inline vec norml(const vec &v){return v*vec(0,1);}
inline vec unormr(const vec &v){return normr(v) / abs(v);}
inline vec unorml(const vec &v){return norml(v) / abs(v);}
int main(){
	while(true){
		int ans = 0;
		int n;
		scanf("%d", &n);
		if(n==0)break;
		if(n==1){
			printf("1\n");
			continue;
		}
		vector< point > points;
		for(int i=0;i<n;++i){
			elem x, y;
			scanf("%Lf%Lf", &x, &y);
			points.push_back( point(x,y) );
		}
		for(unsigned int i = 0; i < points.size(); ++i){
			for(unsigned int j = 0; j < points.size(); ++j){
				if( i != j && dist( points[i], points[j] ) < 2 ){
					int count_p = 0;
					int count_q = 0;
					elem d = sqrt(1-pow(dist(points[i], points[j])/2,2 ) );
					point p,q;
					vec ab = points[j] - points[i];
					point t = 0.5 * ab + points[i];  
					p = d * unormr(ab) + t;
					q = d * unorml(ab) + t;
					for(unsigned int k = 0; k < points.size(); ++k){
						if( dist(p, points[k]) <= 1.0+EPS ){ count_p++; }
						if( dist(q, points[k]) <= 1.0+EPS ){ count_q++; }
					}
					ans = max(ans, max(count_q, count_p) );
				}
			}
		}
		printf("%d\n", ans);
	}
	return 0;
}

C++: #include<iostream>
#include<iomanip>
#include<vector>
#include<cmath>
using namespace std;

const double EPS = 1e-4;

#define x first
#define y second

void normalize(double &a, double &b){
    double norm = sqrt(a*a + b*b);
    a /= norm;
    b /= norm;
}

int main(){
    int n;
    while(cin >> n, n){
        int ans = 1;
        vector<pair<double,double>> v;
        for(int i = 0; i < n; i++){
            double a, b;
            cin >> a >> b;
            v.push_back({a,b});
        }
        for(int i = 0; i < n; i++){
            for(int j = i+1; j < n; j++){
                double dx = v[i].x-v[j].x, dy = v[i].y-v[j].y;
                if(sqrt(dx*dx+dy*dy) > 2.0000)    continue;
                double midx = (v[i].x+v[j].x)/2, midy = (v[i].y+v[j].y)/2;
                double ux = 0, uy = 0, cx, cy;
                if(dx == 0)         ux = 1;
                else if(dy == 0)    uy = 1;
                else                ux = 1, uy = -dx/dy;
                normalize(ux, uy);
                double scale = sqrt(1*1 - (v[i].x-midx)*(v[i].x-midx) - (v[i].y-midy)*(v[i].y-midy));
                for(int k = -1; k <= 1; k++){
                    if(k == 0)  continue;
                    cx = midx + scale * ux * k;
                    cy = midy + scale * uy * k;
                    int tmp = 0;
                    for(int l = 0; l < n; l++){
                        if((v[l].x-cx)*(v[l].x-cx)+(v[l].y-cy)*(v[l].y-cy) <= 1+EPS)    tmp++;
                    }
                    ans = max(ans, tmp);
                }
            }
        }
        cout << ans << endl;
    }
    return 0;
}

C++: #include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;

#define MAX_N 310

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

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

				long double X1, Y1, X2, Y2, X3, Y3, VX2, VY2, Dist, Dist2;

				/*Dist = sqrtl((x[j] - x[i])*(x[j] - x[i]) + (y[j] - y[i])*(y[j] - y[i]));
				Dist2 = sqrtl(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) {
						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) {
						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;
}