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<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.0001) counter++;
  }
  // printf("%d\n",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 = 1;
    for (int i = 0; i < N; i++)
      for (int j = 0; j < i; j++)
        ans = max(ans, check(i, j));
    printf("%d\n", ans);
  }
}

C++: #include <iostream>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <vector>
#include <complex>
using namespace std;

#define REP(i, a, b) for (int i = (a); i < (b); i++)
#define rep(i, n) REP(i, 0, n)

typedef complex<double> P;

#define EPS (1e-10)
#define EQ(a, b) (abs((a)-(b)) < EPS)
#define EQV(a, b) (EQ((a).real(), (b).real) && EQ((a).imag(), (b).imag()))

P unit_vector(P a) { return a/abs(a); }
P normal_vector1(P a) { return a * P(0, 1); }
P normal_vector2(P a) { return a * P(0, -1); }
double distance_p_p(P a, P b) { return abs(a-b); }

int N;

int main()
{
	while (cin >> N) {
		vector<P> ps;

		rep(i, N) {
			double x, y;
			cin >> x >> y;
			ps.push_back(P(x, y));
		}

		int ans = 0;
		rep(i, N) {
			REP(j, i+1, N) {
				P p, q;
				p = ps[i]; q = ps[j];

				double d = distance_p_p(p, q);
				if (d > 2) continue;

				P pq_c = (q - p)/2.0;
				P a = unit_vector(normal_vector1(q-p));
				P b = unit_vector(normal_vector2(q-p));

				P c1 = (p + pq_c) + a*(sqrt(1.0-d*d/4.0));
				P c2 = (p + pq_c) + b*(sqrt(1.0-d*d/4.0));

				int buf1, buf2;
				buf1 = buf2 = 2;
				rep(k, N) {
					if (k == i || k == j) continue;
					if (distance_p_p(ps[k], c1) <= 1)
						buf1++;
					if (distance_p_p(ps[k], c2) <= 1)
						buf2++;
				}
				ans = max(ans, max(buf1, buf2));
			}
		}
		cout << ans << endl;
	}
	return 0;
}

C++: #include <stdio.h>
#include <cmath>
#include <vector>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <string>
using namespace std;

vector<float> getAreaPointx(int x,int y,vector<float> lx,vector<float> ly)
{
	vector<float> result;
	vector<float>::iterator itx=lx.begin();
	vector<float>::iterator ity=ly.begin();
	while (itx != lx.end())
	{
		if ((*itx <= x+1) && (*itx >= x-2) && (*ity <= y+1) && (*ity >= y-2)){
			result.push_back(*itx);
		}
		itx++;
		ity++;
	}
	return result;
}

vector<float> getAreaPointy(int x,int y,vector<float> lx,vector<float> ly)
{
	vector<float> result;
	vector<float>::iterator itx=lx.begin();
	vector<float>::iterator ity=ly.begin();
	while (itx != lx.end())
	{
		if ((*itx <= x+1) && (*itx >= x-2) && (*ity <= y+1) && (*ity >= y-2)){
			result.push_back(*ity);
		}
		itx++;
		ity++;
	}
	return result;
}


float distance(float x1,float y1,float x2,float y2)
{
	return sqrt(pow((x1-x2),2)+pow((y1-y2),2));
}

int countInside(float x,float y,vector<float> lx,vector<float> ly)
{
	int c=0;
	vector<float>::iterator itx=lx.begin();
	vector<float>::iterator ity=ly.begin();
	while (itx != lx.end())
	{
		if (distance(x,y,*itx,*ity) <= 1)
		{
			c++;
		}
		itx++;
		ity++;
	}
	return c;
}

int main(void) {
	vector<int> r;
	while (true)
	{
		int n=0;
		scanf("%d",&n);
		if (n==0){
			break;
		}
		vector<float> x,y;
		int map[12][12]={{0}};
		int maxx=0,maxy=0;
		int maxpoint=0;
		for (int i=0; i<n; i++)
		{
			float tx=0,ty=0;
			scanf("%f %f", &tx, &ty);
			x.push_back(tx);
			y.push_back(ty);
			int ttx=0,tty=0;
			ttx = int(round(tx+0.5));
			tty = int(round(ty+0.5));
			map[ttx][tty]++;
			if (maxpoint < map[ttx][tty]){
				maxpoint = map[ttx][tty];
				maxx=ttx;
				maxy=tty;
			}
		}
		for (float dif=0.5; dif>=0.005; dif/=100)
		{
			for (int i=1; i<=10; i++)
			{
				for (int j=1; j<=10; j++)
				{
					int aroundpoint=0;
					for (int k=-1; k<=1; k++)
					{
						for (int l=-1; l<=1; l++)
						{
							aroundpoint+=map[i+k][j+l];
						}
					}
					if (maxpoint >= aroundpoint){
						continue;
					} else {
						vector<float> lx=getAreaPointx(i,j,x,y);
						vector<float> ly=getAreaPointy(i,j,x,y);
						for (float k=0; k<1; k+=dif)
						{
							for (float l=0; l<1; l+=dif)
							{
								int tmp = countInside(i-1+k,j-1+l,lx,ly);
								if (tmp>maxpoint){
									maxpoint = tmp;
								}
							}
						}
					}
				}
			}
		}
		r.push_back(maxpoint);
	}
	for(int i=0; i<r.size(); i++)
	{
		cout << r[i] << endl;
	}
	return 0;
}

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

complex<double> p[300];

complex<double> core(int a,int b){
  complex<double> res,mid,norm;

  mid = (p[a]+p[b]) / 2.0;

  norm = (p[a] - p[b]) * complex<double>(0,1);
  norm = norm / abs(norm);
  norm = norm * sqrt(1.0 - abs(mid-p[a])*abs(mid-p[a]));

  res = norm + mid;
  return res;
}
  

int main(){
  int n;
  complex<double> o,r;
  int max,tmp;

  while(1){
    cin >> n;
    if(!n)break;

    for(int i=0;i<n;i++)cin >> p[i].real() >> p[i].imag();

    max = 1;
    for(int i=0;i<n;i++){
      for(int j=0;j<n;j++){
	if(i!=j && abs(p[i]-p[j]) < 2.0){
	  o = core(i,j);
	  tmp = 2;
	  for(int k=0;k<n;k++){
	    if(k!=i && k!=j){
	      r = o-p[k];
	      if( r.real()*r.real() + r.imag()*r.imag() <1.0)tmp++;
	    }
	  }
	  if(max < tmp)max = tmp;
	}
      }
    }
    cout << max << endl;
  }
}