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>
#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-5;
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=1;
		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;
}

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

using namespace std;

const double r = 1.0;
const double PI = acos(-1.0);

struct Point {
    double x, y;
    Point() {}
    Point(double x, double y) : x(x), y(y) {}
};
Point operator+(const Point &a, const Point &b) {
    return Point(a.x + b.x, a.y + b.y);
}
Point operator-(const Point &a, const Point &b) {
    return Point(a.x - b.x, a.y - b.y);
}
inline double Distance2(const Point &a, const Point &b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
Point Circle(const Point &p1, const Point &p2, double theta) {
    Point p12 = p2 - p1;
    Point rp(p12.x * cos(theta) - p12.y * sin(theta), p12.x * sin(theta) + p12.y * cos(theta));
    double len = sqrt(p12.x * p12.x + p12.y * p12.y);

    rp.x /= len;
    rp.y /= len;

    return rp + p1;
}

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

    int n;

    while (cin >> n, n) {
        // Input
        vector<Point> p(n);
        for (int i = 0; i < n; ++i)
            cin >> p[i].x >> p[i].y;

        // Solve
        int ans = 1;
        for (int i = 0; i < n; ++i)
            for (int j = i + 1; j < n; ++j) {
                double d = Distance2(p[i], p[j]);

                if (d <= 4.0 * r * r) {
                    double theta = (PI - acos((2.0 - d) * 0.5)) * 0.5;

                    for (int s = -1; s <= 1; s += 2) {
                        Point c = Circle(p[i], p[j], s * theta);

                        int num = 2;
                        for (int k = 0; k < n; ++k) {
                            if (k != i && k != j && Distance2(c, p[k]) <= 1.0)
                                ++num;
                        }

                        ans = max(ans, num);
                    }
                }
            }

        // Output
        cout << ans << endl;
    }

    return 0;
}

C++: #include<cstdio>
#include<algorithm>
#include<math.h>

using namespace std;
typedef pair<double, double> P;

int N;
double X[300], Y[300];
P d[300];

P add(P p1, P p2) {
  return P(p1.first + p2.first, p1.second + p2.second);
}

P minas(P p1, P p2) {
  return P(p1.first - p2.first, p1.second - p2.second);
}

P mul(P p, double c) {
  return P(p.first * c, p.second * c);
}

P ave(P p1, P p2) {
  return P((p1.first + p2.first) / 2, (p1.second + p2.second) / 2);
}

void printP(P p) {
  printf("P:(%f, %f)\n", p.first, p.second);
}

double square(P p) {
  return p.first * p.first + p.second * p.second;
}

double size(P p) {
  return sqrt(square(p));
}

int count(P p) {
  // printP(p);
  int counter = 0;
  for (int i = 0; i < N; i++) {
    if (square(minas(p, d[i])) <= 1.00001) counter++;
  }
  return counter;
}

int check(int i, int j) {
  if (square(minas(d[i], d[j])) > 4.0) return 0;
  int ans = 0;
  P average = ave(d[i], d[j]);
  P dir = mul(minas(d[i], d[j]), 0.5);
  // printP(average);
  double s = sqrt(1.0 - square(dir));
  // printf("!!!!!!%f\n",s);
  P dir1 = P((Y[i] - Y[j]), - (X[i] - X[j]));
  //printf("!!!!!!%f\n",size(dir1));
  double r = s / size(dir1);
  P d1 = add(average, mul(dir1, r));
  P d2 = add(average, mul(dir1, -r));
  return max(count(d1), count(d2));
}

int main() {
  while (1) {
    scanf("%d", &N);
    if (N == 0) break;
    for (int i = 0; i < N; i++) scanf("%lf%lf", &X[i], &Y[i]), d[i] = P(X[i], Y[i]);
    int ans = 0;
    for (int i = 0; i < N; i++)
      for (int j = 0; j < i; j++)
        ans = max(ans, check(i, j));
    printf("%d\n", ans);
  }
}