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
//------------------------------------------
#include <vector>
#include <list>
#include <map>
#include <climits>
#include <set>
#include <deque>
#include <stack>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <string>
#include <cstring>
#include <ctime>
#include <queue>
#include <random>
#include <cctype>
#include <complex>
#include <regex>

using namespace std;

#define C_MAX(a, b) ((a)>(b)?(a):(b))
#define SHOW_VECTOR(v) {std::cerr << #v << "\t:";for(const auto& xxx : v){std::cerr << xxx << " ";}std::cerr << "\n";}
#define SHOW_MAP(v) {std::cerr << #v << endl; for(const auto& xxx: v){std::cerr << xxx.first << " " << xxx.second << "\n";}}

#define EPS 1e-8
#define EQ(a, b) (abs((a) - (b)) < EPS)

typedef complex<double> P;

P getMiddle(P a, P b) {
    return (a + b) / 2.0;
}

P getUnit(P a) {
    return a / abs(a);
}

pair<P, P> getNormals(P a) {
    return make_pair((a / abs(a)) * P(0, 1), (a / abs(a)) * P(0, -1));
}

int main() {

    while (true) {

        int N;
        cin >> N;

        if (N == 0) break;

        vector<P> points(N);
        for (int i = 0; i < N; i++) {
            double x, y;
            cin >> x >> y;
            points[i] = P(x, y);
        }

        int maxV = 1;
        for (int i = 0; i < N; i++) {
            for (int j = i + 1; j < N; j++) {

                P p1 = points[i];
                P p2 = points[j];
                P middle = getMiddle(p1, p2);

                double a, b, c;
                b = abs(middle - p2);
                c = 1.0;
                if (b * b + EPS > 1.0) continue;
                a = sqrt(1.0 - b * b);

                P o1 = middle + ((p2 - p1) / abs(p2 - p1)) * P(0, 1) * a;
                P o2 = middle + ((p2 - p1) / abs(p2 - p1)) * P(0, -1) * a;

                int count = 0;
                for (int k = 0; k < N; k++) {
                    double dx = o1.real() - points[k].real();
                    double dy = o1.imag() - points[k].imag();
                    if (dx * dx + dy * dy < 1.0 + EPS) count++;
                }
                maxV = max(maxV, count);

                count = 0;
                for (int k = 0; k < N; k++) {
                    double dx = o2.real() - points[k].real();
                    double dy = o2.imag() - points[k].imag();
                    if (dx * dx + dy * dy < 1.0 + EPS) count++;
                }
                maxV = max(maxV, count);
            }
        }

        cout << maxV << endl;

    }

    return 0;
}

C++: #include <bits/stdc++.h>
#define For(i, a, b) for(int (i)=(a); (i)<(b); ++(i))
#define rFor(i, a, b) for(int (i)=(a)-1; (i)>=(b); --(i))
#define rep(i, n) For((i), 0, (n))
#define rrep(i, n) rFor((i), (n), 0)
#define fi first
#define se second
using namespace std;
typedef long long lint;
typedef pair<int, int> pii;
typedef pair<lint, int> pli;
typedef complex<double> xy_t;

const lint mod = 1e9 + 7;

int n;
xy_t p[310];
const double pi = acos(-1.0);

bool in_circle(xy_t t, xy_t c){
    return abs(t - c) <= 1;
}

int cnt(int k, double t){
    xy_t c = p[k] + xy_t(cos(t), sin(t));
    int ret = 1;
    rep(i, n)if(i != k){
        ret += in_circle(p[i], c);
    }
    return ret;
}

int main(){
    while(scanf("%d", &n) && n){
        rep(i, n){
            double x, y;
            scanf("%lf%lf", &x, &y);
            p[i] = xy_t(x, y);
        }
        int ans = 0;
        rep(i, n)rep(d, 500){
            ans = max(ans, cnt(i, 2*pi*d/500));
        }
        printf("%d\n", ans);
    }
}

C++: #include <iostream>
#include <iomanip>
#include <sstream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <string>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <map>
#include <algorithm>
#include <functional>
#include <utility>

using namespace std;

typedef unsigned int ui;
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;

const int INF = INT_MAX / 4;

#define REP( i, m, n ) for ( int i = (int)( m ); i < (int)( n ); i++ )
#define UREP( i, m, n ) for ( ui i = (ui)( m ); i < (ui)( n ); i++ )
#define EACH( v, c ) for ( __typeof( c.front() ) &v : c )

#define ITER( c ) __typeof( (c).begin() )
#define IREP( c, it ) for ( ITER(c) it = c.begin(); it != c.end(); ++it )
#define ALL( c ) (c).begin(), (c).end()
#define RALL( c ) (c).rbegin(), (c).rend()

#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

double powdist( double x1, double y1, double x2, double y2 )
{
	return pow( x1 - x2, 2 ) + pow( y1 - y2, 2 );
}

double dist( double x1, double y1, double x2, double y2 )
{
	return sqrt( pow( x1 - x2, 2 ) + pow( y1 - y2, 2 ) );
}

int main()
{
	while ( true )
	{
		int n;
		cin >> n;

		if ( n == 0 )
		{
			break;
		}

		vector<double> x( n ), y( n );

		REP( i, 0, n )
		{
			cin >> x[i] >> y[i];
		}

		int res = 1;

		REP( i, 0, n )
		{
			REP( j, i + 1, n )
			{
				if ( powdist( x[i], y[i], x[j], y[j] ) <= 4 )
				{
					double mx = ( x[i] + x[j] ) / 2, my = ( y[i] + y[j] ) / 2;
					double xl = sqrt( 1 - pow( dist( x[i], y[i], x[j], y[j] ) / 2, 2 ) );
					
					double ox1 = mx, oy1 = my, ox2 = mx, oy2 = my;

					if ( x[i] - x[j] == 0 )
					{
						oy1 += xl;
						oy2 -= xl;
					}
					else if ( y[i] - y[j] == 0 )
					{
						ox1 += xl;
						oy2 -= xl;
					}
					else
					{
						double s = -( ( x[j] - x[i] ) / ( y[j] - y[i] ) );
						double th = atan( s );
						
						ox1 = mx + xl * cos( th );
						oy1 = my + xl * sin( th );

						ox2 = mx - xl * cos( th );
						oy2 = my - xl * sin( th );
					}

					int tmp1 = 2, tmp2 = 2;

					REP( k, 0, n )
					{
						if ( i == k || j == k )
						{
							continue;
						}

						if ( powdist( x[k], y[k], ox1, oy1 ) <= 1 )
						{
							tmp1++;
						}

						if ( powdist( x[k], y[k], ox2, oy2 ) <= 1 )
						{
							tmp2++;
						}
					}

					res = max( res, max( tmp1, tmp2 ) );
				}
			}
		}

		cout << res << endl;
	}

	return 0;
}