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 <algorithm>
#include <iostream>
#include <cstdio>
#include <complex>
#include <vector>
#define X real()
#define Y imag()
#define EPS (1e-10)
using namespace std;
typedef complex<double> P;
typedef pair<P, double> C;

namespace std{
  bool operator<(const P& a, const P& b){
    return a.X != b.X ? a.X < b.X : a.Y < b.Y;
  }
}

P n_vector(P a){
  double l = abs(a);
  return P(-a.Y,a.X)/l;
}

vector<P> cpcc(C a, C b){
  double d = abs(a.first-b.first), ra = a.second, rb = b.second;
  double t = (ra*ra-rb*rb+d*d)/2/d, h = sqrt(ra*ra-t*t);
  P m = t/abs(b.first-a.first)*(b.first-a.first)+a.first;
  vector<P> cp(2,m);
  P n = n_vector(a.first-b.first);
  cp[0] -= h*n;
  cp[1] += h*n;
  return cp;
}


int main(){
  int N;
  double x, y;
  while(cin >> N, N){
    P ps[N];
    for(int i = 0; i < N; ++i){
      cin >> x >> y;
      ps[i] = P(x,y);
    }
    if(N == 1){
      cout << 1 << endl;
      continue;
    }
    vector<P> cps;
    for(int i = 0; i < N; ++i){
      for(int j = i+1; j < N; ++j){
	if(abs(ps[i]-ps[j]) < 2 + EPS){
	  vector<P> cp = cpcc(C(ps[i],1),C(ps[j],1));
	  for(int k = 0; k < cp.size(); ++k)
	    cps.push_back(cp[k]);
	}
      }
    }
    int n = cps.size(), ans = 0;;
    for(int i = 0; i < n; ++i){
      int cnt = 0;
      for(int j = 0; j < N; ++j){
	if(norm(cps[i]-ps[j]) < 1 + EPS)
	  ++cnt;
      }
      ans = max(ans, cnt);
    }
    cout << ans << endl;
  }
  return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;
/*{{{*/  //template
#define rep(i,n) for(int i=0;i<n;i++)
#define INF 1<<29
#define LINF LLONG_MAX/3
#define mp make_pair
#define pb push_back
#define EB emplace_back
#define fi first
#define se second
#define all(v) ALL(v)
#define sz(x) (int)(x).size()
#define debug(x) cerr<<#x<<":"<<x<<endl
#define debug2(x,y) cerr<<#x<<","<<#y":"<<x<<","<<y<<endl
//struct fin{ fin(){ cin.tie(0); ios::sync_with_stdio(false); } } fin_;
struct Double{ double d; explicit Double(double x) : d(x){} };
ostream& operator<<(ostream& os,const Double x){ os << fixed << setprecision(20) << x.d; return os; }
template<typename T> ostream& operator<<(ostream& os,const vector<T>& vec){ os << "["; for(const auto& v : vec){ os << v << ","; } os << "]"; return os; }
template<typename T,typename U> ostream& operator<<(ostream& os,const pair<T,U>& p){ os << "(" << p.first << ","<< p.second <<")"; return os; }
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
typedef vector<int> vi;
typedef vector<vi> vvi;
ll gcd(ll a,ll b){ if(b==0) return a; else return gcd(b,a%b); }
constexpr double eps = 1e-14; 
constexpr ll mod = 1e9+7;
const int dx[]={1,0,-1,0} ,dy[] = {0,1,0,-1};
/*}}}*/

int N;

void solve(){
    vector<double> x(N),y(N);
    rep(i,N){ cin >> x[i] >> y[i]; }

    int ans=1;
    rep(i,N) rep(j,N) if(i!=j){
        double vx = x[i]-x[j];
        double vy = y[i]-y[j];
        double d = sqrt(vx*vx+vy*vy);
        if(d>2) continue;
        double r = sqrt(1-(d/2)*(d/2));

        double rx =  -vy*r/d;
        double ry = vx*r/d;

        double cx = (x[i]+x[j])/2 + rx;
        double cy = (y[i]+y[j])/2 + ry;

        int cnt=0;
        rep(k,N){
            double tx = cx-x[k];
            double ty = cy-y[k];
            double dd = sqrt(tx*tx+ty*ty);
            if(dd<=1) cnt++;
        }
        ans = max(ans,cnt);
    }
    cout << ans << endl;
}

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

C++: #include<iostream>
#include<complex>
#include<vector>
using namespace std;
typedef complex<double> P;
int n;
vector<P> p;
int ret(){
	int cnt=1;
	for(int i=0;i<n;i++){
		for(int j=i+1;j<n;j++){
			if(abs(p[i]-p[j])>2.0)
				continue;
			P c=(p[i]+p[j])*0.5;
			double d=sqrt(1.0-abs(p[i]-c)*abs(p[i]-c));
			P l1=(p[i]-p[j])*P(0,1);
			P l2=(p[i]-p[j])*P(0,-1);
			l1=l1*d/abs(l1)+c;l2=l2*d/abs(l2)+c;
			int x=2,y=2;
			for(int k=0;k<n;k++){
				if(k==i||k==j)
					continue;
				if(abs(l1-p[k])<=1.0){
					x++;
				}
				if(abs(l2-p[k])<=1.0){
					y++;
				}
				
			}
			cnt=max(cnt,max(x,y));
		}
	}
	return cnt;
	
}
int main(){
	while(cin>>n,n){
		p.resize(n);
		for(int i=0;i<n;i++){
			double x,y;
			cin>>x>>y;
			p[i]=P(x,y);
		}
		cout<<ret()<<endl;
	}
	return 0;
}

C++: #include<bits/stdc++.h>
#define f first
#define s second
#define mp make_pair
#define inf 1000000000
#define eps (1e-11)
#define equals(a,b) (fabs((a)-(b))<eps)
using namespace std;

class Point{
public:
  double x,y;
  Point(double x=0,double y=0):x(x),y(y){}

  Point operator+(Point p){ return Point(x+p.x,y+p.y);}
  Point operator-(Point p){ return Point(x-p.x,y-p.y);}
  Point operator*(double k){ return Point(x*k,y*k);}
  Point operator/(double k){ return Point(x/k,y/k);}
  bool operator<(Point p)const{ return (x!=p.x ? x<p.x : y<p.y);}
  bool operator==(Point p)const{ return fabs(x-p.x)<eps && fabs(y-p.y)<eps;}

  double abs(){ return sqrt(norm());}
  double norm(){ return (x*x+y*y);}
};
typedef Point Vector;
typedef vector<Point> Polygon;

class Segment{
public:
  Point p1,p2;
  Segment(Point p1,Point p2):p1(p1),p2(p2){}
};
typedef Segment Line;

class Circle{
public:
  Point c;
  double r;
  Circle(Point c=Point(),double r=0.0):c(c),r(r){}
};

double norm(Vector a){ return (a.x*a.x+a.y*a.y);}
double abs(Vector a){ return sqrt(norm(a));}
double dot(Vector a,Vector b){ return (a.x*b.x+a.y*b.y);}
double cross(Vector a,Vector b){ return (a.x*b.y-a.y*b.x);}

Point median_Point(Point a,Point b){
  Vector v=a-b;
  return b+(v/2);
}

Point PointRotation(Point base,Point a,double r){
  Point b=a-base;
  a.x=b.x*cos((r/180)*M_PI)-b.y*sin((r/180)*M_PI);
  a.y=b.x*sin((r/180)*M_PI)+b.y*cos((r/180)*M_PI);
  a=a+base;
  return a;
}

vector<Point> vp;
int ans;
int n;

pair<Circle,Circle> makeCircle(Point a,Point b){
  Point m=median_Point(a,b);
  double d=sqrt(1-abs(a-b)/2);
  Point A=m+(PointRotation(m,a,90)-m)*d,B=m+(PointRotation(m,a,270)-m)*d;
  return mp(Circle(A,1),Circle(B,1));
}

int PointsInCircle(Circle C)
{
  int num=0;
  for(int i=0;i<n;i++){
    if((abs(vp[i]-C.c)-C.r)<+eps)num++;
  }
  return num;
}

int main()
{
  while(1){
    cin>>n;
    if(n==0)break;
    vp.clear();
    vp.resize(n);
    for(int i=0;i<n;i++)cin>>vp[i].x>>vp[i].y;
    ans=1;
    for(int i=0;i<n;i++){
      for(int j=0;j<n;j++){
	if(i==j)continue;
	if(2<abs(vp[i]-vp[j]))continue;
	pair<Circle,Circle> pc=makeCircle(vp[i],vp[j]);
	ans=max(ans,PointsInCircle(pc.f));
	ans=max(ans,PointsInCircle(pc.s));
      }
    }
    cout<<ans<<endl;
  }
  return 0;
}