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 <cmath>
#include <cstdint>
#include <vector>
#include <array>

using namespace std;

typedef array<double, 2> Point;


uint32_t specialSolve(const vector<Point> v, double centerX, double centerY)
{
    uint32_t cnt = 0;
    for(const auto& p: v)
    {
        const auto r = sqrt(pow(p[0] - centerX, 2) + pow(p[1] - centerY, 2));
        if(r <= 1 + 1.0E-6)
            ++cnt;
    }

    return cnt;
}


uint32_t solve(const vector<Point>& ps)
{
    const auto N = ps.size();

    uint32_t maxCount = 0;
    for(uint32_t i = 0; i < N; ++i)
    {
        auto const x1 = ps[i][0],
                   y1 = ps[i][1];

        for(uint32_t j = i + 1; j < N; ++j)
        {
            auto const x2 = ps[j][0],
                       y2 = ps[j][1];

            auto const a = -(x2 - x1)/(y2 - y1);
            auto const a2 = pow(a, 2.0);
            auto const ds2_4 = pow(x1 - x2, 2) + pow(y1 - y2, 2);
            auto const dt2 = 1 - ds2_4/4;

            // 2点は同じ半径1の円に入ることはない
            if(dt2 < 0)
                continue;

            auto const dx = sqrt(dt2 / (1 + a2));
            auto const dy = dx * a;
            auto const avgPosX = (x2 + x1) / 2,
                       avgPosY = (y2 + y1) / 2;
            auto const m1 = specialSolve(ps, avgPosX + dx, avgPosY + dy),
                       m2 = specialSolve(ps, avgPosX - dx, avgPosY - dy);

            if(maxCount < m1)
                maxCount = m1;
            if(maxCount < m2)
                maxCount = m2;
        }
    }


    return maxCount;
}



int main()
{
    while(1)
    {
        uint32_t N;
        cin >> N;

        if(N == 0)
            break;

        vector<Point> ps;
        for(uint32_t i = 0; i < N; ++i){
            double a, b;
            cin >> a >> b;
            Point p = {a, b};
            ps.push_back(p);
        }


        auto ans = solve(ps);
        cout << ans << endl;
    }

    return 0;
}

C++: #include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <cstring>
#include <map>
#include <utility>
#include <set>
#include <bitset>
#include <iostream>
#include <iomanip>
#include <memory>
#include <string>
#include <vector>
#include <list>
#include <numeric>
#include <algorithm>
#include <functional>
#include <sstream>
#include <complex>
#include <stack>
#include <queue>
#include <deque>
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
const int INF = 1e9;
const ll LINF = 1e16;

/*
<url:>
問題文============================================================

=================================================================

解説=============================================================

================================================================
*/
typedef complex<double> Point;
const double eps = 1e-9;
Point P[305];
int main(void) {
	cin.tie(0); ios::sync_with_stdio(false);
    int N;
    while(true){
        scanf("%d",&N);
        if(N == 0) break;
        for(int i = 0; i < N;i++){
            double x,y; scanf("%lf%lf",&x,&y);
            P[i] = Point(x,y);
        }
        if(N == 1){
            printf("1\n");
            continue;
        }
        int ans = 1;
        for(int i = 0; i < N;i++){
            for(int j = i + 1; j < N;j++){
                double dist = abs(P[i] - P[j]);
                if(dist > 2) continue;
                double r = sqrt(1 - dist*dist/4);
                double dx = -(P[i] - P[j]).imag()*r/dist;
                double dy = (P[i] - P[j]).real()*r/dist;
                Point np(dx,dy);
                Point cp(np + (P[i] + P[j])/2.0);
                
                int cnt = 0;
                for(int k = 0; k < N;k++){
                    if(abs(cp - P[k]) < 1.0 + eps) cnt++;
                }
                ans = max(ans,cnt);
            }
        }
        printf("%d\n",ans);
    }
	return 0;
}

C++: #include <bits/stdc++.h>

using namespace std;

#define X first
#define Y second

using P=pair<double, double>;

const double EPS = 1e-14;

double dist(P a, P b) {
	double dx = a.X - b.X;
	double dy = a.Y - b.Y;
	return sqrt(dx*dx+dy*dy);
}

int main() {
	int n;
	while(cin >> n, n) {
		vector<P> p(n), pl;
		for(int i = 0; i < n; i++) cin >> p[i].X >> p[i].Y;
		for(int i = 0; i < p.size(); i++) {
			P p1 = p[i];
			for(int j = 0; j < p.size(); j++) {
				if(i==j) continue;
				P p2 = p[j];
				if(dist(p1, p2) > 2) continue;

				double dxy = dist(p1, p2);
				double l = sqrt(1-(dxy/2));
				double cx = (p1.X+p2.X)/2;
				double cy = (p1.Y+p2.Y)/2;
				double dx = (p2.X - p1.X);
				double dy = (p2.Y - p1.Y);

				double dlx = l/sqrt(dx*dx+dy*dy) * -dy;
				double dly = l/sqrt(dx*dx+dy*dy) * dx;

				pl.push_back(P(cx+dlx, cy+dly));
				pl.push_back(P(cx-dlx, cy-dly));
			}
		}

		int res = 0;
		for(P p1:pl) {
			int tmp = 0;
			for(P p2:p) {
				if(dist(p1, p2) <= 1+EPS) tmp++;
			}
			res = max(res, tmp);
		}

		cout << res << endl;
	}

	return 0;
}

C++: #include <algorithm>
#include <cassert>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <deque>
#include <functional>
#include <iomanip>
#include <iostream>
#include <list>
#include <map>
#include <queue>
#include <set>
#include <sstream>
#include <string>
#include <utility>
#include <vector>

#define REP(i,s,n) for(int i=(int)(s);i<(int)(n);i++)
#define DEBUGP(val) cerr << #val << "=" << val << "\n"

using namespace std;
typedef long long int ll;
typedef vector<int> VI;
typedef vector<ll> VL;
typedef pair<int, int> PI;
typedef pair<double, double> PD;

const double EPS = 1e-11;
double sq(double x) { return x * x; }

vector<PD> solve(double ax, double ay, double bx, double by) {
  double sqdist = sq(ax - bx) + sq(ay - by);
  double dist = sqrt(sqdist);
  double ndist = sqrt(1 - sqdist / 4);
  double ex = (bx - ax) / dist, ey = (by - ay) / dist;
  vector<PD> ans(2);
  REP(i, 0, 2) {
    double coef = i == 0 ? 1 : -1;
    ans[i].first = (ax + bx) / 2 + ey * coef * ndist;
    ans[i].second = (ay + by) / 2 - ex * coef * ndist;
  }
  return ans;
}

int main(void) {
  ios::sync_with_stdio(false);
  cin.tie(0);
  int n;
  while(cin >> n && n) {
    vector<double> x(n), y(n);
    REP(i, 0, n) cin >> x[i] >> y[i];
    vector<PD> cand;
    REP(i, 0, n) {
      REP(j, 0, i) {
	for (PD e: solve(x[i], y[i], x[j], y[j])) {
	  cand.push_back(e);
	}
      }
    }
    int ma = 0;
    for (PD c: cand) {
      int cnt = 0;
      REP(i, 0, n) {
	double sqdist = sq(x[i] - c.first) + sq(y[i] - c.second);
	if (sqdist <= 1 + EPS) {
	  cnt += 1;
	}
      }
      ma = max(ma, cnt);
    }
    cout << ma << "\n";
  }
}