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 <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.0049; 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 <algorithm>
#include <complex>
#include <cstdio>
#include <vector>

using namespace std;
#define reps(i,f,n) for(int i=f; i<int(n); ++i)
#define rep(i,n) reps(i,0,n)

typedef vector<int> vi;
typedef vector<vi> vvi;
typedef complex<double> P;

const double EPS = 1e-8;

int main()
{
	int n;
	while(scanf("%d", &n), n){
		P p[300];
		rep(i, n)
			scanf("%lf%lf", &p[i].real(), &p[i].imag());

		int ans = 1;
		rep(i, n) rep(j, i){
			double l = 1.0 - norm(p[i]-p[j])/4.0;
			if(l < 0.0)
				continue;
			l = sqrt(l);

			P center = (p[i] + p[j]) / 2.0;
			P ll = l * (p[i] - p[j])/abs(p[i] - p[j]) * P(0.0, 1.0);
			P c[] = {center+ll, center-ll};
			rep(k, 2){
				int cnt = 0;
				rep(m, n){
					if(norm(c[k] - p[m]) <= 1.0 + EPS)
						++cnt;
				}
				ans = max(ans, cnt);
			}
		}
		printf("%d\n", ans);
	}
	return 0;
}

C++: #include<bits/stdc++.h>
using namespace std;
typedef complex<double> xy_t;
int n, temp, ans;
xy_t p[300];
xy_t c[2];
// 内積: a.x*b.x +a.y*b.y
double dot_product(xy_t a, xy_t b) { return (conj(a)*b).real(); }

bool center(xy_t a, xy_t b){
    xy_t dir = b - a;
    if(dot_product(dir, dir) > 4){
        return false;
    }
    xy_t ver = xy_t(-dir.imag(), dir.real());
    ver /= sqrt(dot_product(ver, ver));
    double len = sqrt(1 - dot_product(dir, dir) / 4);
    c[0] = a + dir / 2.0 - ver * len;
    c[1] =  a + dir / 2.0 + ver * len;
    return true;
}
int main(){
    while(cin >> n && n){
        ans = 1; // 最低でも一つは囲める
        for(int i = 0; i < n; i++){
            double x, y;
            cin >> x >> y;
            p[i] = xy_t(x, y);
        }
        for(int i = 0; i < n; i++){
            for(int j = i + 1; j < n; j++){
                if(center(p[i], p[j])){
                    for(int l = 0; l < 2; l++){
                        temp = 2;
                        for(int k = 0; k < n; k++){
                            if(k == i || k == j){
                                continue;
                            }
                            xy_t vec = c[l] - p[k];
                            if(dot_product(vec, vec) <= 1){
                                temp++;
                            }
                        }
                        ans = max(ans, temp);
                    }
                }
            }
        }
        cout << ans << endl;
    }
}