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 <cstdio>
#include <string>
#include <cstring>
#include <deque>
#include <list>
#include <queue>
#include <stack>
#include <vector>
#include <utility>
#include <algorithm>
#include <map>
#include <set>
#include <complex>
#include <cmath>
#include <limits>
#include <cfloat>
#include <climits>
#include <ctime>
#include <cassert>
#include <numeric>
#include <fstream>
#include <functional>
#include <bitset>

using namespace std;

#define INF 1e10
#define EPS (1e-10)    //誤差
#define PI 3.1415926535897932384626433832795028841971
typedef complex<double> P;
 
struct star{
    P point[5];    
};
 
//実数同士の比較
bool EQ(double a, double b){
    return (abs(a - b) < EPS);
}
 
//二つのベクトルが等しいかどうか
bool EQV(P a, P b){
    return (EQ(a.real(), b.real()) && EQ(a.imag(), b.imag()));
}
 
//内積
double dot(P a, P b){
    return (a.real() * b.real() + a.imag() * b.imag());
}
 
//外積
double cross(P a, P b){
    return (a.real() * b.imag() - a.imag() * b.real());
}
 
// 点a,bを通る直線と点cとの距離
double distance_l_p(P a, P b, P c) {
  return abs(cross(b-a, c-a)) / abs(b-a);
}
 
// 点cが線分a,b上にあるかないか(3)
int is_point_on_line(P a, P b, P c) {
  // |a-c| + |c-b| <= |a-b| なら線分上
  return (abs(a-c) + abs(c-b) < abs(a-b) + EPS);
}
 
// 点a,bを端点とする線分と点cとの距離
double distance_ls_p(P a, P b, P c) {
  if ( dot(b-a, c-a) < EPS ) return abs(c-a);
  if ( dot(a-b, c-b) < EPS ) return abs(c-b);
  return abs(cross(b-a, c-a)) / abs(b-a);
}

//図を思い出す
int ccw(P a, P b, P c) {
    b -= a; c -= a;
    if( cross(b,c) > EPS ) return +1;
    if( cross(b,c) < -EPS ) return -1;
    if( dot(b,c) < 0 ) return +2;
    if( norm(b) < norm(c) ) return -2;
    return 0;
}
 
// a1,a2を端点とする線分とb1,b2を端点とする線分の交差判定
int is_intersected_ls(P a1, P a2, P b1, P b2) {
  return ( ccw(a1,a2,b1) * ccw(a1,a2,b2) <= 0 ) &&
           ( ccw(b1,b2,a1) * ccw(b1,b2,a2) <= 0 );
}
 
// a1,a2を通る直線とb1,b2を通る直線の交点計算
P intersection_l(P a1, P a2, P b1, P b2) {
  P a = a2 - a1; P b = b2 - b1;
  return a1 + a * cross(b, b1-a1) / cross(b, a);
}


// 点a1,a2を短点とする線分と、点b1,b2を端点とする線分との距離
double distance_ls_ls(P a1, P a2, P b1, P b2){
    if(is_intersected_ls(a1, a2, b1, b2)) return 0.0;
    double ret = INF;
    ret = min(ret, distance_ls_p(a1, a2, b1));
    ret = min(ret, distance_ls_p(a1, a2, b2));
    ret = min(ret, distance_ls_p(b1, b2, a1));
    ret = min(ret, distance_ls_p(b1, b2, a2));
    return ret;
}
 


signed main(){
    
    while(1){

        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 ans = 0;
        for(int i = 0; i < n; i++){
            for(int j = i + 1; j < n; j++){
                P a = points[i];
                P b = points[j];
                if(abs(a - b) >= 2.0 + EPS) continue;
                P m = (a + b) * 0.5;
                P norm = (b - a) * P(0, 1);
                norm = norm / abs(norm);
                double l = abs(b - a) / 2;
                double dist = sqrt(1.0 - l * l);
                P O = m + norm * dist;
                
                int tmp = 0;
                for(int k = 0; k < n; k++){
                    if(abs(O - points[k]) <= 1.0 + EPS) tmp++;
                }

                ans = max(tmp, ans);
            }
        }

        cout << max(ans, 1) << endl;
    }
    return 0;
}

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-8;

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+eps) t++;
			chmax(res,t);
			t = 0;
			rep(k,n) if (sqrt((xj-x[k])*(xj-x[k]) + (yj-y[k])*(yj-y[k]))<=1+eps) t++;
			chmax(res,t);
		}
		cout << res << ln;
	}
}