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;
#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;
                int sign[]={-1,1};
                rep(s,2){
                    Point C;
                    C.x = p[i].x + cos(atan2(p[j].y-p[i].y,p[j].x-p[i].x) + sign[s]*acos(d/2.0));
                    C.y = p[i].y + sin(atan2(p[j].y-p[i].y,p[j].x-p[i].x) + sign[s]*acos(d/2.0));
                    int cnt=0;
                    rep(k,N) cnt += (pow(C.x-p[k].x,2) + pow(C.y-p[k].y,2) <= 1.0+eps);
                    if(maxcnt<cnt)maxcnt=cnt;
                }
            }
        }
        cout<<maxcnt<<endl;
    }

    return 0;
}

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

        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>
using namespace std;
#define FOR(i,a,b) for(int i=(a);i<(b);++i)
#define rep(i,n)   FOR(i,0,n)
#define pb emplace_back
typedef long long ll;
typedef pair<int,int> pint;

#define eps (1e-10)
struct Point{
    double x,y;
    Point(){}
    Point(double x,double y):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*(Point p) {return Point(x*p.x-y*p.y,x*p.y+y*p.x);}
    Point operator*(double k){return Point(x*k,y*k);}
    double norm(){return x*x+y*y;}
    double abs(){return sqrt(norm());}
    bool operator == (const Point &p) const{return fabs(x-p.x)<eps && fabs(y-p.y)<eps;}
    double arg(){return atan2(y,x);}
    double dot(Point p){return x*p.x+y*p.y;}
    double det(Point p){return x*p.y-y*p.x;}
};
bool cmp_x(const Point& p,const Point& q){
  if(p.x!=q.x) return p.x<q.x;
  return p.y<q.y;
}
struct Circle{
    double r;
    Point p;
    Circle(){}
    Circle(Point p,double r):p(p),r(r){}
};
bool isIntersectCC(Circle c1,Circle c2){
    return (c1.p-c2.p).abs()<=c1.r+c2.r+eps;
}
pair<Point,Point> CrossPointsCC(Circle c1,Circle c2){
    assert(isIntersectCC(c1,c2));
    double d=(c1.p-c2.p).abs();
    double k=acos((d*d+c1.r*c1.r-c2.r*c2.r)/(c1.r*d*2));
    return make_pair(c1.p+(c2.p-c1.p)*Point(cos(k),sin(k))*(c1.r/d),c1.p+(c2.p-c1.p)*Point(cos(-k),sin(-k))*(c1.r/d));
}
Point p[301];
int main(){
    int n;
    double x,y;
    while(cin>>n,n){
        rep(i,n){
            cin>>x>>y;
            p[i]=Point(x,y);
        }
        int mx=1;
        rep(i,n)FOR(j,i+1,n){
            if((p[i]-p[j]).abs()>2) continue;
            pair<Point,Point> pp=CrossPointsCC(Circle(p[i],1.0),Circle(p[j],1.0));
            int cnt=0,cnt2=0;
            rep(k,n){
                if((pp.first-p[k]).abs()<=1.0) ++cnt;
                if((pp.second-p[k]).abs()<=1.0) ++cnt2;
            }
            mx=max(mx,max(cnt,cnt2));
        }
        cout<<mx<<endl;
    }
    return 0;
}