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;

typedef long long ll;
typedef pair<int,int> P;
const double EPS = 1e-12;
const int INF = numeric_limits<int>::max()/2;
const int MOD = 1e9+7;

typedef long double ld;
typedef complex<ld> Point;

const ld eps = 1e-9, pi = acos(-1.0);

namespace std {
	bool operator<(const Point &lhs, const Point &rhs) {
		if (lhs.real() < rhs.real() - eps) return true;
		if (lhs.real() > rhs.real() + eps) return false;
		return lhs.imag() < rhs.imag();
	}
}

Point input_point() {
	ld x, y;
	cin >> x >> y;
	return Point(x, y);
}

bool eq(ld a, ld b) {
	return (abs(a - b) < eps);
}

class Circle {
public:
	Point p;
	ld r;
	Circle() : p(Point(0, 0)), r(0) {}
	Circle(Point p, ld r) : p(p), r(r) {}
};

vector<Point> is_cc(Circle c1, Circle c2) {
	vector<Point> res;
	ld d = abs(c1.p - c2.p);
	ld rc = (d * d + c1.r * c1.r - c2.r * c2.r) / (2 * d);
	ld dfr = c1.r * c1.r - rc * rc;
	if (abs(dfr) < eps) dfr = 0.0;
	else if (dfr < 0.0) return res; // no intersection
	ld rs = sqrt(dfr);
	Point diff = (c2.p - c1.p) / d;
	res.push_back(c1.p + diff * Point(rc, rs));
	if (dfr != 0.0) res.push_back(c1.p + diff * Point(rc, -rs));
	return res;
}

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

    int n;
    while(cin>>n,n){
        vector<Point> ps(n);
        for(int i=0;i<n;i++){
            ps[i] = input_point();
        }
        if(n<=2){
            cout<<n<<endl;
            continue;
        }
        int res=0;
        for(int i=0;i<n;i++){
            for(int j=i+1;j<n;j++){
                Circle c1=Circle(ps[i],1.0),c2=Circle(ps[j],1.0);
                vector<Point> vs=is_cc(c1,c2);
                int sz=vs.size();
                if(sz==0) continue;
                for(int k=0;k<sz;k++){
                    Point p=vs[k];
                    int tmp=0;
                    for(int idx=0;idx<n;idx++){
                        if(abs(ps[idx]-p)<=1+eps) tmp++;
                    }
                    res=max(res,tmp);
                }
            }
        }
        cout<<res<<endl;
    }
}

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

    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 <cstdio>
#include <complex>
#include <vector>
using namespace std;
typedef complex<double> P;
const double PI=3.1415926535897932;

bool inC(P p, P c, double r){
    return abs(p-c)<r+1e-8;
}

int main(){
    int n,i,j,k,ans,cnt;
    double x, y;
    while(scanf("%d",&n)){
        if(n==0) break;
        ans=0;
        vector<P> ps(n);
        for(i=0; i<n; i++){
            scanf("%lf%lf",&x,&y);
            ps[i] = P(x,y);
        }
        for(i=0;i<n;i++)
            for(j=i+1;j<n;j++)if(abs(ps[i]-ps[j])<2){
                cnt=0;
                double d=sqrt(1-abs(ps[j]-ps[i])*abs(ps[j]-ps[i])*0.25)
                /abs(ps[j]-ps[i]);
                P c1=(ps[i]+ps[j])/2.0+(ps[j]-ps[i])*polar(d,PI/2);
                P c2=(ps[i]+ps[j])/2.0+(ps[j]-ps[i])*polar(d,-PI/2);
                for(k=0;k<n;k++){
                    cnt += inC(ps[k],c1,1.0);
                }
                if(cnt>ans)ans=cnt;
                cnt=0;
                for(k=0;k<n;k++){
                    cnt += inC(ps[k],c2,1.0);
                }
                if(cnt>ans)ans=cnt;
            }
        printf("%d\n",ans);
    }
    return 0;
}

C++: #include <bits/stdc++.h>
#define N 301
using namespace std;
typedef complex <double> P;
const double EPS = 1e-8;
P p[N];
int n;

int cnt(P c){
  int res=0;
  for(int i=0;i<n;i++)if(abs(c-p[i])<1+EPS) res++;
  return res;
}

P get_C(P a,P b){
  double x=abs(b-a);
  return P(x/2,sqrt(1-x*x/4))*(b-a)/abs(b-a)+a;  
}


int main(){
  while(1){
  cin>>n;
  if(!n)break;
  double a,b;
  for(int i=0;i<n;i++)cin>>a>>b,p[i]=P(a,b);
  
  int ans=1;
  for(int i=0;i<n;i++)
    for(int j=i+1;j<n;j++)
      if(abs(p[i]-p[j])<=2){
	ans=max(ans,cnt(get_C(p[i],p[j])));
	ans=max(ans,cnt(get_C(p[j],p[i])));
      }
  cout <<ans<<endl;  
  }
  return 0;
}