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 <algorithm>/*{{{*/
#include <cassert>
#include <cctype>
#include <climits>
#include <cmath>
#include <complex>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <list>
#include <map>
#include <numeric>
#include <queue>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <vector>
using namespace std;

static const double PI (3.14159265358979323846);
static const double EPS (1e-10);

typedef long long int ll;
typedef unsigned long long int ull;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef vector<string> vs;
typedef pair<int, int> pii;
typedef complex<double> P;

inline int toInt(string s) {int v; istringstream sin(s);sin>>v;return v;}
template<class T> inline string toString(T x) {ostringstream sout;sout<<x;return sout.str();}

#define FOR(i, b, e) for (typeof(b) i = (b); assert((i) <= (e)), i != (e); ++i)
#define REP(i, n) FOR(i, 0, n)

#define OPOVER(_op, _type) inline bool operator _op (const _type &t) const

#define arrsz(a) ( sizeof(a) / sizeof(a[0]) )

#define F first
#define S second
#define MP(a, b) make_pair(a, b)

#define SZ(a) ((ll)a.size())
#define PB(e) push_back(e)
#define SORT(v) sort((v).begin(), (v).end())
#define RSORT(v) sort((v).rbegin(), (v).rend())

#define ALL(a) (a).begin(), (a).end()
#define RALL(a) (a).rbegin(), (a).rend()
#define EACH(t,i,c) for(t::iterator i=(c).begin(); i!=(c).end(); ++i)
#define EXIST(s,e) ((s).find(e)!=(s).end())

#define BIT(n) (1ULL << (n))
#define BITOF(n, m) ((n) >> (m) & 1)

#define RANGE(a, b, c) ((a) <= (b) && (b) <= (c))

#define DEBUG_ON

#ifdef DEBUG_ON
	#define dprt(fmt, ...) fprintf(stderr, fmt, __VA_ARGS__)
	#define darr(a) copy( (a), (a) + arrsz(a), ostream_iterator<int>(cerr, " ") ); cerr << endl
	#define darr_range(a, f, t) copy( (a) + (f), (a) + (t), ostream_iterator<int>(cerr, " ") ); cerr << endl
	#define dvec(v) copy( ALL(v), ostream_iterator<int>(cerr, " ") ); cerr << endl
	#define darr2(a, n, m) FOR(i, 0, (n)){ darr_range( (a)[i], 0, (m) ); }
	#define dvec2(v) FOR(i, 0, SZ(v)){ dvec( (v)[i] ); }
	#define WAIT() { string _wait_; cerr << "(hit return to continue)" << endl; getline(cin, _wait_); }
#else
	#define dprt(fmt, ...)
	#define darr(a)
	#define darr_range(a, f, t)
	#define dvec(v)
	#define darr2(a, n, m)
	#define dvec2(v)
	#define WAIT()
#endif
/*}}}*/

int N;
vector<P> points;

P middle_point(P a, P b)
{
	return a + (b-a) / 2.0;
}

// 点pを中心としてr(radian)回転 p(0,0)で原点を中心として回転
P rotate(P t, P p ,double r){
	// double r=radians(angle);
	return (t - p) * P(cos(r), sin(r)) + p;
}

// 2円 |x-a|=raと|x-b|=rbの交点計算
bool intersection_c_c(P a, double ra, P b, double rb, P ans[]){
	double di = abs(a-b);
	if(di > ra+rb || di < abs(ra-rb)) return false;
	double t=(ra*ra - rb*rb + di*di) / (di+di);
	double rd = acos(t/ra);

	P dv = (b-a) / abs(b-a);
	P g1 = rotate(dv, P(0,0), rd);
	P g2 = rotate(dv, P(0,0), -rd);
	ans[0] = a + g1*ra;
	ans[1] = a + g2*ra;
	return true;
}

int main(int argc, char const* argv[])
{
	while (cin >> N, N) {
		int ans = 1;
		points.clear();
		REP(i, N) {
			double x, y;
			cin >> x >> y;
			points.PB(P(x, y));
		}
		REP(i, N) {
			FOR(j, i+1, N) {
				P p1 = points[i], p2 = points[j];
				dprt("(%lf, %lf) and (%lf, %lf)\n", p1.real(), p1.imag(), p2.real(), p2.imag());
				P centers[2];
				if (intersection_c_c(p1, 1, p2, 1, centers)) {
					REP(k, 2) {
						dprt(" center: (%lf, %lf)\n", centers[k].real(), centers[k].imag());
						int cnt = 0;
						REP(l, N) {
							if (abs(centers[k] - points[l]) < 1.0 + EPS) {
								++cnt;
							}
						}
						ans = max(ans, cnt);
					}
				}
			}
		}
		cout << ans << endl;
	}
	return 0;
}

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

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;
        }
        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 <bits/stdc++.h>
using namespace std;
#define rep(i,n) for(int i=0;i<n;++i)
struct Point{double x,y;};
double dist(Point a,Point b) {return sqrt(pow(a.x-b.x,2)+pow(a.y-b.y,2));}
#define eps 0.0001

int main(void){
    int N;
    while(cin>>N,N){
        vector<Point> p(N);
        int maxcnt=1;
        rep(i,N)cin>>p[i].x>>p[i].y;
        rep(i,N){
            for(int j=i+1;j<N;j++){
                double d = dist(p[i],p[j]);
                if(d>2.0)continue;
                Point C1,C2;
                C1.x = p[i].x + cos(atan2(p[j].y-p[i].y,p[j].x-p[i].x) + acos(d/2.0));
                C1.y = p[i].y + sin(atan2(p[j].y-p[i].y,p[j].x-p[i].x) + acos(d/2.0));
                C2.x = p[i].x + cos(atan2(p[j].y-p[i].y,p[j].x-p[i].x) - acos(d/2.0));
                C2.y = p[i].y + sin(atan2(p[j].y-p[i].y,p[j].x-p[i].x) - acos(d/2.0));
                int cnt=0;
                rep(k,N){
                    if(pow(C1.x-p[k].x,2) + pow(C1.y-p[k].y,2) <= 1.0+eps)cnt++;
                }
                if(maxcnt<cnt)maxcnt=cnt;
                cnt=0;
                rep(k,N){
                    if(pow(C2.x-p[k].x,2) + pow(C2.y-p[k].y,2) <= 1.0+eps)cnt++;
                }
                if(maxcnt<cnt)maxcnt=cnt;
            }
        }
        cout<<maxcnt<<endl;
    }

    return 0;
}