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 <queue>
#include <vector>
#include <algorithm>
#include <cmath>
#include <complex>
using namespace std;
#define rep(i,n) for(int i=0;i<n;i++)
static const double EPS = 1e-8;
static const double PI = 3.14159265;

typedef complex<double> P;
struct C {P p; double r;C(const P &p, double r) : p(p), r(r) { }};
#define toRad(n) (n*(PI/180.0))
vector<P> data;
pair<C,C> get(P a,P b){
	P diff = a-b;
	if(abs(diff) > 2)return make_pair(C(P(-1,-1),-1),C(P(-1,-1),-1));
	double theta = atan2(diff.imag(),diff.real()) + toRad(90);
	P mid = (a+b) / 2.0;
	double d = sqrt( 1 - (abs(diff)*abs(diff))/4.0 );
	double X = cos(theta)*d , Y = sin(theta)*d;
	return make_pair(C(mid+P(X,Y),1),C(mid-P(X,Y),1));
}
int check(C t){
	int ret = 0;
	rep(i,data.size()){
		if(abs(data[i]-t.p) < 1+EPS)ret++;
	}
	return ret;
}
int main(){
	int n; 
	while(cin >> n , n){
		data.resize(n);
		rep(i,n){
			cin >> data[i].real() >> data[i].imag();
		}
		int ret = 0;
		rep(i,n){
			for(int j=i+1;j<n;j++){
				pair<C,C> g = get(data[i],data[j]);
				ret = max(ret,check(g.first));
				ret = max(ret,check(g.second));
			}
		}
		cout << ret << endl;
	}
}

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

using namespace std;

int ans;

complex<double> MaxCoverCircle(vector<complex<double>> &p, double radius) {
  const int n = p.size();
  const double eps = 1e-6;
  int num_points = 1;
  complex<double> center_max(p[0]);

  for (int i = 0; i < n; ++i) {
    for (int j = i + 1; j < n; ++j) {
      auto diff = p[i] - p[j];
      double dist = abs(diff);

      if (dist - 2*radius > eps) continue;

      auto mid = p[j] + (diff/2.0);
      auto mid_unit_normal = (diff*complex<double>(0, 1))/dist;
      auto mid_to_center = sqrt(radius - dist*dist/4.0)*mid_unit_normal;

      for (int k = 0; k < 2; ++k) {
        mid_to_center *= -1;
        auto center = mid + mid_to_center;

        int cnt = 0;
        for (int l = 0; l < n; ++l) {
          auto q = p[l] - center;
          if (q.real()*q.real() + q.imag()*q.imag() - radius*radius < eps) cnt++;
        }

        if (cnt > num_points) {
          num_points = cnt;
          center_max = center;
        }
      }
    }
  }
  ans = num_points;
  return center_max;
}

int main() {
  int n;
  while (cin >> n, n) {
    vector<complex<double>> p;
    p.reserve(n);
    for (int i = 0; i < n; ++i) {
      double x, y;
      cin >> x >> y;
      p.emplace_back(x, y);
    }
    MaxCoverCircle(p, 1);
    cout << ans << endl;
  }
}

C++: //include
//------------------------------------------
#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 <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <string>
#include <cstring>
#include <ctime>
#include <climits>
#include <queue>
#include <complex>

using namespace std;

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

//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())

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

//constant
//--------------------------------------------
const double EPS = 1e-6;
const double PI  = acos(-1.0);


typedef complex<double> CP;
int N;
CP point[310];

int main(){
  cin.tie(0);
  ios_base::sync_with_stdio(false);

  while(true){
	cin >> N; if(N == 0) break;
	REP(i,N){
	  double x,y; cin >> x >> y;
	  point[i] = CP(x, y);
	}

	int ans = 1;
	for(int i1=0;i1<N;++i1){
	  for(int i2=i1+1;i2<N;++i2){
		double d = norm(point[i1] - point[i2]);
		if(sqrt(d) > 2.00){
		  continue;
		}
		CP m = (point[i1] + point[i2]) / 2.0;
		CP n = point[i2] - point[i1];
		n /= abs(n);
		n = CP(-imag(n), real(n));

		CP c = m + n * sqrt(1 - d / 4.0);
		int res = 2;
		for(int i=0;i<N;++i){
		  if(i==i1 || i==i2) continue;
		  if(abs(c-point[i]) <= 1.0)
			++res;
		}
		ans = max(ans, res);

		c = m - n * sqrt(1 - d / 4.0);
		res = 2;
		for(int i=0;i<N;++i){
		  if(i==i1 || i==i2) continue;
		  if(abs(c-point[i]) <= 1.0)
			++res;
		}
		ans = max(ans, res);
		
	  }
	}
	
	cout << ans << endl;
  }

  return 0;
}

C++: #include <cstdio>
#include <complex>
#include <algorithm>
#include <utility>
#include <iostream>
#define rep(i,n) for(int i=0;i<n;i++)
#define re real()
#define im imag()
using namespace std;
int main(){
	while(true){
		int n,mx=1;
		scanf("%d",&n);
		if(n==0) break;
		complex<double> p[300];
		rep(i,n) scanf("%lf%lf",&p[i].re,&p[i].im);
		rep(i,n) rep(j,i){
			double d=abs(p[i]-p[j]);
			if(d>2) continue;
			complex<double> c=(p[i]+p[j])/2.0,dis=(p[i]-p[j])/d*sqrt(1-d*d/4);
			swap(dis.re,dis.im);
			dis.im=-dis.im;
			complex<double> cen=c+dis;
			int cnt=0;
			rep(k,n) if(abs(cen-p[k])<=1.0000000001) cnt++;
			mx=max(mx,cnt);
			cen=c-dis;
			cnt=0;
			rep(k,n) if(abs(cen-p[k])<=1.0000000001) cnt++;
			mx=max(mx,cnt);
		}
		printf("%d\n",mx);
	}
	return 0;
}