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++: 
//Name: Circle and Points
//Level: 3
//Category: 幾何,円
//Note:

/*
 * 半径が固定なので，円周上に2点が乗っているときが境界状態．
 * すべての2点の組み合わせについて中心の来る場所（2点）を求め，内包判定を行う．
 * N=1のケースと，全ての点が互いに距離2以上離れているケースに注意．
 *
 * オーダーは O(N^3)．
 */
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <complex>
#include <utility>

using namespace std;

typedef complex<double> Point;

bool inCircle(const Point &p, const Point &center, double radius) {
    double dist2 = norm(p-center);
    double rad2 = radius*radius;
    return dist2 < rad2+1e-9;
}

int check(const Point &center, const vector<Point> &points) {
    const int N = points.size();
    int cnt = 0;
    for(int i = 0; i < N; ++i) {
        if(inCircle(points[i], center, 1)) ++cnt;
    }
    return cnt;
}

int test(int a, int b, const vector<Point> &points) {
    const Point v(points[a]-points[b]);
    if(abs(v) >= 2) return 1;
    const double hlen = abs(v)/2;

    const double plen = sqrt(1 - hlen*hlen);
    Point perp = v*Point(0, 1);
    perp /= abs(perp);
    const Point mid = points[b] + v/2.0;
    return max(check(mid + perp*plen, points), check(mid - perp*plen, points));
}

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

    while(true) {
        int N;
        cin >> N;
        if(!N) break;

        vector<Point> points(N);
        for(int i = 0; i < N; ++i) {
            cin >> points[i].real() >> points[i].imag();
        }
        int ans = 0;
        for(int i = 0; i < N; ++i) {
            for(int j = 0; j < i; ++j) {
                ans = max(ans, test(i, j, points));
            }
        }
        if(N == 1) ans = 1;
        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 = 1;
		REP( i, 0, n )
		{
			REP( j, 0, i )
			{
				const Point &a = ps[i], &b = ps[j];
				if ( 2 + EPS < abs( a - b ) )
				{
					continue;
				}
				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 <iostream>
#include <vector>
#include <string>
#include <map>
#include <set>
#include <queue>
#include <list>
#include <stack>
#include <algorithm>
#include <functional>
#include <numeric>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <climits>
#include <sstream>
#include <complex>
using namespace std;

#define REP(i,a,n) for(i=a; i<n; i++)
#define rep(i,n) REP(i,0,n)
#define all(x) x.begin(), x.end()
#define rall(x) x.rbegin(), x.rend()
#define foreach(it,x) for(typeof(x.begin()) it=x.begin(); it!=x.end(); it++)

typedef complex<double> Point;
Point p[300];

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

bool isok(Point& e, Point& q){
    return (sqr(e.real() - q.real()) + sqr(e.imag() - q.imag()) < 1.0);
}

/*
 Point(0,1)@¶É90xñ]
 Point(0,-1)@EÉ90xñ]
 */
int main(){
    int n,i,j,k, res;

    while( scanf("%d",&n), n ){
        res = 1;
        rep(i,n) cin >> p[i].real() >> p[i].imag();

        rep(i,n){
            REP(j,i+1,n){
                double dist = sqr(p[i].real() - p[j].real()) + sqr(p[i].imag() - p[j].imag());
                if( dist < 4.0 ){
                    dist /= 4.0; // (d/2)^2
                    double x = sqrt(1.0 - dist); // x
                    Point M((p[i].real() + p[j].real())/2.0, (p[i].imag() + p[j].imag())/2.0);

                    Point d = (p[i] - p[j])/abs(p[i] - p[j]);
                    Point e1 = M + d*Point(0,-1)*x;
                    Point e2 = M + d*Point(0,1)*x;

                    int cnt1 = 2, cnt2 = 2;

                    rep(k,n) if( i != k && j != k ){ // e1
                        if( isok(e1, p[k]) ) cnt1++;
                        if( isok(e2, p[k]) ) cnt2++;
                    }
                    res = max(res, max(cnt1,cnt2));
                }
            }
        }

        printf("%d\n",res);
    }
    return 0;
}