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 <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 -= 0.047;
    if(y_min > 1)y_min -= 0.047;
    if(x_max < 9)x_max += 0.047;
    if(y_max < 9)y_max += 0.047;

    int count,max_count = 0;

    for(double y = y_min; y <= y_max; y += 0.01){
        for(double x = x_min; x <= x_max; x += 0.01){
            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 <iostream>
#include <utility>
#include <cmath>
#include <algorithm>

using namespace std;
typedef pair<double, double> pt;
const double EPS = 1e-5;
#define x first
#define y second

double sq(double a) { return a * a; }

double dist(pt p, pt q) {
  return sqrt(sq(p.x - q.x) + sq(p.y - q.y));
}

bool lessthan(double a, double b) {
    if (a - b < EPS)
        return true;
    return false;
}

pt midpoint(pt a, pt b) {
    return pt((a.x + b.x)/2, (a.y + b.y)/2);
}

pt operator+(pt a, pt b) {
    return pt(a.x + b.x, a.y + b.y);
}

pt operator-(pt a, pt b) {
  return pt(a.x - b.x, a.y - b.y);
}

pt operator*(pt a, double b) {
    return pt(a.x * b, a.y * b);
}

pt rotate(pt p, double a) {
  pt res;
  res.x = p.x * cos(a) - p.y * sin(a);
  res.y = p.x * sin(a) + p.y * cos(a);
  return res;
}



const int maxn = 300;
pt pts[maxn + 1];

int main() {
    int N;
    
    cin >> N;
    while (N != 0) {
        for (int i = 0; i < N; ++i) {
            double x, y;
            cin >> x >> y;
            pts[i] = pt(x, y);
        }
        int maxcount = 0;
        for (int i = 0; i < N ; i++){
            for (int j = 0; j < N; j++) {
                if (j == i) continue;
                // for every pair of points
                if (dist(pts[i], pts[j]) > 2) continue;
                pt A = pts[i];
                pt B = pts[j];
                //cout << "for points " << A.x << " " << A.y << " b: " << B.x << " "<< B.y << endl;
                pt M = midpoint(A, B);
                //cout << " M is " << M.x << " " << M.y << endl;
                pt BM = M - B;
                double distBM = dist(B,M);
                //cout << "dist BM " << distBM << endl; 
                double alpha = acos(distBM);
                pt BK = rotate(BM, alpha);
                pt BC = BK * (1 / distBM);
                pt C = B + BC;
                //cout << "C is  "<< C.x << " " << C.y << endl;
                int count = 2;
                for (int k = 0; k < N; k++) {
                    if (k == i || k == j) continue;
                    if (lessthan(dist(C, pts[k]), 1)) count++; 
                }
                maxcount = max(count, maxcount);
            }
        }
        cout << maxcount << endl;
        cin >> N;
    }
}

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

    int count,max_count = 0;

    for(double y = y_min; y <= y_max; y += 0.01){
        for(double x = x_min; x <= x_max; x += 0.01){
            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)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 <iostream>
#include <complex>
using namespace std;

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

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

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

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

vector<P> center(P a, P b)
{
    auto t = atan2(b.first - a.first, b.second - a.second);
    auto x = dist(a, b) / 2;
    auto y = sqrt(1.0 - x * x);
    auto cs = cos(t), sn = sin(t);
    return { 
        { cs * x - sn * y + a.first, sn * x + cs * y + a.second },
        { cs * x + sn * y + a.first, sn * x - cs * y + a.second }
    };
}

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) {
            auto ps = center(pnt[i], pnt[j]);
            for (auto o : ps) {
                int cnt = 0;
                loop (n, l) {
                    //if (dist(o, pnt[l]) < 1 + EPS) cnt++;
                }
                ans = max(ans, cnt);
            }
        }
        cout << ans << endl;
    }
    return 0;
}