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 <algorithm>
#include <cmath>
#define COS45 0.707107
using namespace std;

class Point {
public:
  double x, y;
  Point(double x = 0, 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);
  }
  Point operator *(double a) {
    return Point(x*a, y*a);
  }
  Point operator /(double a) {
    return Point(x/a, y/a);
  }
  bool operator ==(const Point &p) const {
    return x == p.x && y == p.y;
  }
  bool operator <(const Point &p) const {
    double l1 = x+y;
    double l2 = p.x+p.y;
    if(l1 != l2) return l1 < l2;
    if(y != p.y) return y < p.y;
    return x < p.x;
  }
};

typedef Point Vector;

double norm(Point p1, Point p2) {
  return (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y);
}

int n;
Point p[300];

int calc(Point p1, Point p2) {
  Point m = Point((p1.x+p2.x)/2.0, (p1.y+p2.y)/2.0);
  Point o;
  Vector moto = Vector(p2.y-p1.y, p1.x-p2.x);
  int cnt = 0;
  int newCnt;
  double p1p2 = norm(p1,p2);
  if(p1p2 > 4) return 1;
  double l = sqrt(4.0/p1p2-1)/2.0;

  for(int t = 0; t < 2; t++) {
    newCnt = 0;
    if(t == 1) {
      moto = moto * (-1);
    }
    o = m + moto * l;

    for(int i = 0; i < n; i++) {
      //if(p[i].x > o.x + COS45 && p[i].y > o.y + COS45) break;
      if(newCnt + n-i <= cnt) break;
      if(p[i] == p1 || p[i] == p2) continue;
      if((p[i].x-o.x)*(p[i].x-o.x) + (p[i].y-o.y)*(p[i].y-o.y) <= 1.0) {
	newCnt++;
      }
    }
    if(newCnt > cnt) cnt = newCnt;
  }
  return cnt + 2;
}

main() {
  while(1) {
    int ans;
    cin >> n;
    if(n == 0) break;
    for(int i = 0; i < n; i++) {
      cin >> p[i].x >> p[i].y;
    }
    sort(p, &p[n]);
    ans = -1;
    for(int i = 0; i < n; i++) {
      for(int j = i+1; j < n; j++) {
	ans = max(ans, calc(p[i],p[j]));
      }
    }
    cout << ans << endl;
  }
}

C++: #include <bits/stdc++.h>

#define rep(n) for(int i=0;i<n;i++)
#define repp(j, n) for(int j=0;j<n;j++)
#define reppp(i, m, n) for(int i=m;i<n;i++)
#define all(c) c.begin(), c.end()
#define rall(c) c.rbegin(), c.rend()
#define debug(x) cerr << #x << ": " << x << endl

using namespace std;

typedef long long ll;
typedef pair<ll, ll> Pll;
typedef pair<int, int> Pii;

const long double EPS = 10e-10;

class Segment;

class Point{
    public: 
    long double x,y;
    
    Point(long double x = 0.0, long double y = 0.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);
    }
    Point operator * (long double a) {
        return Point(x*a, y*a);
    }
    Point operator / (long double a) {
        return Point(x/a, y/a);
    }
    bool operator < (const Point &p) const {
        return (x != p.x ? x < p.x : y < p.y);
    }
    bool operator == (const Point &p) const {
        return fabsl(x - p.x) < EPS && fabsl(y - p.y) < EPS;
    }
    
    long double abs() {
        return sqrt(norm());
    }
    
    long double norm() {
        return x*x + y*y;
    }
    
    long double dot(Point p) {
        return x*p.x + y*p.y;
    }
    
    long double cross(Point p) {
        return x*p.y - y*p.x;
    }
};

class Segment{
    public: 
    Point p,q;
    long double grad;
    
    Segment(Point p, Point q) {
        this->p = p;
        this->q = q;
        this->grad = (q.y-p.y)/(q.x-p.x);
    }
    
};

class Circle{
    public:
    Point c;
    long double r;
    
    Circle(Point c = Point(), long double r = 0.0): c(c), r(r) {}
};

Point calc_center(Point p, Point mid, Point midv){
    midv = midv / midv.abs() * (1.0 - (p - mid).abs());
    return mid+midv;
}

int count(vector<Point> points, Point center){
    int tmp = 0;
    for(Point p: points){
        if((p-center).norm() <= 1.0) tmp++;
    }
    return tmp;
}

int main() {
    int n;
    while(scanf("%d", &n) && n){
        vector<Point> points(n);
        rep(n) cin >> points[i].x >> points[i].y;

        int ans = 1;
        rep(n){
            repp(j, n){
                if(j == i) continue;
                Point mid = (points[i] + points[j]) / 2;
                Point center = calc_center(points[i], mid, Point(mid.y, -mid.x));
                ans = max(count(points, center), ans);

                center = calc_center(points[i], mid, Point(-mid.y, mid.x));
                ans = max(count(points, center), ans);
            }
        }
        cout << ans << endl;
    }
}

C++: #include <bits/stdc++.h>
#define PB push_back
#define MP make_pair
#define REP(i,n) for (int i=0;i<(n);i++)
#define FOR(i,a,b) for(int i=(a);i<(b);i++)
#define ALL(a) (a).begin(),(a).end()
int INF=1e9;
int MOD=1000000007;
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
double EPS = 1e-10;
double add(double a,double b){
	if(abs(a+b)<EPS*(abs(a)+abs(b)))return 0;
	return a+b;
}
struct P{
	double x,y;
	P(){}
	P(double x,double y):x(x),y(y){}
	P operator + (P p){
		return P(add(x,p.x),add(y,p.y));
	}
	P operator - (P p){
		return P(add(x,p.x),add(y,-p.y));
	}
	P operator * (P p){
		return P(x*p.x,y*p.y);
	}
};
double dist(P a,P b){
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int cnt(P pos,vector<P> dat){
	int res=0;
	REP(i,dat.size()){
		if(dist(pos,dat[i])<=1.0+EPS)res++;
	}
	return res;
}
int main(){
	int n;
	double x,y;
	while(scanf("%d",&n)&&n){
		vector<P> dat;
		REP(i,n){
			cin>>x>>y;
			dat.PB(P(x,y));
		}
		int ans=0;
		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				P a=dat[i],b=dat[j];
				P mid=P((a.x+b.x)/2.0,(a.y+b.y)/2.0);
				double dx=a.y-b.y,dy=b.x-a.x;
				double d=dist(a,b);
				if(d>2)continue;
				double dd=sqrt(dx*dx+dy*dy);
				double ddd=sqrt(1-d*d/4);
				dx=dx/dd*ddd;dy=dy/dd*ddd;
				ans=max(ans,cnt(P(mid.x+dx,mid.y+dy),dat));
				ans=max(ans,cnt(P(mid.x-dx,mid.y-dy),dat));
			}
		}
		printf("%d\n",ans);
	}
	return 0;
}