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>
#define N 300
using namespace std;
typedef complex<double> P;
const double EPS = 1e-5;

struct Circle{
  P c;
  double r;
  Circle(P c=P(0,0),double r=0):c(c),r(r){}
};

vector<P> getIntersectCC(Circle c1, Circle c2){
  vector<P> res;
  double r1=c1.r,r2=c2.r,d;
  P p1=c1.c,p2=c2.c;
  d=abs(p1-p2);
  if(d<EPS&&abs(r1-r2)<EPS) return res;
  if(r1+r2<d-EPS||d+EPS<abs(r1-r2)) return res; 
  double a=(r1*r1-r2*r2+d*d)/(2*d);
  double h=sqrt(max(r1*r1-a*a,0.0));
  P tmp1=p1+a/d*(p2-p1);
  P tmp2=h/d*(p2-p1);
  if(abs(tmp2)<EPS)res.push_back(tmp1);
  else{
    res.push_back(P(real(tmp1)-imag(tmp2),imag(tmp1)+real(tmp2)));
    res.push_back(P(real(tmp1)+imag(tmp2),imag(tmp1)-real(tmp2)));
  }
  return res;
}

double getdistancepp(P a,P b){
  return sqrt((a.real()-b.real())*(a.real()-b.real())+(a.imag()-b.imag())*(a.imag()-b.imag()));
}

int main(){
  int n,ans;
  double a,b;
  P p[N];
  while(1){
    ans=0;
    cin>>n;
    if(!n)break;
    for(int i=0;i<n;i++){
      cin>>a>>b;
      p[i]=P(a,b);
    }
    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
	Circle c1,c2;
	c1.c=p[i];
	c2.c=p[j];
	c1.r=c2.r=1;
	vector<P> r=getIntersectCC(c1,c2);
	if(!r.size())continue;
	for(int k=0;k<r.size();k++){
	  int cnt=0;
	  for(int l=0;l<n;l++)
	    if(getdistancepp(r[k],p[l])<=1+EPS)cnt++;	    
	  ans=max(ans,cnt);
	}
      }
    }
    // if(!ans&&n>=1)ans=1;
    cout<<ans<<endl;
  }
  return 0;
}

C++: #include <complex>
#include <iostream>
#include <cmath>
#include <vector>

constexpr double pi = 3.14159265358979323846264338327950288419716939937510582097494459;

int main() {
	for (;;){
		int n;
        std::cin >> n;

		if (n == 0) {
			return 0;
		}

        std::vector<std::complex<double>> arr;
        arr.reserve(n);

		for (int i = 0; i < n; i++) {
			double x, y;

            std::cin >> x >> y;

			arr.emplace_back(x, y);
		}

		int max = 1;
		for(int i = 0; i < n; i++) {
		    for(int j = i + 1; j < n; j++) {
				auto sub = arr[j] - arr[i];

				if (std::abs(sub) >= 2) {
					continue;
				}
				auto sum = arr[j] + arr[i];

				{
					auto center = sum/2. + sub*std::polar(std::sqrt(1-std::abs(sub/2.) * std::abs(sub/2.))/std::abs(sub), pi/2);

					int counter = 0;
					for(int k = 0; k < n; k++) {
						if(i == k || j == k || std::abs(arr[k]-center) <= 1) {
							counter++;
						}
					}
					if (max < counter) {
						max = counter;
					}
				}

				{
					auto center = sum/2. + sub*std::polar(std::sqrt(1-std::abs(sub/2.) * std::abs(sub/2.))/std::abs(sub), -pi/2);

					int counter = 0;
                    for(int k = 0; k < n; k++) {
						if(i == k || j == k || std::abs(arr[k]-center) <= 1) {
							counter++;
						}
					}
					if (max < counter) {
						max = counter;
					}
				}
			}
		}

        std::cout << max << std::endl;
	}
}

C++: #include <iostream>
#include <cmath>
#include <algorithm>
#include <queue>
#include <map>
#define LP(i,n) for(int i = 0; i < n; i++)

using namespace std;
typedef pair<double,double> P;

double Dist(double x1[], double x2[]){
	return pow((x1[0]-x2[0])*(x1[0]-x2[0])+(x1[1]-x2[1])*(x1[1]-x2[1]),0.5);

}

int main(){
	int n;
	cin >> n;
	while(n){
		P p[300];
		LP(i,n){
			cin >> p[i].first >> p[i].second;
		}
		sort(p,p+n);
		double x[300][2];
		LP(i,n){
			x[i][0] = p[i].first;
			x[i][1] = p[i].second;
		}
		double dist[300][300];
		LP(i,n){
			LP(j,n){
				dist[i][j] = Dist(x[i],x[j]);
			}
		}
		int mx = 0;
		LP(i,n){
			for(int j = i+1; j < n; j++){
				if(x[j][0]-x[i][0] > 2){
					break;
				}else{
					if(dist[i][j] < 2){
						double m[] = {(x[j][0]+x[i][0])/2,(x[j][1]+x[i][1])/2};
						double d[] = {x[j][1]-x[i][1],-x[j][0]+x[i][0]};
						double O[] = {0,0};
						double _d = Dist(d,O);
						LP(k,2){
							d[k] *= (pow(_d,-2)-0.25);
						}
						double o[2][2];
						LP(k,2){
							LP(l,2){
								o[k][l] = m[l];
								if(k) o[k][l] += d[l];
								else o[k][l] -= d[l];
							}
						}
						LP(k,2){
							int s = 2;
							for(int l = i+1; l < j; l++){
								if(Dist(x[l],o[k]) < 1) s++;
							}
							for(int l = i-1; 0<=l && x[i][0]-x[l][0] < 1;l--){
								if(Dist(x[l],o[k]) < 1) s++;
							}
							for(int l = j+1; l < n && x[l][0]-x[j][0] < 1; l++){
								if(Dist(x[l],o[k]) < 1) s++;
							}
							mx = max(mx,s);
						}
					}
				}
			}
		}
		cout << mx << endl;
		cin >> n;
	}
}

C++: #include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <stack>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <queue>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <string>
#include <cstring>
#include <ctime>
#include <fstream>
#include <stdio.h>
#include <complex>

using namespace std;

//conversion
//------------------------------------------
inline int toInt(string s) {int v; istringstream sin(s);sin>>v;return v;}
template<class T> inline string toString(T x) {ostringstream sout;sout<<x;return sout.str();}

//typedef
//------------------------------------------
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<string> VS;
typedef pair<int, int> PII;
typedef pair<int, PII> TIII;
typedef long long LL;
typedef vector<LL> VLL;

//container util

//------------------------------------------
#define ALL(a)  (a).begin(),(a).end()
#define RALL(a) (a).rbegin(), (a).rend()
#define PB push_back
#define MP make_pair
#define SZ(a) int((a).size())
#define EACH(i,c) for(typeof((c).begin()) i=(c).begin(); i!=(c).end(); ++i)
#define EXIST(s,e) ((s).find(e)!=(s).end())
#define SORT(c) sort((c).begin(),(c).end())
#define MT(a,b,c) MP(a, MP(b, c))

//repetition
//------------------------------------------
#define FOR(i,a,b) for(int i=(a);i<(b);++i)
#define REP(i,n)  FOR(i,0,n)

double x[300], y[300];
int n;

int count(double x0, double y0, int p, int q){
	int ret = 2;
	REP(i,n){
		if(i==p||i==q)continue;
		double dx = x0-x[i];
		double dy = y0-y[i];
		if(pow(dx, 2.0)+pow(dy, 2.0)<1.0)ret++;
	}
	return ret;
}

int main(){
	while(1){
		int ret = 1;
		cin>>n;
		if(n==0)break;
		REP(i,n)scanf("%f %f", &x[i],&y[i]);
		REP(i,n){
			FOR(j, i+1,n){
				double dx = x[j]-x[i];
				double dy = y[j]-y[i];
				double d = sqrt(pow(dx, 2.0)+pow(dy, 2.0));
				if(d < 2.0){
					double mx = (x[j]+x[i])*0.5;
					double my = (y[j]+y[i])*0.5;
					{
						double x1 = mx-dy/d*sqrt(1.0-d*d/4.0);
						double y1 = my+dx/d*sqrt(1.0-d*d/4.0);
						ret = max(ret, count(x1, y1, i, j));
					}
					{
						double x1 = mx+dy/d*sqrt(1.0-d*d/4.0);
						double y1 = my-dx/d*sqrt(1.0-d*d/4.0);
						ret = max(ret, count(x1, y1, i, j));
					}
				}
			}
		}
		cout << ret << endl;
	}
	return 0;
}