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 rep2(i, b, n) for(int i=b; i < n; ++i)
#define rep(i, n) rep2(i, 0, n)
#define X first
#define Y second

using namespace std;

typedef pair<double, double> pdd;

double distance_sq( pdd a, pdd b )
{
	double cx = a.X - b.X;
	double cy = a.Y - b.Y;
	return ((cx*cx) + (cy*cy));
}

double norm( pdd &v )
{
	double len = sqrt((v.X*v.X) + (v.Y*v.Y));
	v.X /= len;
	v.Y /= len;
	return len;
}

void rotate( pdd &p )
{
	swap(p.X, p.Y);
	p.Y = -p.Y;
}

int main( void )
{
	int n;
	pdd p[300];
	while(cin >> n, n)
	{
		rep(i, n) cin >> p[i].X >> p[i].Y;
		int Max = 0;
		rep(i, n) {
			rep2(j, i+1, n) {
				if(distance_sq(p[i], p[j]) < 4.0) {
					pdd v = make_pair(p[i].X-p[j].X, p[i].Y-p[j].Y);
					pdd c = make_pair(v.X / 2.0000, v.Y / 2.0000);
					pdd m = c;
					double len_c = norm(c);
					double sin_x = sqrt(1.0000 - (len_c*len_c));
					rotate(c);
					c.X *= sin_x; c.Y *= sin_x;
					//p[j].X += m.X; p[j].Y += m.Y;
					pdd p3 = make_pair(p[j].X+m.X+c.X, p[j].Y+m.Y+c.Y);
					pdd p4 = make_pair(p[j].X+m.X-c.X, p[j].Y+m.Y-c.Y);
					int count_p3, count_p4;
					count_p3 = count_p4 = 0;
					rep(k, n) {
						if(distance_sq(p3, p[k]) < 1.000001) ++count_p3;
						if(distance_sq(p4, p[k]) < 1.000001) ++count_p4;
					}
					Max = max(Max, max(count_p3, count_p4));
				}
			}
		}
		cout << Max << endl;
	}
	return 0;
}

C++: /**
 *  Circle and Points (PKU 1981)
 *  by Kenji Inoue, 2007-01-05
 **/

#include <iostream>
#include <math.h>

using namespace std;

#define EPS (1e-6)

double x[300], y[300];

int sign[] = {+1, -1};

int main(void) {
  int n;
  while (cin >> n, n) {
    for (int i=0; i<n; i++) {
      cin >> x[i] >> y[i];
    }

    int ans = 1;
    for (int i=0; i<n; i++) {
      for (int j=i+1; j<n; j++) {
        double d = hypot(x[i]-x[j], y[i]-y[j]);
        if ( d > 2.0 + EPS ) { continue; }

        double mx = (x[i] + x[j]) / 2;
        double my = (y[i] + y[j]) / 2;
        double vx = (x[j] - x[i]) / d;
        double vy = (y[j] - y[i]) / d;
        double e = sqrt(1.0 - d*d/4.0);

        for (int s=0; s<2; s++) {
          double px = mx + sign[s]*e*vy;
          double py = my - sign[s]*e*vx;

		  printf("%f,%f\n",px,py);

          int cnt = 2;
          for (int k=0; k<n; k++) {
            if (k==i || k==j) { continue; }
            if ( (px-x[k])*(px-x[k]) + (py-y[k])*(py-y[k]) < 1.0 + EPS ) { cnt++; }
          }
          if (cnt > ans) { ans = cnt; }
        }
      }
    }

    cout << ans << endl;
  }
  return 0;
}

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

using namespace std;

class Point
{
private:

	double x_;
	double y_;

public:

	Point() : x_(0), y_(0) { }

	explicit Point(double x, double y) : x_(x), y_(y) { }

	explicit Point(double x1, double y1, double x2, double y2) : x_(x2 - x1), y_(y2 - y1) { }

	double x() { return x_; }
	double y() { return y_; }

	double x() const { return x_; }
	double y() const { return y_; }

	Point& operator+= (const Point& p) { x_ += p.x_; y_ += p.y_; return *this; }
	Point& operator-= (const Point& p) { x_ -= p.x_; y_ -= p.y_; return *this; }

	Point& operator*= (double s) { x_ *= s; y_ *= s; return *this; }
	Point& operator/= (double s) { x_ /= s; y_ /= s; return *this; }

	friend bool operator== (const Point& p1, const Point& p2) { return p1.x_ == p2.x_ && p1.y_ == p2.y_; }
	friend bool operator!= (const Point& p1, const Point& p2) { return !(p1 == p2); }

	friend Point operator+ (const Point& p1, const Point& p2) { return Point(p1.x_ + p2.x_, p1.y_ + p2.y_); }
	friend Point operator- (const Point& p1, const Point& p2) { return Point(p1.x_ - p2.x_, p1.y_ - p2.y_); }

	friend Point operator* (const Point& p, double s) { return Point(p.x_ * s, p.y_ * s); }
	friend Point operator/ (const Point& p, double s) { return Point(p.x_ / s, p.y_ / s); }

	friend Point operator* (double s, const Point& p) { return Point(p.x_ * s, p.y_ * s); }
};

double distance(const Point& p1, const Point& p2)
{
	return sqrt(pow(p2.x() - p1.x(), 2) + pow(p2.y() - p1.y(), 2));
}

int main()
{
	int N; double X, Y;

	vector<Point> Plane;

	while (true)
	{
		cin >> N;

		if (N == 0) { break; }

		Plane.clear();

		for (int i = 0; i < N; i++)
		{
			cin >> X >> Y;

			Plane.push_back(Point(X, Y));
		}

		int max_ = 0;

		int count_;

		for (int i = 0; i < N; i++)
		{
			for (int j = i + 1; j < N; j++)
			{
				count_ = 1;

				if (distance(Plane[i], Plane[j]) < 2)
				{
					double vx = Plane[j].x() - Plane[i].x();
					double vy = Plane[j].y() - Plane[i].y();

					swap(vx, vy); vx = -vx;

					Point R(vx, vy);

					count_ = 2;

					for (int k = 0; k < N; k++)
					{
						if (k != i && k != j)
						{
							if (distance(Plane[i], R) < 1)
							{
								count_++;
							}
						}
					}
				}

				max_ = max(count_, max_);

				count_ = 0;

				if (distance(Plane[i], Plane[j]) < 2)
				{
					double vx = Plane[j].x() - Plane[i].x();
					double vy = Plane[j].y() - Plane[i].y();

					swap(vx, vy); vy = -vy;

					Point R(vx, vy);

					count_ = 2;

					for (int k = 0; k < N; k++)
					{
						if (k != i && k != j)
						{
							if (distance(Plane[i], R) < 1)
							{
								count_++;
							}
						}
					}
				}

				max_ = max(count_, max_);
			}
		}

		cout << max_ << endl;
	}

	return 0;
}

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

#define MAX_N 310

long double x[MAX_N], y[MAX_N];
int n, maxn;

int main() {
	while (true) {
		cin >> n;
		if (n == 0) { break; }
		for (int i = 0; i < n; i++) {
			cin >> x[i] >> y[i];
		}
		maxn = 0;
		for (int i = 0; i < n; i++) {
			for (int j = i + 1; j < n; j++) {
				int sum = 0;
				long double X = (x[i] + x[j]) / 2.0;
				long double Y = (y[i] + y[j]) / 2.0;
				for (int k = 0; k < n; k++) {
					if (((x[k] - X)*(x[k] - X) + (y[k] - Y)*(y[k] - Y)) <= 1) {
						sum++;
					}
				}
				maxn = max(maxn, sum);
			}
		}
		for (int i = 0; i < n; i++) {
			int sum = 0;
			long double X = x[i];
			long double Y = y[i];
			for (int k = 0; k < n; k++) {
				if (((x[k] - X)*(x[k] - X) + (y[k] - Y)*(y[k] - Y)) <= 1) {
					sum++;
				}
			}
			maxn = max(maxn, sum);
		}
		cout << maxn << endl;
	}
	return 0;
}