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;

const double r = 1.0;
const double PI = acos(-1.0);
const double sq2 = sqrt(2.0);
const double dx[] = {0.0, 1.0, 0.0, 1.0};
const double dy[] = {0.0, 0.0, 1.0, 1.0};

struct Point {
    double x, y;
    Point() {}
    Point(double x, double y) : x(x), y(y) {}
};
struct Rect {
    double x, y, len; // ??????(x, y)??§?????????len?????£?????¢
    int ub;
    Rect() { ub = 0; }
    Rect(double x, double y, double len, int ub) : x(x), y(y), len(len), ub(ub) {}
};
bool operator<(const Rect &a, const Rect &b) {
    return a.ub < b.ub;
}
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);
}

int n;
vector<Point> p;

Rect MakeRect(double x, double y, double l, int &ans) {
    int num = 0, u_bound = 0;
    const double len = l * l * 0.5 + sq2 * l + 1.0;
    Point c(x + l * 0.5, y + l * 0.5);

    for (int i = 0; i < n; ++i) {
        double d = Distance2(c, p[i]);

        if (d <= 1.0)
            ++num;
        if (d <= len)
            ++u_bound;
    }

    ans = max(ans, num);

    return Rect(x, y, l, u_bound);
}

int Solve() {
    int ans = 1;
    priority_queue<Rect> que;

    que.push(MakeRect(0.0, 0.0, 10.0, ans));

    while (!que.empty()) {
        Rect now = que.top();
        que.pop();

        if (ans < now.ub) {
            double half = now.len * 0.5;

            for (int d = 0; d < 4; ++d) {
                Rect r = MakeRect(now.x + half * dx[d], now.y + half * dy[d], half, ans);
                if (ans < r.ub)
                    que.push(r);
            }
        }
    }

    return ans;
}

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

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

        cout << Solve() << endl;
    }

    return 0;
}

C++: #include "bits/stdc++.h"
using namespace std;
typedef complex<double> Point;
const double eps = 1e-9;
int main(void) {
	cin.tie(0); ios::sync_with_stdio(false);
	int N; 
	while (cin >> N, N) {
		vector<Point> P(N);
		for (int i = 0; i < N;i++) {
			double x, y; cin >> x >> y;
			P[i] = Point(x, y);
		}
		if (N == 1) { cout << 1 << endl; continue; }
		int ans = 1;
		for (int i = 0; i < N;i++) for (int j = i + 1; j < N;j++) {
			double dist = abs(P[i] - P[j]);
			if (dist > 2) continue;
			double r = sqrt(1 - dist*dist / 4);
			double difx = -(P[i] - P[j]).imag()*r / dist;
			double dify = (P[i] - P[j]).real()*r / dist;
			Point new_P(difx, dify);
			Point check_P(new_P + (P[i] + P[j]) / 2.0);
			
			int cnt = 0;
			for (int k = 0; k < N;k++) {
				if (abs(check_P - P[k]) < 1.0 + eps)cnt++;
			}
			ans = max(ans, cnt);
		}
		cout << ans << endl;
	}
	return 0;
}

C++: #include <iostream>
#include <algorithm>
#include <cmath>
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;
  }
};

typedef Point Vector;

double dist(Point p1, Point p2) {
  return sqrt(pow(p1.x-p2.x, 2.0) + pow(p1.y-p2.y, 2.0));
}

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

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;
  double l = sqrt(4.0/norm(p1,p2)-1)/2.0;
  int cnt = 0;
  int newCnt;
  for(int t = 0; t < 2; t++) {
    newCnt = 0;
    moto = Vector(p2.y-p1.y, p1.x-p2.x);
    if(t == 0) {
      moto = moto * (-1);
    }
    o = m + moto * l;

    for(int i = 0; i < n; i++) {
      if(p[i] == p1 || p[i] == p2) continue;
      if(pow(p[i].x-o.x, 2.0) + pow(p[i].y-o.y, 2.0) <= 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;
    }
    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;
  }
}