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<double> getAreaPointx(int x,int y,vector<double> lx,vector<double> ly)
{
	vector<double> result;
	vector<double>::iterator itx=lx.begin();
	vector<double>::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<double> getAreaPointy(int x,int y,vector<double> lx,vector<double> ly)
{
	vector<double> result;
	vector<double>::iterator itx=lx.begin();
	vector<double>::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;
}


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

int countInside(double x,double y,vector<double> lx,vector<double> ly)
{
	int c=0;
	vector<double>::iterator itx=lx.begin();
	vector<double>::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<double> x,y;
		int map[12][12]={{0}};
		int maxx=0,maxy=0;
		int maxpoint=0;
		for (int i=0; i<n; i++)
		{
			double tx=0,ty=0;
			scanf("%lf %lf", &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 (double dif=0.8; dif>=0.008; dif/=10)
		{
			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<double> lx=getAreaPointx(i,j,x,y);
						vector<double> ly=getAreaPointy(i,j,x,y);
						for (double k=0; k<1; k+=dif)
						{
							for (double 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<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
//#include<cctype>
#include<climits>
#include<iostream>
#include<string>
#include<vector>
#include<map>
//#include<list>
#include<queue>
#include<deque>
#include<algorithm>
//#include<numeric>
#include<utility>
//#include<memory>
#include<functional>
#include<cassert>
#include<set>
#include<stack>
#include<random>

const int dx[] = {1, 0, -1, 0};
const int dy[] = {0, 1, 0, -1};
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef vector<int> vi;
typedef vector<ll> vll;
typedef pair<int, int> pii;
typedef long double Real;

Real eps = 1e-9;

Real add(Real a, Real b) {
    if (abs(a+b) < eps * (abs(a)+abs(b))) return 0;
    return a+b;
}

bool equal(Real a, Real b) {
    return add(a, -b) == 0;
}

struct P {
    Real x, y;
    P() {}
    P(Real x, Real y) : x(x), y(y) {}
    P operator+(P p) const {return P(add(x, p.x), add(y, p.y));}
    P operator-(P p) const {return P(add(x, -p.x), add(y, -p.y));}
    P operator*(Real d) const {return P(x*d, y*d);}
    Real dot(P p) const {return add(x*p.x, y*p.y);} // ??????
    Real det(P p) const {return add(x*p.y, -y*p.x);} // ??????
    Real dist(P p) const {return sqrt((x-p.x)*(x-p.x) + (y-p.y)*(y-p.y));} // ?????¢
    void normalize() {Real d = sqrt(x*x+y*y); x /= d; y /= d;} // ??£??????
    bool operator<(const P& rhs) const {
        if (x != rhs.x) return x < rhs.x;
        return y < rhs.y;
    }
    bool operator==(const P& rhs) const {
        return equal(x, rhs.x) && equal(y, rhs.y);
    }
};
const int MAXN = 333;
P pnt[MAXN];
int N;

int calc(P center) {
	int ret = 0;
	for (int i = 0; i < N; i++) {
		P vec = pnt[i]-center;
		if (vec.dot(vec) < 1+eps) {
			ret++;
		}
	}
	return ret;
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    while (cin >> N) {
    	if (N==0) break;
    	for (int i = 0; i < N; i++) 
    		cin >> pnt[i].x >> pnt[i].y;
    	int ans = 1;
    	for (int i = 0; i < N; i++) for (int j = i+1; j < N; j++) {
    		if (pnt[i].dist(pnt[j]) > 2) continue;
    		// diff ?????????
    		P vec = pnt[j]-pnt[i];
    		// vec ????????´????????????????????????
    		P n = P(vec.y, -vec.x);
    		n.normalize();
    		vec = vec*0.5;
    		Real len = sqrt(1-vec.dot(vec));
    		P cand1 = pnt[i] + vec + n*len;
    		P cand2 = pnt[i] + vec - n*len;
    		ans = max(ans, calc(cand1));
    		ans = max(ans, calc(cand2));
    	}
    	cout << ans << endl;
    }
    return 0;
}

C++: #include <stdio.h>
#include <algorithm>
#include <iostream>
#include <stack>
#include <vector>
#include <math.h>
#include <queue>

using namespace std;

struct Point{
    double x,y;
};

void func(int N){

    Point point[N];
    double x_min = 100,y_min = 100,x_max = -1,y_max = -1;

    for(int i = 0; i < N; i++){
        scanf("%lf %lf",&point[i].x,&point[i].y);
        x_min = min(point[i].x,x_min);
        y_min = min(point[i].y,y_min);
        x_max = max(point[i].x,x_max);
        y_max = max(point[i].y,y_max);
    }

    if(x_min > 1)x_min -= 0.04;
    if(y_min > 1)y_min -= 0.04;
    if(x_max < 9)x_max += 0.04;
    if(y_max < 9)y_max += 0.04;

    int count,max_count = 0;

    for(double y = y_min; y <= y_max; y += 0.01){
        for(double x = x_min; x <= x_max; x += 0.01){
            count = 0;
            for(int i = 0; i < N; i++){
                if(sqrt((x-point[i].x)*(x-point[i].x)+(y-point[i].y)*(y-point[i].y)) <= 1.0)count++;
            }
            max_count = max(max_count,count);
        }
    }

    printf("%d\n",max_count);
}


int main(){

    int N;

    while(true){
        scanf("%d",&N);
        if(N == 0)break;

        func(N);
    }


    return 0;
}