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<(int)(N);i++)
#define REP(i,k,N) for(int i=k;i<(int)(N);i++)

#define vi vector<int>
#define pb push_back

#define pdd pair<double,double>
#define fi first
#define se second
const double  EPS = 1e-6;

double dist(double x,double y){
    return x*x+y*y;
};

int main(){
    ifstream in("d_in.txt");
    cin.rdbuf(in.rdbuf());
    int n;
    while(cin >> n,n){
        int ma = 0;
        pdd p[330];
        rep(i,n) cin >> p[i].fi >> p[i].se;
        rep(i,n-1){
            REP(j,i+1,n){

                double x1,y1,x2,y2,x3,y3;
                tie(x1,y1) = p[i];
                tie(x2,y2) = p[j];
                double cx1,cy1,cx2,cy2;
                if(dist(x1-x2,y1-y2)  > 4+EPS  ) continue;

                //はりまコーナー
//                double ax,ay,bx1,by1,bx2,by2,mx,my,dx,dy;
//                ax=x1-x2;ay=y1-y2;
//                mx=(x1+x2)/2;my=(y1+y2)/2;
//                by1 = (-0.5)*pow(ax,4)+pow(ax,2)+(-0.5)*pow(ax,2)*pow(ay,2)
//                        /dist(ax,ay);
//                by1 = sqrt(abs(by1));
//                by2 = -1*by1;
//                cy1= my+by1;
//                cy2= my+by2;
//                bx1= -1*(ay*by1)/ax;
//                bx2= -1*(ay*by2)/ax;
//                cx1=bx1+mx;
//                cx2=bx2+mx;


                //あじのりコーナー
                double mx=(x1+x2)/2,my=(y1+y2)/2;

                double ld = sqrt(1-dist((x1-x2),(y1-y2))/4);
                double theta = atan((y2-y1)/(x2-x1));
                cx1 = mx - sin(theta)*ld;
                cy1 = my + cos(theta)*ld;
                cx2 = mx + sin(theta)*ld;
                cy2 = my - cos(theta)*ld;

                int cnt1 = 0,cnt2=0;
                for(int k = 0;k<n;k++) {
                    tie(x3, y3) = p[k];
                    if (EPS + 1> dist((cx1 - x3), (cy1 - y3))) cnt1++;
                    if (EPS + 1> dist((cx2 - x3), (cy2 - y3))) cnt2++;
                }
                ma = max({cnt1,cnt2,ma});
            }
        }
        cout<< ma << endl;
    }
    return 0;
}

C++: #define M_PI 3.14159265359
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace std;
class P{
public:
	double x,y;
	P(){}
	P(double xx,double yy){
		x=xx;
		y=yy;
	}
	P operator +(P p){
		return P(x+p.x,y+p.y);
	}
	P operator -(P p){
		return P(x-p.x,y-p.y);
	}
	P operator /(double d){
		return P(x/d,y/d);
	}
};

int n;
P p[301];

int main(void){
	while(1){
		scanf("%d",&n);
		if(n==0)break;
		for(int i=0;i<n;i++){
			scanf("%lf%lf",&p[i].x,&p[i].y);
		}
		int res=1;
		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				P np=p[j]-p[i];
				//printf("%lf %lf\n",np.x,np.y);
				double blen=np.x*np.x/4.0+np.y*np.y/4.0;
				if(blen>1.0)continue;
				double len=sqrt(1-blen);
				double dir=atan2(np.y,np.x);
				P cp=(p[i]+p[j])/2+P(len*cos(dir+M_PI/2),len*sin(dir+M_PI/2));
				int cnt=0;
				for(int k=0;k<n;k++){
					double len2=sqrt((cp.x-p[k].x)*(cp.x-p[k].x)+(cp.y-p[k].y)*(cp.y-p[k].y));
					if(len2<=1.0001)cnt++;
				}
				res=max(res,cnt);
				cp=(p[i]+p[j])/2+P(len*cos(dir-M_PI/2),len*sin(dir-M_PI/2));
				cnt=0;
				for(int k=0;k<n;k++){
					double len2=sqrt((cp.x-p[k].x)*(cp.x-p[k].x)+(cp.y-p[k].y)*(cp.y-p[k].y));
					if(len2<=1.0001)cnt++;
				}
				res=max(res,cnt);
			}
		}
		printf("%d\n",res);
	}
	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 <iostream>
#include <complex>
#include <set>
#include <vector>
#include <algorithm>
#include <cmath>

#define EPS (1e-10)
#define EQ(a,b) (abs((a)-(b)) < EPS)
#define EQV(a,b) ( EQ((a).real(), (b).real()) && EQ((a).imag(), (b).imag()) )

using namespace std;

typedef complex<double> P;

P center1(P a ,P b ,double r){
  return (a+b)/(P)2+(b-a)*P(0,1)*(P)(sqrt(r*r-norm((b-a)/(P)2))/abs(b-a));
}

P center2(P a ,P b ,double r){
  return (a+b)/(P)2+(b-a)*P(0,-1)*(P)(sqrt(r*r-norm((b-a)/(P)2))/abs(b-a));
}

int main(){
  int n;
  while(true){
    cin >>n;
    if (!n)
      break;
    P dot[n];
    double x,y;
    for(int i=0;i<n;i++){
      cin >>x >>y;
      dot[i] = P(x,y);
    }
    set<int> near_dot[n];
    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
	if (norm(dot[i]-dot[j])<4)
	  near_dot[i].insert(j);
      }
    }
    
    set<int>::iterator it;
    int cnt,max=2;
    P c;
    for(int i=0;i<n;i++){
      it = near_dot[i].begin();
      while(it != near_dot[i].end()){
	c = center1(dot[i],dot[*it],1);
	cnt=2;
	for(int j=0;j<n;j++)
	  if (norm(dot[j]-c)<=1 && j!=i && j!=*it)
	    cnt++;
	if (max<cnt){
	  max=cnt;
	}

	c = center2(dot[i],dot[*it],1);
	cnt=2;
	for(int j=0;j<n;j++)
	  if (norm(dot[j]-c)<=1 && j!=i && j!=*it)
	    cnt++;
	if (max<cnt){
	  max=cnt;
	}
	it++;
      }
    }
    cout << max << endl;
  }
  return 0;
}