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 rep(i,n) for(long long i = 0; i < (long long)(n); i++)
#define pb push_back
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define mt make_tuple
template<class T1, class T2> bool chmin(T1 &a, T2 b) { return b < a && (a = b, true); }
template<class T1, class T2> bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }

using ll = long long; using ld = long double; using vll = vector<ll>; using vvll = vector<vll>; using vld = vector<ld>; 
using vi = vector<int>; using vvi = vector<vi>;
using P = pair<ll, ll>;

template <typename T, typename U> ostream &operator<<(ostream &o, const pair<T, U> &v) {  o << "(" << v.first << ", " << v.second << ")"; return o; }
template<size_t...> struct seq{}; template<size_t N, size_t... Is> struct gen_seq : gen_seq<N-1, N-1, Is...>{}; template<size_t... Is> struct gen_seq<0, Is...> : seq<Is...>{};
template<class Ch, class Tr, class Tuple, size_t... Is>
void print_tuple(basic_ostream<Ch,Tr>& os, Tuple const& t, seq<Is...>){ using s = int[]; (void)s{0, (void(os << (Is == 0? "" : ", ") << get<Is>(t)), 0)...}; }
template<class Ch, class Tr, class... Args> 
auto operator<<(basic_ostream<Ch, Tr>& os, tuple<Args...> const& t) -> basic_ostream<Ch, Tr>& { os << "("; print_tuple(os, t, gen_seq<sizeof...(Args)>()); return os << ")"; }
ostream &operator<<(ostream &o, const vvll &v) { rep(i, v.size()) { rep(j, v[i].size()) o << v[i][j] << " "; cout << endl; } return o; }
template <typename T> ostream &operator<<(ostream &o, const vector<T> &v) { o << '['; rep(i, v.size()) o << v[i] << (i != v.size()-1 ? ", " : ""); o << "]";  return o; }
template <typename T>  ostream &operator<<(ostream &o, const set<T> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]";  return o; }
template <typename T, typename U>  ostream &operator<<(ostream &o, const map<T, U> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it << (next(it) != m.end() ? ", " : ""); o << "]";  return o; }
template <typename T, typename U>  ostream &operator<<(ostream &o, const unordered_map<T, U> &m) { o << '['; for (auto it = m.begin(); it != m.end(); it++) o << *it; o << "]";  return o; }
#define VN(v) # v
#define print(a) cout << a << "#" << VN(a) << endl;
#define ldout fixed << setprecision(40) 

static const double EPS = 1e-14;
static const long long INF = 1e18;
static const long long mo = 1e9+7;

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

    ll n; 
    while (cin >> n && n) {
        vld x(n), y(n); rep(i, n) cin >> x[i] >> y[i];
        auto d = [&](ll i, ll j) { return (x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j]); };
        ll ret = 1;
        for (int i = 0; i < n-1; i++) {
            for (int j = i+1; j < n; j++) {
                ld dij = d(i, j);
                if (dij >= 4) continue;
                ld x1 = x[i], y1 = y[i], x0 = x[j], y0 = y[j];
                for (auto dir : {1, -1}) {
                    ld xc = x0 + (x1 - x0) / 2.0l + (-(y1 - y0)) * dir 
                        * sqrtl(1.0l - ((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) / 4.0l) 
                        / sqrtl((y1 - y0) * (y1 - y0)  + (x1 - x0) * (x1 - x0));
                    ld yc = y0 + (y1 - y0) / 2.0l + (+(x1 - x0)) * dir 
                        * sqrtl(1.0l - ((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)) / 4.0l) 
                        / sqrtl((y1 - y0) * (y1 - y0)  + (x1 - x0) * (x1 - x0));
//                    cout << x1 << " " << y1 << " " << x0 << " " << y0 << " : " << xc << " " << yc << endl;
                    ll m = 2;
                    rep(k, n) {
                        if (k == i || k == j) continue;
                        if ((xc - x[k]) * (xc - x[k]) + (yc - y[k]) * (yc - y[k]) < 1) 
                            m++;
                    }
                    chmax(ret, m);
                }
            }
        }
        cout << ret << endl;
    }
    return 0;
}

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

double pi = 3.14159265358979;

signed main(){
    ios::sync_with_stdio(false);
	cin.tie(0);
    cout << fixed << setprecision(20);

    int n;
    double eps = 1e-6;
    double b = 2*pi/3000;
    while(cin>>n){
        if(n==0) break;
        double x[n],y[n];
        for(int i=0;i<n;i++){
            cin>>x[i]>>y[i];
        }
        int ans = 0;
        for(int i=0;i<n;i++){
            int ret = 0;
            for(int j=0;j<n;j++){
                double d = (x[i] - x[j]) * (x[i] - x[j]);
               d += (y[i] - y[j]) * (y[i]-y[j]);
                if(d <= 1.0) ret++;
            }
            ans = max(ans,ret);
        }
        //cerr << " a " << ans << endl;
        for(int i=0;i<n;i++){
            for(int cnt=0;cnt<3000;cnt++){
                double theta = b*cnt;
                double ox = x[i] + cos(theta);
                double oy = y[i] + sin(theta);
                int ret = 0;
                for(int j=0;j<n;j++){
                    double d = (ox - x[j]) * (ox - x[j]);
                    d += (oy - y[j]) * (oy - y[j]);
                    if(d <= 1.0 + eps) ret++;
                }
                ans = max(ans,ret);
            }
        }
        cout << ans << "\n";
    }
}

C++: #include <iostream>
#include <algorithm>
#include <cmath>
#include <complex>

using namespace std;

typedef complex<double> Point;

const double EPS = 1e-8;

int n;
Point data[300];

int calc(Point a, Point b){
  int res = 0;
  double d = abs(b - a) / 2.0;
  Point ab = (b - a) / (2.0 * d);
  if(d < 1){
    Point o = ab * d;
    d = sqrt(1.0 - d * d);
    o += Point(d * ab.imag() + a.real(), -d * ab.real() + a.imag());
    for(int i=0;i<n;i++){
      if(abs(o - data[i]) < 1.0 + EPS) res++;
    }
  }
  return res;
}

main(){
  while(cin >> n && n){
    for(int i=0;i<n;i++){
      cin >> data[i].imag() >> data[i].real();
    }
    int ans = 1;
    for(int i=0;i<n;i++){
      for(int j=0;j<n;j++){
        ans = max(ans, calc(data[i], data[j]));
        ans = max(ans, calc(data[j], data[i]));
      }
    }
    cout << ans << endl;
  }
}