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 <complex>
using namespace std;

#define P complex<double>
#define EPS 1.0e-4

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 m  = (p[i]+p[j]) / 2.0; // p[i]Æp[j]Ì_
	P nn = (p[i]-p[j]) * P(0,1);       // @üxNg
	nn = nn / abs(nn);

	double x = abs( p[i]-p[j] )/2.0;
	P c = m + nn*sqrt( 1.0-x*x );        // ~ÌS

	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 <stdio.h>
#include <algorithm>
#include <iostream>
#include <stack>
#include <vector>
#include <math.h>
#include <queue>

using namespace std;

struct Point{
	double x,y;
};

void func(int N){

	Point point[N];
	double x_min = 100,y_min = 100,x_max = -1,y_max = -1;

	for(int i = 0; i < N; i++){
		scanf("%lf %lf",&point[i].x,&point[i].y);
		x_min = min(point[i].x,x_min);
		y_min = min(point[i].y,y_min);
		x_max = max(point[i].x,x_max);
		y_max = max(point[i].y,y_max);
	}

	if(x_min > 1)x_min -= 1;
	if(y_min > 1)y_min -= 1;
	if(x_max < 9)x_max += 1;
	if(y_max < 9)y_max += 1;

	int count,max_count = 0;

	for(double y = y_min; y <= y_max; y += 0.02){
		for(double x = x_min; x <= x_max; x += 0.02){
			count = 0;
			for(int i = 0; i < N; i++){
				if(sqrt((x-point[i].x)*(x-point[i].x)+(y-point[i].y)*(y-point[i].y)) <= 1.0)count++;
			}
			max_count = max(max_count,count);
		}
	}

	printf("%d\n",max_count);
}


int main(){

	int N;

	while(true){
		scanf("%d",&N);
		if(N == 0)break;

		func(N);
	}


    return 0;
}

C++: #include <bits/stdc++.h>
 
#define loop(n, i) for(int i=0;i<n;i++)
#define all(v) v.begin(),v.end()

using namespace std;
typedef pair<double, double> P;
const double EPS = 1e-6;

double dist(P a, P b)
{
    return hypot(a.first - b.first, a.second - b.second);
}

P center(P a, P b, P c)
{
    double ax = a.first, ay = a.second;
    double bx = b.first, by = b.second;
    double cx = c.first, cy = c.second;
    double d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by));
    auto norm = [](P p) { return p.first * p.first + p.second * p.second; };
    return P(
        (norm(a) * (by - cy) + norm(b) * (cy - ay) + norm(c) * (ay - by)) / d,
        (norm(a) * (cx - bx) + norm(b) * (ax - cx) + norm(c) * (bx - ax)) / d
    );
}

int main()
{
    while (1) {
        int n; cin >> n;
        if (!n) break;
        vector<P> pnt;
        loop (n, i) {
            double x, y; cin >> x >> y;
            pnt.emplace_back(x, y);
        }
        int ans = 1;
        loop (n, i) loop (n, j) {
            if (dist(pnt[i], pnt[j]) < 2) ans = 2;
        }
        loop (n, i) loop (n, j) loop (n, k) {
            int cnt = 0;
            P o = center(pnt[i], pnt[j], pnt[k]);
            loop (n, l) {
                if (dist(o, pnt[l]) < 1 + EPS) cnt++;
            }
            ans = max(ans, cnt);
        }
        cout << ans << endl;
    }
    return 0;
}

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

#define MAX_N 310

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

				long double X = (x[i] + x[j]) / 2.0;
				long 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);

				long double X1, Y1, X2, Y2, X3, Y3, VX2, VY2, Dist, Dist2;

				Dist = sqrtl((x[j] - x[i])*(x[j] - x[i]) + (y[j] - y[i])*(y[j] - y[i]));
				Dist2 = sqrtl(4 - Dist*Dist);

				X2 = x[j] - x[i]; Y2 = y[j] - y[i];
				VX2 = (X2 / Dist)*Dist2 / 2; VY2 = (Y2 / Dist)*Dist2 / 2;

				X1 = X + VY2; Y1 = Y - VX2;
				X3 = X - VY2; Y3 = Y + VX2;

				sum = 0;
				for (int k = 0; k < n; k++) {
					if (((x[k] - X1)*(x[k] - X1) + (y[k] - Y1)*(y[k] - Y1)) <= 1.0001) {
						sum++;
					}
				}
				maxn = max(maxn, sum);
				sum = 0;
				for (int k = 0; k < n; k++) {
					if (((x[k] - X3)*(x[k] - X3) + (y[k] - Y3)*(y[k] - Y3)) <= 1) {
						sum++;
					}
				}
				maxn = max(maxn, sum);
			}
		}
		for (int i = 0; i < n; i++) {
			int sum = 0;
			long double X = x[i];
			long double Y = y[i];
			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;
}