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-12;
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 = 0;
        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 <iostream>
#include <iomanip>
#include <sstream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <climits>
#include <cfloat>
#include <vector>
#include <string>
#include <stack>
#include <queue>
#include <set>
#include <map>
#include <algorithm>
#include <functional>
#include <utility>
#include <numeric>
#include <iterator>

using namespace std;

typedef long long LL;
typedef unsigned long long ULL;
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<string> VS;
typedef pair<int,int> PII;
typedef vector<PII> VPII;
typedef istringstream ISS;
typedef ostringstream OSS;

#define REP( i, m, n ) for ( int i = (int)( m ); i < (int)( n ); ++i )
#define FOR( v, c ) for ( auto &v : c )
#define EACH( it, c ) for ( auto it = c.begin(); it != c.end(); ++it )
#define ALL( c ) (c).begin(), (c).end()
#define DRANGE( c, p ) (c).begin(), (c).begin() + p, (c).end()

#define PB( n ) push_back( n )
#define MP( a, b ) make_pair( ( a ), ( b ) )
#define EXIST( c, e ) ( (c).find( e ) != (c).end() )

#define fst first
#define snd second

#define DUMP( x ) cerr << #x << " = " << ( x ) << endl
#define DEBUG( x ) cerr << __FILE__ << ":" << __LINE__ << ": " << #x << " = " << ( x ) << endl

#include <complex>
typedef complex<double> Point;
const double EPS = 1e-8;

// 入力ストリームから実数二つをとって Point へ
istream& operator >> ( istream &s, Point &a )
{
	double r, i;
	s >> r >> i;
	a = Point( r, i );
	return s;
}

// 単位ベクトル
Point unitVector( const Point &a )
{
	return a / abs( a );
}

// 法線ベクトル
Point normalVector1( const Point &a )
{
	return a * Point( 0, 1 );
}
Point normalVector2( const Point &a )
{
	return a * Point( 0, -1 );
}
Point ( *normalVectors[] )( const Point& ) = { normalVector1, normalVector2 };

int main()
{
	cin.tie( 0 );
	ios::sync_with_stdio( false );

	for ( int n; cin >> n, n; )
	{
		vector<Point> ps( n );
		REP( i, 0, n )
		{
			cin >> ps[i];
		}

		int res = 0;
		REP( i, 0, n )
		{
			REP( j, 0, i )
			{
				const Point &a = ps[i], &b = ps[j];
				REP( k, 0, 2 )
				{
					const Point c = ( a + b ) / Point( 2, 0 ) + unitVector( normalVectors[k]( a - b ) ) * sqrt( 1 - pow( abs( a - b ) / 2., 2. ) );

					int tmp = 2;
					REP( l, 0, n )
					{
						if ( l == i || l == j )
						{
							continue;
						}

						tmp += abs( c - ps[l] ) <= 1 + EPS;
					}
					res = max( res, tmp );
				}
			}
		}
		cout << res << endl;
	}

	return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;
using ld = long double;
using P = complex<ld>;
const ld eps = 1e-6;

int main()
{
	int N;
	while (cin >> N, N) {
		vector<ld> x(N), y(N);
		vector<P> p(N);
		for (int i = 0; i < N; i++) {
			cin >> x[i] >> y[i];
			p[i] = P(x[i], y[i]);
		}
		vector<P> ko;
		for (int i = 0; i < N; i++) {
			for (int j = i + 1; j < N; j++) {
				if (norm(p[i] - p[j]) < 4.0 + eps) {
					auto d = p[i] - p[j];
					ko.push_back(p[j] + (d * (ld)0.5) + d * P(0, 1) / abs(d) * sqrtl(1 - norm(d * (ld)0.5)));
					ko.push_back(p[j] + (d * (ld)0.5) - d * P(0, 1) / abs(d) * sqrtl(1 - norm(d * (ld)0.5)));
				}
			}
		}
		int res = 1;
		for (auto c : ko) {
			int cnt = 0;
			for (int i = 0; i < N; i++) {
				if (norm(c - p[i]) < 1.0001) {
					cnt++;
				}
			}
			res = max(res, cnt);
		}
		cout << res << endl;
	}
	return 0;
}