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;
const int dx[]={0,1,0,-1,1,-1,-1,1};
const int dy[]={-1,0,1,0,1,1,-1,-1};
const int INF = 1<<30;
const double EPS = 1e-15;
#define PB push_back
#define mk make_pair
#define fi first
#define se second
#define ll long long
#define reps(i,j,k) for(int i = (j); i < (k); i++)
#define rep(i,j) reps(i,0,j)
#define MOD 1000000007
typedef pair<int,int> Pii;
typedef pair<int,Pii> P;
typedef vector<int> vi;
typedef vector<vi> vvi;
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 * (double k){
        return Point(k*x,k*y);
    }
     
    Point operator * (Point p){
        return Point(p.x*x-p.y*y,p.x*y+p.y*x);
    }

    Point operator / (double k){
        return Point(x/k,y/k);
    }
 	
    Point operator = (Point p){
    	x = p.x;
    	y = p.y;
    	return *this;
    }

    bool operator<(const Point &p)const{
        return x!=p.x?x<p.x:y<p.y;
    }

};
double dot(Point p, Point q){
    return p.x*q.x+p.y*q.y;
}
  
double cross(Point p,Point q){
    return p.x*q.y-q.x*p.y;
}
  
double norm(Point p){
	//ベクトルの大きさの2乗
    return p.x*p.x+p.y*p.y;
}
  
double pabs(Point p){
	//ベクトルの大きさ
    return sqrt(norm(p));
}
Point UnitVector(Point a){
	//単位ベクトル
	return a/pabs(a);
}
typedef vector<Point> Polygon;
struct Line : Polygon{
    Line(){}
    Line(Point p,Point q){
        PB(p);
        PB(q);  
    }
};
struct Circle{
    Point p; //中心点
    double r; //半径
    Circle(){}
    Circle(Point _p,double _r){
        p = _p;
        r = _r;
    }
};
int ccw(Point a,Point b,Point c){
    b = b-a;
    c = c-a;
    if(cross(b,c) > 0)return +1; //a→b で反時計に折れてc ccw
    if(cross(b,c) < 0)return -1; //a→b で時計に折れてc ccw
    if(dot(b,c) < 0)return +2; // c--a--b;
    if(pabs(b)<pabs(c))return -2; //a--b--c
    return 0; //a--c--b
}
bool intersectSS(Line l,Line m){
    //lとmが交差してるかどうか
    return (ccw(l[0],l[1],m[0])*ccw(l[0],l[1],m[1])<EPS)&&
            (ccw(m[0],m[1],l[0])*ccw(m[0],m[1],l[1])<EPS);
}

bool intersectSP(Line l, Point p){
	//線分lと点pが交差してるか
	return !ccw(l[0],l[1],p);
}

bool intersectLL(Line l,Line m){
	return fabs(cross(l[1]-l[0],m[1]-m[0])) > EPS || //平行でない
			fabs(cross(l[1]-l[0],m[0]-l[0])) < EPS; //同一の線
}

int intercectCC(Circle a,Circle b,Point *p,Point *q){
	double dist = pabs(a.p-b.p);
	if(dist < fabs(a.r-b.r) || dist > a.r+b.r){
		//円の中心が一致しているとき
		//円の交点が存在しないとき
		*p = Point(1e15,1e15);
		*q = Point(1e15,1e15);
		return 0;
	}
	else if(fabs(dist-a.r-b.r) < EPS){
		//円が1点で交わっているとき
		Point way;
		way = b.p-a.p;
		*p = a.p+UnitVector(way)*a.r;
		*q = *p;
		return 1;
	}
	else{
		//円が2点で交わっているとき
        if(a.r < b.r){
            swap(a.r,b.r);
            swap(a.p,b.p);
        }
        double s = (a.r + b.r + dist) / 2.0 ;
        double area = sqrt(s * (s - a.r) * (s - b.r) * (s - dist));
        double h = 2.0 * area / dist;

        Point  v = b.p-a.p;
        v = UnitVector(v);
        Point  m = a.p + (v*sqrt(a.r*a.r-h*h));
        Point  n = v * Point(0.0,1.0);
        
        *p = m + (n * h);
        *q = m - (n * h);
        return 2;
	}
}

double TriangleArea(Point a,Point b,Point c){
    b = b - a;
    c = c - a;
    return cross(b,c)*0.5; //負になることもあるので注意
}
 
double VertexArea(Polygon v){
    double subArea = 0.0;
    rep(i,v.size()-1){
        subArea += TriangleArea(v[0],v[i],v[i+1]);
    }
    return subArea; //負になることもあるので注意
}

Point crosspoint(Line l,Line m){
	//線分lと線分mの交点
    if(intersectSS(l,m) == false)return Point(INF,INF);
    return m[0]+(m[1]-m[0])*cross(l[1]-l[0],l[1]-m[0])/cross(l[1]-l[0],m[1]-m[0]);
}
 
Polygon ConvexHull(Polygon v){
	int n = v.size();
	int k = 0;
	sort(v.begin(),v.end());
	Polygon e(2*n);
	for(int i=0;i<n;e[k++]=v[i++])
		while(k>=2&&ccw(e[k-2],e[k-1],v[i])<=0)k--;
	for(int i=n-2,t=k+1;i>=0;e[k++]=v[i--])
		while(k>=t&&ccw(e[k-2],e[k-1],v[i])<=0)k--;
	e.resize(k-1);
	return e;
}

int main(){
    int N;
    while(scanf("%d",&N),N){
        Circle C[128];
        vector< Point > points;
        rep(i,N){
            double x,y;
            scanf("%lf %lf",&x,&y);
            C[i] = Circle(Point(x,y),1.);
        }
        int len = 0;
        rep(i,N){
            reps(j,i+1,N){
                Point p,q;
                int tmp = 0;
                tmp = intercectCC(C[i],C[j],&p,&q);
                if(tmp > 0){
                    points.PB(p);
                    points.PB(q);
                    /*
                    cout << i << " " << j << "->\n";
                    cout << "1:"<< p.x << " " << p.y << "\n";
                    cout << "2:"<< q.x << " " << q.y << "\n";
                    */
                }
            }
        }
        int ans = -INF;
        rep(i,points.size()){
            Point tar = points[i];
            int cnt = 0;
            rep(j,N){
                Point base = C[j].p;
                //cout << pabs(base-tar) << " ";
                if(pabs(base-tar)-1.<EPS){
                    cnt++;
                }
            }
            ans = max(ans,cnt);
        }

        printf("%d\n",ans);

    }
    return 0;
}

C++: #include<iostream>
#include<complex>
#define rep(i,n) for(int i=0;i<n;i++)
#define EPS 1e-8
using namespace std;
typedef complex<double> P;
int main(){
	int n; double x,y;
	while(cin>>n,n){
		P p[n]; rep(i,n)cin>>x>>y,p[i]=P(x,y);
		int mx=0,m1,m2;
		rep(i,n)rep(j,i)if(norm(p[i]-p[j])<4+EPS){
			P M=(p[i]+p[j])/P(2,0),C1,C2; double l=sqrt(1-norm(M-p[i]));
			C1=M+polar(l,arg(p[j]-p[i])+M_PI/2);
			C2=M+polar(l,arg(p[j]-p[i])-M_PI/2);
			m1=m2=0;
			rep(k,n){
				if(norm(p[k]-C1)<1+EPS)m1++; if(norm(p[k]-C2)<1+EPS)m2++;
			}
			if(m1>mx)mx=m1; if(m2>mx)mx=m2;
		}
		cout<<mx<<endl;
	}
	return 0;
}