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 <cmath>
#define EPS 1e-9
using namespace std;
typedef complex<double> xy_t;
int n;
xy_t points[301];
int count_in(xy_t center)
{
  int ans=0;
  for(int i=0;i<n;i++)
  {
    if (abs(center-points[i])<1.0+EPS)
      ans++;
  }
  return ans;
}

xy_t ppdcl(xy_t a)
{
  return xy_t(a.imag(),-a.real());
}
xy_t delta_v(xy_t C1,xy_t C2)
{
  xy_t half=(C1-C2)*0.5;
  double l=abs(half);
  return ppdcl(half)*sqrt(1-l*l)*(1.0/l);
}

int main()
{
  while(1)
  {
    cin>>n;
    if (n==0) break;
    double x,y;
    for(int i=0;i<n;i++)
    {
      cin>>x>>y;
      points[i]=xy_t(x,y);
    }
    int maxn=1;
    for(int i=0;i<n;i++)
    {
      for(int j=i;j<n;j++)
      {
        if(abs(points[i]-points[j])<2.0+EPS)
        {
          xy_t mid_p=(points[i]+points[j])*0.5;
          maxn=max(maxn,count_in(mid_p+delta_v(points[i],points[j])));
          maxn=max(maxn,count_in(mid_p-delta_v(points[i],points[j])));
        }
        /*double d=abs(points[j]-points[i])/2;
        xy_t v=(points[j]-points[i])/(2*d);
        if(d<1)
        {
            maxn=max(maxn,count_in(points[i]+d*v+sqrt(1-d*d)*v*xy_t(0,+1)));
            maxn=max(maxn,count_in(points[i]+d*v+sqrt(1-d*d)*v*xy_t(0,-1)));
        }*/
      }
    }
    cout<<maxn<<endl;
  }

}

C++: #include <iostream>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <cmath>

using namespace std;

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

struct Vector2{
	double x;
	double y;
};

Vector2 P[500];

int main(){
	int n;
	while(true){
		cin >> n;
		if(n == 0){
			break;
		}
		for(int i = 0; i < n; i++){
			cin >> P[i].x >> P[i].y;
		}
		double k, L_2, a, b;
		
		int max = 0;
		
		for(int num1 = 0; num1 < n; num1++){
			for(int num2 = num1 + 1; num2 < n; num2++){
				L_2 = pow2(P[num1].x - P[num2].x) + pow2(P[num1].y - P[num2].y);
				if(L_2 >= 4){
					continue;
				}
				k = sqrt((4 - L_2) / L_2);
				
				a = (k * (P[num1].y - P[num2].y) + P[num1].x + P[num2].x) / 2;
				b = (k * (P[num2].x - P[num1].x) + P[num1].y + P[num2].y) / 2;
				
				int count;
				
				count = 0;
				
				for(int i = 0; i < n; i++){
					if(i == num1 || i == num2){
						count++;
						continue;
					}
					if(pow2(P[i].x - a) + pow2(P[i].y - b) <= 1){
						count++;
					}
				}
				
				if(max < count){
					max = count;
				}
				
				a = (k * (P[num2].y - P[num1].y) + P[num1].x + P[num2].x) / 2;
				b = (k * (P[num1].x - P[num2].x) + P[num1].y + P[num2].y) / 2;
				
				count = 0;
				
				for(int i = 0; i < n; i++){
					if(i == num1 || i == num2){
						count++;
						continue;
					}
					if(pow2(P[i].x - a) + pow2(P[i].y - b) <= 1){
						count++;
					}
				}
				
				if(max < count){
					max = count;
				}
			}
		}
		printf("%d\n", max);
	}
	return 0;
}

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

using namespace std;
typedef complex<double> P;

int N;
P dots[300];
#define EPS (1e-6)

// 2_ðÊé¼arÌ~Ì_ðßé
pair<P,P> calcCircleCenterPoint(P dot1,P dot2,double r){
    double v=abs(dot1-dot2);
    double x = sqrt(r*r - (v/2)*(v/2));
    // ñgÌPÊ@üxNg
    P hose[2];
    //double y1 = -(1.0)/(dot1.imag() - dot2.imag());
    //double x1 = 1.0/(dot1.real() - dot2.real());
    double x1=-(dot1.imag()-dot2.imag());
    double y1=(dot1.real() - dot2.real());
    hose[0]=P(x1,y1);
    //double y2 = 1.0/(dot1.imag() - dot2.imag());
    //double x2 = -1.0/(dot1.real() - dot2.real());
    double y2=-y1;
    double x2=-x1;
    hose[1]=P(x2,y2);
    double dd1=abs(hose[0]);
    double dd2=abs(hose[1]);
    hose[0]/=dd1;hose[0]*=x;
    hose[1]/=dd2;hose[1]*=x;
    P tmp=(dot1+dot2);
    tmp/=2;
    hose[0]+=tmp;
    hose[1]+=tmp;
    return make_pair(hose[0],hose[1]);
}

int main(){

	while(cin >> N && N != 0){
        for(int i = 0; i < N; i++){
            P p;
            double x,y;
            cin>>x>>y;
            dots[i]=P(x,y);
            //cin >> dots[i].second >> dots[i].first;
		}
		int maxCnt = 1;
		// ~ðìé½ßÌQ_ðßé
		for(int i = 0; i < N; i++){
			for(int j = i+1; j < N; j++){
                //pair<double,double> dot1,dot2;
                P dot1,dot2;
				dot1 = dots[i];
				dot2 = dots[j];
                //double v = sqrt((dot1.first - dot2.first)*(dot1.first - dot2.first)
                //+(dot1.second - dot2.second)*(dot1.second - dot2.second));
                double v=abs(dot1-dot2);
				// ñ_Ô£ª2æè¨¨«¯êÎÂ­êÈ¢
				if(v > 2.0 + EPS){
					continue;
				}

    //            double x = sqrt(1 - (v/2)*(v/2));
				//// ñgÌPÊ@üxNg
                P hose[2];
                pair<P,P> tmp=calcCircleCenterPoint(dot1,dot2,1);
                hose[0]=tmp.first;
                hose[1]=tmp.second;
    //            double y1 = -1.0/(dots[i].imag() - dots[j].imag());
    //            double x1 = 1.0/(dots[i].real() - dots[j].real());
    //            hose[0]=P(x1,y1);
    //            double y2 = 1.0/(dots[i].imag() - dots[j].imag());
    //            double x2 = -1.0/(dots[i].real() - dots[j].real());
    //            hose[1]=P(x2,y2);

				////double dd1 = sqrt((hose[0].second)*(hose[0].second)+(hose[0].first)*(hose[0].first));
    //            //double dd2 = sqrt((hose[1].second)*(hose[1].second)+(hose[1].first)*(hose[1].first));
    //            double dd1=abs(hose[0]);
    //            double dd2=abs(hose[1]);
    //            hose[0]/=dd1;hose[0]*=x;
    //            hose[1]/=dd2;hose[1]*=x;
    //            P tmp=(dot1+dot2);
    //            tmp/=2;
    //            hose[0]+=tmp;
    //            hose[1]+=tmp;

    //            hose[0].first /= dd1; hose[0].second /= dd1;
				//hose[1].first /= dd2; hose[1].second /= dd2;
				// x{·é
				//hose[0].first *= x;hose[1].first *= x;
				//hose[0].second *= x;hose[1].second *= x;

				// mð«·
				//hose[0].first += (dot1.first + dot2.first)/2;
				//hose[0].second += (dot1.second + dot2.second)/2;
				//hose[1].first += (dot1.first + dot2.first)/2;
				//hose[1].second += (dot1.second + dot2.second)/2;

				// hose[0]Æhose[1]Í~ÌSÀW
				// e_ÉÖµÄ~ÌÉ¶Ý·é©Ç¤©ðßé
				for(int k = 0; k < 2; k++){
                    int cnt = 0;
					for(int l = 0; l < N; l++){
						//double dist = sqrt((dots[l].first -hose[k].first)*(dots[l].first -hose[k].first)
						//	+ (dots[l].second -hose[k].second)*(dots[l].second -hose[k].second));
                        double dist=abs(dots[l]-hose[k]);
						if(dist <= 1.0 + EPS){
							cnt++;
						}
					}
					maxCnt = max(maxCnt,cnt);
				}
			}
		}
		cout << maxCnt << endl;
	}

	return 0;
}