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;
struct Point {
	long double x;
	long double y;
	Point(long double _x = 0, long double _y = 0) :x(_x), y(_y) {}
	Point operator+(Point p) {
		return Point(x + p.x, y + p.y);
	}
	Point operator-(Point p)  {
		return Point(x - p.x, y - p.y);
	}
	long double abs() {
		return sqrt(norm());
	}
	long double norm() {
		return x*x + y*y;
	}
};
typedef Point Vector;
struct Circle {
	Point c;
	long double r;
	Circle(Point _c = Point(), long double _r = 0.0) :c(_c), r(_r) {}
};
long double arg(Vector p) {
	return atan2(p.y, p.x);
}
Vector polar(long double a, long double r) {
	return Point(cos(r)*a, sin(r)*a);
}
pair<Point, Point> getCrossPoints(Circle c1, Circle c2) {
	long double d = (c1.c - c2.c).abs();
	long double a = acos((c1.r*c1.r + d*d - c2.r*c2.r) / (2 * c1.r*d));
	long double t = arg(c2.c - c1.c);
	return make_pair(c1.c + polar(c1.r, t + a), c1.c + polar(c1.r, t - a));
}
long double X[300], Y[300];
int main() {
	int N;
	while (cin >> N, N) {
		for (int i = 0; i < N; i++) {
			cin >> X[i] >> Y[i];
		}
		int ans = 0;
		for (int i = 0; i < N; i++) {
			for (int j = i + 1; j < N; j++) {
				Circle c1(Point(X[i], Y[i]), 1.0);
				Circle c2(Point(X[j], Y[j]), 1.0);
				if ((c1.c - c2.c).abs() > 2.0) continue;
				pair<Point, Point> c = getCrossPoints(c1, c2);
				int cnt1 = 0, cnt2 = 0;
				for (int k = 0; k < N; k++) {
					Point p(X[k], Y[k]);
					if ((p - c.first).abs() <= 1.0) cnt1++;
					if ((p - c.second).abs() <= 1.0) cnt2++;
				}
				ans = max(ans, cnt1);
				ans = max(ans, cnt2);
			}
		}
		cout << ans << endl;
	}
}

C++: #include <iostream>
#include <deque>
#include <string>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <string>
#include <cmath>
#include <complex>

#define loop(i,a,b) for(int i=a; i<int(b); i++)
#define rep(i,n) loop(i,0,n)

using namespace std;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef complex<double> vec;

double const pi = acos(-1);
double const eps = 1e-5;

int main(){
    int n;
    while (cin >> n && n){
        vector<vec> v(n);
        rep(i, n){
            double x, y; cin >> x >> y;
            vec p(x, y);
            v[i] = p;
        }

        if (n == 1) {
            puts("1");
            continue;
        }
        int ans = 1;
        rep(i, n)rep(j, i){
            double d = abs(v[i] - v[j]);
            if (d - eps > 2) continue;
            ans = max(ans, 2);
            double h = 1 - d*d / 4;
            vec c1 = polar(h, arg(v[i] - v[j]) + pi / 2) + (v[i] + v[j]) / 2.;
            vec c2 = polar(h, arg(v[i] - v[j]) - pi / 2) + (v[i] + v[j]) / 2.;
            int a1 = 0, a2 = 0;
            rep(k, n){
                if (abs(c1 - v[k]) <= 1 - eps) a1++;
                if (abs(c2 - v[k]) <= 1 - eps) a2++;
            }
            ans = max(ans, a1);
            ans = max(ans, a2);
        }
        cout << ans << endl;
    }
}

C++: #include<cstdio>  
#include<cstdlib>  
#include<cstring>  
#include<cmath>  
#include<algorithm>  
#define eps 1e-8  
using namespace std;  
struct point{  
    double x,y;  
}A[310];  
double dist(point a,point b){  
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));  
}  
void getmid(point p1,point p2,point &center){//????????????????????????????????????  
    point mid;  
    mid.x=(p1.x+p2.x)/2.0;mid.y=(p1.y+p2.y)/2.0;  
    double angle=atan2(p1.x-p2.x,p2.y-p1.y);  
    double dcm=sqrt(1-dist(p1,mid)*dist(p1,mid));  
    center.x=mid.x+dcm*cos(angle);center.y=mid.y+dcm*sin(angle);  
}  
int main()  
{  
    int n,i,j,k;  
    while(scanf("%d",&n),n){  
        for(i=0;i<n;++i){  
            scanf("%lf%lf",&A[i].x,&A[i].y);  
        }  
        int ans=1;  
        for(i=0;i<n;++i){  
            for(j=i+1;j<n;++j){  
                if(dist(A[i],A[j])>2.0)continue;  
                point center;  
                getmid(A[i],A[j],center);  
                int cnt=0;  
                for(k=0;k<n;++k){  
                    if(dist(A[k],center)<1.0+eps)cnt++;  
                }  
                ans=max(ans,cnt);   
            }  
        }  
        printf("%d\n",ans);  
    }  
    return 0;  
}

C++: #include <cmath>  
#include <algorithm>
#include <iostream>
#include <utility>

using namespace std;

const double eps=1e-8; 
struct point{
	double x,y;	
}P[310];
int x,y,N;
int ans=1;

double cal_distance(point a, point b){
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

point cal_center(point a, point b){
	point mid,center;
	mid.x = (a.x+b.x)/2.0;
	mid.y = (a.y+b.y)/2.0;
	double angle = atan2(a.x-b.x,b.y-a.y);
	double line = sqrt(1-(cal_distance(a,mid)*cal_distance(a,mid)));
	center.x = mid.x+line*cos(angle);
	center.y = mid.y+line*sin(angle);
	return center;
}

int main(){
	while(cin >> N && N>0){
		ans = 1;
		for (int i=0; i<N; i++){
				cin >> P[i].x >> P[i].y;
			}
			for (int i=0; i<N; i++){
				for (int j=i+1; j<N; j++){
					if (cal_distance(P[i],P[j]) > 2.0) continue;
					point center = cal_center(P[i], P[j]);
					int count = 0;
					for (int k=0; k<N;k++){
						if (cal_distance(P[k],center) < (1.0+eps)) count++;
					}
					ans = max(ans,count);
				}
			}
		cout << ans << endl;
	}
}

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

int main() {
  const complex<double> I(0,1);
  complex<double> d2(2,0);
  for(int tci = 0; ; tci++) {
    int n; scanf("%d", &n);
    if(!n) break;
    complex<double> *pts = new complex<double>[n];
    for(int i = 0; i < n; i++) {
      scanf("%lf%lf", &pts[i].real(), &pts[i].imag());
    }
    int maxcnt = 1;
    for(int i = 0; i < n; i++) {
      for(int j = 0; j < n; j++) {
        if(i==j) continue;
        complex<double> d = (pts[i] - pts[j]) / d2;
        double absd = abs(d);
        complex<double> k = d*I/absd * sqrt(1-absd*absd);
        complex<double> c = (pts[i] + pts[j]) / d2 + k;
        int cnt = 0;
        for(int k = 0; k < n; k++) {
          cnt += (abs(pts[k]-c) <= 1 + 1e-9);
        }
        maxcnt = max(maxcnt, cnt);
      }
    }
    delete[] pts;
    printf("%d\n", maxcnt);
  }
  return 0;
}