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>
#define eps 1e-9
using namespace std;
typedef complex<double> point;

pair<point, point> norm(point p){
  double unit = 1 / sqrt(p.real()*p.real() + p.imag()*p.imag());
  return pair<point, point>(p*point(0, 1)/unit,
                            p*point(0, -1)/unit);
}

vector<point> calcC(point a, point b){
  point p = (b - a)/2.0;
  pair<point, point> n = norm(p);
  double unit = sqrt(1 - abs(p)*abs(p));
  n.first  =  n.first*unit;
  n.second = n.second*unit;
  vector<point> res(2);
  res[0] = a + p + n.first;
  res[1] = a + p + n.second;
  return res;
}
int main(int argc, char *argv[]){
  int n;
  while(cin >> n, n){;
    point p[n];
    for (int i = 0; i < n; i++) {
      double real, imag;
      cin >> real >> imag;
      p[i] = complex<double>(real, imag);
    }
    int ans = 1, cnt = 0;
    for (int i = 0; i < n; i++) {
      for (int j = i + 1; j < n; j++) {
        point a = p[i], b = p[j];
        if(abs(b - a) > 2 - eps)continue;
        vector<point> center = calcC(a, b);
        for (int c = 0; c < 2; c++) {
          cnt = 0;
          // std::cout << center[c].real() << " " <<
          //   center[c].imag()<< std::endl;
          for (int k = 0; k < n; k++) {
            if(abs(center[c] - p[k]) <= 1 + eps)cnt++;
          }
          ans = max(ans, cnt);
        }
      }
    }
    std::cout << ans << std::endl;
  }
  return 0;
}

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

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,double lx[],double ly[],int n)
{
	int c=0;
	for (int i=0; i<n; i++)
	{
		if (distance(x,y,lx[i],ly[i]) <= 1)
		{
			c++;
		}
	}
	return c;
}

int main(void) {
	vector<int> r;
	while (true)
	{
		int n=0;
		scanf("%d",&n);
		if (n==0){
			break;
		}
		double x[300]={0},y[300]={0};
		int map[12][12]={{0}};
		int maxpoint=0;
		for (int i=0; i<n; i++)
		{
			scanf("%lf %lf", &x[i], &y[i]);
			int tx=0,ty=0;
			tx = int(round(x[i]+0.5));
			ty = int(round(y[i]+0.5));
			map[tx][ty]++;
			if (maxpoint < map[tx][ty]){
				maxpoint = map[tx][ty];
			}
		}
		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 {
					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,x,y,n);
							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>
using namespace std;
#include <math.h>
int main(){
	int n=0;
	while(cin>>n && n!=0){
		double x[300]={},y[300]={};
		int ans=1;
		for(int i=0;i<n;i++){
			cin>>x[i]>>y[i];
		}
		for(int i=0;i<n;i++){
			for(int j=0;j<n;j++){
				if(i!=j){
					double dx=x[i]-x[j],dy=y[i]-y[j];
					double dd=dx*dx+dy*dy;
					if(dd<4.0){
						double h=pow(1.0-dd/4.0,1.0/2.0),d=pow(dd,1.0/2.0);
						double px=(x[i]+x[j])/2.0+h*dy/d,py=(y[i]+y[j])/2.0-h*dx/d;
						int a=0;
						for(int k=0;k<n;k++){
							if((x[k]-px)*(x[k]-px)+(y[k]-py)*(y[k]-py)<1.0001) a++;
						}
						if(a>ans)ans=a;
					}
				}
			}
		}
		cout<<ans<<endl;
	}
}