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 <complex>
#include <cmath>
#define EPS 10e-6
using namespace std;
typedef complex<double> xy_t;
int n;
xy_t points[301];
int count_in(xy_t center)
{
  int ans=0;
  for(int i=0;i<n;i++)
  {
    if (abs(center-points[i])<1+EPS)
      ans++;
  }
  return ans;
}
xy_t ppdcl(xy_t a)
{
  return xy_t(a.imag(),-a.real());
}
xy_t delta_v(xy_t C1,xy_t C2)
{
  xy_t half=(C1-C2)*0.5;
  double l=abs(half);
  return ppdcl(half)*sqrt(1-l*l)*(1/l);
}
int main()
{
  while(1)
  {
    cin>>n;
    if (n==0) break;
    double x,y;
    for(int i=0;i<n;i++)
    {
      cin>>x>>y;
      points[i]=xy_t(x,y);
    }
    int maxn=-1;
    for(int i=0;i<n;i++)
    {
      for(int j=i+1;j<n;j++)
      {
        if(abs(points[i]-points[j])<4+EPS)
        {
          xy_t mid_p=(points[i]+points[j])*0.5;
          maxn=max(maxn,count_in(mid_p+delta_v(points[i],points[j])));
          maxn=max(maxn,count_in(mid_p-delta_v(points[i],points[j])));
        }
      }
    }
    cout<<maxn<<endl;
  }

}

C++: #include <iostream>
#include <complex>
#include <vector>
#include <utility>
using namespace std;

using Float = double;
using Point = complex<Float>;

vector<Point> ps;

const Float radius = 1.0;
const Float EPS = 1e-8;

int calc(Point c){
	int ret = 0;
	for(const Point &p : ps){
		if(abs(p - c) <= radius + EPS) ret++;
	}
	return ret;
}

signed main(){
	int n;
	while(cin >> n,n){
		ps.clear();
		for(int i = 0;i < n;i++){
			Float x,y;
			cin >> x >> y;
			ps.emplace_back(x,y);
		}
		int ans = 1;
		for(int i = 0;i < n;i++){
			for(int j = i + 1;j < n;j++){
				if(abs(ps[i] - ps[j]) >= radius * 2) continue;
				Point mid = (ps[i] + ps[j]) / 2.0;
				Point vect = Point(imag(ps[i] - ps[j]),-real(ps[i] - ps[j]));
				Float dist = sqrt(1.0 - norm(mid - ps[i]));
				ans = max(ans,max(calc(mid + vect / abs(vect) * dist),calc(mid - vect / abs(vect) * dist)));
			}
		}
		cout << ans << endl;
	}
}

C++: #include <iostream>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <string>
#include <vector>
#include <utility>
#include <complex>
#include <set>
#include <map>
#include <queue>
#include <stack>
#include <deque>
#include <tuple>
#include <bitset>
#include <algorithm>
using namespace std;
typedef long double ld;
typedef long long ll;
typedef vector<int> vint;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
typedef pair<double, double> pdd;
typedef complex<double> comd;
#define rep(i,n)	for(int i=0;i<n;i++)
#define srep(i,a,n)	for(int i=a;i<n;i++)
#define REP(i,n)	for(int i=0;i<=n;i++)
#define SREP(i,a,n)	for(int i=a;i<=n;i++)
#define rrep(i,n)	for(int i=n-1;i>=0;i--)
#define RREP(i,n)	for(int i=n;i>=0;i--)
#define all(a)	(a).begin(),(a).end()
#define mp(a,b)	make_pair(a,b)
#define mt	make_tuple
#define fst	first
#define scn second
#define bucnt(x)	__buildin__popcount(x)
#define debug(x)	cout<<"debug: "<<x<<endl

const ll inf = (ll)1e9;
const ll mod = (ll)1e9 + 7;
const ld eps = 1e-9;
const int dx[] = { 0,1,0,-1 };
const int dy[] = { 1,0,-1,0 };

double Dist2(pdd a, pdd b) {
	return (a.fst - b.fst)*(a.fst - b.fst) + (a.scn - b.scn)*(a.scn - b.scn);
}

vector<pdd> getCenter(pdd a, pdd b) {
	vector<pdd> ret;
	pdd c = mp((a.fst + b.fst) / 2.0, (a.scn + b.scn) / 2.0);
	pdd v = mp(a.scn - b.scn, b.fst - a.fst);
	double D = sqrt(v.fst*v.fst + v.scn*v.scn);
	v.fst /= D;	v.scn /= D;
	D = sqrt(1 - Dist2(a, b)/4.0);
	v.fst *= D;	v.scn *= D;
	ret.push_back(mp(c.fst + v.fst, c.scn + v.scn));
	ret.push_back(mp(c.fst - v.fst, c.scn - v.scn));
	return ret;
}

int main() {
	while (true) {
		int n;	cin >> n;
		if (n == 0)	break;
		vector<pdd> p(n,mp(0,0));
		rep(i, n) cin >> p[i].fst >> p[i].scn;
		int ret = 1;
		rep(i, n) {
			rep(j, i) {
				if (Dist2(p[i], p[j]) <= 4.0) {
					vector<pdd> c = getCenter(p[i], p[j]);
					int tmp = 0;
					rep(k, n) {
						if (Dist2(c[0], p[k]) <= 1.0 + eps)	tmp++;
					}
					ret = max(ret, tmp);
					tmp = 0;
					rep(k, n) {
						if (Dist2(c[1], p[k]) <= 1.0 + eps)	tmp++;
					}
					ret = max(ret, tmp);
				}
			}
		}
		cout << ret << endl;
	}
	return 0;
}

C++: #include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;

//eÍÍÌë·
//1 * 10^(-5)
#define EPS (1e-8)

int n;
pair<double, double> p[301];

int main() {
	while (1) {
		cin >> n;
		if (n == 0) break;
		for (int i = 0; i < n; i++) {
			cin >> p[i].first >> p[i].second;
		}
		int maxPoints = 1;
		//2_ð~üãÉÂ~ðè
		for (int i = 0; i < n; i++) {
			for (int j = i + 1; j < n; j++) {
				double d = hypot(p[i].first - p[j].first, p[i].second - p[j].second);
				if (d > 2.0 + EPS) continue;	//£ê·¬
				
				//_(»ê¼êÌÀWÌ½Ï)
				pair<double, double> m((p[i].first + p[j].first) / 2.0, (p[i].second + p[j].second) / 2.0);
				//PÊ@üxNg
				//	@üxNgÌX«Íüªp[i]p[j]Æ|¯Ä-1C·Èí¿
				//		v = (p[j].second - p[i].second, p[j].first - p[i].first)
				//	±êÌ·³Ívx^2 + vy^2 = d^2
				//	æÁÄPÊxNgÉ·éÉÍ
				//		(vx/d)^2 + (vy/d)^2 = 1
				pair<double, double> ve((p[j].second - p[i].second) / d, (p[j].first - p[i].first) / d);
				//_©çSÜÅÌ£
				double l = sqrt(1.0 - d * d / 4.0);
				
				//_2Â©çèÜé~Í2Â
				for (int s = -1; s < 2; s += 2) {//
					pair<double, double> c(m.first + s * l * ve.first, m.second - s * l * ve.second);
					int count = 2;
					for (int k = 0; k < n; k++) {
						if (k==i || k==j) continue;//¯êÌ_ÍÜßÈ¢
						if (hypot(c.first - p[k].first, c.second - p[k].second) < 1.0 + EPS) count++;
					}
					maxPoints = max(maxPoints, count);
				}
			}
		}
		cout << maxPoints << endl;
	}
	return 0;
}