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;

void calc(double x1,double y1,double x2,double y2,double *c_x1,double *c_y1,double *c_x2,double *c_y2){
  double midx,midy,hox,hoy,dx,dy,dis;
  midx = (x1+x2) / 2 , midy = (y1+y2) / 2;
  dx = x1-x2 , dy = y1-y2;
  dis = sqrt(1 - pow(x1-midx,2) - pow(y1-midy,2));
  hox = sqrt(pow(dy,2) / (pow(dy,2) + pow(dx,2)));
  hoy = hox * (-1) * (dx / dy);
  *c_x1 = midx + hox*dis,*c_y1 = midy + hoy*dis;
  *c_x2 = midx - hox*dis,*c_y2 = midy - hoy*dis;
}

int main(){
  int N;

  while(cin >> N){
    double x[300],y[300],c_x1,c_y1,c_x2,c_y2;
    int count,m;
    m = 0;
    if(!(N)) break;
    for(int i=0;i<N;i++) cin >> x[i] >> y[i];

    for(int i=0;i<N;i++){
      for(int j=0;j<N;j++){
        if(i == j) continue;
        if((x[i] - x[j])*(x[i] - x[j]) + (y[i] - y[j])*(y[i] - y[j]) > 4.0 ) continue;
        count = 2;
        calc(x[i],y[i],x[j],y[j],&c_x1,&c_y1,&c_x2,&c_y2);
        for(int k=0;k<N;k++){
          if(k == i || k == j) continue;
          if((x[k] - c_x1)*(x[k] - c_x1) + (y[k] - c_y1)*(y[k] - c_y1) <= 1.0  ) count++;
        }
        if(m < count) m = count;
        count = 2;
        for(int k=0;k<N;k++){
          if(k == i || k == j) continue;
          if((x[k] - c_x2)*(x[k] - c_x2) + (y[k] - c_y2)*(y[k] - c_y2) <= 1.0  ) count++;
        }
        if(m < count) m = count;
      }
    }

    printf("%d\n",m);
  }

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

//ツ点
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);
	if (c1.r + c2.r - d < -EPS || fabs(c1.r - c2.r) - d > EPS) return vector<P>();
	double x = (d*d + c1.r * c1.r - c2.r * c2.r) / 2 / d;
	double y = sqrt(c1.r * c1.r - x * x);
	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 + EPS) cnt++;
				res = max(res, cnt);
			}
		}
		if (res == 0) res++;
		cout << res <<endl;
	}
	
	return 0;
}

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

using namespace std;

#define x first
#define y second
const int N = 300 + 1;

typedef pair<double, double> point;
point points[N];
double distances[N][N];

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

double dist(point a, point b) {
  return sqrt(sq(a.x - b.x) + sq(a.y - b.y));
}

int main() {
  int n;
  cin >> n;
  while (n) {
    for (int i = 0; i < n; ++i) {
      double x, y;
      cin >> x >> y;
      points[i] = {x, y};
    }

    for (int i = 0; i < n; ++i)
      for (int j = 0; j < n; ++j)
        distances[i][j] = dist(points[i], points[j]);

    // for each point,
    //  for every other point
    //    for every point currently selected
    //      check if it goes past the bounds
    //    add the point if its contained
    //
    int m = 0;
    for (int i = 0; i < n; ++i) {
      vector<int> selected = {i};
      for (int j = 0; j < n; ++j) {
        if (i == j) continue;
        bool contained = true;
        for (auto k : selected) {
          if (distances[j][k] > 2.0) {
            contained = false;
            break;
          }
        }
        if (contained) selected.push_back(j);
        /*
        cout << i << ' ' << j << ": " << endl;
        for (auto p : selected) cout << p << ' ';
        cout << endl;
        */
      }

      m = max(m, (int)selected.size());
    }

    cout << m << endl;

    cin >> n;
  }

}

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

using namespace std;

//conversion
//------------------------------------------
inline int toInt(string s) {int v; istringstream sin(s);sin>>v;return v;}
template<class T> inline string toString(T x) {ostringstream sout;sout<<x;return sout.str();}

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

//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())
#define MT(a,b,c) MP(a, MP(b, c))

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

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

int count(double x0, double y0, int p, int q){
	int ret = 2;
	REP(i,n){
		if(i==p||i==q)continue;
		double dx = x0-x[i];
		double dy = y0-y[i];
		if(pow(dx, 2.0)+pow(dy, 2.0)<1.0)ret++;
	}
	return ret;
}

int main(){
	while(1){
		int ret = 1;
		cin>>n;
		if(n==0)break;
		REP(i,n)scanf("%f %f", &x[i],&y[i]);
		REP(i,n){
			FOR(j, i+1,n){
				double dx = x[j]-x[i];
				double dy = y[j]-y[i];
				double d = sqrt(pow(dx, 2.0)+pow(dy, 2.0));
				if(d < 1.0){
					double mx = (x[j]+x[i])*0.5;
					double my = (y[j]+y[i])*0.5;
					{
						double x = mx-dy/d*sqrt(1-d*d/4.0);
						double y = my+dx/d*sqrt(1-d*d/4.0);
						ret = max(ret, count(x, y, i, j));
					}
					{
						double x = mx+dy/d*sqrt(1-d*d/4.0);
						double y = my-dx/d*sqrt(1-d*d/4.0);
						ret = max(ret, count(x, y, i, j));
					}
				}
			}
		}
		cout << ret << endl;
	}
	return 0;
}