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 <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>
#define REP(i,n) for(int i=0; i<(int)(n); ++i)

using namespace std;

#define X real()
#define Y imag()

typedef double D;
D EPS = 1e-8;

typedef complex<D> P;
struct L{
    P a, b;
    P vec(){ return b - a; }
};
typedef vector<P> Pol;

int sig(D a) { return a < -EPS ? -1 : a > EPS ? 1 : 0; }
D dot(P a, P b){ return a.X * b.X + a.Y * b.Y; }
D cross(P a, P b){ return a.X * b.Y - a.Y * b.X; }

D area(Pol pol) {
    int n = pol.size();
    D res = 0;
    REP(i, n) {
        res += cross(pol[i], pol[ (i + 1) % n ]);
    }
    return abs(res) / 2.0;
}

// 2点a, bを通る半径rの円の中心
vector<P> tc2(P a, P b, D r){
    vector<P> res;

    D d = abs(b - a);
    if(d > 2 * r) return res;

    P mid = 0.5 * (a + b);
    P dir = polar(sqrt(r * r - d * d / 4), arg(b - a) + M_PI / 2);
    res.push_back(mid + dir);
    res.push_back(mid - dir);
    return res;
}

int main(){
    int N;
    while(cin >> N && N > 0) {
        vector<P> ps(N);
        REP(i, N) {
            D x, y;
            cin >> x >> y;
            ps[i] = P(x, y);
        }

        int ans = 0;
        REP(i, N) REP(j, N) if(i != j) {
            vector<P> cv = tc2(ps[i], ps[j], 1.0);
            for(auto c : cv) {
                int sum = 0;
                REP(k, N) {
                    if(sig(abs(c - ps[k]) - 1.0) <= 0) {
                        sum++;
                    }
                }
                ans = max(ans, sum);
            }
        }
        cout << ans << endl;
    }
    return 0;
}

C++: #include <stdio.h>
#include <algorithm>
#include <iostream>
#include <stack>
#include <vector>
#include <math.h>
#include <queue>

using namespace std;

struct Point{
    double x,y;
};

void func(int N){

    Point point[N];
    double x_min = 100,y_min = 100,x_max = -1,y_max = -1;

    for(int i = 0; i < N; i++){
        scanf("%lf %lf",&point[i].x,&point[i].y);
        x_min = min(point[i].x,x_min);
        y_min = min(point[i].y,y_min);
        x_max = max(point[i].x,x_max);
        y_max = max(point[i].y,y_max);
    }

    int count,max_count = 0;

    for(double y = y_min; y <= y_max; y += 0.01){
        for(double x = x_min; x <= x_max; x += 0.01){
            count = 0;
            for(int i = 0; i < N; i++){
                if(sqrt((x-point[i].x)*(x-point[i].x)+(y-point[i].y)*(y-point[i].y)) <= 1.0)count++;
            }
            max_count = max(max_count,count);
        }
    }

    printf("%d\n",max_count);
}


int main(){

    int N;

    while(true){
        scanf("%d",&N);
        if(N == 0)break;

        func(N);
    }


    return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;

#define REP(i,n) for(int i=0;i<(int)n;++i)
#define FOR(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();++i)
#define ALL(c) (c).begin(), (c).end()

int n;
double x[310], y[310];
double eps = 1e-9;

double pow(double a) {
    return a * a;
}

double dist2(double ax, double ay, double bx, double by) {
    return pow(ax - bx) + pow(ay - by);
}

double dist2(int i, int j) {
    return dist2(x[i], y[i], x[j], y[j]);
}

int count(double cx, double cy) {
    int ret = 0;
    for (int i = 0; i < n; ++i) {
        if (dist2(x[i], y[i], cx, cy) < 1 + eps)
            ++ret;
    }
    return ret;
}

int main(){
    while (true) {
        cin >> n;
        if (!n) return 0;
        for (int i = 0; i < n; ++i) {
            cin >> x[i] >> y[i];
        }
        int ans = 1;
        for (int i = 0; i < n; ++i) {
            for (int j = i + 1; j < n; ++j) {
                double d2 = dist2(i, j);
                double d = sqrt(d2);
                if (d > 2)
                    continue;
                
                double dis2 = 1 - d2 / 4;
                double dis = sqrt(dis2);
                double cx = (x[i] + x[j]) / 2;
                double cy = (y[i] + y[j]) / 2;
                double dx = (y[i] - y[j]) / d * dis;
                double dy = (x[i] - x[j]) / d * dis;
                ans = max(ans, max(count(cx + dx, cy - dy), count(cx - dx, cy + dy)));
            }
        }
        cout << ans << "\n";
    }
}