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>
using namespace std;

#define int long long
typedef vector<int>vint;
typedef pair<int,int>pint;
typedef vector<pint>vpint;
#define rep(i,n) for(int i=0;i<(n);i++)
#define reps(i,f,n) for(int i=(f);i<(n);i++)
#define all(v) (v).begin(),(v).end()
#define each(it,v) for(__typeof((v).begin()) it=(v).begin();it!=(v).end();it++)
#define pb push_back
#define fi first
#define se second
template<typename A,typename B>inline void chmin(A &a,B b){if(a>b)a=b;}
template<typename A,typename B>inline void chmax(A &a,B b){if(a<b)a=b;}

const double PI=acos(-1);
const double EPS=1e-6;
typedef complex<double>Point;
Point in(){
    double x,y;
    scanf("%lf%lf",&x,&y);
    return Point(x,y);
}

int N;
Point P[333];

signed main(){

    Point r=polar(1.0,PI/2);
    while(scanf("%lld",&N),N){
        rep(i,N)P[i]=in();

        int ma=1;
        rep(i,N)rep(j,N){
            Point b=(P[i]+P[j])*0.5;
            Point v=(P[i]-P[j])*0.5*r;
            for(double d:{-1.0,1.0}){
                Point u=v*d*sqrt(1.0-abs(v)*abs(v));
                Point p=b+u;
                int cnt=0;
                rep(k,N)if(abs(P[k]-p)<1.0+EPS)cnt++;
                chmax(ma,cnt);
            }
        }
        printf("%lld\n",ma);
    }
    return 0;
}

C++: #include <iostream>
#include <deque>
#include <string>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <string>
#include <cmath>
#include <complex>

#define loop(i,a,b) for(int i=a; i<int(b); i++)
#define rep(i,n) loop(i,0,n)

using namespace std;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef complex<double> vec;

double const pi = acos(-1);
double const eps = 1e-6;

int main(){
    int n;
    while (cin >> n && n){
        vector<vec> v(n);
        rep(i, n){
            double x, y; cin >> x >> y;
            vec p(x, y);
            v[i] = p;
        }

        if (n == 1) {
            puts("1");
            break;
        }
        if (n == 2){
            if (abs(v[0] - v[1]) < 2 - eps) {
                puts("2"); break;
            }
            else{
                puts("1"); break;
            }
        }

        int ans = 0;
        rep(i, n)rep(j, i){
            double d = abs(v[i] - v[j]);
            if (d - eps > 2) continue;
            double h = 1 - d*d / 4;
            vec c1 = polar(h, arg(v[i] - v[j]) + pi / 2) + (v[i] + v[j]) / 2.;
            vec c2 = polar(h, arg(v[i] - v[j]) - pi / 2) + (v[i] + v[j]) / 2.;
            int a1 = 0, a2 = 0;
            rep(k, n){
                if (abs(c1 - v[k]) <= 1 - eps) a1++;
                if (abs(c2 - v[k]) <= 1 - eps) a2++;
            }
            ans = max(ans, a1);
            ans = max(ans, a2);
        }
        cout << ans << endl;
    }
}

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

#define FOR(i, a, b) for (int i = (a); i < int(b); ++i)
#define REP(i, n) FOR(i, 0, n)

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

    int N;
    while (cin >> N, N) {
        static const complex<double> I = complex<double>(0, -1);
        using Point = complex<double>;
        vector<Point> points(N);
        REP(i, N) {
            double x, y; cin >> x >> y;
            points[i] = Point(x, y);
        }

        int ans = 0;
        REP(i, N) FOR(j, i+1, N) {
            auto x = points[i];
            auto y = points[j];
            double l_sq = norm(x - y);
            if (l_sq > 4) continue;
            for(double d : {-1, 1}) {
                Point middle = 0.5 * (x + y);
                Point dir = (x - y) * I;
                dir /= abs(dir);
                Point center = middle + d*sqrt(1 - l_sq/4)*dir;

                int count = 0;
                REP(k, N) {
                    if (norm(center - points[k]) <= 1 + 1e-10)
                        ++count;
                }
                ans = max(ans, count);
            }
        }
        cout << ans << '\n';
    }
}

C++: #include<iostream>
#include<iomanip>
#include<cstdio>
#include<string>
#include<cstring>
#include<climits>
#include<algorithm>
#include<vector>
#include<cmath>
#include<map>
#include<set>
#include<stack>
#include<queue>
#include<deque>
#include<tuple>
#include<list>
#include<unordered_map>
#define ll long long
#define pii pair<int,int>
#define pll pair<ll,ll>
#define vii vector<int>
#define vll vector<ll>
#define lb lower_bound
#define pb push_back
#define mp make_pair
#define rep(i,n) for(ll i=0;i<n;i++)
#define rep2(i,a,b) for(ll i=a;i<b;i++)
#define repr(i,n) for(ll i=n-1;i>=0;i--)
#define all(x) x.begin(),x.end()
// #define int ll
using namespace std;
const int INF (1 << 30);
const ll LLINF (1LL << 55LL);
const int MOD = 1000000007;
const int MAX = 510000;
const double pi = acos(-1);
const double eps = 1e-9;

int main(){
    cin.tie(0);
    ios::sync_with_stdio(false);
    cout << fixed << setprecision(15);
    while(true){
        ll n,ans=1;
        cin>>n;
        if(n==0) break;
        double x[310],y[310];
        rep(i,n){
            cin>>x[i]>>y[i];
        }
        rep(i,n){
            rep(j,n){
                if(i==j) continue;
                double dx=x[i]-x[j],dy=y[i]-y[j];
                double d=sqrt(dx*dx+dy*dy);
                if(d>2+eps) continue;
                double s=sqrt((1.0-d*d/4));
                double vx=-dy/d*s;
                double vy=dx/d*s;
                double ox=(x[i]+x[j])/2+vx;
                double oy=(y[i]+y[j])/2+vy;
                ll cnt=0;
                rep(k,n){
                    if((x[k]-ox)*(x[k]-ox)+(y[k]-oy)*(y[k]-oy)<1.0+eps)cnt++;
                }
                ans=max(ans,cnt);
            }   
        }
        cout<<ans<<endl;
    }
    return 0;
}