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;

using ld = long double;
using Point = complex<ld>;

struct Line{
    Point a, b;
    Line (Point a, Point b) : a(a), b(b) {}
    Line () : Line(Point(), Point()) {}
};

struct Circle{
    Point p;
    ld r;
    Circle (Point p, ld r) : p(p), r(r) {}
    Circle () : Circle(Point(),0.0) {}
};

constexpr 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();
  }
}

namespace Geometry{
    bool eq(ld a, ld b){
        return abs(a-b) < eps;
    }
    //内積
    ld dot(Point a, Point b){
        return real(conj(a) * b);
    }
    //外積
    ld cross(Point a, Point b){
        return imag(conj(a) * b);
    }
    //3点の位置関係
    int ccw (Point a, Point b, Point c) {
        b -= a; c -= a;
        if (cross(b, c) > eps) return 1;//a,b,cで反時計周り
        if (cross(b, c) < -eps) return -1;//a,b,cで時計周り
        if (dot(b, c) < 0) return 2;//c,a,bで直線
        if (norm(b) < norm(c)) return -2;//a,b,cで直線
        return 0;//a,c,bで直線
    }

//====================================================
    Point inputPoint(){
        ld x, y;
        cin>>x>>y;
        return Point(x, y);
    }
    //2直線の交差判定
    bool isCrossed_ll(Line l, Line m){
        return !eq(cross(l.b-l.a, m.b-m.a), 0);
    }
    //直線と線分の交差判定
    bool isCrossed_ls(Line l, Line s){
        return isCrossed_ll(l, s) && 
            cross(l.b - l.a, s.a - l.a) * cross(l.b - l.a, s.b - l.a) < eps;
    }
    //線分と線分の交差判定
    bool isCrossed_ss(Line s, Line t){
        return ccw(s.a, s.b, t.a) * ccw(s.a, s.b, t.b) <= 0 &&
            ccw(t.a, t.b, s.a) * ccw(t.a, t.b, s.b) <= 0;
    }
    //点が直線上にあるか
    bool isON_l(Line l, Point p){
        return abs(cross(l.b-p, l.a-p)) < eps;
    }
    //点が線分上にあるか
    bool isON_s(Line s, Point p){
        return abs(s.a-p)+abs(s.b-p)-abs(s.b-s.a) < eps;
    }

    //点から直線への垂線の足
    Point foot(Line l, Point p){
        ld t = dot(p-l.a,l.a-l.b)/norm(l.a-l.b);
        return l.a+t*(l.a-l.b);
    }
    //直線と直線の交点
    Point intersection_ll(Line l, Line m){
        Point lv=l.b-l.a, mv=m.b-m.a;
        assert(cross(lv,mv)!=0);//平行
        return l.a+lv*cross(mv,m.a-l.a)/cross(mv,lv);
    }
    //線分と線分の交点
    Point intersection_ss(Line s, Line t){
        assert(isCrossed_ll(s, t));
        return intersection_ll(s,t);
    }
    //点と点の距離
    ld dist_pp(Point p, Point q){
        ld x=p.real()-q.real(), y=p.imag()-q.imag();
        return sqrt(x*x+y*y);
    }
    //点と直線の距離
    ld dist_lp(Line l, Point p) {
        return abs(p - foot(l, p));
    }
    //直線と直線の距離
    ld dist_ll(Line l, Line m){
        return isCrossed_ll(l, m)?0:dist_lp(l, m.a);
    }
    //直線と線分の距離
    ld dist_ls(Line l, Line s){
        return isCrossed_ls(l, s)?0:min(dist_lp(l, s.a), dist_lp(l, s.b));
    }
    //線分と点の距離
    ld dist_sp (Line s, Point p) {
        Point r = foot(s, p);
        return isON_s(s, r) ? abs(r - p) : min(abs(s.a - p), abs(s.b - p));
    }
    //線分と線分の距離
    ld dist_ss (Line s, Line t) {
        if (isCrossed_ss(s, t)) return 0;
        return min({dist_sp(s, t.a), dist_sp(s, t.b), dist_sp(t, s.a), dist_sp(t, s.b)});
    }
    //円と円の交点
    vector<Point>intersection_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;
        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(){
    int n;
    using namespace Geometry;
    while(cin>>n,n){
        vector<Point>g(n);
        vector<Circle>Cs(n);
        for(int i=0;i<n;++i){
            g[i]=inputPoint();
            Cs[i]={g[i],1.0};
        }
        int ans=1;
        for(int i=0;i<n;++i){
            for(int j=i+1;j<n;++j){
                vector<Point>intersections(intersection_cc(Cs[i],Cs[j]));
                for(auto intersection:intersections){
                    int ret=0;
                    for(int k=0;k<n;++k){
                        if(dist_pp(intersection,g[k]) < 1.0){
                            if(i == k || j==k)continue;
                            ++ret;
                        }
                    }
                    ans=max(ans,ret+2);
                }
            }
        }
        cout<<ans<<endl;
    }
}

C++: #include<iostream>
#include<complex>
#include<cmath>
#include<algorithm>
#define REP(i,s,n) for(int i=s;i<n;i++)
#define rep(i,n) REP(i,0,n)
#define EPS 1e-10
#define PI 3.141592653589793238
// sqrtは重い
using namespace std;
typedef pair<double,double> P;

void display(complex<double> a,complex<double> b){
 cout << "(" << a.real() <<"," << a.imag() <<")"<< " " <<"("<< b.real() <<"," <<b.imag() << ")" << endl;
}

int main(){
  int n,ans;
  double xin,yin,X;
  complex<double> m;
  while(cin >> n && n){
    ans = 1;
    complex<double> point[n];
    for(int i=0;i<n;i++){
      cin >> xin >> yin;
      point[i] = complex<double>(xin,yin);
    }
    for(int i=0;i<n;i++){
      for(int j=0;j<n;j++){
	if(i == j)continue;
	complex<double> Rr = point[i]-point[j];
	if(abs(Rr) <= 2+EPS){

	  complex<double> UKPC,U,P1,P2;
	  UKPC = (point[i]+point[j])*0.5;
       
	  double T = sqrt(1-norm(point[i]-point[j])*0.25);

	  UKPC = (point[i]+point[j])*0.5;
	               /* unit vector(単位べく利ゅ)             *//* T倍に拡大 ９０度回転 */ 
  	  P1 = UKPC + (point[i]-point[j])/abs(point[i]-point[j])*polar(T,PI/2.0);
	  P2 = UKPC + (point[i]-point[j])/abs(point[i]-point[j])*polar(T,-PI/2.0);
	  //cout << "P1 = " << P1 << ", P2 = " << P2 << endl;
	  int cnt=0;
	  complex<double> CalcP;
	  double Calc;
	  for(int k=0;k<n;k++){
	    Calc = abs(P1-point[k]);
	    if(Calc < 1.0+EPS){
	      cnt++;
	    }
	  }	

	  ans = max(ans,cnt);
	  cnt = 0;
	  for(int k=0;k<n;k++){
	    Calc = abs(P2-point[k]);
	    if(fabs(Calc) < 1.0+EPS){
	      cnt++;
	    }
	  }	
  
	  ans = max(ans,cnt);
	  

	}

      }
    }

    cout << ans << endl;

  }
  return 0;
}