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;

double distance(pair<double,double> p1,pair<double,double> p2){
	return sqrt((p1.first-p2.first)*(p1.first-p2.first)+(p1.second-p2.second)*(p1.second-p2.second));
}

int main(){
	int n;
	while(cin>>n,n){
		vector<double> x(n);
		vector<double> y(n);
		int res=1;
		for(int i=0;i<n;i++) cin>>x[i]>>y[i];
		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				pair<double,double> midp=make_pair((x[i]+x[j])/2,(y[i]+y[j])/2);
				double a=-(x[i]-x[j]),b=(y[i]-y[j]);
				if(a<0.000001) a=0;
				if(b<0.000001) b=0;
				//a(x-mid[i].first)+b(y-mid.second)=0
					for(int k=0;k<2;k++){
						pair<double,double> lef=midp;
						pair<double,double> rig;
						if(k==0) rig=make_pair(-b*1e7+midp.first,a*1e7+midp.second);
						if(k==1) rig=make_pair(b*1e6+midp.first,-a*1e6+midp.second);
						while(distance(lef,rig)>=0.00000001){
								pair<double,double> mid=make_pair((lef.first+rig.first)/2,(lef.second+rig.second)/2);
								if(distance(mid,make_pair(x[i],y[i]))>1){
									rig=mid;
								}
								else{
									lef=mid;
								}

						}
						int count=0;
						for(int r=0;r<n;r++){
							if(distance(make_pair(x[r],y[r]),lef)<1+0.000001){
								count++;
							}
						}
						res=max(res,count);
					}
			}

		}
		cout<<res<<endl;

	}

	return 0;
}

C++: #include <algorithm>
#include <iostream>
#include <cstdio>
#include <complex>
#include <vector>
#define X real()
#define Y imag()
#define EPS (1e-10)
using namespace std;
typedef complex<double> P;
typedef pair<P, double> C;

namespace std{
  bool operator<(const P& a, const P& b){
    return a.X != b.X ? a.X < b.X : a.Y < b.Y;
  }
}

P n_vector(P a){
  double l = abs(a);
  return P(-a.Y,a.X)/l;
}

vector<P> cpcc(C a, C b){
  double d = abs(a.first-b.first), ra = a.second, rb = b.second;
  double t = (ra*ra-rb*rb+d*d)/2/d, h = sqrt(ra*ra-t*t);
  P m = t/abs(b.first-a.first)*(b.first-a.first)+a.first;
  vector<P> cp(2,m);
  P n = n_vector(a.first-b.first);
  cp[0] -= h*n;
  cp[1] += h*n;
  return cp;
}


int main(){
  int N;
  double x, y;
  while(cin >> N, N){
    P ps[N];
    for(int i = 0; i < N; ++i){
      cin >> x >> y;
      ps[i] = P(x,y);
    }
    vector<P> cps;
    for(int i = 0; i < N; ++i){
      for(int j = i+1; j < N; ++j){
	if(abs(ps[i]-ps[j]) <= 2 + EPS){
	  vector<P> cp = cpcc(C(ps[i],1),C(ps[j],1));
	  for(int k = 0; k < cp.size(); ++k)
	    cps.push_back(cp[k]);
	}
      }
    }
    int n = cps.size(), ans = 1;
    for(int i = 0; i < n; ++i){
      int cnt = 0;
      for(int j = 0; j < N; ++j){
	if(norm(cps[i]-ps[j]) <= 1 + EPS)
	  ++cnt;
      }
      ans = max(ans, cnt);
    }
    cout << ans << endl;
  }
  return 0;
}

C++: #include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef complex<double> P;
typedef pair<int,int> pii;
#define REP(i,n) for(ll i=0;i<n;++i)
#define REPR(i,n) for(ll i=1;i<n;++i)
#define FOR(i,a,b) for(ll i=a;i<b;++i)

#define DEBUG(x) cout<<#x<<": "<<x<<endl
#define DEBUG_VEC(v) cout<<#v<<":";REP(i,v.size())cout<<" "<<v[i];cout<<endl
#define ALL(a) (a).begin(),(a).end()

#define MOD (ll)(1e9+7)
#define ADD(a,b) a=((a)+(b))%MOD
#define FIX(a) ((a)%MOD+MOD)%MOD

int n;
P pts[353];

int check(P c){
  int ret = 0;
  REP(i,n){
    if(abs(c-pts[i])<=1.00001)ret++;
  }
  return ret;
}

int main(){
  while(true){
    scanf("%d",&n);
    if(n==0)break;
    REP(i,n){
      double x,y;
      scanf("%lf%lf",&x,&y);
      pts[i] = P(x,y);
    }
    int ans = 1;
    REP(i,n)REP(j,i){
      P a,b;
      a = pts[i];
      b = pts[j];
      double l = abs(b-a);
      if(l>2.0)continue;
      double h = sqrt(1.0-l*l/4.0);
      P c;
      c = a + (b-a)/2.0 + (b-a)*P(0,1)*h/l;
      ans = max(ans,check(c));
      c = a + (b-a)/2.0 + (b-a)*P(0,-1)*h/l;
      ans = max(ans,check(c));
    }
    printf("%d\n",ans);
  }
  return 0;
}

C++: #include "bits/stdc++.h"

#define REP(i,n) for(ll i=0;i<n;++i)
#define RREP(i,n) for(ll i=n-1;i>=0;--i)
#define FOR(i,m,n) for(ll i=m;i<n;++i)
#define RFOR(i,m,n) for(ll i=n-1;i>=m;--i)
#define ALL(v) (v).begin(),(v).end()
#define UNIQUE(v) v.erase(unique(ALL(v)),v.end());
#define DUMP(v) REP(aa, (v).size()) { cout << v[aa]; if (aa != v.size() - 1)cout << " "; else cout << endl; }
#define INF 1000000001ll
#define MOD 1000000007ll
#define EPS 1e-9

const int dx[8] = { 1,1,0,-1,-1,-1,0,1 };
const int dy[8] = { 0,1,1,1,0,-1,-1,-1 };


using namespace std;
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef vector<vi> vvi;
typedef vector<vl> vvl;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
ll max(ll a, int b) { return max(a, ll(b)); }
ll max(int a, ll b) { return max(ll(a), b); }
ll min(ll a, int b) { return min(a, ll(b)); }
ll min(int a, ll b) { return min(ll(a), b); }
///(?´????????`)(?´????????`)(?´????????`)(?´????????`)(?´????????`)(?´????????`)///

int main() {
	cin.tie(0);
	ios::sync_with_stdio(false);
	int n;
	while (cin >> n, n) {
		vector<double> x(n), y(n);
		int ans = -1;
		REP(i, n)cin >> x[i] >> y[i];
		REP(i, n) {
			REP(j, n) {
				if (i == j)continue;
				double ax = x[j] - x[i], ay = y[j] - y[i];
				if (ax*ax + ay*ay+EPS > 4)continue;
				double px = x[i] + ax / 2 + sqrt((4 - (ax*ax + ay*ay)) / (4 * (ax*ax + ay*ay)))*ay,
				py = y[i] + ay / 2 - sqrt((4 - (ax*ax + ay*ay)) / (4 * (ax*ax + ay*ay)))*ax;
				int cnt = 0;
				REP(k, n) {
					if ((x[k] - px)*(x[k] - px) + (y[k] - py)*(y[k] - py) < 1+EPS)cnt++;
				}
				ans = max(ans, cnt);
			}
		}
		REP(i, n) {
			REP(j, n) {
				if (i == j)continue;
				double ax = x[j] - x[i], ay = y[j] - y[i];
				if (ax*ax + ay*ay+EPS > 4)continue;
				double px = x[i] + ax / 2 - sqrt((4 - (ax*ax + ay*ay)) / (4 * (ax*ax + ay*ay)))*ay,
					py = y[i] + ay / 2 + sqrt((4 - (ax*ax + ay*ay)) / (4 * (ax*ax + ay*ay)))*ax;
				int cnt = 0;
				REP(k, n) {
					if ((x[k] - px)*(x[k] - px) + (y[k] - py)*(y[k] - py) < 1+EPS)cnt++;
				}
			
				ans = max(ans, cnt);
			}
		}
		cout <<ans << endl;
	}
}