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 <vector>
#include <math.h>

using namespace std;

const double EPS = 1e-8;
typedef pair<double, double> pt;
#define x first
#define y second

pt operator+(pt p, pt q) {
  return (pt){p.x + q.x, p.y + q.y};
}
pt operator-(pt p, pt q) {
  return (pt){p.x - q.x, p.y - q.y};
}
pt operator*(pt p, double c) {
  return (pt){p.x * c, p.y * c};
}
double sq(double a) { return a * a; }
double dist(pt p, pt q) {
  return sqrt(sq(p.x - q.x) + sq(p.y - q.y));
}

pt normal(pt a, pt b) {
	pt AB = b - a;
	double d = dist(a, b);
	pt n = pt(AB.y/d, -AB.x/d);
	return n;
}

void solve(int N) {
	vector<pt> point;
	for(int i=0; i < N; i++) {
		double a, b;
		cin >> a >> b;
		point.push_back(pt(a, b));
	}

	int ans = 1;
	for(int i=0; i < N; i++) {
		for(int j=i+1; j < N; j++) {
			pt A = point[i];
			pt B = point[j];
			pt M = pt((A.x + B.x)/2, (A.y + B.y)/2);

			if(dist(A, B) > 2.0 + EPS) continue;

			double l = sqrt(1 - sq(dist(A, M)));
			pt V = normal(A, B);

			pt C[2];
			C[0] = M + V*l;
			C[1] = M - V*l;

			for(int t=0; t < 2; t++) {
				int count = 0;
				for(int k=0; k < N; k++) {
					if(dist(C[t], point[k]) < 1.0 + EPS) count++;
				}
				ans = max(ans, count);
			}
		}
	}
	cout << ans << endl;
}

int main() {

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

	return 0;
}

C++: #include<bits/stdc++.h>
#define N 300
using namespace std;
typedef complex<double> P;
const double EPS = 1e-5;

struct Circle{
  P c;
  double r;
  Circle(P c=P(0,0),double r=0):c(c),r(r){}
};

vector<P> getIntersectCC(Circle c1, Circle c2){
  vector<P> res;
  double r1=c1.r,r2=c2.r,d;
  P p1=c1.c,p2=c2.c;
  d=abs(p1-p2);
  if(d<EPS&&abs(r1-r2)<EPS) return res;
  if(r1+r2<d-EPS||d+EPS<abs(r1-r2)) return res; 
  double a=(r1*r1-r2*r2+d*d)/(2*d);
  double h=sqrt(max(r1*r1-a*a,0.0));
  P tmp1=p1+a/d*(p2-p1);
  P tmp2=h/d*(p2-p1);
  if(abs(tmp2)<EPS)res.push_back(tmp1);
  else{
    res.push_back(P(real(tmp1)-imag(tmp2),imag(tmp1)+real(tmp2)));
    res.push_back(P(real(tmp1)+imag(tmp2),imag(tmp1)-real(tmp2)));
  }
  return res;
}

double getdistancepp(P a,P b){
  return sqrt((a.real()-b.real())*(a.real()-b.real())+(a.imag()-b.imag())*(a.imag()-b.imag()));
}

int main(){
  int n,ans=0;
  double a,b;
  P p[N];
  while(1){
    cin>>n;
    if(!n)break;
    for(int i=0;i<n;i++){
      cin>>a>>b;
      p[i]=P(a,b);
    }
    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
	Circle c1,c2;
	c1.c=p[i];
	c2.c=p[j];
	c1.r=c2.r=1;
	vector<P> r=getIntersectCC(c1,c2);
	if(!r.size())continue;
	for(int k=0;k<r.size();k++){
	  int cnt=0;
	  for(int l=0;l<n;l++)
	    if(getdistancepp(r[k],p[l])<=1+EPS)cnt++;	    
	  ans=max(ans,cnt);
	}
      }
    }
    cout<<ans<<endl;
  }
  return 0;
}

C++: #include<stdio.h>
#include<math.h>
#include<algorithm>
using namespace std;
double x[300],y[300];
double px[90000],py[90000];
int main(){
    int a;
    while(scanf("%d",&a),a){
        for(int i=0;i<a;i++)scanf("%lf%lf",x+i,y+i);
        int now=0;
        for(int i=0;i<a;i++){
            for(int j=i+1;j<a;j++){
                if((x[i]-x[j])*(x[i]-x[j])+(y[i]-y[j])*(y[i]-y[j])<4.0000001){
                    double d=sqrt((x[i]-x[j])*(x[i]-x[j])+(y[i]-y[j])*(y[i]-y[j]));
                    double mx=(x[i]+x[j])/2;
                    double my=(y[i]+y[j])/2;
                    double s=sqrt(1-d*d/4);
                    double theta=atan((y[i]-y[j])/(x[i]-x[j]));
                    px[now]=mx-sin(theta)*s;
                    py[now++]=my+cos(theta)*s;
                    px[now]=mx+sin(theta)*s;
                    py[now++]=my-cos(theta)*s;
                }
            }
        }
        int ret=1;
        for(int i=0;i<now;i++){
            int k=0;
            for(int j=0;j<a;j++)if((px[i]-x[j])*(px[i]-x[j])+(py[i]-y[j])*(py[i]-y[j])<1.0000001)k++;
            ret=max(ret,k);
        }
        printf("%d\n",ret);
    }
}

C++: #define _USE_MATH_DEFINES
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <climits>
#include <cfloat>
#include <ctime>
#include <cassert>
#include <map>
#include <utility>
#include <set>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>
#include <sstream>
#include <complex>
#include <stack>
#include <queue>
#include <numeric>
#include <list>
#include <iomanip>


#ifdef _DEBUG
#define typeof(X) std::identity<decltype(X)>::type //C++0x (for vs2010)
#else
#define typeof(X) __typeof__(X) // for gcc
#endif

#define sz(a)  int((a).size())
#define FOREACH(it, c) for (typeof((c).begin()) it=(c).begin(); it != (c).end(); ++it)
#define FOR(i,count) for (int i = 0; i < (count); i++)

template <class T> void max_swap(T& a, const T& b) { a = max(a, b); }
template <class T> void min_swap(T& a, const T& b) { a = min(a, b); }

using namespace std;
static const double EPS = 1e-6;
typedef long long ll;
const int MODULO = 1000000007;
const int INF = 100000000; //1e8

typedef long long ll;
typedef pair<int,int> Pii;
typedef pair<ll,ll> Pll;
typedef complex<double> Cd;

int n;
Cd pt[300];

double sqr(double x) { return x * x; }

int C(Cd &o){
	int cnt = 0;
	for (int i = 0; i < n; i++){
		Cd c = pt[i] - o;
		if(sqr(c.real()) + sqr(c.imag()) <= 1.0 + EPS)
			cnt++;
	}
	return cnt;
}

int solve()
{
	for (int i = 0; i < n; i++){
		double x,y; cin>>x>>y;
		pt[i] = Cd(x,y);
	}

	int ans = 1;

	for (int i = 0; i < n; i++){
		for (int j = 0; j < n; j++){
			if(i == j) continue;
			Cd e = (pt[j] - pt[i]);
			double len = abs(e);
			if(len > 2.0 + EPS) continue;
			Cd md = pt[i] + e / 2.0;
			e /= len;
			e *= Cd(0,1);
			double o_dis = sqrt(1.0 - sqrt(len / 2));
			static int dir[] = {1,-1};
			for (int k = 0; k < 2; k++){
				Cd o = md + dir[k] * o_dis * e; // 中心を出す
				//printf("%d %d %lf %lf\n",i,j,o.real(),o.imag());
				ans = max(ans,C(o));
			}
		}
	};

	return ans;

}


int main(){
	while(cin>>n,n){
		cout << solve() << endl;
	}
	return 0;
}