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"

using namespace std;

typedef long long ll;
typedef pair<int,int> P;
const double EPS = 1e-12;
const int INF = numeric_limits<int>::max()/2;
const int MOD = 1e9+7;

typedef long double ld;
typedef complex<ld> Point;

const ld eps = 1e-9, pi = acos(-1.0);

namespace std {
	bool operator<(const Point &lhs, const Point &rhs) {
		if (lhs.real() < rhs.real() - eps) return true;
		if (lhs.real() > rhs.real() + eps) return false;
		return lhs.imag() < rhs.imag();
	}
}

Point input_point() {
	ld x, y;
	cin >> x >> y;
	return Point(x, y);
}

bool eq(ld a, ld b) {
	return (abs(a - b) < eps);
}

class Circle {
public:
	Point p;
	ld r;
	Circle() : p(Point(0, 0)), r(0) {}
	Circle(Point p, ld r) : p(p), r(r) {}
};

vector<Point> is_cc(Circle c1, Circle c2) {
	vector<Point> res;
	ld d = abs(c1.p - c2.p);
	ld rc = (d * d + c1.r * c1.r - c2.r * c2.r) / (2 * d);
	ld dfr = c1.r * c1.r - rc * rc;
	if (abs(dfr) < eps) dfr = 0.0;
	else if (dfr < 0.0) return res; // no intersection
	ld rs = sqrt(dfr);
	Point diff = (c2.p - c1.p) / d;
	res.push_back(c1.p + diff * Point(rc, rs));
	if (dfr != 0.0) res.push_back(c1.p + diff * Point(rc, -rs));
	return res;
}

int main(){
    cin.tie(0);
    ios::sync_with_stdio(false);

    int n;
    while(cin>>n,n){
        vector<Point> ps(n);
        for(int i=0;i<n;i++){
            ps[i] = input_point();
        }
        if(n==1){
            cout<<n<<endl;
            continue;
        }
        int res=1;
        for(int i=0;i<n;i++){
            for(int j=i+1;j<n;j++){
                Circle c1=Circle(ps[i],1.0),c2=Circle(ps[j],1.0);
                vector<Point> vs=is_cc(c1,c2);
                int sz=vs.size();
                if(sz==0) continue;
                for(int k=0;k<sz;k++){
                    Point p=vs[k];
                    int tmp=0;
                    for(int idx=0;idx<n;idx++){
                        if(abs(ps[idx]-p)<=1+eps) tmp++;
                    }
                    res=max(res,tmp);
                }
            }
        }
        cout<<res<<endl;
    }
}

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)


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];
                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;
}

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;
  double x,y;
  while(true){
    cin >>n;
    if (!n)
      break;
    P dot[n];
    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=1;
    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=0;
	for(int j=0;j<n;j++)
	  if (norm(dot[j]-c)<=1.0001)
	    cnt++;
	if (max<cnt){
	  max=cnt;
	}

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