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;

const double r = 1.0;
const double PI = acos(-1.0);

struct Point {
    double x, y;
    Point() {}
    Point(double x, double y) : x(x), y(y) {}
};
inline double Distance2(const Point &a, const Point &b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

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

    int n;
    std::random_device rnd;
    std::mt19937 mt(rnd());
    std::uniform_real_distribution<> rand100(0.0, 10.0);

    while (cin >> n, n) {
        // Input
        vector<Point> p(n);
        for (int i = 0; i < n; ++i)
            cin >> p[i].x >> p[i].y;

        // Solve
        int ans = 1;
        auto start = std::chrono::high_resolution_clock::now();
        auto end = std::chrono::high_resolution_clock::now();
        auto take_time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

        while (true) {
            end = std::chrono::high_resolution_clock::now();
            take_time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
            if (400 < take_time.count())
                break;

            for (int i = 0; i < 10; ++i) {
                Point c(rand100(mt), rand100(mt));
                int num = 0;
                for (int j = 0; j < n; ++j) {
                    double d = Distance2(p[j], c);
                    if (d <= 1.0)
                        ++num;
                }
                ans = max(ans, num);
            }
        }

        // Output
        cout << ans << endl;
    }

    return 0;
}

C++: #include <iostream>
#include <iomanip>
#include <complex>

using namespace std;

complex<double> points[300];

constexpr double epsilon = 1.0e-6;

int count_points(const complex<double> &center,int N){
  int num=0;
  for(int i=0;i<N;++i){
    double dist = abs(points[i]-center);
    if(dist-1.0<=epsilon){
      num++;
    }
  }
  return num;
}

int main(){
  int n;
  while(cin >> n&&n>0){
    for(int i=0;i<n;++i){
      double x,y;
      cin >> x >> y;
      points[i] = complex<double>(x,y);
    }
    
    int maxnum = 1;
    for(int i=0;i<n-1;++i){
      for(int j=i+1;j<n;++j){
        complex<double> v = points[j]-points[i];
        double dist = abs(v);
        int circnum=0;
        complex<double> center[2];
        if(abs(dist-2.0)<=epsilon){
          circnum = 1;
          center[0] = (points[i]+points[j])/2.0;
        }else if(dist-2.0<-epsilon){
          circnum = 2;
          complex<double> w = v/dist;
          w *= complex<double>(0,1);
          w *= sqrt(1.0-dist*dist/4.0);
          center[0] = (points[i]+points[j])/2.0 + w;
          center[1] = (points[i]+points[j])/2.0 - w;
        }
        for(int c=0;c<circnum;++c){
          maxnum = max(maxnum,count_points(center[c],n));
        }
      }
    }
    cout << maxnum << endl;
  }
  return 0;
}

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

#define EPS 1.0e-8
typedef complex<double> P;

P center(P a, P b){
  P m  = (a+b) / 2.0;         // p[i]Æp[j]Ì_
  P n = (a-b) * P(0.0,1.0);  // @üxNg
  n = n / abs(n);
  
  double x = abs( a-b )/2.0;
  return m + n * sqrt( 1.0-x*x );     // ~ÌS
}

int main(){
  int n;
  P p[301];

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

    int ans=0;
    for( int i=0;i<n;i++ ){
      for( int j=0;j<n;j++ ){
	if( i==j ) 	  continue;
	if( abs( p[i]-p[j] )+EPS > 2.0 ) continue;
	P c = center( p[i],p[j] );

	int num=0;
	for( int k=0;k<n;k++ )
	  if( abs( c-p[k] ) < 1.0 + EPS  )
	    num++;
	ans = max( ans,num );
      }
    }
    cout << ans << endl;
  }

  return 0;
}

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*1e5+midp.first,a*1e5+midp.second);
						if(k==1) rig=make_pair(b*1e5+midp.first,-a*1e5+midp.second);
						while(distance(lef,rig)>=0.00001){
								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.00001){
								count++;
							}
						}
						res=max(res,count);
					}
			}

		}
		cout<<res<<endl;

	}

	return 0;
}

C++: #include <iostream>
#include <cmath>
#include <algorithm>
#include <queue>
#include <map>
#define LP(i,n) for(int i = 0; i < n; i++)

using namespace std;
typedef pair<double,double> P;

double Dist(double x1[], double x2[]){
	return pow((x1[0]-x2[0])*(x1[0]-x2[0])+(x1[1]-x2[1])*(x1[1]-x2[1]),0.5);

}

int main(){
	int n;
	cin >> n;
	while(n){
		int mx = 0;
		P p[300];
		LP(i,n){
			cin >> p[i].first >> p[i].second;
		}
		sort(p,p+n);
		double x[300][2];
		LP(i,n){
			x[i][0] = p[i].first;
			x[i][1] = p[i].second;
		}
		double dist[300][300];
		LP(i,n){
			LP(j,n){
				dist[i][j] = Dist(x[i],x[j]);
			}
		}
		LP(i,n){
			for(int j = i+1; j < n; j++){
				if(x[j][0]-x[i][0] > 2){
					break;
				}else{
					if(dist[i][j] <= 2){
						double m[] = {(x[j][0]+x[i][0])/2,(x[j][1]+x[i][1])/2};
						double d[] = {x[j][1]-x[i][1],-x[j][0]+x[i][0]};
						double O[] = {0,0};
						double _d = Dist(d,O);
						LP(k,2){
							d[k] *= (pow(_d,-2)-0.25);
						}
						double o[2][2];
						LP(k,2){
							LP(l,2){
								o[k][l] = m[l];
								if(k) o[k][l] += d[l];
								else o[k][l] -= d[l];
							}
						}
						LP(k,2){
							int s = 0;
							for(int l = i; l < j; l++){
								if(Dist(x[l],o[k]) <= 1) s++;
							}
							for(int l = i-1; 0<=l && x[i][0]-x[l][0] <= 1;l--){
								if(Dist(x[l],o[k]) <= 1) s++;
							}
							for(int l = j; l < n && x[l][0]-x[j][0] <= 1; l++){
								if(Dist(x[l],o[k]) <= 1) s++;
							}
							mx = max(mx,s);
						}
					}
				}
			}
		}
		if(mx == 0) mx = 1;
		cout << mx << endl;
		cin >> n;
	}
}