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;
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define each(itr,c) for(__typeof(c.begin()) itr=c.begin(); itr!=c.end(); ++itr)
#define all(x) (x).begin(),(x).end()
#define mp make_pair
#define pb push_back
#define fi first
#define se second

const double eps=1e-8;
typedef pair<double,double> p;

double dist(p a, p b)
{
    return sqrt((a.fi-b.fi)*(a.fi-b.fi)+(a.se-b.se)*(a.se-b.se));
}

double size(p a)
{
    return sqrt(a.fi*a.fi+a.se*a.se);
}

int main()
{
    int n;
    while(cin >>n,n)
    {
        vector<p> v(n);
        rep(i,n) scanf(" %lf %lf", &v[i].fi, &v[i].se);

        if(n==1)
        {
            printf("1\n");
            continue;
        }

        int ans=0;
        //select 2 points
        rep(i,n)rep(j,i)
        {
            double d=dist(v[i],v[j]);
            if(d>2.0) continue;

            double r=sqrt(1.0 - (d*d/4.0));

            double X=v[j].fi-v[i].fi, Y=v[j].se-v[i].se;

            p v1=p(X,Y);
            double v1size=size(v1);
            v1.fi/=v1size; v1.se/=v1size;
            v1.fi*=d/2; v1.se*=d/2;

            p v2;
            double v2size;
            p c;

            //2 candidate for center of circle
            //first
            int ct=0;
            v2=p(-Y,X);
            v2size=size(v2);
            v2.fi/=v2size; v2.se/=v2size;
            v2.fi*=r; v2.se*=r;

            c.fi=v[i].fi+v1.fi+v2.fi; c.se=v[i].se+v1.se+v2.se;
            assert(dist(c,v[i]));
            rep(k,n)
            {
                if(k==i || k==j) ++ct;
                else
                {
                    if(dist(c,v[k])<1.0) ++ct;
                }
            }
            ans=max(ans,ct);

            //second
            ct=0;
            v2=p(Y,-X);
            v2size=size(v2);
            v2.fi/=v2size; v2.se/=v2size;
            v2.fi*=r; v2.se*=r;

            c.fi=v[i].fi+v1.fi+v2.fi; c.se=v[i].se+v1.se+v2.se;
            assert(dist(c,v[i]));
            rep(k,n)
            {
                if(k==i || k==j) ++ct;
                else
                {
                    if(dist(c,v[k])<1.0) ++ct;
                }
            }
            ans=max(ans,ct);
        }

        cout << ans << endl;
    }
    return 0;
}

C++: #include <bits/stdc++.h>
 
#define loop(n, i) for(int i=0;i<n;i++)
#define all(v) v.begin(),v.end()

using namespace std;
typedef pair<double, double> P;
const double EPS = 1e-4;

#include <bits/stdc++.h>
 
#define loop(n, i) for(int i=0;i<n;i++)
#define all(v) v.begin(),v.end()

using namespace std;

vector<P> center(P a, P b)
{
    auto t = atan2(b.second - a.second, b.first - a.first);
    auto x = hypot(a.first - b.first, a.second - b.second) / 2;
    auto y = sqrt(1.0 - x * x);
    auto cs = cos(t), sn = sin(t);
    return { 
        { cs * x - sn * y + a.first, sn * x + cs * y + a.second },
        { cs * x + sn * y + a.first, sn * x - cs * y + a.second }
    };
}

int main()
{
    while (1) {
        int n; cin >> n;
        if (!n) break;
        vector<P> pnt;
        loop (n, i) {
            double x, y; cin >> x >> y;
            pnt.emplace_back(x, y);
        }
        int ans = 1;
        loop (n, i) loop (n, j) {
            auto ps = center(pnt[i], pnt[j]);
            for (auto o : ps) {
                int cnt = 0;
                loop (n, l) {
                    auto x = o.first - pnt[l].first;
                    auto y = o.second - pnt[l].second;
                    if (x * x + y * y < 1 + EPS) cnt++;
                }
                ans = max(ans, cnt);
            }
        }
        cout << ans << endl;
    }
    return 0;
}