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 <cstdio>
#include <cstdlib>
#include <cstring>
#include <climits>
#include <cmath>
#include <ctime>
#include <cassert>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <complex>
#include <string>
#include <vector>
#include <list>
#include <deque>
#include <stack>
#include <queue>
#include <set>
#include <map>
#include <bitset>
#include <iterator>
#include <functional>
#include <utility>
#include <algorithm>
#include <numeric>
#include <typeinfo>

using namespace std;

#define dump(n) cout<<"# "<<#n<<"="<<(n)<<endl
#define repi(i,a,b) for(int i=int(a);i<int(b);i++)
#define rep(i,n) repi(i,0,n)
#define iter(c) __typeof((c).begin())
#define foreach(i,c) for(iter(c) i=(c).begin();i!=(c).end();++i)
#define allof(c) (c).begin(),(c).end()
#define mp make_pair

typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef vector<double> vd;
typedef vector<vd> vvd;
typedef vector<string> vs;
typedef pair<int,int> pii;

const double EPS=1e-10;

struct Square{
	double x,y,d,n;
	Square(){}
	Square(double x,double y,double d,int n):x(x),y(y),d(d),n(n){}
	bool operator<(const Square& s)const{return n<s.n;}
	bool operator>(const Square& s)const{return n>s.n;}
};

inline double myhypot(double x,double y)
{
	return sqrt(x*x+y*y);
}

int main()
{
	for(int n;cin>>n,n;){
		vd xs(n),ys(n);
		rep(i,n)
			cin>>xs[i]>>ys[i];
		double xmin=*min_element(allof(xs));
		double xmax=*max_element(allof(xs));
		double ymin=*min_element(allof(ys));
		double ymax=*max_element(allof(ys));
		priority_queue<Square> pq;
		pq.push(Square((xmin+xmax)/2,(ymin+ymax)/2,max(xmax-xmin,ymax-ymin)/2,n));
		int res=1;
		while(!pq.empty()){
			Square s=pq.top(); pq.pop();
			if(s.n<=res)
				continue;
			
			int temp=0;
			rep(i,n)
				if(myhypot(xs[i]-s.x,ys[i]-s.y)<1+EPS)
					temp++;
			res=max(res,temp);
			
			if(s.d<1e-5)
				continue;
			
			double nd=s.d/2;
			double nx[]={s.x-nd,s.x-nd,s.x+nd,s.x+nd};
			double ny[]={s.y-nd,s.y+nd,s.y-nd,s.y+nd};
			rep(k,4){
				double x=nx[k],y=ny[k];
				temp=0;
				rep(i,n){
					if(myhypot(x-xs[i],y-ys[i])<nd*sqrt(2)+1+EPS)
						temp++;
				}
				if(temp>res)
					pq.push(Square(x,y,nd,temp));
			}
		}
		cout<<res<<endl;
	}
	
	return 0;
}

C++: #define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include <complex>
#include <cstdio>
#include <cstring>
#include <ctime>
#include <string>
#include <vector>
#include <algorithm>
#include <cmath>
#include <sstream>
#include <utility>
using namespace std;
typedef long long ll;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef pair<int,int> pii;
#define all(c) (c).begin(), (c).end()
#define loop(i,a,b) for(ll i=a; i<ll(b); i++)
#define rep(i,b) loop(i,0,b)
#define each(e,c) for(auto&e:c)
#define pb push_back
#define eb emplace_back
#define mp make_pair
#define mt make_tuple
#define lb lower_bound
#define ub upper_bound
#ifdef DEBUG
#define dump(...) (cerr<<#__VA_ARGS__<<" = "<<(DUMP(),__VA_ARGS__).str()<<" ["<<__LINE__<<"]"<<endl)
struct DUMP:ostringstream{template<class T>DUMP &operator,(const T&t){if(this->tellp())*this<<", ";*this<<t;return *this;}};
#else
#define dump(...)
#endif
template<class T> ostream& operator<<(ostream& os, vector<T> const& v){
    rep(i,v.size()) os << v[i] << (i+1==v.size()?"":" ");
    return os;
}

typedef double R;
typedef complex<R> P;
R const eps = 1e-8;
int main(){
    int n;
    while(cin>>n && n){
        vector<P> ps(n);
        rep(i,n){
            R x,y; cin >> x >> y;
            ps[i] = P(x,y);
        }

        int ans = 0;
        rep(i,n)rep(j,n)if(i!=j){
            P & p = ps[i];
            P & q = ps[j];
            R d = abs(p-q);
            if(d > 2) continue;
            P m = (p+q)/2.;
            R l = sqrt(1-(d/2)*(d/2));
            P c = m + (p-m)*(d/2)*l * P(0,1);
            int cnt = 0;
            rep(k,n){
                if(abs(c-ps[k]) < 1+eps) cnt++;
            }
            ans = max(ans,cnt);
        }
        cout << ans << endl;
    }
}

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

#define rep(i,N) for(int i=0;i<(int)(N);i++)
#define REP(i,k,N) for(int i=k;i<(int)(N);i++)

#define vi vector<int>
#define pb push_back

#define pdd pair<double,double>
#define fi first
#define se second
const double  EPS = 1e-8;

double dist(double x,double y){
    return x*x+y*y;
};

int main(){

    int n;
    while(cin >> n,n){
        int ma = 0;
        pdd p[330];
        rep(i,n) cin >> p[i].fi >> p[i].se;
        rep(i,n-1){
            REP(j,i+1,n){
                int cnt1 = 0,cnt2=0;
                double x1,y1,x2,y2,x3,y3;
                tie(x1,y1) = p[i];
                tie(x2,y2) = p[j];
                double cx1,cy1,cx2,cy2;
                if(dist(x1-x2,y1-y2)  > 4+EPS  ) continue;

                //はりまコーナー
//                double ax,ay,bx1,by1,bx2,by2,mx,my,dx,dy;
//                ax=x1-x2;ay=y1-y2;
//                mx=(x1+x2)/2;my=(y1+y2)/2;
//                by1 = (-0.5)*pow(ax,4)+pow(ax,2)+(-0.5)*pow(ax,2)*pow(ay,2)
//                        /dist(ax,ay);
//                by1 = sqrt(abs(by1));
//                by2 = -1*by1;
//                cy1= my+by1;
//                cy2= my+by2;
//                bx1= -1*(ay*by1)/ax;
//                bx2= -1*(ay*by2)/ax;
//                cx1=bx1+mx;
//                cx2=bx2+mx;


                //あじのりコーナー
                double mx=(x1+x2)/2,my=(y1+y2)/2;

                double ld = 1-dist((x1-x2),(y1-y2))/4;
                ld = sqrt(ld);
                double theta = atan((y2-y1)/(x2-x1));
                cx1 = mx - sin(theta)*ld;
                cy1 = my + cos(theta)*ld;
                cx2 = mx + sin(theta)*ld;
                cy2 = my - cos(theta)*ld;

                for(int k = 0;k<n;k++) {
                    tie(x3, y3) = p[k];
                    if (EPS + 1> dist((cx1 - x3), (cy1 - y3))) cnt1++;
                    if (EPS + 1> dist((cx2 - x3), (cy2 - y3))) cnt2++;
                }
                ma = max({cnt1,cnt2,ma});
            }
        }
        cout<< ma << endl;
    }
    return 0;
}