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 <set>
#include <vector>
#include <algorithm>
#include <cmath>

#define EPS (1e-10)
#define EQ(a,b) (abs((a)-(b)) < EPS)
#define EQV(a,b) ( EQ((a).real(), (b).real()) && EQ((a).imag(), (b).imag()) )

using namespace std;

typedef complex<double> P;

P center1(P a ,P b ,double r){
  return (a+b)/(P)2+(b-a)*P(0,1)*(P)(sqrt(r*r-norm((b-a)/(P)2))/abs(b-a));
}

P center2(P a ,P b ,double r){
  return (a+b)/(P)2+(b-a)*P(0,-1)*(P)(sqrt(r*r-norm((b-a)/(P)2))/abs(b-a));
}

int main(){
  int n;
  while(true){
    cin >>n;
    if (!n)
      break;
    P dot[n];
    double x,y;
    for(int i=0;i<n;i++){
      cin >>x >>y;
      dot[i] = P(x,y);
    }
    set<int> near_dot[n];
    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
	if (norm(dot[i]-dot[j])<4)
	  near_dot[i].insert(j);
      }
    }
    
    set<int>::iterator it;
    int cnt,max=2;
    P c;
    for(int i=0;i<n;i++){
      it = near_dot[i].begin();
      while(it != near_dot[i].end()){
	c = center1(dot[i],dot[*it],1);
	cnt=2;
	for(int j=0;j<n;j++)
	  if (norm(dot[j]-c)<1 && j!=i && j!=*it)
	    cnt++;
	if (max<cnt){
	  max=cnt;
	}

	c = center2(dot[i],dot[*it],1);
	cnt=2;
	for(int j=0;j<n;j++)
	  if (norm(dot[j]-c)<1 && j!=i && j!=*it)
	    cnt++;
	if (max<cnt){
	  max=cnt;
	}
	it++;
      }
    }
    cout << max << endl;
  }
  return 0;
}

C++: #include <iostream>
#include <cmath>
#include <algorithm>
using namespace std;

struct point{
    double x;
    double y;
};

int n;
point points[300];

double getDist2(point a, point b){
    return pow(a.x - b.x ,2) + pow(a.y - b.y, 2);
}

double getDist(point a, point b){
    return sqrt(getDist2(a, b));
}

int findMaxCover(){
    int ans = 1;
    int count;
    int sign[2] = {-1, 1};
    for (int i = 0; i < n - 1; i++) {
        for (int j = i + 1; j < n; j++) {
            point a, b;
            a = points[i];
            b = points[j];
            if (getDist2(a, b) > 4.0) {
                continue;
            }
            point m;
            m.x = (a.x + b.x) / 2.0;
            m.y = (a.y + b.y) / 2.0;
            double v, u, d, e;
            d = getDist(a, b);
            v = (a.y - b.y) / d;
            u = (a.x - b.x) / d;
            d = getDist(a, b);
            e = sqrt(1 - d * d / 4);
            for (int k = 0; k < 2; k++) {
                count = 0;
                double p = m.x - e * v * sign[k];
                double q = m.y + e * u * sign[k];
                for (int l = 0; l < n; l++) {
                    point t;
                    t.x = p;
                    t.y = q;
                    if (l == i || l == j) {
                        count++;
                        continue;
                    }
                    if (getDist2(t, points[l]) <= 1.0) {
                        count++;
                    }
                }
                ans = max(ans, count);
            }
        }
    }
    return ans;
}

int main()
{
    while (true) {
        cin >> n;
        if (n == 0) {
            break;
        }
        for (int i = 0; i < n; i++) {
            cin >> points[i].x >> points[i].y;
        }
        cout << findMaxCover() << endl;
        
    }

    return 0;
}

C++: #include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <queue>
#include <stack>
#include <set>

using namespace std;


int main(){
	int n;
	double x1,y1,xc,yc,dis,vx,vy,ex,ey;
	while(1){
		cin >> n;
		if(n==0)break;
		vector<double> x(n);
		vector<double> y(n);
		for(int i=0;i<n;i++){
			cin >> x[i] >> y[i];
		}
		int ans =0;
		for(int i=0;i<n;i++){
			for(int j=0;j<n;j++){
				if(i==j)break;
				xc = (x[i]+x[j])/2;
				yc = (y[i]+y[j])/2;
				dis = sqrt(1.0-(xc-x[i])*(xc-x[i])-(yc-y[i])*(yc-y[i]));
				vx = x[i]-x[j];
				vy = y[i]-y[j];
				ex = dis * vy * sqrt(vx*vx+vy*vy);
				ey = -dis* vx * sqrt(vx*vx+vy*vy);
				x1 = xc + ex;
				y1 = yc + ey;
				int count=0;
				for(int k=0;k<n;k++){
					if(k==i||k==j)count +=1;
					else{
						if((x1-x[k])*(x1-x[k])+(y1-y[k])*(y1-y[k])<=1){
							count +=1;
						}
					}
				}
				ans = max(ans,count);
				x1 = xc - ex;
				y1 = yc - ey;
				count = 0;
				for(int k=0;k<n;k++){
					if(k==i||k==j)count +=1;
					else{
						if((x1-x[k])*(x1-x[k])+(y1-y[k])*(y1-y[k])<=1){
							count +=1;
						}
					}
				}
				ans = max(ans,count);
			}
		}
		cout << ans << endl;
	}
	return 0;
}

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

#define int long long
typedef vector<int>vint;
typedef pair<int,int>pint;
typedef vector<pint>vpint;
#define rep(i,n) for(int i=0;i<(n);i++)
#define reps(i,f,n) for(int i=(f);i<(n);i++)
#define all(v) (v).begin(),(v).end()
#define each(it,v) for(__typeof((v).begin()) it=(v).begin();it!=(v).end();it++)
#define pb push_back
#define fi first
#define se second
template<typename A,typename B>inline void chmin(A &a,B b){if(a>b)a=b;}
template<typename A,typename B>inline void chmax(A &a,B b){if(a<b)a=b;}

const double PI=acos(-1);
const double EPS=1e-6;

int N;
double x[333],y[333];

signed main(){
    while(cin>>N,N){
        rep(i,N)cin>>x[i]>>y[i];

        int ma=1;
        rep(i,N)rep(j,N){
            if(i==j)continue;
            double dx=x[i]-x[j],dy=y[i]-y[j];
            double d=sqrt(dx*dx+dy*dy);
            if(d>2+EPS)continue;
            double r=sqrt(1.0-d*d/4);
            double vx=-dy/d*r;
            double vy=dx/d*r;

            double xx=(x[i]+x[j])/2+vx;
            double yy=(y[i]+y[j])/2+vy;
            int cnt=0;
            rep(k,N){
                if((x[k]-xx)*(x[k]-xx)+(y[k]-yy)*(y[k]-yy)<1.0+EPS)cnt++;
            }
            chmax(ma,cnt);
        }
        printf("%lld\n",ma);
    }
}