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;

double Dis(double x1, double y1, double x2, double y2);
int main(void){
	while(1){
		int n;
		int ans = 1;
		vector< pair<double, double> > p;

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

		for(int i=0; i<n; i++){
			pair<double, double> t;
			cin >> t.first >> t.second;
			p.push_back(t);
		}

		for(int i=0; i<n; i++){
			for(int j=0; j<n; j++){
				//中点
				double mpx = (p[i].first + p[j].first) / 2.0;
				double mpy = (p[i].second + p[j].second) / 2.0;
				//2点間の距離
				double d1 = Dis(p[i].first, p[i].second, p[j].first, p[j].second);
				//p[i]からp[j]へ向かうベクトル
				double v1x = (p[i].first - p[j].first) / d1;
				double v1y = (p[i].second - p[j].second) / d1;
				//中点から円の中心へ向かう単位ベクトル
				double v2x = -v1y;
				double v2y = v1x;
				//中心から円の中心へ向かうベクトルの長さ
				double d2 = sqrt(1.0 - (d1/2.0)*(d1/2.0));

				//円の中心の座標
				double cx = mpx + v2x * d2;
				double cy = mpy + v2y * d2;

				int cnt = 0;
				for(int k=0; k<n; k++){
					if(k == i || k == j || Dis(cx, cy, p[k].first, p[k].second) <= 1.0)	cnt++;
				}
				if(cnt > ans){
					ans = cnt;
				}

				cx = mpx - v2x * d2;
				cy = mpy - v2y * d2;

				cnt = 0;
				for(int k=0; k<n; k++){
					if(Dis(cx, cy, p[k].first, p[k].second) <= 1.0)	cnt++;
				}
				if(cnt > ans){
					ans = cnt;
				}
			}
		}

		cout << ans << endl;


	}
	
	return 0;
}

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

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

double distance(pair<double,double> p1,pair<double,double> p2){
	return sqrt((p1.first-p2.first)*(p1.first-p2.first)+(p1.second-p2.second)*(p1.second-p2.second));
}

int main(){
	int n;
	while(cin>>n,n){
		vector<double> x(n);
		vector<double> y(n);
		int res=1;
		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++){
				pair<double,double> midp=make_pair((x[i]+x[j])/2,(y[i]+y[j])/2);
				double a=-(x[i]-x[j]),b=(y[i]-y[j]);
				//a(x-mid[i].first)+b(y-mid.second)=0
					for(int k=0;k<2;k++){
						pair<double,double> lef=midp;
						pair<double,double> rig;
						if(k==0) rig=make_pair(-b*1e6+midp.first,a*1e6+midp.second);
						if(k==1) rig=make_pair(b*1e6+midp.first,-a*1e6+midp.second);
						while(distance(lef,rig)>=0.00000001){
								pair<double,double> mid=make_pair((lef.first+rig.first)/2,(lef.second+rig.second)/2);
								if(distance(mid,make_pair(x[i],y[i]))>1){
									rig=mid;
								}
								else{
									lef=mid;
								}

						}
						int count=0;
						for(int r=0;r<n;r++){
							if(distance(make_pair(x[r],y[r]),lef)<1+0.000001){
								count++;
							}
						}
						res=max(res,count);
					}
			}

		}
		cout<<res<<endl;

	}

	return 0;
}

C++: #include <iostream>
#include <math.h>

using namespace std;

int n;

double xi[300];
double yi[300];

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

int search(int i,int j)
{
	double mx = (xi[i] + xi[j]) / 2,my = (yi[i] + yi[j]) / 2;
	double diff2 = distant2(mx,my,xi[i],yi[i]);

	if(diff2 > 1.0)
		return 1; //1ツ古つつセツつッ
	double len = sqrt(1 - diff2);
	double diff = sqrt(diff2);
	double dx = len * (yi[i] - my) / diff,dy = len * (xi[i] - mx) / diff;
	
	int sign[] = {1,-1};
	int counter[2] = {0};
	for(int s = 0; s < 2; s++)
	{
		double x = mx - sign[s] * dx,y = my + sign[s] * dy;

		for(int k = 0; k < n; k++){
			if(k == i || k == j){counter[s]++; continue;}
			if((x-xi[k])*(x-xi[k]) + (y-yi[k])*(y-yi[k])  <= 1.0)
				counter[s]++;
		}
	}

	return max(counter[0],counter[1]);
}



int solve()
{
	for(int i = 0; i < n; i++){
		cin >> xi[i] >> yi[i];
	}

	int count = 0;

	for(int i = 0; i < n; i++){
		for(int j = i+1; j < n; j++){
			count = max(count,search(i,j));
		}
	}

	return count;
}

int main() {
  while (cin >>n,n) {
    cout << solve() << endl;
  }

  return 0;
}

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

class Solve {
private:
	using Point = std::complex<long double>;
	long double epsilon{1e-6};

	int N;
	std::vector<Point> points;

	int countPoints(int pi_1, int pi_2)
	{
		Point &p1{points[pi_1]}, &p2{points[pi_2]};
		if (std::abs(p1 - p2) > 2 + epsilon) return 0;
		Point center{p1 + std::polar((long double)1.0, std::acos(std::abs(p2 - p1) / 2) + std::arg(p2 - p1))};
		int count{};
		for (auto& e: points)
			if (std::abs(e - center) <= 1 + epsilon)
				count++;
		return count;
	}
public:
	bool is_last_query_{};
	Solve()
	{
		scanf("%d", &N);
		if (N == 0)
		{
			is_last_query_ = true;
			return;
		}
		points.resize(N);
		for (auto& e: points)
		{
			long double r, i;
			std::cin >> r >> i;
			e.real(r);
			e.imag(i);
		}
		int max{1};
		for (int i{}; i < N; i++)
			for (int j{}; j < N; j++)
				if (i != j)
					max = std::max(max, countPoints(i, j));
		printf("%d\n", max);
		return;
	}
};

int main()
{
	while (!Solve().is_last_query_);
	return 0;
}

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

double Dis(double x1, double y1, double x2, double y2);
int main(void){
	while(1){
		int n;
		int ans = 1;
		vector< pair<double, double> > p;

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

		for(int i=0; i<n; i++){
			pair<double, double> t;
			cin >> t.first >> t.second;
			p.push_back(t);
		}

		for(int i=0; i<n; i++){
			for(int j=0; j<n; j++){
				//中点
				double mpx = (p[i].first + p[j].first) / 2.0;
				double mpy = (p[i].second + p[j].second) / 2.0;
				//2点間の距離
				double d1 = Dis(p[i].first, p[i].second, p[j].first, p[j].second);
				//p[i]からp[j]へ向かうベクトル
				double v1x = (p[i].first - p[j].first) / d1;
				double v1y = (p[i].second - p[j].second) / d1;
				//中点から円の中心へ向かう単位ベクトル
				double v2x = -v1y;
				double v2y = v1x;
				//中心から円の中心へ向かうベクトルの長さ
				double d2 = sqrt(1.0 - (d1*d1/4.0));

				//円の中心の座標
				double cx = mpx + v2x * d2;
				double cy = mpy + v2y * d2;

				int cnt = 0;
				for(int k=0; k<n; k++){
					if(k == i || k == j || Dis(cx, cy, p[k].first, p[k].second) <= 1.0)	cnt++;
				}
				if(cnt > ans){
					ans = cnt;
				}

				cx = mpx - v2x * d2;
				cy = mpy - v2y * d2;

				cnt = 0;
				for(int k=0; k<n; k++){
					if(Dis(cx, cy, p[k].first, p[k].second) <= 1.0)	cnt++;
				}
				if(cnt > ans){
					ans = cnt;
				}
			}
		}

		cout << ans << endl;


	}
	
	return 0;
}

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