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 <deque>
#include <string>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <string>
#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 vector<int> vi;
typedef vector<vi> vvi;
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;
        }

        if (n == 1) {
            puts("1");
            continue;
        }
        if (n == 2){
            if (abs(v[0] - v[1]) < 2 - eps) {
                puts("2"); 

            }
            else{
                puts("1"); 
            }
            continue;
        }

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

C++: #include<cstdio>
#include<complex>

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

using namespace std;

typedef complex<double> Point;

const double EPS=1e-9;

int n;
Point p[300];

int f(const Point &c){
	int res=0;
	rep(i,n) if(abs(p[i]-c)<1+EPS) res++;
	return res;
}

int main(){
	for(;scanf("%d",&n),n;){
		rep(i,n){ double x,y; scanf("%lf%lf",&x,&y); p[i]=Point(x,y); }

		int ans=1;
		rep(j,n) rep(i,j) {
			double d=abs(p[j]-p[i])/2;
			Point v=(p[j]-p[i])/(2*d);
			if(d<1){
				ans=max(ans,f(p[i]+d*v+sqrt(1-d*d)*v*Point(0,+1)));
				ans=max(ans,f(p[i]+d*v+sqrt(1-d*d)*v*Point(0,-1)));
			}
		}
		printf("%d\n",ans);
	}

	return 0;
}

C++: #include <iostream>
#include <sstream>
#include <iomanip>
#include <algorithm>
#include <cmath>
#include <string>
#include <vector>
#include <list>
#include <queue>
#include <stack>
#include <set>
#include <map>
#include <bitset>
#include <numeric>
#include <cfloat>
using namespace std;

const double EPS = 1.0e-9;

class Point
{
public:
    double y, x;
    Point(){
        y = x = 0.0;
    }
    Point(double y0, double x0){
        y = y0;
        x = x0;
    }
    const Point operator+(const Point& p) const{
        return Point(y + p.y, x + p.x);
    }
    const Point operator-(const Point& p) const{
        return Point(y - p.y, x - p.x);
    }
    const Point operator*(double a) const{
        return Point(y * a, x * a);
    }
    double dist(const Point& p) const{
        return sqrt(pow(y-p.y, 2.0) + pow(x-p.x, 2.0));
    }
    double dist2(const Point& p) const{
        return pow(y-p.y, 2.0) + pow(x-p.x, 2.0);
    }
};

int solve(vector<Point>& vp)
{
    int n = vp.size();
    int ret = 1;

    for(int i=0; i<n; ++i){
        for(int j=i+1; j<n; ++j){
            double dist = vp[i].dist(vp[j]);
            if(dist > 2.0)
                continue;

            for(int k=0; k<2; ++k){
                Point p = vp[j] - vp[i];
                Point p1(p.x, p.y);
                if(k == 0)
                    p1.x *= -1;
                else
                    p1.y *= -1;
                p = vp[i] + p * 0.5 + p1 * (sqrt(pow(dist, -2) - 0.25));
                int num = 0;
                for(int l=0; l<n; ++l){
                    if(p.dist2(vp[l]) < 1.0 + EPS)
                        ++ num;
                }
                ret = max(ret, num);
            }
        }
    }

    return ret;
}

int main()
{
    for(;;){
        int n;
        cin >> n;
        if(n == 0)
            break;

        vector<Point> p(n);
        for(int i=0; i<n; ++i)
            cin >> p[i].x >> p[i].y;
        cout << solve(p) << endl;
    }
}