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 int long long
#define reps(i,s,n) for(int (i)=(s);(i)<(n);++(i))
#define rep(i,n) reps(i,0,n)
#define rept(i,n) rep(i,(n)+1)
#define repst(i,s,n) reps(i,s,(n)+1)
#define reprt(i,n,t) for(int (i)=(n);(i)>=(t);--(i))
#define repr(i,n) reprt(i,n,0)
#define each(itr,v) for(auto &(itr):(v))
#define all(c) (c).begin(),(c).end()
#define rall(c) (c).rbegin(),(c).rend()
#define pb push_back
#define mp make_pair
#define fi first
#define se second
#define tmax(x,y,z) max(x,max(y,z))
#define tmin(x,y,z) min(x,min(y,z))
#define chmin(x,y) x=min(x,y)
#define chmax(x,y) x=max(x,y)
#define ln '\n'
#define bln(i,n) (((i)==(n)-1)?'\n':' ')
#define dbg(x) cout<<#x" = "<<(x)<<ln<<flush
#define dbga(x,n) {cout<<#x" : ";for(int (i)=0;i<(n);++i){cout<<((x)[i])<<(i==((n)-1)?'\n':' ')<<flush;}}
#define zero(a) memset(a,0,sizeof(a))
#define unq(a) sort(all(a)),a.erase(unique(all(a)),a.end())

typedef complex<double> P;
typedef long double ld;
typedef long long ll;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef vector<string> vst;
typedef vector<bool> vb;
typedef vector<ld> vld;
typedef vector<pii> vpii;
typedef vector<pll> vpll;
typedef vector<vector<int> > mat;

const ll inf = (ll)1e9+10;
const ll linf = (ll)1e18+10;
const ll mod = (ll)(1e9+7);
const int dx[] = {0, 1, 0, -1};
const int dy[] = {1, 0, -1, 0};
const int ddx[] = {0, 1, 1, 1, 0, -1, -1, -1};
const int ddy[] = {1, 1, 0, -1, -1, -1, 0, 1};
const double eps = 1e-10;

ll mop(ll a,ll b,ll m=mod) {ll r=1;a%=m;for(;b;b>>=1){if(b&1)r=r*a%m;a=a*a%m;}return r;}
ll gcd(ll a,ll b) {return b?gcd(b,a%b):a;}
ll lcm(ll a,ll b) {return a*b/gcd(a,b);}
bool ool(int x,int y,int h,int w) {return((x<0)||(h<=x)||(y<0)||(w<=y));}
bool deq(double a,double b) {return fabs(a-b)<eps;}

struct oreno_initializer {
	oreno_initializer() {
		cin.tie(0);
		ios::sync_with_stdio(0);
	}
} oreno_initializer;

// ━━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…
// .｡.:( ^ω^)・ﾟ＋.｡.:( ^ω^)・ﾟ＋.｡.:( ^ω^)・ﾟ＋.｡.:( ^ω^)・ﾟ＋.｡.:( ^ω^)・ﾟ＋
// ・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・‥…━━━☆・

int n;
double x[300], y[300];

signed main() {
	while (1) {
		cin >> n;
		if (n==0) break;
		rep(i,n) cin >> x[i] >> y[i];
		
		int res = 1;
		// 2点に接する円のみ考える
		rep(i,n) reps(j,i+1,n) {
			double d = sqrt((x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]));
			if (d>2) continue;
			P I = P(x[i],y[i]), J = P(x[j],y[j]), C = (I+J)/2.0;
			double l = sqrt(1 - abs(I-C)*abs(I-C));
			P HA = P(I-J) * P(0,1), HB = P(I-J) * P(0,-1);
			HA /= abs(HA), HA *= l;
			HB /= abs(HB), HB *= l;
			double xi = (C+HA).real(), yi = (C+HA).imag(), xj = (C+HB).real(), yj = (C+HB).imag();
			int t = 0;
			rep(k,n) if (sqrt((xi-x[k])*(xi-x[k]) + (yi-y[k])*(yi-y[k]))<=1) t++;
			chmax(res,t);
			t = 0;
			rep(k,n) if (sqrt((xj-x[k])*(xj-x[k]) + (yj-y[k])*(yj-y[k]))<=1) t++;
			chmax(res,t);
		}
		cout << res << ln;
	}
}

C++: #include<cstdio>
#include<cmath>
#include<iostream>
using namespace std;
#define EPS 1.0e-6
#define F(I,J,N) for(int I=J;I<N;I++)

class point{
public:
	double x,y;
	point(){}
	point(double x,double y){
		this->x=x;
		this->y=y;
	}
	double len(point &p){
		return sqrt((x-p.x)*(x-p.x)+(y-p.y)*(y-p.y));
	}
	double len(){
		return sqrt(x*x+y*y);
	}
	point operator*(double d){
		point p(x*d,y*d);
		return p;
	}
	point operator/(double d){
		point p(x/d,y/d);
		return p;
	}
	point operator+(const point &o){
		point p(x+o.x,y+o.y);
		return p;
	}
	point operator-(const point &o){
		point p(x-o.x,y-o.y);
		return p;
	}
	point& operator=(const point &o){
		this->x=o.x;
		this->y=o.y;
		return *this;
	}
	bool operator==(const point &o){
		return (x==o.x && y==o.y);
	}
	bool operator!=(const point &o){
		return !(x==o.x && y==o.y);
	}
};

int main(){
	int n;
	point *p;
	int ans;
	point c,ce,h;
	while(cin>>n){
		c=ce;
		if(n==0) break;
		p=new point[n];
		ans=0;
		F(i,0,n) cin>>p[i].x>>p[i].y;
		F(i,0,n-1){
			F(j,i+1,n){
				if(p[i].len(p[j])<2+EPS){
					c=(p[i]+p[j])/2;
					h.x=-(p[i]-p[j]).y;
					h.y=(p[i]-p[j]).x;
					double len=p[i].len(p[j]);
					len=sqrt(1-len*len/4);
					h=h*(len/h.len());

					ce=c+h;
					int cnt=0;
					F(k,0,n) if(ce.len(p[k])<1+EPS) cnt++;
					ans=max(ans,cnt);

					ce=c-h;
					cnt=0;
					F(k,0,n) if(ce.len(p[k])<1+EPS) cnt++;
					ans=max(ans,cnt);
				}
			}
		}
		cout<<ans<<endl;
		delete []p;
	}
	return 0;
}

C++: #include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>

#include <iostream>
#include <complex>
#include <string>
#include <algorithm>
#include <vector>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <unordered_map>
#include <unordered_set>

#include <functional>
#include <cassert>

typedef long long ll;
using namespace std;

#define debug(x) cerr << __LINE__ << " : " << #x << " = " << (x) << endl;

#define mod 1000000007 //1e9+7(prime number)
#define INF 1000000000 //1e9
#define LLINF 2000000000000000000LL //2e18
#define SIZE 300

typedef double P_type;
typedef complex<P_type> P;
const P_type P_eps = 1e-8; //整数の時はゼロ

namespace std{
  template<class T> bool operator<(const complex<T> &a, const complex<T> &b){
    return a.real() == b.real() ? a.imag() < b.imag() : a.real() < b.real();
  }
};

/* 円の点包含判定 */
double isContainedCP(P c, double r, P p){
  //return abs(c-p) < r - P_eps; //円周上を含まない
  return abs(c-p) < r + P_eps; //円周上を含む
}

int solve(){
  int n;
  double x[SIZE], y[SIZE];

  scanf("%d", &n);
  if(n == 0) return false;

  for(int i=0;i<n;i++){
    scanf("%lf%lf", x+i, y+i);
  }

  int ans = 0;
  
  for(int i=0;i<n;i++){
    for(int j=0;j<n;j++){
      if(i == j) continue;
      P vec = P(x[i]-x[j], y[i]-y[j]);

      if(abs(vec) > 2.0 + P_eps) continue;
      
      P p = P(-vec.imag(), vec.real());
      p = p / abs(p) * sqrt(1 - norm(vec) / 4);
      //debug(p);

      p += vec / 2.0 + P(x[j], y[j]); 

      int counter = 0;
      // debug(abs(p - P(x[i], y[i])));
      
      for(int k=0;k<n;k++){
        counter += isContainedCP(p, 1, P(x[k], y[k]));
      }

      ans = max(ans, counter);
    }
  }

  printf("%d\n", ans);
  return true;
}


int main(){
  while(solve());
  return 0;
}