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;

#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) continue;
                ld x1 = x[i], y1 = y[i], x0 = x[j], y0 = y[j];
                for (auto dir : {1, -1}) {
                    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));
//                    cout << x1 << " " << y1 << " " << x0 << " " << y0 << " : " << xc << " " << yc << endl;
                    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) 
                            m++;
                    }
                    chmax(ret, m);
                }
            }
        }
        cout << ret << endl;
    }
    return 0;
}

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

typedef complex<double> P;
#define EPS (1e-10)

int main()
{
	int n;
	while(cin>>n && n!=0){
		int ans=1;
		vector<P> points;
		for(int i=0;i<n;i++){
			double x,y; cin>>x>>y;
			points.push_back(P(x,y));
		}
		for(int a=0;a<n;a++){
			for(int b=a+1;b<n;b++){
				if(abs(points[b]-points[a])+EPS>2.0) continue;
				P sa=P(points[b]-points[a])/2.0;
				P rev[2];
				rev[0]=sa*P(0,1.0)/abs(sa);
				rev[1]=sa*P(0,-1.0)/abs(sa);
				double l=sqrt(1-abs(sa)*abs(sa));

				for(int k=0;k<2;k++){
					P center=points[a]+sa+rev[k]*l;
					int t_ans=2;
					for(int id=0;id<n;id++){
						if(id==a || id==b) continue;
						if(abs(center-points[id])+EPS<1.0){
							t_ans++;
//							cout<<abs(center-points[id])<<endl;
						}
					}
					
//					cout<<t_ans<<endl;
					ans=max(ans,t_ans);
				}
			}
		}

		cout<<ans<<endl;
	}
	return 0;
}

C++: #include <bits/stdc++.h>
#define N 301
using namespace std;
typedef complex <double> P;
const double EPS = 1e-8;
P p[N];
int n;

int cnt(P c){
  int res=0;
  for(int i=0;i<n;i++)if(abs(c-p[i])<1+EPS) res++;
  return res;
}

P get_C(P a,P b){
  double x=abs(b-a);
  if(x>2||a==b) return a;
  return P(x/2,sqrt(1-x*x/4))*(b-a)/abs(b-a)+a;  
}


int main(){
  while(1){
  cin>>n;
  if(!n)break;
  double a,b;
  for(int i=0;i<n;i++)cin>>a>>b,p[i]=P(a,b);
  
  int ans=0;
  for(int i=0;i<n;i++)
    for(int j=0;j<n;j++)
      ans=max(ans,cnt(get_C(p[i],p[j])));
  cout <<ans<<endl;  
  }
  return 0;
}

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

using namespace std;

typedef complex<double> P;

#define EPS 1e-12

int main() {
	int N;
	while (cin >> N, N) {
		if (N == 1) {
			cout << 1 << endl;
			continue;
		}

		vector<P> v;
		for (int i = 0; i < N; ++i) {
			double x, y; cin >> x >> y;
			v.push_back( P(x, y) );
		}

		vector<P> c;
		for (int i = 0; i < N; ++i) {
			for (int j = i+1; j < N; ++j) {
				P r = (v[j] - v[i])/2.0;
				double l = abs(r);
				if (l > 1 + EPS) {
					continue;
				}
				else if (l >= 1 - EPS) {
					c.push_back( v[i] + r );
				}
				else {
					P t = r / l * sqrt(1 - l*l);
					c.push_back(v[i] + r + t * P(0, 1));
					c.push_back(v[i] + r + t * P(0, -1));
				}
			}
		}

		int ans = 0;
		for (int i = 0; i < c.size(); ++i) {
			int cnt = 0;
			for (int j = 0; j < v.size(); ++j) {
				if ( abs(c[i] - v[j]) <= 1 + EPS ) {
					++cnt;
				}
			}
			ans = max(ans, cnt);
		}
		cout << ans << endl;
	}
}