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;
typedef complex<double> P;
#define F real()
#define S imag()
#define EPS (1e-6)
int N;
P D[333];
struct Circle{
  double r;
  P c;
  Circle(){ r = 1.0; }
  Circle(P c1){
    c = c1;
    r = 1.0;
  }
};
double getK(P a,P b ){
  return sqrt((a.F-b.F)*(a.F-b.F)+(a.S-b.S)*(a.S-b.S));
}

int getC(P a){
  // cout << a.F << " " << a.S << endl;
  int ret = 0;
  for(int i=0;i<N;i++){
    //cout << " k " << getK(a,D[i]) << endl;
    if( 1.0+EPS >= getK(a,D[i]) ) ret++;
  }
  return ret;
}

double getAngle( P a ,P b ){
  return atan((b.S-a.S)/(b.F-a.F));
}

vector<P> getIntersectCC(Circle C1,Circle C2){
  vector<P> res;
  double r1 = C1.r, r2 = C2.r, d;
  P p1 = C1.c, p2 = C2.c;
  d = abs(p1-p2);
  
  if( d < EPS && abs( r1 - r2 ) < EPS ) return res;
  if( r1 + r2 < d - EPS || d + EPS < abs(r1-r2) ) return res;
  
  double a = (r1*r2 - r2*r2 + d*d ) / (2*d);
  double h = sqrt( max( r1 * r1 - a*a , 0.0 ) );
  P tmp1 = p1 + a/d * (p2-p1);
  P tmp2 = h/d * (p2-p1);
  
  if( abs(tmp2) < EPS) res.push_back(tmp1);
  else {
    res.push_back( P(real(tmp1) - imag(tmp2), imag(tmp1) + real(tmp2)));
    res.push_back( P(real(tmp1) + imag(tmp2), imag(tmp1) - real(tmp2)));
  }
}

vector<P> getP(P a, P b){
  return getIntersectCC( Circle(a), Circle(b) );
}

int main(){
  while( cin >> N && N ){
    for(int i=0;i<N;i++) {
      cin >> D[i].real() >> D[i].imag();
    }
    int res = 1;
    for(int i=0;i<N;i++){
      for(int j=0;j<N;j++){
	if( i==j ) continue;
	P r,l;
	vector<P> p = getP( D[i],D[j]);
	for(int i=0;i<(int)p.size();i++){
	  res = max( res, getC(p[i]) );
	  res = max( res, getC(p[i]) );
	}
	//res = max( res, solve(i) );
      }
    }
    cout << res << endl;
  }
}

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

				long double X1, Y1, X2, Y2, X3, Y3, VX2, VY2, Dist, Dist2;

				Dist = sqrtl((x[j] - x[i])*(x[j] - x[i]) + (y[j] - y[i])*(y[j] - y[i]));
				Dist2 = sqrtl(4 - Dist*Dist);
				X2 = x[j] - x[i]; Y2 = y[j] - y[i];
				VX2 = (X2 / Dist)*Dist2 / 2; VY2 = (Y2 / Dist)*Dist2 / 2;
				X1 = X + VY2; Y1 = Y - VX2;
				X3 = X - VY2; Y3 = Y + VX2;

				sum = 0;
				for (int k = 0; k < n; k++) {
					if (((x[k] - X1)*(x[k] - X1) + (y[k] - Y1)*(y[k] - Y1)) <= 1) {
						sum++;
					}
				}
				maxn = max(maxn, sum);
				sum = 0;
				for (int k = 0; k < n; k++) {
					if (((x[k] - X3)*(x[k] - X3) + (y[k] - Y3)*(y[k] - Y3)) <= 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;
}

C++: #include <bits/stdc++.h>
using namespace std;

int main(void){
	int n;
	while(cin >> n, n){
		vector<double> x(n), y(n);
		for (int i = 0; i < n; ++i){
			cin >> x[i] >> y[i];
		}
		int ans = 0;
		for (int i = 0; i < n; ++i){
			for(int j = i+1; j < n; ++j){
				double rx[2], ry[2], a, ax, ay, bx, by;
				ax = (x[j] - x[i]) / 2.0;
				ay = (y[j] - y[i]) / 2.0;
				a = sqrtl(ax*ax + ay*ay);
				by = ax * sqrtl(1.0-a*a) / a;
				bx = ay * sqrtl(1.0-a*a) / a;
				rx[0] = x[i] + ax + bx; rx[1] = x[i] + ax - bx;
				ry[0] = y[i] + ay - by; ry[1] = y[i] + ay + by;
				for(int k = 0; k < 2; ++k){
					int sum = 0;
					//cout << "[" << rx[k] << " " << ry[k] << "]" << endl;
					for(int l = 0; l < n; ++l){
						//cout << (x[l]-rx[k])*(x[l]-rx[k])+(y[l]-ry[k])*(y[l]-ry[k]) << endl;
						if((x[l]-rx[k])*(x[l]-rx[k])+(y[l]-ry[k])*(y[l]-ry[k]) <= 1.0001){
							sum++;
							//cout << x[l] << " " << y[l] << endl;
						}
					}
					ans = max(ans, sum);
				}
			}
		}
		cout << (ans?ans:1) << endl;
	}
	return 0;
}

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 = 1;
		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;
					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 = 2;
					rep(k, n) {
						if(i==k || j==k) continue;
						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;
}