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(int i=0;i<n;++i)
struct Point{double x,y;};
double dist(Point a,Point b) {return sqrt(pow(a.x-b.x,2)+pow(a.y-b.y,2));}

int main(void){
    int N;
    while(cin>>N,N){
        vector<Point> p(N);
        int maxcnt=1;
        rep(i,N)cin>>p[i].x>>p[i].y;
        rep(i,N){
            for(int j=i+1;j<N;j++){
                double d = dist(p[i],p[j]);
                if(d>2.0)continue;
                int sign[]={-1,1};
                rep(s,2){
                    Point C;
                    C.x = p[i].x + cos(atan2(p[j].y-p[i].y,p[j].x-p[i].x) + sign[s]*acos(d/2.0));
                    C.y = p[i].y + sin(atan2(p[j].y-p[i].y,p[j].x-p[i].x) + sign[s]*acos(d/2.0));
                    int cnt=0;
                    rep(k,N) cnt += (pow(C.x-p[k].x,2) + pow(C.y-p[k].y,2) <= 1.0001);
                    if(maxcnt<cnt)maxcnt=cnt;
                }
            }
        }
        cout<<maxcnt<<endl;
    }

    return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;

typedef pair<double,double> P;

bool cmp_x (const P& firstEof, const P& secondEof) {
    if (firstEof.first == secondEof.first) return firstEof.second < secondEof.second;
    return firstEof.first < secondEof.first;
}

int main () {
    int n;
    while (cin >> n, n) {
        vector<P> v;

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

        sort(v.begin(), v.end(), cmp_x);

        int res = 0;
        for (int i = 0; i < n; i++) {
            vector<P> tmp;
            tmp.push_back(v[i]);
            int k = i + 1;
            while (k < n && fabs(v[k].first - v[i].first) < 2.0) {
                bool flag = true;
                for (int j = 0; j < tmp.size(); j++) {
                    int dis = pow(v[k].first - tmp[j].first, 2.) + pow(v[k].second - tmp[j].second, 2.);
                    if (dis > 4.) {
                        flag = false;
                        break;
                    }
                }
                if (flag) {
                    tmp.push_back(v[k]);
                }
                k++;
            }
            res = max(res, (int)tmp.size());
        }

        cout << res << endl;
    }

    return 0;
}

C++: #include <cstdio>
#include <cmath>
#include <cstring>
#include <cassert>
#include <iostream>
#include <string>
#include <vector>
#include <queue>
#include <stack>
#include <algorithm>
#include <map>

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

void init() {
  for(int i = 0; i < (int)512; ++i) {
    x[i] = y[i] = 0.0;
  }
} 

inline double Heron(double a, double b, double c) {
  double s = (a + b + c) / 2;
  return sqrt(s*(s-a)*(s-b)*(s-c));
}

class solve_ {
public:
  int operator() (void) {
    int res = 1;
    for(int i = 0; i < (int)n; ++i) {
      for(int j = 0; j < (int)i; ++j) {
        res = std::max(res, CountPoints(i, j));
      }
    }
    return res;
  }
  int CountPoints(int i, int j) {
    double dx = x[i] - x[j];
    double dy = y[i] - y[j];
    double dist = sqrt(dx * dx + dy * dy);
    if( dist > 2.0 ) return 0;
    auto t1 = std::pair<double, double>((x[i] + x[j]) / 2, (y[i] + y[j]) / 2);
    auto t2 = std::pair<double, double>(dx / dist, dy / dist);
    double h = 2 * Heron(1, 1, dist) / dist;
    auto t3 = std::pair<double, double>(t2.first * h, t2.second * h);
    auto c1 = std::pair<double, double>(t1.first + t3.second, t1.second - t3.first);
    auto c2 = std::pair<double, double>(t1.first - t3.second, t1.second + t3.first);
    int count1 = 0;
    int count2 = 0;
    {
      double cx = c1.first, cy = c1.second;
      for(int k = 0; k < (int)n; ++k) {
        if( k == i or k == j ) continue;
        if( (cx - x[k]) * (cx - x[k]) + (cy - y[k]) * (cy - y[k]) <= 1 ) {
          count1 += 1;
        }
      }
    }
    {
      double cx = c2.first, cy = c2.second;
      for(int k = 0; k < (int)n; ++k) {
        if( k == i or k == j ) continue;
        if( (cx - x[k]) * (cx - x[k]) + (cy - y[k]) * (cy - y[k]) <= 1 ) {
          count2 += 1;
        }
      }
    }
    int res = std::max(count1, count2) + 2;
    //printf("(%d, %d) -> %d\n", i, j, res);
    return res;
  }
};
solve_ solve;

int main() {
  for(;;) {
    init();
    scanf("%d", &n);
    if( n == 0 ) break;
    for(int i = 0; i < (int)n; ++i) {
      scanf("%lf %lf", &x[i], &y[i]);
    }
    printf("%d\n", solve());
  }
  return 0;
}

C++: #include"bits/stdc++.h"
using namespace std;

#define REP(i,n) for(int i=0;i<(n);i++)
#define FOR(i,m,n) for(int i=m;i<(n);i++)
#define ALL(x) x.begin(),x.end()
#define ll long long

typedef pair<int,int> pii;

typedef complex<double> point;
const double R = 1.0;
const double EPS = 1e-6;
double dist(point &p1,point &p2){
	return abs(p1-p2);
}
struct Square{
	double d,u,l,r;
	int up_bound;
	
	Square(double _d,double _u, double _l,double _r,const vector<point>&ps):d(_d),u(_u),l(_l),r(_r){
		int res = 0;
		point c = {(l+r)/2,(d+u)/2};
		double ra = sqrt(2)*(r-l)/2+1;
		for(auto p:ps){
			if(dist(c,p)<ra+EPS)res++;
		}
		up_bound=res;
	}
	int count(const vector<point>&ps,const double &R){
		point c = {(l+r)/2,(d+u)/2};
		int res = 0;
		for(auto p:ps){
			if(dist(c,p)<R+EPS)res++;
		}
		return res;
	}
};
namespace std{
	bool operator<(const Square &S1, const Square &S2){
		return S1.up_bound<S2.up_bound;
	}
	bool operator>(const Square &S1, const Square &S2){
		return S1.up_bound>S2.up_bound;
	}
}
void solve() {
	int N;cin>>N;
	if(N==0)exit(0);

	int ans = 1;
	double l = 0,r = 10;
	vector<point>ps(N);
	REP(i,N){
		double x,y;cin>>x>>y;
		ps[i]={x,y};
	}
	Square s(-1e9,1e9,-1e9,1e9,ps);
	priority_queue<Square,vector<Square>,greater<Square>>que;
	que.push(s);
	REP(_,500){
		vector<Square>nxt;
		while(!que.empty()){
			Square s = que.top();que.pop();
			if(s.up_bound<=ans)continue;
			ans = max(ans,s.count(ps,R));
			double xl = s.l , xm = (s.l+s.r)/2.0,xr = s.r;
			double yl =s.d,ym = (s.d+s.u)/2.0,yr = s.u;
			Square s1(yl,ym,xl,xm,ps);
			if(s1.up_bound>ans)nxt.push_back(s1);

			Square s2(yl,ym,xm,xr,ps);
			if(s2.up_bound>ans)nxt.push_back(s2);
			
			Square s3(ym,yr,xl,xm,ps);
			if(s3.up_bound>ans)nxt.push_back(s3);
			
			Square s4(ym,yr,xm,xr,ps);
			if(s4.up_bound>ans)nxt.push_back(s4);
		}
		for(auto s:nxt)que.push(s);
	}
	cout<<ans<<endl;
}
signed main() {
	while (true) {
		solve();
	}
}