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<complex>
#include<vector>
using namespace std;
typedef complex<double> P;
int n;
vector<P> p;
int ret(){
	int cnt=0;
	for(int i=0;i<n;i++){
		for(int j=i+1;j<n;j++){
			if(abs(p[i]-p[j])>2.0)
				continue;
			P c=(p[i]+p[j])*0.5;
			double d=sqrt(1.0-abs(p[i]-c)*abs(p[i]-c));
			P l1=(p[i]-p[j])*P(0,1);
			P l2=(p[i]-p[j])*P(0,-1);
			l1=l1*d/abs(l1)+c;l2=l2*d/abs(l2)+c;
			int x=0,y=0;
			for(int k=0;k<n;k++){
				if(abs(l1-p[k])<=1.0){
					x++;
				}
				if(abs(l2-p[k])<=1.0){
					y++;
				}
				cnt=max(cnt,max(x,y));
			}
		}
	}
	return cnt;
	
}
int main(){
	while(cin>>n,n){
		p.resize(n);
		for(int i=0;i<n;i++){
			double x,y;
			cin>>x>>y;
			p[i]=P(x,y);
		}
		cout<<ret()<<endl;
	}
	return 0;
}

C++: #include <iostream>
#include <complex>
#include <vector>

using namespace std;


class Circle {
public:
	complex<double> p;
	double r;
	Circle (complex<double> p0, double r0) {
		p = p0;
		r = r0;
	}
};

vector<complex<double> > intersectCC(Circle a, Circle b) {
	vector<complex<double> > res(2);

	auto d = abs(a.p - b.p);
	auto c = (d*d + a.r*a.r - b.r*b.r) / (2*d);
	auto s = sqrt(a.r*a.r - c*c);
	auto diff = (b.p - a.p) / d;

	res[0] = a.p + diff * complex<double>(c, s);
	res[1] = a.p + diff * complex<double>(c, -s);

	return res;
}


bool containsCP(Circle c, complex<double> p) {
	const double EPS = 10e-7;
	auto d = abs(c.p - p);
	return d <= c.r + EPS;
}


int main()
{
	int n;

	while (cin >> n, n) {
		vector<complex<double> > ps(n);

		for (int i = 0; i < n; i++) {
			double x[2];
			cin >> x[0] >> x[1];
			ps[i] = complex<double>(x[0], x[1]);
		}

		auto res = 1;

		for (int i = 0; i < n; i++) {
			for (int j = i+1; j < n; j++) {
				if (abs(ps[i]-ps[j]) > 2) continue;
				auto cp = intersectCC(Circle(ps[i], 1.0), Circle(ps[j], 1.0));
				int cnt = 0;
				
				for (int k = 0; k < n; k++) {
					if (containsCP(Circle(cp[0], 1.0), ps[k])) cnt++;
				}

				res = max(cnt, res);
				cnt = 0;

				for (int k = 0; k < n; k++) {
					if (containsCP(Circle(cp[1], 1.0), ps[k])) cnt++;
				}
				res = max(cnt, res);
			}
		}
		cout << res << endl;
	}

	return 0;
}

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 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];
			}
		}
		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+=0.01)
					{
						for (double l=0; l<1; l+=0.01)
						{
							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 <vector>

using namespace std;

#define x first
#define y second
const int N = 300;
const int SCALE = 1e5;
const int radius = 1;
const double e = 1e-6;
typedef pair<double, double> point;
point points[N];
double distances[N][N];

double sq(double d) {
  return d * d;
}

double dist(point a, point b) {
  return sqrt(sq(a.x - b.x) + sq(a.y - b.y));
}

bool le(double a, double b){
  return a < b or fabs(a - b) < e;
}

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

    for (int i = 0; i < n; ++i)
      for (int j = 0; j < n; ++j)
        distances[i][j] = dist(points[i], points[j]);

    // for each point,
    //  for every other point
    //    for every point currently selected
    //      check if it goes past the bounds
    //    add the point if its contained
    //
    int m = 0;
    for (int i = 0; i < n; ++i) {
      vector<int> selected = {i};
      for (int j = 0; j < n; ++j) {
        if (i == j) continue;
        bool contained = true;
        for (auto k : selected) {
          if (distances[j][k] > 2.0) {
            contained = false;
            break;
          }
        }
        if (contained) selected.push_back(j);
        /*
        cout << i << ' ' << j << ": " << endl;
        for (auto p : selected) cout << p << ' ';
        cout << endl;
        */
      }

      m = max(m, (int)selected.size());
    }

    cout << m << endl;

    cin >> n;
  }

}