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<stdio.h>  
#include<algorithm>  
#include<math.h>  
using namespace std;  
const double eps=1e-8;  
struct Point {  
    double x,y;  
};  
struct Node{  
    double angle;  
    int in;  
};  
Node arc[1005];  
Point p[305];  
  
double dist(Point p,Point q){  
    return sqrt((p.x-q.x)*(p.x-q.x)+(p.y-q.y)*(p.y-q.y));  
}  
int cmp(Node a,Node b){  
    if(a.angle!=b.angle) return a.angle<b.angle;  
    else return a.in>b.in;  
}  
int main(){  
    #ifndef ONLINE_JUDGE  
    freopen("in.txt","r",stdin);  
    #endif // ONLINE_JUDGE  
    int n;  
    while(scanf("%d",&n),n){  
        int ans=1;  
        for(int i=0;i<n;i++) scanf("%lf%lf",&p[i].x,&p[i].y);  
        for(int i=0;i<n;i++){  
            int cnt=0;  
            for(int j=0;j<n;j++){  
                if(i==j) continue;  
                if(dist(p[i],p[j])>2.0+eps) continue;  
                double b=atan2(p[j].y-p[i].y,p[j].x-p[i].x);  
                double a=acos(dist(p[i],p[j])/2);  
                arc[cnt].angle=b-a; arc[cnt++].in=1; //把角度小的点作为入点  
                arc[cnt].angle=b+a; arc[cnt++].in=-1;  
            }  
            sort(arc,arc+cnt,cmp);  
            int tmp=1;  
            for(int j=0;j<cnt;j++){  
                //if(arc[j].in) tmp++;值为负时if也成立（只要非0都成立）  
                tmp+=arc[j].in;  
                ans=max(ans,tmp);  
            }  
        }  
        printf("%d\n",ans);  
    }  
    return 0;  
}

C++: #include<bits/stdc++.h>
#define _ ios_base::sync_with_stdio(0);cin.tie(0);
#define REP(i,n) for(int i=0;i<(int)(n);i++)
#define REPP(i,j,n) for(int i=(j);i<(int)(n);i++)

using namespace std;
typedef long double ld;
typedef complex<ld> Point;
const ld eps = 1e-5, pi = acos(-1.0);
ld dot(Point a, Point b){return real(conj(a)*b);}
ld cross(Point a, Point b){ return imag(conj(a)*b);}

Point p[310];
int n;

void solve(){
  int result=1;
  REP(i,n){
    REPP(j,i+1,n){
      Point ab=p[j]-p[i];
      if(norm(ab)>4)continue;
      Point c=(p[i]+p[j])/2.L+Point{imag(ab),-real(ab)}*sqrt(1-norm(ab)/4.L)/abs(ab)/2.L;
      Point d=(p[i]+p[j])/2.L+Point{-imag(ab),real(ab)}*sqrt(1-norm(ab)/4.L)/abs(ab)/2.L;
      int cntc=0,cntd=0;
      REP(k,n){
	if(abs(p[k]-c)<=1.0001L)cntc++;
	if(abs(p[k]-d)<=1.0001L)cntd++;
      }
      result=max({result,2,cntc,cntd});
    }
  }
  cout<<result<<endl;
}

int main(){ _;
  ld r,i;
  while(cin>>n,n!=0){
    REP(j,n){
      cin>>r>>i;
      p[j]=Point{r,i};
    }
    solve();
  }
}

C++: #include <iostream>
#include <algorithm>
#include <cmath>

#define rep2(i, b, n) for(int i=b; i < n; ++i)
#define rep(i, n) rep2(i, 0, n)
#define X first
#define Y second

using namespace std;

typedef pair<double, double> pdd;

double distance_sq( pdd a, pdd b )
{
	double cx = a.X - b.X;
	double cy = a.Y - b.Y;
	return ((cx*cx) + (cy*cy));
}

double norm( pdd &v )
{
	double len = sqrt((v.X*v.X) + (v.Y*v.Y));
	v.X /= len;
	v.Y /= len;
	return len;
}

void rotate( pdd &p )
{
	swap(p.X, p.Y);
	p.Y *= -1.0;
}

int main( void )
{
	int n;
	pdd p[300];
	while(cin >> n, n)
	{
		rep(i, n) cin >> p[i].X >> p[i].Y;
		int Max = 0;
		rep(i, n) {
			rep2(j, i+1, n) {
				if(distance_sq(p[i], p[j]) <= 4.0) {
					pdd v = make_pair(p[i].X-p[j].X, p[i].Y-p[j].Y);
					pdd c = make_pair(v.X / 2.0, v.Y / 2.0);
					pdd m = c;
					double len_c = norm(c);
					double sin_x = sqrt(1.0 - (len_c * len_c));
					rotate(c);
					c.X *= sin_x; c.Y *= sin_x;
					pdd p3 = make_pair(p[j].X+m.X+c.X, p[j].Y+m.Y+c.Y);
					pdd p4 = make_pair(p[j].X+m.X-c.X, p[j].Y+m.Y-c.Y);
					int count_p3, count_p4;
					count_p3 = count_p4 = 0;
					rep(k, n) {
						if(distance_sq(p3, p[k]) <= 1.0) ++count_p3;
						if(distance_sq(p4, p[k]) <= 1.0) ++count_p4;
					}
					Max = max(Max, max(count_p3, count_p4));
				}
			}
		}
		cout << Max << endl;
	}
	return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;
#define REP(i, n) FOR(i, 0, n)
#define FOR(i, a, b) for(int i=a; i<b; i++)
typedef long double ld;
typedef complex<ld> pt;
#define LAR(a, b) ((a)=max((a), (b)))

int n;
pt p[305];

int cnt(pt c, int x, int y){
	int ret=0;
	REP(i, n){
		if(i==x || i == y) ret++;
		else ret+=(abs(c-p[i])<1);
	}
	return ret;
}

int main(){
	while(scanf("%d", &n), n){
		REP(i, n){
			ld x, y; scanf("%Lf%Lf", &x, &y);
			p[i]=pt(x, y);
		}
		int ans=1;
		REP(i, n){
			FOR(j, i+1, n){
				if(abs(p[i]-p[j])>2) continue;
				pt r = (p[j]-p[i])/2.0L;
				pt s = r * (pt)1il;
				s *= sqrt(1 - abs(r)*abs(r)) / abs(s);
				LAR(ans, cnt(p[i]+r+s, i, j));
				LAR(ans, cnt(p[i]+r-s, i, j));
			}
		}
		printf("%d\n", ans);
	}
}