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 <cmath>
#include <algorithm>
#include <queue>
#include <map>
#define LP(i,n) for(int i = 0; i < n; i++)

using namespace std;
typedef pair<double,double> P;

double Dist(double x1[], double x2[]){
	return pow((x1[0]-x2[0])*(x1[0]-x2[0])+(x1[1]-x2[1])*(x1[1]-x2[1]),0.5);

}

int main(){
	int n;
	cin >> n;
	while(n){
		int mx = 0;
		P p[300];
		LP(i,n){
			cin >> p[i].first >> p[i].second;
		}
		sort(p,p+n);
		double x[300][2];
		LP(i,n){
			x[i][0] = p[i].first;
			x[i][1] = p[i].second;
		}
		double dist[300][300];
		LP(i,n){
			LP(j,n){
				dist[i][j] = Dist(x[i],x[j]);
			}
		}
		LP(i,n){
			for(int j = i+1; j < n; j++){
				if(x[j][0]-x[i][0] > 2){
					break;
				}else{
					if(dist[i][j] < 2){
						double m[] = {(x[j][0]+x[i][0])/2,(x[j][1]+x[i][1])/2};
						double d[] = {x[j][1]-x[i][1],-x[j][0]+x[i][0]};
						double O[] = {0,0};
						double _d = Dist(d,O);
						LP(k,2){
							d[k] *= (pow(_d,-2)-0.25);
						}
						double o[2][2];
						LP(k,2){
							LP(l,2){
								o[k][l] = m[l];
								if(k) o[k][l] += d[l];
								else o[k][l] -= d[l];
							}
						}
						LP(k,2){
							int s = 2;
							for(int l = i+1; l < j; l++){
								if(Dist(x[l],o[k]) < 1) s++;
							}
							for(int l = i-1; 0<=l && x[i][0]-x[l][0] < 2;l--){
								if(Dist(x[l],o[k]) < 1) s++;
							}
							for(int l = j+1; l < n && x[l][0]-x[j][0] < 2; l++){
								if(Dist(x[l],o[k]) < 1) s++;
							}
							mx = max(mx,s);
						}
					}
				}
			}
		}
		if(mx == 0) mx = 1;
		cout << mx << endl;
		cin >> n;
	}
}

C++: #include <iostream>
#include <cmath>

using namespace std;

void calc(double x1,double y1,double x2,double y2,double *c_x1,double *c_y1,double *c_x2,double *c_y2){
  double midx,midy,hox,hoy,dx,dy,dis;
  midx = (x1+x2) / 2.0 , midy = (y1+y2) / 2.0;
  dx = x1-x2 , dy = y1-y2;
  dis = sqrt(1 - pow(x1-midx,2.0) - pow(y1-midy,2.0));
  hox = sqrt(pow(dy,2.0) / (pow(dy,2.0) + pow(dx,2.0)));
  hoy = hox * (-1.0) * (dx / dy);
  *c_x1 = midx + hox*dis,*c_y1 = midy + hoy*dis;
  *c_x2 = midx - hox*dis,*c_y2 = midy - hoy*dis;
}

int main(){
  int N;

  while(cin >> N){
    double x[300],y[300],c_x1,c_y1,c_x2,c_y2;
    int count,m;
    m = 0;
    if(!(N)) break;
    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) continue;
        if((x[i] - x[j])*(x[i] - x[j]) + (y[i] - y[j])*(y[i] - y[j]) >= 4.0 ) continue;
        count = 2;
        calc(x[i],y[i],x[j],y[j],&c_x1,&c_y1,&c_x2,&c_y2);
        for(int k=0;k<N;k++){
          if(k == i || k == j) continue;
          if((x[k] - c_x1)*(x[k] - c_x1) + (y[k] - c_y1)*(y[k] - c_y1) <= 1.0  ) count++;
        }
        if(m < count) m = count;
        count = 2;
        for(int k=0;k<N;k++){
          if(k == i || k == j) continue;
          if((x[k] - c_x2)*(x[k] - c_x2) + (y[k] - c_y2)*(y[k] - c_y2) <= 1.0  ) count++;
        }
        if(m < count) m = count;
      }
    }

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

  return 0;
}

C++: #include<bits/stdc++.h>
using namespace std;

double distance(pair<double,double> p1,pair<double,double> p2){
	return sqrt((p1.first-p2.first)*(p1.first-p2.first)+(p1.second-p2.second)*(p1.second-p2.second));
}

int main(){
	int n;
	while(cin>>n,n){
		vector<double> x(n);
		vector<double> y(n);
		int res=1;
		for(int i=0;i<n;i++) cin>>x[i]>>y[i];
		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				pair<double,double> midp=make_pair((x[i]+x[j])/2,(y[i]+y[j])/2);
				double a=-(x[i]-x[j]),b=(y[i]-y[j]);
				if(a<0.000001) a=0;
				if(b<0.000001) b=0;
				//a(x-mid[i].first)+b(y-mid.second)=0
					for(int k=0;k<2;k++){
						pair<double,double> lef=midp;
						pair<double,double> rig;
						if(k==0) rig=make_pair(-b*1e8+midp.first,a*1e8+midp.second);
						if(k==1) rig=make_pair(b*1e8+midp.first,-a*1e8+midp.second);
						while(distance(lef,rig)>=0.00000001){
								pair<double,double> mid=make_pair((lef.first+rig.first)/2,(lef.second+rig.second)/2);
								if(distance(mid,make_pair(x[i],y[i]))>1){
									rig=mid;
								}
								else{
									lef=mid;
								}

						}
						int count=0;
						for(int r=0;r<n;r++){
							if(distance(make_pair(x[r],y[r]),lef)<1+0.0001){
								count++;
							}
						}
						res=max(res,count);
					}
			}

		}
		cout<<res<<endl;

	}

	return 0;
}

C++: #include <iostream>
#include <iomanip>
#include <complex>
#include <vector>
#include <algorithm>
#include <cmath>
#define EPS 1e-8
using namespace std;
typedef complex<double> P;
struct C{
    P p;
    double r;
    C(const P& p, const double& r):
        p(p), r(r) {}
};

P center(const P& a, const P& b, int lr){
    P mid = (a+b)/P(2,0);
    P len(sqrt(1-abs(b-mid)*abs(b-mid)),0);
    P rot = P(imag(b-a),-real(b-a)) / P(abs(a-b),0);
    if(lr==0){
        return len*rot+mid;
    }else{
        return len*rot*P(-1,0)+mid;
    }
}

int main(){
    while(1){
        int n;
        cin >> n;
        if(n==0) break;

        vector<P> pset(n);
        for(int i=0; i<n; i++){
            double x,y;
            cin >> x >> y;
            pset[i] = P(x,y);
        }

        int mnop=0;
        for(int i=0; i<n-1; i++){
            for(int j=i+1; j<n; j++){
                if(abs(pset[i]-pset[j]) < 2+EPS){
                    for(int lr=0; lr<2; lr++){
                        int count=0;
                        P mid = center(pset[i], pset[j], lr);
                        for(int k=0; k<n; k++){
                            if(abs(pset[k]-mid) < 1+EPS){
                                count++;
                            }
                        }
                        mnop = max(mnop, count);
                    }
                }
            }
        }

        if(mnop==0) mnop=1;
        cout << mnop << endl;
    }
            
    return 0;
}