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 <iostream>
#include <algorithm>
#include <cmath>

using namespace std;

const double EPS=1e-6;

int main(void){
	while(1){
		int n; cin >> n;
		if(!n) break;

		double x[333], y[333];
		for(int i=0; i<n; i++){
			cin>>x[i]>>y[i];
		}

		int ans=1;
		for(int i=0; i<n; i++){
			for(int j=0; j<n; j++){
				if(i==j) continue;
				double dx = x[i]-x[j];
				double dy = y[i]-y[j];
				double dist = sqrt(dx*dx+dy*dy);
				if(dist > 2.0+EPS) continue;

				double r = sqrt(1.0-dist*dist/4);
				double vx = -dy/dist*r;
				double vy = dx/dist*r;
				double xx = (x[i]+x[j])/2+vx;
				double yy = (y[i]+y[j])/2+vy;
				int cnt=0;
				for(int k=0; k<n; k++){
					if((x[k]-xx)*(x[k]-xx)+(y[k]-yy)*(y[k]-yy) < 1.0+EPS)cnt++;
				}
				ans = max(ans, cnt);
			}
		}
		cout << ans << endl;
	}
	return 0;
}

C++: #include<iostream>
#include<string>
#include<algorithm>
#include<vector>
#include<iomanip>
#include<math.h>
#include<complex>
#include<queue>
#include<deque>
#include<map>
#include<set>
#include<bitset>
using namespace std;
#define REP(i,m,n) for(int i=(int)m ; i < (int) n ; i++ )
#define rep(i,n) REP(i,0,n)
typedef long long ll;
typedef pair<int,int> pint;
const int inf=1e9+7;
const ll longinf=1LL<<60 ;
const ll mod=1e9+7 ;
int dx[4]={1,0,-1,0} , dy[4]={0,1,0,-1} ;

const double eps=1e-10;
double add(double a,double b){
    if(abs(a+b)<eps*(abs(a)+abs(b)))return 0;
    return a+b;
}

typedef pair<double,double> P;
P operator+(const P& x,const P& y){
    return {add(x.first,y.first),add(x.second,y.second)};
}

P operator-(const P& x,const P& y){
    return {add(x.first,-y.first),add(x.second,-y.second)};

}

P operator*(const P& x,const double& y){
    return {x.first*y,x.second*y};
}
double dot(P x,P y){
    return add(x.first*y.first,x.second*y.second);


}

double det(P x,P y){
    return add(x.first*y.second,-x.second*y.first);
}

bool onseg(P p1,P p2,P q){
    return det(p1-q,p2-q)==0 &&dot(p1-q,p2-q)<=0; 
}

P intersection(P p1,P p2,P q1,P q2){
    return p1+(p2-p1)*(det(q2-q1,q1-p1)/det(q2-q1,p2-p1));
}

bool cross(P p1,P p2,P q1,P q2){
    if(det(p1-p2,q1-q2)==0)return onseg(p1,p2,q1)||onseg(p1,p2,q2)||onseg(q1,q2,p1)||onseg(q1,q2,p2);
    P ret=intersection(p1,p2,q1,q2);
    return onseg(p1,p2,ret)&&onseg(q1,q2,ret);
}

double dist(P p1,P p2,P q){
    if(dot(q-p1,p2-p1)<=0)return sqrt(dot(q-p1,q-p1));
    if(dot(q-p2,p1-p2)<=0)return sqrt(dot(q-p2,q-p2));
    P d=p1-p2;
    P ret=q+P{d.second,-d.first};
    P h=intersection(p1, p2, q, ret);
    return sqrt(dot(q-h,q-h));
}
double d(P p,P q){
    return sqrt(dot(p-q,p-q));
}

int main(){
    int n;
    while(cin>>n,n!=0){
        P a[n];
        rep(i,n)cin>>a[i].first>>a[i].second;
        int ans=0;
        rep(i,n)rep(j,n){
            if(i==j)continue;
            if(d(a[i],a[j])>2)continue;
            P v=(a[i]-a[j])*(1.0/d(a[i],a[j]));
            P h=P{v.second,-v.first};
            P c=(a[i]+a[j])*(1.0/2)+h*sqrt(1.0-dot(a[i]-a[j],a[i]-a[j])/4);
            int ret=2;
            rep(k,n){
                if(k==i||k==j)continue;
                if(d(a[k],c)<=1)ret++;
            }
            ans=max(ans,ret);
            c=(a[i]+a[j])*(1.0/2)-h*sqrt(1.0-dot(a[i]-a[j],a[i]-a[j])/4);
            ret=2;
            rep(k,n){
                if(k==i||k==j)continue;
                if(d(a[k],c)<=1)ret++;
            }
            ans=max(ans,ret);
        }
        cout<<ans<<endl;
    }
    return 0;
}

C++: #include <cstdio>
#include <complex>
#include <vector>
using namespace std;
typedef complex<double> P;
const double PI=3.141592653589793238462643383279;

bool inC(P p, P c, double r){return abs(p-c)<r+1e-9;}

P getE(P p){return p/abs(p);}

double sq(double x){ return  x*x;}

int main(){
    int n,i,j,k,ans,cnt;
    double x, y;
    while(scanf("%d",&n)){
        if(n==0) break;
        ans=0;
        vector<P> ps(n);
        for(i=0; i<n; i++){
            scanf("%lf%lf",&x,&y);
            ps[i] = P(x,y);
        }
        for(i=0;i<n;i++)
            for(j=i+1;j<n;j++)if(abs(ps[i]-ps[j])<2){
                double d=sqrt(1-sq(abs(ps[j]-ps[i])*0.5));
                P c1=0.5*(ps[i]+ps[j])+getE(getE(ps[j]-ps[i]))*polar(d,PI/2);
                P c2=0.5*(ps[i]+ps[j])+getE(getE(ps[j]-ps[i]))*polar(d,-PI/2);
                cnt=0;
                for(k=0;k<n;k++)if(inC(ps[k],c1,1.0))cnt++;
                if(cnt>ans)ans=cnt;
                cnt=0;
                for(k=0;k<n;k++)if(inC(ps[k],c2,1.0))cnt++;
                if(cnt>ans)ans=cnt;
            }
        if(ans==0) ans++;
        printf("%d\n",ans);
    }
    return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;
#define all(c) (c).begin(),(c).end()
#define rrep(i,n) for(int i=(int)(n)-1;i>=0;i--)
#define REP(i,m,n) for(int i=(int)(m);i<(int)(n);i++)
#define rep(i,n) REP(i,0,n)
#define iter(c) __typeof((c).begin())
#define tr(it,c) for(iter(c) it=(c).begin();it!=(c).end();it++)
#define mem(a) memset(a,0,sizeof(a))
#define pd(a) printf("%.10f\n",a)
#define pb(a) push_back(a)
#define in(a) insert(a)
#define pi M_PI
#define R cin>>
#define F first
#define S second
#define C class
#define ll long long
#define ln cout<<'\n'
template<C T>void pr(T a){cout<<a;ln;}
template<C T,C T2>void pr(T a,T2 b){cout<<a<<' '<<b;ln;}
template<C T,C T2,C T3>void pr(T a,T2 b,T3 c){cout<<a<<' '<<b<<' '<<c;ln;}
template<C T>void PR(T a,int n){rep(i,n){if(i)cout<<' ';cout<<a[i];}ln;}
bool check(int n,int m,int x,int y){return x>=0&&x<n&&y>=0&&y<m;}
const ll MAX=1000000007,MAXL=1LL<<60,dx[4]={-1,0,1,0},dy[4]={0,1,0,-1};
const double EPS = 1e-10;
typedef complex<double> P;

double D(P a, P b) {
  return sqrt((a.real()-b.real())*(a.real()-b.real())+(a.imag()-b.imag())*(a.imag()-b.imag()));
}
double toRad(double agl){return agl*M_PI/180.0;}
P rotate(P a, double r){return P(a.real()*cos(r)-a.imag()*sin(r),a.real()*sin(r)+a.imag()*cos(r));}

void Main() {
  int n;
  while(R n && n) {
    P a[n];
    rep(i,n) {
      double x,y;
      cin >> x >> y;
      a[i]=P(x,y);
    }
    int ans=1;
    rep(i,n) {
      REP(j,i+1,n) {
        if(D(a[i],a[j])>2) continue;
        int cnt=0;
        P q=P((a[i].real()+a[j].real())/2,(a[i].imag()+a[j].imag())/2);
        P p=a[i]-q;
        double r=sqrt(1.0-D(q,a[i])*D(q,a[i]));
        if(r) p*=r/D(q,a[i]);
        p=rotate(p,toRad(90));
        p+=q;
        rep(k,n) if(D(p,a[k])<=1+EPS) cnt++;
        ans=max(ans,cnt);
        cnt=0;
        p-=q;
        p=rotate(p,toRad(180));
        p+=q;
        rep(k,n) if(D(p,a[k])<=1+EPS) cnt++;
        ans=max(ans,cnt);
      }
    }
    pr(ans);
  }
}

int main() {
  ios::sync_with_stdio(0);cin.tie(0);
  Main();return 0;
}