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-3;

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=0;
  double a,b;
  P p[N];
  while(1){
    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);
	}
      }
    }
    cout<<ans<<endl;
  }
  return 0;
}

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

#define EPS (1e-10)

using namespace std;

int count( vector<complex<double>> &vpos );

int main(int argc, char const *argv[]) {
    int n, result;
    double x, y;
    vector<complex<double>> vpos;
    vector<int> results;

    while(1){
        cin>>n;
        if( n == 0 ) break;
        for(int i = 0; i < n; i++){
            cin>>x>>y;
            vpos.push_back( complex<double> (x,y) );
        }
        results.push_back( count( vpos ) );
        vpos.clear();
    }

    for (auto it = results.begin(); it != results.end(); it++) {
        cout<<*it<<endl;
    }

    return 0;
}

int count( vector<complex<double>> &vpos ){
    int count = 0;
    int tmp_count1, tmp_count2;
    complex<double> m, n1, n2, center1, center2;

    for(auto fi = vpos.begin(); fi != vpos.end(); fi++){
        for(auto sc = vpos.begin(); sc != vpos.end(); sc++){
            if( (*fi).real() == (*sc).real() ) continue;
            tmp_count1 = 0;
            tmp_count2 = 0;

            m = ( (*fi) + (*sc) ) / 2.0;
            n1 = ( (*fi) -  (*sc) ) * complex<double> ( 0.0, 1.0 );
            n1 = n1 / abs(n1);
            n2 = ( (*fi) -  (*sc) ) * complex<double> ( 0.0, -1.0 );
            n2 = n2 / abs(n2);

            double x = abs( (*fi) -  (*sc) ) / 2.0;
            center1 =  m + n1 * sqrt( 1.0 - x * x );
            center2 =  m + n2 * sqrt( 1.0 - x * x );

            for(auto it = vpos.begin(); it != vpos.end(); it++){
                double dist1 = abs( center1 - (*it) );
                double dist2 = abs( center2 - (*it) );
                if( dist1 <= 1.0 + EPS ){
                    tmp_count1++;
                }
                if( dist2 <= 1.0 + EPS ){
                    tmp_count2++;
                }
            }
            if( count < max(tmp_count1, tmp_count2) ){
                count = max(tmp_count1, tmp_count2);
            }
        }
    }
    return count;
}

C++: #include <bits/stdc++.h>
using namespace std;

#define rep(i,n) for(long long i = 0; i < (long long)(n); i++)
template<class T1, class T2> bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }

using ll = long long; using ld = long double; using vld = vector<ld>; 

int main(void) {
    ll n; 
    while (cin >> n && n) {
        vld x(n), y(n); rep(i, n) cin >> x[i] >> y[i];
        auto d = [&](ll i, ll j) { return (x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j]); };
        ll ret = 1;
        for (int i = 0; i < n-1; i++) {
            for (int j = i+1; j < n; j++) {
                ld dij = d(i, j);
                if (dij >= 4.0l) continue;
                ld x1 = x[i], y1 = y[i], x0 = x[j], y0 = y[j];
                for (auto dir : {1.0l, -1.0l}) {
                    ld xc = x0 + (x1 - x0) / 2.0l + (-(y1 - y0)) * dir 
                        * sqrtl(1.0l - ((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) / 4.0l) 
                        / sqrtl((y1 - y0) * (y1 - y0)  + (x1 - x0) * (x1 - x0));
                    ld yc = y0 + (y1 - y0) / 2.0l + (+(x1 - x0)) * dir 
                        * sqrtl(1.0l - ((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) / 4.0l) 
                        / sqrtl((y1 - y0) * (y1 - y0)  + (x1 - x0) * (x1 - x0));
                    ll m = 2;
                    rep(k, n) {
                        if (k == i || k == j) continue;
                        if ((xc - x[k]) * (xc - x[k]) + (yc - y[k]) * (yc - y[k]) < 1.0l) 
                            m++;
                    }
                    chmax(ret, m);
                }
            }
        }
        cout << ret << endl;
    }
    return 0;
}