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 <cstdio>
#include <cassert>
#include <cstring>
#include <vector>
#include <valarray>
#include <array>
#include <queue>
#include <set>
#include <unordered_set>
#include <map>
#include <unordered_map>
#include <algorithm>
#include <cmath>
#include <complex>
#include <random>

using namespace std;
typedef long long ll;
typedef unsigned long long ull;

typedef long double R;
typedef complex<R> P;

const R EPS = 1e-10;
const R PI = acos((R)(-1));

R radNorP(R x) {
    return fmod(fmod(x, 2*PI)+2*PI, 2*PI);
}

const int MN = 330;

P p[MN];
typedef pair<R, int> Pi;

bool solve() {
    int n;
    cin >> n;
    if (!n) return false;

    for (int i = 0; i < n; i++) {
        R x, y;
        cin >> x >> y;
        p[i] = P(x, y);
    }

    int res = 1;
    for (int i = 0; i < n; i++) {
        vector<Pi> v;
        for (int j = 0; j < n; j++) {
            if (i == j) continue;
            if (2.0 < abs(p[j]-p[i])) continue;
            R th = acos(abs(p[j]-p[i])/2.0);
            R l = radNorP(arg(p[j]-p[i])-th-EPS);
            R r = radNorP(arg(p[j]-p[i])+th+EPS);
            if (r < l) {
                v.push_back(Pi(l-EPS, 1));
                v.push_back(Pi(2*PI+EPS, -1));
                v.push_back(Pi(0-EPS, 1));
                v.push_back(Pi(r+EPS, -1));
            } else {
                v.push_back(Pi(l-EPS, 1));
                v.push_back(Pi(r+EPS, -1));
            }

        }

        sort(v.begin(), v.end());
        int sm = 1;
        for (Pi p: v) {            
            sm += p.second;
            res = max(res, sm);
        }
    }
    cout << res << endl;
    return true;
}

int main() {
    while (solve()) {}
    return 0;
}

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

#define rep(i,x,y) for(int i=(x);i<(y);++i)
#define debug(x) #x << "=" << (x)

#ifdef DEBUG
#define _GLIBCXX_DEBUG
#define dump(x) std::cerr << debug(x) << " (L:" << __LINE__ << ")" << std::endl
#else
#define dump(x)
#endif

typedef long long int ll;
typedef pair<int,int> pii;
//template<typename T> using vec=std::vector<T>;

const int inf=1<<30;
const long long int infll=1LL<<58;
const double eps=1e-9;
const int dx[]={1,0,-1,0},dy[]={0,1,0,-1};

template <typename T> ostream &operator<<(ostream &os, const vector<T> &vec){
    os << "[";
    for (const auto &v : vec) {
    	os << v << ",";
    }
    os << "]";
    return os;
}

void solve(){
    while(true){
        int n;
        cin >> n;
        if(n==0) break;
        vector<long double> xs(n),ys(n);
        rep(i,0,n) cin >> xs[i] >> ys[i];

        auto count=[&](const long double x,const long double y){
            int res=0;
            rep(i,0,n) if((xs[i]-x)*(xs[i]-x)+(ys[i]-y)*(ys[i]-y)<=1+eps) ++res;
            return res;
        };

        int ans=0;
        rep(i,0,401){
            const long double y=10.*i/400;
            rep(j,0,n){
                const long double tmp=1-(ys[j]-y)*(ys[j]-y);
                if(tmp<0) continue;
                const long double x1=xs[j]-sqrtl(tmp),x2=xs[j]+sqrt(tmp);
                ans=max({ans,count(x1,y),count(x2,y)});
            }
        }
        cout << ans << endl;
    }
}

int main(){
    std::ios::sync_with_stdio(false);
    std::cin.tie(0);
    cout << fixed << setprecision(8);
    solve();
    return 0;
}

C++: #include <iostream>
#include <vector>
#include <algorithm>
#include <math.h>
using namespace std;

int main() {

  while( true ) {

    long long int n;
    cin >> n;
    if ( n == 0 ) break;

    vector< double > x, y;

    for ( long long int i = 0; i < n; i++ ) {

      double in_x, in_y;
      cin >> in_x >> in_y;
      x.push_back( in_x );
      y.push_back( in_y );

    }

    long long int ans = 1;

    for ( long long int i = 0; i < n; i++ ) {

      vector< pair< double, long long int > > v;

      for ( long long int j = 0; j < n; j++ ) {

	if ( i == j ) continue;

	double dx = x[i] - x[j];
	double dy = y[i] - y[j];
	double ang = atan2( dy, dx );
	double d = sqrt( dx * dx + dy * dy );
	if ( d > 2.0 ) continue;
	double a = acos( d / 2 );

	v.push_back( make_pair( ang - a, 0 ) );
	v.push_back( make_pair( ang + a, 1 ) );
	v.push_back( make_pair( ang - a - M_PI * 2, 0 ) );
	v.push_back( make_pair( ang + a - M_PI * 2, 1 ) );
	v.push_back( make_pair( ang - a + M_PI * 2, 0 ) );
	v.push_back( make_pair( ang + a + M_PI * 2, 1 ) );

      }

      sort( v.begin(), v.end() );

      long long int cnt = 1;

      for ( long long int j = 0; j < v.size(); j++ ) {

	if ( v[j].second == 0 ) {
	  cnt++;
	  ans = max( ans, cnt );
	}else {
	  cnt--;
	}
      }

    }

    cout << ans << endl;

  }

  return 0;

}

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

#define REP(i, n) FOR(i, 0, n)
#define FOR(i, a, b) for(int i = int(a); i < int(b); ++i)
#define RREP(i, n) RFOR(i, 0, n)
#define RFOR(i, a, b) for(int i = int(b) - 1; i >= int(a); --i)

constexpr double EPS = 1e-7;

signed main() {

	while (true) {
		int n;
		std::cin >> n;
		if (n == 0) break;

		std::vector<double> x(n), y(n);
		REP(i, n) std::cin >> x[i] >> y[i];

		int ans = 0;

		REP(i, n) FOR(j, i + 1, n) {
			double dist = std::hypot(x[i] - x[j], y[i] - y[j]);
			if (dist > 2 + EPS) continue;

			// iとjの中点m
			double mx = (x[i] + x[j]) / 2;
			double my = (y[i] + y[j]) / 2;

			// 線分ijから中心への距離
			double r = std::sqrt(1 - (dist / 2) * (dist / 2));

			// 中心c
			for (int sign = -1; sign <= 1; sign += 2) {
				// 移動のベクトルv
				double vx = sign * r * (y[j] - y[i]) / dist;
				double vy = sign * r * -(x[j] - x[i]) / dist;

				double cx = mx + vx;
				double cy = my + vy;

				int cnt = 0;
				REP(k, n) {
					if (std::hypot(cx - x[k], cy - y[k]) < 1 + EPS) ++cnt;
				}
				ans = std::max(ans, cnt);
			}
		}
		std::cout << ans << "\n";
	}

	return 0;
}