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>

#define rep(n) for(int i=0;i<n;i++)
#define repp(j, n) for(int j=0;j<n;j++)
#define reppp(i, m, n) for(int i=m;i<n;i++)
#define all(c) c.begin(), c.end()
#define rall(c) c.rbegin(), c.rend()
#define debug(x) cerr << #x << ": " << x << endl

using namespace std;

typedef long long ll;
typedef pair<ll, ll> Pll;
typedef pair<int, int> Pii;

const long double EPS = 10e-10;

class Segment;

class Point{
    public: 
    long double x,y;
    
    Point(long double x = 0.0, long double y = 0.0): x(x), y(y) {}
    
    Point operator + (Point p) {
        return Point(x+p.x, y+p.y);
    }
    Point operator - (Point p) {
        return Point(x-p.x, y-p.y);
    }
    Point operator * (long double a) {
        return Point(x*a, y*a);
    }
    Point operator / (long double a) {
        return Point(x/a, y/a);
    }
    bool operator < (const Point &p) const {
        return (x != p.x ? x < p.x : y < p.y);
    }
    bool operator == (const Point &p) const {
        return fabsl(x - p.x) < EPS && fabsl(y - p.y) < EPS;
    }
    
    long double abs() {
        return sqrt(norm());
    }
    
    long double norm() {
        return x*x + y*y;
    }
    
    long double dot(Point p) {
        return x*p.x + y*p.y;
    }
    
    long double cross(Point p) {
        return x*p.y - y*p.x;
    }
};

class Segment{
    public: 
    Point p,q;
    long double grad;
    
    Segment(Point p, Point q) {
        this->p = p;
        this->q = q;
        this->grad = (q.y-p.y)/(q.x-p.x);
    }
    
};

class Circle{
    public:
    Point c;
    long double r;
    
    Circle(Point c = Point(), long double r = 0.0): c(c), r(r) {}
};



int main() {
    int n;
    while(scanf("%d", &n) && n){
        vector<Point> points(n);
        rep(n) cin >> points[i].x >> points[i].y;

        int ans = 1;
        rep(n){
            repp(j, n){
                if(j == i) continue;
                Point mid = (points[i] + points[j]) / 2;
                Point midv = Point(mid.y, -mid.x);
                midv = midv / midv.abs() * (1.0 - (points[i] - mid).abs());
                Point center = mid + midv;
                int tmp = 0;
                repp(k, n){
                    if((points[k]-center).norm() <= 1.0) tmp++;
                }
                ans = max(tmp, ans);
            }
        }
        cout << ans << endl;
    }
}

C++: #include <iostream>
#include <cmath>
#include <algorithm>
#include <queue>
#include <map>
#define LP(i,n) for(int i = 0; i < n; i++)

using namespace std;
typedef pair<double,double> P;

double Dist(double x1[], double x2[]){
	return pow((x1[0]-x2[0])*(x1[0]-x2[0])+(x1[1]-x2[1])*(x1[1]-x2[1]),0.5);

}

int main(){
	int n;
	cin >> n;
	while(n){
		P p[300];
		LP(i,n){
			cin >> p[i].first >> p[i].second;
		}
		sort(p,p+n);
		double x[300][2];
		LP(i,n){
			x[i][0] = p[i].first;
			x[i][1] = p[i].second;
		}
		double dist[300][300];
		LP(i,n){
			LP(j,n){
				dist[i][j] = Dist(x[i],x[j]);
			}
		}
		int mx = 0;
		LP(i,n){
			for(int j = i+1; j < n; j++){
				if(x[j][0]-x[i][0] > 2){
					break;
				}else{
					if(dist[i][j] <= 2){
						double m[] = {(x[j][0]+x[i][0])/2,(x[j][1]+x[i][1])/2};
						double d[] = {x[j][1]-x[i][1],-x[j][0]+x[i][0]};
						double O[] = {0,0};
						double _d = Dist(d,O);
						LP(k,2){
							d[k] *= (pow(_d,-2)-0.25);
						}
						double o[2][2];
						LP(k,2){
							LP(l,2){
								o[k][l] = m[l];
								if(k) o[k][l] += d[l];
								else o[k][l] -= d[l];
							}
						}
						LP(k,2){
							int s = 2;
							for(int l = i+1; l < j; l++){
								if(Dist(x[l],o[k]) <= 1) s++;
							}
							for(int l = i-1; 0<=l && x[i][0]-x[l][0] <= 1;l--){
								if(Dist(x[l],o[k]) <= 1) s++;
							}
							for(int l = j+1; l < n && x[l][0]-x[j][0] <= 1; l++){
								if(Dist(x[l],o[k]) <= 1) s++;
							}
							mx = max(mx,s);
						}
					}
				}
			}
		}
		cout << mx << endl;
		cin >> n;
	}
}

C++: #include <bits/stdc++.h>
using namespace std;

#define rep(i,N) for(int i=0;i<(int)(N);i++)
#define REP(i,k,N) for(int i=k;i<(int)(N);i++)

#define vi vector<int>
#define pb push_back

#define pdd pair<double,double>
#define fi first
#define se second
const double  EPS = 1e-9;

double dist(double x,double y){
    return x*x+y*y;
};

int main(){
    int n;
    while(cin >> n,n){
        int ma = 0;
        pdd p[330];
        rep(i,n) cin >> p[i].fi >> p[i].se;
        rep(i,n-1){
            REP(j,i+1,n){
                int cnt1 = 0,cnt2=0;
                double x1,y1,x2,y2,x3,y3;
                tie(x1,y1) = p[i];
                tie(x2,y2) = p[j];
                double cx1,cy1,cx2,cy2;
                if(dist(x1-x2,y1-y2)>4) continue;

                //はりまコーナー
//                double ax,ay,bx1,by1,bx2,by2,mx,my,dx,dy;
//                ax=x1-x2;ay=y1-y2;
//                mx=(x1+x2)/2;my=(y1+y2)/2;
//                by1 = (-0.5)*pow(ax,4)+pow(ax,2)+(-0.5)*pow(ax,2)*pow(ay,2)
//                        /dist(ax,ay);
//                by1 = sqrt(abs(by1));
//                by2 = -1*by1;
//                cy1= my+by1;
//                cy2= my+by2;
//                bx1= -1*(ay*by1)/ax;
//                bx2= -1*(ay*by2)/ax;
//                cx1=bx1+mx;
//                cx2=bx2+mx;


                //あじのりコーナー
                double mx=(x1+x2)/2,my=(y1+y2)/2;

                double ld = 1-dist((x1-x2)/2,(y1-y2)/2);
                ld = sqrt(ld);
                cx1 = mx - sin(atan2(y2-y1,x2-x1))*ld;
                cy1 = my + cos(atan2(y2-y1,x2-x1))*ld;
                cx2 = mx + sin(atan2(y2-y1,x2-x1))*ld;
                cy2 = my - cos(atan2(y2-y1,x2-x1))*ld;

                for(int k = 0;k<n;k++) {
                    tie(x3, y3) = p[k];
                    if (EPS > dist((cx1 - x3), (cy1 - y3))-1) cnt1++;
                    if (EPS > dist((cx2 - x3), (cy2 - y3))-1) cnt2++;
                }
                ma = max({cnt1,cnt2,ma});
            }
        }
        cout<< ma << endl;
    }
    return 0;
}