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++: #if 0
かかった時間：約1時間半
最初はわかりませんが、aojで他の人のソースコードを参考して、
先生のヒント「2点を通る円だけを考えれば十分である」がわかりました。
簡単に作成し、でも何度もwrong answerでした。
困っている時に、先生のヒントの「例外ケースにも注意」を思い出しました。
そして点が一つしかない場合にプログラムが何も処理していないと発見しました。
初条件を変えればacceptedを得ました。
#endif
#include <iostream>
#include <complex>
#include <cmath>
#define EPS 1e-9
using namespace std;
typedef complex<double> xy_t;
int n;
xy_t points[301];
int count_in(xy_t center)
{
  int ans=0;
  for(int i=0;i<n;i++)
  {
    if (abs(center-points[i])<1.0+EPS)
      ans++;
  }
  return ans;
}

xy_t ppdcl(xy_t a)
{
  return xy_t(a.imag(),-a.real());
}
xy_t delta_v(xy_t C1,xy_t C2)
{
  xy_t half=(C1-C2)*0.5;
  double l=abs(half);
  return ppdcl(half)*sqrt(1-l*l)*(1.0/l);
}

int main()
{
  while(1)
  {
    cin>>n;
    if (n==0) break;
    double x,y;
    for(int i=0;i<n;i++)
    {
      cin>>x>>y;
      points[i]=xy_t(x,y);
    }
    int maxn=1;
    for(int i=0;i<n;i++)
    {
      for(int j=i;j<n;j++)
      {
        if(abs(points[i]-points[j])<2.0+EPS)
        {
          xy_t mid_p=(points[i]+points[j])*0.5;
          maxn=max(maxn,count_in(mid_p+delta_v(points[i],points[j])));
          maxn=max(maxn,count_in(mid_p-delta_v(points[i],points[j])));
        }
      }
    }
    cout<<maxn<<endl;
  }

}

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 = 1;

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

C++: #include<iostream>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<queue>
#include<algorithm>
#include<cstring>
#include<string>
#include<vector>
#define mp make_pair
#define pb push_back
#define rep(i,n) for(int i = 0;i < n;i++)
using namespace std;

double dist(double x1,double y1,double x2,double y2){
    return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}

double dist(double x1,double y1){
    return dist(x1,y1,0,0);
}

void calccenter(double c[][2],double u1[],double u2[]){
    double v[2]={u2[0]-u1[0],u2[1]-u1[1]};
    rep(i,2) {
        rep(j,2){
            c[i][j]=(u1[j]+u2[j])/2.0;
        }
    }
    double lenv = dist(v[0],v[1]);
    double lent = sqrt(1-pow((lenv/2.0),2));
    rep(i,2){
        double dg = i?1.0:-1.0;
        c[i][0]+=dg*(v[1]*lent/lenv);
        c[i][1]-=dg*(v[0]*lent/lenv);
    }
}

int main(){
    int N;
    while(1){
        cin>>N;
        if(!N)break;
        double a[300][2];
        rep(i,N){
            cin>>a[i][0]>>a[i][1];
        }
        int ans = 1;
        int temp = 0;
        if(N==1){
            cout<<1<<endl;
            continue;
        }
        rep(u1,N){
            for(int u2 = u1+1;u2<N;u2++){
                double c[2][2];
                if(dist(a[u1][0],a[u1][1],a[u2][0],a[u2][1])>2){
                    continue;
                }
                calccenter(c,a[u1],a[u2]);
                rep(k,2){
                    temp = 0;
                    rep(l,N){
                        if(l==u1||l==u2){
                            temp++;
                        }else if(dist(c[k][0],c[k][1],a[l][0],a[l][1])<=1.0){
                            temp++;
                        }
                    }
                    ans = max(ans,temp);
                }
            }
        }
        cout<<ans<<endl;
    }
}

C++: #include <bits/stdc++.h>
using namespace std;
using ld = long double;
using P = complex<ld>;
const ld eps = 1e-8;

int main()
{
	int N;
	while (cin >> N, N) {
		vector<ld> x(N), y(N);
		vector<P> p(N);
		for (int i = 0; i < N; i++) {
			cin >> x[i] >> y[i];
			p[i] = P(x[i], y[i]);
		}
		vector<P> ko;
		for (int i = 0; i < N; i++) {
			for (int j = i + 1; j < N; j++) {
				if (norm(p[i] - p[j]) < 4.0 + eps) {
					auto d = p[i] - p[j];
					ko.push_back(p[j] + (d * (ld)0.5) + d * P(0, 1) / abs(d) * sqrtl(1 - norm(d * (ld)0.5)));
					ko.push_back(p[j] + (d * (ld)0.5) - d * P(0, 1) / abs(d) * sqrtl(1 - norm(d * (ld)0.5)));
				}
			}
		}
		int res = 0;
		for (auto c : ko) {
			int cnt = 0;
			for (int i = 0; i < N; i++) {
				if (norm(c - p[i]) < 1.0001 + eps) {
					cnt++;
				}
			}
			res = max(res, cnt);
		}
		cout << res << endl;
	}
	return 0;
}

C++: #include <iostream> 
#include <string> 
#include <vector> 
#include <cmath> 
#include <algorithm> 
#include <cstdlib> 
#include <ctime> 
#include <cstdio> 
#include <functional> 
#include <set> 
#include <sstream> 


using namespace std; 

int main(){

	int n;
	const double eps=1e-8;
	while(true){

		cin>>n;

		if(n==0) return 0;

		int ans=0;
		vector< double > x(n),y(n);

		for(int i=0;i<n;i++) cin>>x[i]>>y[i];


		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				const double d=hypot(x[i]-x[j],y[i]-y[j]);
				if(d>2-eps) continue;

				const double mx=(x[i]+x[j])/2.0;
				const double my=(y[i]+y[j])/2.0;
				const double dvx=(x[i]-x[j])/d;
				const double dvy=(y[i]-y[j])/d;

				const double cx1=mx+dvy*sqrt(1.0-d*d/4.0);
				const double cy1=my-dvx*sqrt(1.0-d*d/4.0);

				const double cx2=mx-dvy*sqrt(1.0-d*d/4.0);
				const double cy2=my+dvx*sqrt(1.0-d*d/4.0);

				int res=2;
				for(int k=0;k<n;k++){
					if(k==i || k==j) continue;
					if(hypot(x[k]-cx1,y[k]-cy1)<1.0+eps) res++;
				}
				
				ans=max(ans,res);

				res=2;
				for(int k=0;k<n;k++){
					if(k==i || k==j) continue;
					if(hypot(x[k]-cx2,y[k]-cy2)<1.0+eps) res++;
				}

				ans=max(ans,res);

			}
		}

		cout<<ans<<endl;

	}

	return 0;
}