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>
#define EPS 1.0e-5
using namespace std;

typedef struct {
  double x,y;
} P;

double dis(P p ){
  return sqrt( p.x*p.x + p.y*p.y );
}
double dis(P p1,P p2 ){
  P p;
  p.x = p1.x-p2.x;
  p.y = p1.y-p2.y;
  return dis( p );
}


int main(){
  int n, mx;
  P p[301], m, p1,p2, mv,dmv, c[2];
  double x;

  while( cin>>n && n ){
    for( int i=0;i<n;i++ )
      cin >> p[i].x >> p[i].y ;

    mx=0;

    for( int i=0;i<n;i++ ){
      for( int j=i+1;j<n;j++ ){
	if( dis(p[i],p[j])<2.0 ){
	  p1=p[i]; p2=p[j];
	  m.x = (p1.x+p2.x)/2.0;
	  m.y = (p1.y+p2.y)/2.0;

	  x = sqrt( 1-dis(p1,m) );
	  // ツ陳心ツ点 cツづーツ仰づ淞づゥ
	  if( fabs( p1.x-p2.x ) < EPS ){
	    mv.x=1.0;  mv.y=0.0;
	  }else{
	    double a1 = ( p2.y-p1.y )/( p2.x-p1.x );
	    mv.x=1.0; mv.y=-1.0/a1;
	  }
	  double scale = dis( mv );
	  dmv.x = x*mv.x / scale;
	  dmv.y = x*mv.y / scale;

	  c[0].x=m.x+dmv.x; c[0].y=m.y+dmv.y;
	  c[1].x=m.x-dmv.x; c[2].y=m.y-dmv.y;

	  for( int cn=0;cn<2;cn++ ){
	    int sum=0;
	    for( int k=0;k<n;k++ ){
	      if( dis( c[cn],p[k] ) < 1.0+EPS )
		sum++;
	    }
	    if( mx<sum ) mx=sum;
	  }
	}
      }
    }
    cout << mx << 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 <cmath>
#include <vector>
#include <iostream>

#define EPS	1.0e-10
#define rep(idx, m)	for(int idx = 0, idx##Max = m; idx < idx##Max; idx ++)

template<class LHS, class OP, class RHS>
class CExp
{
	const LHS	&m_lhs;
	const RHS	&m_rhs;
	
public:
	CExp(const LHS &lhs, const RHS &rhs) : m_lhs(lhs), m_rhs(rhs) {}
};

typedef struct _POINT
{
	double x, y;
	
/*	template<class T>
	_POINT operator =(const &T rhs)
	{
		return(*this);
	}*/
	_POINT operator *=(double s)
	{
		x *= s, y *= s;
		return(*this);
	}
	_POINT operator /=(double s)
	{
		x /= s, y /= s;
		return(*this);
	}
	double Abs()	{	return(::sqrt(x * x + y * y));	}
	double Abs2()	{	return(x * x + y * y);	}
	_POINT(double _x, double _y) : x(_x), y(_y) {}
	_POINT() {}
} POINT;

int main()
{
	int count;
	int pointNum;
	std::vector<POINT> points;
	
	while(std::cin >> pointNum, pointNum)
	{
		points.resize(pointNum);
		rep(i, points.size())
		{
			std::cin >> points[i].x >> points[i].y;
		}
		
		int maxCount = 0;
		rep(i, points.size())
		{
			rep(j, points.size())
			{
				if(i == j)	{	continue;	}
				POINT diff(points[i].x - points[j].x, points[i].y - points[j].y);
				if(diff.Abs() > 2.0 + EPS)	{	continue;	}
				POINT normal(diff.y, -diff.x);
				double distToCenter = ::sqrt(1 - normal.Abs2() / 4);
				
				normal *= distToCenter / normal.Abs();
				
				POINT center(points[j].x + diff.x / 2 + normal.x, points[j].y + diff.y / 2 + normal.y);
				
				count = 0;
				rep(k, points.size())
				{
					double x = points[k].x - center.x;
					double y = points[k].y - center.y;
					
					if(x * x + y * y <= 1 + EPS)	{	count ++;	}
				}
				if(maxCount < count)	{	maxCount = count;	}
			}
		}
		std::cout << maxCount << std::endl;
	}
}

C++: #include<iostream>
#include<algorithm>
using namespace std;

#define MAX_N 310

double x[MAX_N], y[MAX_N];
int n, maxn;

int main() {
	while (true) {
		cin >> n;
		if (n == 0) { break; }
		for (int i = 0; i < n; i++) {
			cin >> x[i] >> y[i];
		}
		maxn = 0;
		for (int i = 0; i < n; i++) {
			for (int j = i + 1; j < n; j++) {
				int sum = 0;
				double X = (x[i] + x[j]) / 2.0;
				double Y = (y[i] + y[j]) / 2.0;
				for (int k = 0; k < n; k++) {
					if (((x[k] - X)*(x[k] - X) + (y[k] - Y)*(y[k] - Y)) <= 1) {
						sum++;
					}
				}
				maxn = max(maxn, sum);
			}
		}
		cout << maxn << endl;
	}
	return 0;
}

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

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