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>

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

using namespace std;
namespace geom{
#define at(i) ((*this)[i])
#define pb push_back
#define X real()
#define Y imag()
#define SELF (*this)
typedef long double R;
typedef complex<R> P;
const R EPS = 1e-8;
const R PI=3.14159265358979323846264338327950288;
enum { TURE = 1, FALSE = 0, BORDER= -1};
inline int sig(const R&x){return (abs(x)<EPS?0:x>0?1:-1);}
inline int less(const R &x,const R &y){return sig(x-y)?x<y:BORDER;}
inline R inp(const P &a,const P &b){return (conj(a)*b).X;}
inline R outp(const P &a,const P & b){return (conj(a)*b).Y;}
inline R norm(const P &p){return p.X*p.X+p.Y*p.Y;}
inline P unit(const P & p){return p/abs(p);}
inline P proj(const P &s,const P &t){return t*inp(s,t)/norm(t);}
struct L :public vector<P>{
  L(const P &p1,const P &p2){this->pb(p1);this->pb(p2);}
  L(){}
  P dir()const {return at(1)-at(0);}
  int online(const P &p) const {return !sig(outp(p-at(0),dir()));}
};
struct S: public L {
    S(const P &p1, const P &p2): L(p1,p2){}
    S(){}
    int online(const P &p) const{
        if(!sig(norm( p-at(0) )) || !sig(norm(p - at(1))) ) return BORDER;
        return !sig(outp(p-at(0), dir())) && inp(p-at(0), dir()) > EPS && inp(p-at(1),-dir()) > -EPS;
    }
};
struct C :public P{
  C(){}
  C(const P &p,const R r):P(p),r(r){}
  R r;
  int inside(const P &p)const { return less(norm(p-SELF),r*r);}
};
//inline P proj(const P &s,const L &t){return t[0] + proj(s-t[0], t[1]-t[0]);}

inline int intersect (const C &a,const C &b){
  return less((a.r-b.r)*(a.r-b.r),norm(a-b))+less(norm(a-b),(a.r+b.r)*(a.r+b.r))-1;
}

inline S crosspoint(const C &c1,const C &c2){
  if(!intersect(c1,c2)) return S();
  R d=abs(c1-c2);
  R x=(c1.r*c1.r-c2.r*c2.r+d*d)/(2*d);
  R h=sqrt(max<R>(0.,c1.r*c1.r-x*x));
  P u=unit(c2-c1);
  return S(c1+u*x+u*P(0,-1)*h,c1+u*x+u*P(0,1)*h);
}

// inline S crosspoint(const C &c,const L &l){
//   R d2=dist2(l,c);
//   if(c.r*c.r+EPS < d2) return S();
//   P m= proj(c,l);
//   P u = unit(l[1]-l[0]);
//   R d=sqrt(max<R>(.0,c.r*c.r-d2));
//   return S(m+u*d,m-u*d);
// }

S circlePPR(const P& a, const P& b, R r){
    return crosspoint( C(a,r), C(b,r) );
}
}
using namespace geom;

int main(){
    int N;
    while(cin>>N, N){
        int ans = 1;
        P po[305];
        rep(i,N){
            R x, y;
            cin >> x >> y;
            po[i] = P(x, y);
        }
        rep(i,N) rep(j,i){
            //if( i==j ) continue;
            S s = circlePPR( po[i], po[j], 1.0 );
            //if( s.empty() ) continue;
            //for(auto t: s) cout << t << " "; cout << endl;
            for(auto k: s){
                int count = 0;
                rep(h,N){
                    if( abs(k-po[h]) < 1.0 + EPS ){
                        count++;
                    }
                }
                ans = max(ans, count);
            }
        }
        cout << ans << endl;
    }
    return 0;
}

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

const double EPS = 1e-9;

typedef struct
{
  double x, y;
}Point;

int main()
{
  int N;
  while (~scanf("%d", &N)) {
    if (N == 0)
      break;

    Point c;
    vector<Point> p(N);
    for (int i = 0; i < N; ++i)
      scanf("%lf %lf", &p[i].x, &p[i].y);

    int ans = 1;
    for (int i = 0; i < N; ++i) {
      for (int j = i+1; j < N; ++j) {
	double H = hypot(p[i].x-p[j].x, p[i].y-p[j].y);
	if (H - 2.0 > EPS)
	  continue;
	double D = sqrt(1.0 - (H*H)/4.0);
	for (int pm = 0; pm < 2; ++pm) {
	  c.x = (p[i].x+p[j].x)/2.0 + (p[j].y-p[i].y)*D/H;
	  c.y = (p[i].y+p[j].y)/2.0 - (p[j].x-p[i].x)*D/H;
	  int cnt = 2;
	  for (int k = 0; k < N; ++k) {
	    if (k == i || k == j)
	      continue;
	    if ((c.x-p[k].x)*(c.x-p[k].x)+(c.y-p[k].y)*(c.y-p[k].y)-1.0 <= EPS)
	      ++cnt;
	  }
	  ans = max(ans, cnt);
	  D *= -1.0;
	}
     }
    }
    printf("%d\n", ans);
  }
  return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int,int> P;
typedef pair<int,P> P1;
typedef pair<P,P> P2;
#define pu push
#define pb push_back
#define mp make_pair
#define INF 1000000000
#define fi first
#define sc second
#define rep(i,x) for(int i=0;i<x;i++)
#define SORT(x) sort(x.begin(),x.end())
#define ERASE(x) x.erase(unique(x.begin(),x.end()),x.end())
#define POSL(x,v) (lower_bound(x.begin(),x.end(),v)-x.begin())
#define POSU(x,v) (upper_bound(x.begin(),x.end(),v)-x.begin())
typedef complex<double> pt;
typedef pair<pt,pt> L;
typedef vector<P> poly;
const double EPS = 1e-11;
#define x real()
#define y imag()
 
double dot(pt a,pt b){
	return (conj(a)*b).x;
}
double cross(pt a,pt b){
	return (conj(a)*b).y;
}
int ccw(pt a,pt b,pt c){
	b -= a; c -= a;
	if(cross(b,c) > EPS) return 1; // counter clockwise
	if(cross(b,c) < -EPS) return -1; // clockwise
	if(dot(b,c) < -EPS) return 2; //c-a-b
	if(norm(b) < norm(c)) return -2; //a-b-c
	return 0; //a-c-b
}
bool cmp(const pt& a,const pt& b){
	if(-EPS < a.x-b.x && a.x-b.x < EPS) return a.y < b.y;
	else return a.x < b.x;
}
vector<pt>convex_hull(vector<pt>ps)
{
	sort(ps.begin(),ps.end(),cmp);
	int k=0,n = ps.size();
	vector<pt>qs(n*2);
	for(int i=0;i<n;i++)
	{
		while(k>1 && ccw(qs[k-2],qs[k-1],ps[i]) == -1) k--;
		qs[k++] = ps[i];
	}
	for(int i=n-2,t=k;i>=0;i--)
	{
		while(k>t && ccw(qs[k-2],qs[k-1],ps[i]) == -1) k--;
		qs[k++] = ps[i];
	}
	qs.resize(k-1);
	return qs;
}
pt intersect(pt a,pt b,pt c,pt d){
	double dev = (b.y-a.y)*(d.x-c.x)-(b.x-a.x)*(d.y-c.y);
	
	double d1, d2;
 
	d1 = (c.y*d.x-c.x*d.y);
	d2 = (a.y*b.x-a.x*b.y);
 
	pt ap1;
	ap1.x = d1*(b.x-a.x) - d2*(d.x-c.x);
	ap1.x /= dev;
	ap1.y = d1*(b.y-a.y) - d2*(d.y-c.y);
	ap1.y /= dev;
 
	return ap1;
}
bool par(pt a,pt b,pt c,pt d){
	b.x -= a.x; b.y -= a.y;
	d.x -= c.x; d.y -= c.y;
	if(abs(b.y*d.x - b.x*d.y) < EPS){
		return 1;
	}
	return 0;
}
double PI = 3.1415926535897932;
int main(){
	while(1){
		int n; scanf("%d",&n); if(!n) return 0; if(n==1){ puts("1"); continue;}
		pt p[305];
		for(int i=1;i<=n;i++){
			double a,b; scanf("%lf%lf",&a,&b);
			p[i] = pt(a,b);
		}
		int ans = 0;
		for(int i=1;i<=n;i++){
			for(int j=i+1;j<=n;j++){
				pt q = (p[i]+p[j])/2.0;
				pt r = (p[i]-p[j]); if(abs(r) > 2.0+EPS) continue;
				pt p1 = q + r / abs(r) * polar(sqrt(1.0-abs(r)*abs(r)/4.0),PI/2);
				pt p2 = q - r / abs(r) * polar(sqrt(1.0-abs(r)*abs(r)/4.0),-PI/2);
				int c=0,d=0;
				for(int k=1;k<=n;k++){
					if(k==i || k==j || abs(p1-p[k]) <= 1.0+EPS) c++;
					if(k==i || k==j || abs(p2-p[k]) <= 1.0+EPS) d++;
				}
				ans = max(ans,max(c,d));
			}
		}
		printf("%d\n",ans);
	}
}