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<cmath>

using namespace std;

int n;
double x[300];
double y[300];

double dist(double x1, double y1, double x2, double y2) {
    double dx = abs(x1 - x2);
    double dy = abs(y1 - y2);
    return sqrt(dx * dx + dy * dy);
}

int count(double cx, double cy, int i, int j) {
    int ans = 0;
    for (int k = 0; k < n; ++k) {
        if (k == i || k == j || dist(cx, cy, x[k], y[k]) < 1) ans++;
    }
    return ans;
}

void solve() {
    for (int i = 0; i < n; ++i) {
        cin >> x[i] >> y[i];
    }

    int ans = 1;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            if (i == j) continue;
            double d = dist(x[i], y[i], x[j], y[j]);
            if (d > 2.) continue;
            double dx = x[j] - x[i];
            double dy = y[j] - y[i];
            double r = sqrt(1 - (d * d / 4));
            // cout << r << endl;
            double cx = x[i] + dx / 2 - dy * r / d;
            // cout << d << endl;
            double cy = y[i] + dy / 2 + dx * r / d;
            // cout << x[i] << ", " << y[i] << " - " << x[j] << ", " << y[j] << " : " << cx << ", "<< cy << endl;
            ans = max(ans, count(cx, cy, i, j));
        }
    }
    cout << ans << endl;
}

int main() {
    while (1) {
        cin >> n;
        if (n == 0) break;
        solve();
    }
}

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

#define rep(i, n)       rep2(i, 0, n)
#define rep2(i, m, n)   for(int i = m; i < (int)(n); ++i)
#define X real()
#define Y imag()

typedef complex<double> Pt;

const double EPS = 1e-9;

int main()
{
    int n;
    while (cin >> n, n) {
        double x, y;
        vector<Pt> p;
        rep (i, n) {
            cin >> x >> y;
            p.push_back(Pt(x, y));
        }
        
        int ans = 1;
        rep (i, n) rep2 (j, i + 1, n) {
            double d = abs(p[i] - p[j]);
            if (d > 2.0) continue;
            
            Pt m = (p[i] + p[j]) / 2.0;
            Pt v = Pt(p[i].Y - p[j].Y, p[j].X - p[i].X);
            v /= abs(v);
            v *= sqrt(1.0 - d / 2.0);
            
            int t1 = 0, t2 = 0;
            Pt c1 = m + v, c2 = m - v;
            rep (k, n) {
                if (abs(p[k] - c1) < 1.0 + EPS) ++t1;
                if (abs(p[k] - c2) < 1.0 + EPS) ++t2;
            }
            ans = max(max(t1, t2), ans);
        }
        
        cout << ans << endl;        
    }
}

C++: #include <cstdio>
#include <complex>
#include <vector>
using namespace std;
typedef complex<double> P;
const double PI=3.141592653589793238462643383279;

bool inC(P p, P c, double r){return abs(p-c)<r+1e-9;}

P getE(P p){return p/abs(p);}

double sq(double x){ return  x*x;}

int main(){
    int n,i,j,k,ans,cnt;
    double x, y;
    while(scanf("%d",&n)){
        if(n==0) break;
        ans=0;
        vector<P> ps(n);
        for(i=0; i<n; i++){
            scanf("%lf%lf",&x,&y);
            ps[i] = P(x,y);
        }
        for(i=0;i<n;i++)
            for(j=i+1;j<n;j++)if(abs(ps[i]-ps[j])<2){
                double d=sqrt(1-sq(abs(ps[j]-ps[i])*0.5));
                P c1=0.5*(ps[i]+ps[j])+getE(getE(ps[j]-ps[i]))*polar(d,PI/2);
                P c2=0.5*(ps[i]+ps[j])+getE(getE(ps[j]-ps[i]))*polar(d,-PI/2);
                cnt=0;
                for(k=0;k<n;k++)if(inC(ps[k],c1,1.0))cnt++;
                if(cnt>ans)ans=cnt;
                cnt=0;
                for(k=0;k<n;k++)if(inC(ps[k],c2,1.0))cnt++;
                if(cnt>ans)ans=cnt;
            }
        if(ans==0) ans++;
        printf("%d\n",ans);
    }
    return 0;
}

C++: #include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <complex>

using namespace std;

#define all(c) ((c).begin()),((c).end())
#define debug(c) cerr<<"> "<<#c<<" = "<<(c)<<endl;
#define iter(c) __typeof((c).begin())
#define tr(i,c) for(iter(c) i=(c).begin();i!=(c).end();i++)
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define REP(i,a,b) for(int i=(int)(a);i<=(int)(b);i++)

#define mp make_pair
#define fst first
#define snd second
#define pb push_back

typedef long long ll;
typedef vector<int> vi;
typedef pair<int,int> pii;

#define EPS 1.0e-10

//ツ点
typedef complex<double> P;

struct C{
	P p;
	double r;
};

vector<P> isCC(const C &c1, const C &c2) {
	double d = abs(c1.p - c2.p);
	double x = (d*d + c1.r * c1.r - c2.r * c2.r) / 2 / d;
	double y = c1.r * c1.r - x * x;
	if(y < -EPS) return vector<P>();
	if(y < 0) y = 0;
	P diff = (c2.p - c1.p) / d;
	vector<P> res(2);
	res[0] = c1.p + diff * P(x, y);
	res[1] = c1.p + diff * P(x, -y);
	return res;
}

int main() {
	for (int N; cin >> N, N != 0; ) {
		vector<P> ps;
		rep (i, N) {
			double x, y; cin >> x >> y;
			ps.pb(P(x, y));
		}
		int res = 0;
		rep (i, N) REP (j, i + 1, N-1) {
			vector<P> isp = isCC((C){ps[i],1},(C){ps[j],1});
			rep (k, isp.size()) {
				int cnt = 0;
				rep (l, N) if (abs(isp[k] - ps[l]) <= 1) cnt++;
				res = max(res, cnt);
			}
		}
		cout << res <<endl;
	}
	
	return 0;
}

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

using namespace std;

#define EPS 1e-8

struct Point{
  double x,y;
  Point(){}
  Point(double x,double y) : x(x),y(y) {}
  Point operator + (const Point &p)const{
    return Point(x+p.x,y+p.y);
  }
  Point operator / (double a)const{
    return Point(x/a,y/a);
  }
  double d(const Point &p){
    return abs(sqrt(pow(x-p.x,2)+pow(y-p.y,2)));
  }
};

int main(){
  int N;
  while(cin >> N,N){
    Point p[N];
    for(int i = 0 ; i < N ; i++){
      cin >> p[i].x >> p[i].y;
    }
    int res = 1;
    for(int i = 0 ; i < N ; i++){
      for(int j = i+1 ; j < N ; j++){
        if(p[i].d(p[j]) + EPS > 2.0){
          continue;
        }
        Point cp = (p[i]+p[j])/2.0;
        double y = cp.d(p[i]);
        double x = sqrt(1-pow(y,2));
        Point c[2] = {cp,cp};
        c[0].x += x; c[1].x -= x;
        for(int k = 0 ; k < 2 ; k++){
          int cnt = 2;
          for(int l = 0 ; l < N ; l++){
            if(i == l || j == l) continue;
            if(c[k].d(p[l]) + EPS < 1.0){
              cnt++;
            }
          }
          res = max(res,cnt);
        }
      }
    }
    cout << res << endl;
  }
  return 0;
}