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

complex<double> p[300];

complex<double> core(int a,int b){
  complex<double> res,mid,norm;

  mid = (p[a]+p[b]) / 2.0;

  norm = (p[a] - p[b]) * complex<double>(0,1);
  norm = norm / abs(norm);
  norm = norm * sqrt(1.0 - abs(mid-p[a])*abs(mid-p[a]));

  res = norm + mid;
  return res;
}
  

int main(){
  int n;
  complex<double> o,r;
  int max,tmp;

  while(1){
    cin >> n;
    if(!n)break;

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

    max = 1;
    for(int i=0;i<n;i++){
      for(int j=0;j<n;j++){
	r = p[i]-p[j];
	if(i!=j && r.real()*r.real() + r.imag()*r.imag() < 4.0){
	  o = core(i,j);
	  tmp = 2;
	  for(int k=0;k<n;k++){
	    if(k!=i && k!=j){
	      r = o-p[k];
	      if( r.real()*r.real() + r.imag()*r.imag() <1.0)tmp++;
	    }
	  }
	  if(max < tmp)max = tmp;
	}
      }
    }
    cout << max << endl;
  }
}

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

using namespace std;

typedef complex<double> Point;

bool eqv(double a, double b) {
    return abs(a-b) < 1e-10;
}

bool inCircle(const Point &p, const Point &center, double radius) {
    double dist2 = norm(p-center);
    double rad2 = radius*radius;
    return dist2 < radius*radius || eqv(dist2, rad2);
}

int check(const Point &center, const vector<Point> &points) {
    const int N = points.size();
    int cnt = 0;
    for(int i = 0; i < N; ++i) {
        if(inCircle(points[i], center, 1)) ++cnt;
    }
    return cnt;
}

int test(int a, int b, const vector<Point> &points) {
    const Point v(points[a]-points[b]);
    const double hlen = abs(v)/2;

    if(hlen >= 1) return 0;
    const double plen = sqrt(1 - hlen*hlen);
    Point perp = v*Point(0, 1);
    perp /= abs(perp);
    const Point mid = points[b] + v/2.0;
    return max(check(mid + perp*plen, points), check(mid - perp*plen, points));
}

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

    while(true) {
        int N;
        cin >> N;
        if(!N) break;

        vector<Point> points(N);
        for(int i = 0; i < N; ++i) {
            cin >> points[i].real() >> points[i].imag();
        }
        int ans = 0;
        for(int i = 0; i < N; ++i) {
            for(int j = 0; j < i; ++j) {
                ans = max(ans, test(i, j, points));
            }
        }
        if(N == 1) ans = 1;
        cout << ans << endl;
    }
    return 0;
}

C++: #include <iostream>
#include <vector>
#include <cmath>

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

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

bool zero(double x) {
  return fabs(x) <= EPS;
}

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

pt midpoint(pt p1, pt p2){
    double midx = (p1.x+p2.x)/2;
    double midy = (p1.y+p2.y)/2;
    return make_pair(midx, midy);
}

int get_max_points(pt centre, vector<pt> points){

    int count = 0;
    for (int i = 0; i < points.size(); i++){
        pt current = points[i];
        double num = sq(current.x - centre.x) + sq(current.y - centre.y);
        // cout << num << endl;
        // cout << fabs(num - 1) << endl;
        if (sqrt(num) < 1+EPS) {
            count++;
        }
    }
    return count;
}

int find_max_points(vector<pt> points){
    // for every pair of two points, create a cricle and check num points in circle

    int max_points = 0;
    for (int i = 0; i < points.size(); ++i){
        pt p1 = points[i];
        for (int j = i+1; j < points.size(); ++j){
            pt p2 = points[j];

            pt mid = midpoint(p1, p2);
            double q = dist(p1, p2);

            double x1 = mid.x + sqrt(1-sq(q/2))*(p1.y-p2.y)/q;
            double y1 = mid.y + sqrt(1-sq(q/2))*(p2.x-p1.x)/q;
            pt c1(x1, y1);
            int m = get_max_points(c1, points);
            max_points = max(m, max_points);
            // cout << mid.x << " " << mid.y << endl;
            // int m = get_max_points(mid, points);
            double x2 = mid.x - sqrt(1-sq(q/2))*(p1.y-p2.y)/q;
            double y2 = mid.y - sqrt(1-sq(q/2))*(p2.x-p1.x)/q;
            pt c2(x1, x2);
            m = get_max_points(c2, points);

            max_points = max(max_points, m);
        }
    }
    // cout << max_points 
    return max_points;
}

int main(){

    int n;

    while (1){
        cin >> n;
        if (n == 0) break;
        vector<pt> points;
        for (int i = 0; i < n; ++i){
            double x, y;
            cin >> x >> y;
            points.push_back(make_pair(x, y));
        }
        cout << find_max_points(points) << endl;
    }


    return 0;
}

C++: #include <iostream>
#include <complex>
#include <vector>
#include <cmath>
#define EPS 1.0e-10
using namespace std;

#define X real()
#define Y imag()
typedef complex<double> P;
 
vector<P> p;
int check(P a){
	int ret=0;
	for(int i=0;i<p.size();i++) if(abs(a-p[i]) < 1+EPS) ret ++;
	return ret;
}

int adjcircle(P a, P b){
	P dir = a-b, pos=(a+b)*0.5;
	dir = P(dir.Y, -dir.X);
	dir = dir*(1/abs(dir))*sqrt(1-abs(pos-a)*abs(pos-a));
	return max(check(pos+dir), check(pos-dir));
}

main(){
	int i,j,N;
	while(cin >> N, N){
		p = vector<P>(N);
		int ans = 0;
		for(i=0;i<N;i++) cin >> p[i].X >> p[i].Y;
		for(i=1;i<N;i++)
			for(j=0;j<i;j++)
				if(abs(p[i]-p[j])<=2.0) ans = max(ans, adjcircle(p[i], p[j]));
		cout << ans << endl;
	}
}