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>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
using namespace std;
using ll = long long ;
using P = pair<int,int> ;
using pll = pair<long long,long long>;
constexpr int INF = 1e9;
constexpr long long LINF = 1e17;
constexpr int MOD = 1000000007;
constexpr double PI = 3.14159265358979323846;

using d = double;
constexpr double EPS = 1e-6;
constexpr bool eq(const d& x,const d& y){return abs(x-y)<EPS;}
struct point{
    d x=0,y=0;
    constexpr point(const d& x=0,const d& y=0):x(x),y(y){}
    constexpr bool operator==(const point& r)const{return eq(x,r.x)&&eq(y,r.y);};
    constexpr bool operator!=(const point& r)const{return !(*this==r);};
    constexpr bool operator<(const point& r)const{return (eq(x,r.x))?(y<r.y):(x<r.x);};
    constexpr bool operator>(const point& r)const{return (eq(x,r.x))?(y>r.y):(x>r.x);};
    constexpr point operator+(const point& r)const{return point(*this)+=r;}
    constexpr point operator+=(const point& r){
        x+=r.x;y+=r.y;
        return *this;
    }
    constexpr point operator-(const point& r)const{return point(*this)-=r;}
    constexpr point operator-=(const point& r){
        x-=r.x;y-=r.y;
        return *this;
    }
    constexpr point operator*(const d& r)const{return point(*this)*=r;}
    constexpr point operator*=(const d& r){
        x*=r;y*=r;
        return *this;
    }
    constexpr point operator*(const point& r)const{return point(*this)*=r;}
    constexpr point operator*=(const point& r){
        d tmp=x;
        x=x*r.x-y*r.y;
        y=tmp*r.y+y*r.x;
        return *this;
    }
    constexpr point operator/(const double& r)const{return point(*this)/=r;}
    constexpr point operator/=(const double& r){
        x/=r;y/=r;
        return *this;
    }
};
d norm(const point& a){return a.x*a.x+a.y*a.y;}
double abs(const point& a){return sqrt(a.x*a.x+a.y*a.y);}
d dot(const point& a,const point& b){return a.x*b.x+a.y*b.y;}
d cross(const point& a,const point& b){return a.x*b.y-a.y*b.x;}
point projection(const point& p,const point& b){return b*dot(p,b)/norm(b);}
point unit_vec(const point& p){return p/abs(p);}
point normal_vec(const point& p){return p*point(0,1);}
point rotation(const point& center,const point& p,const point& rot){return rot*(p-center)+center;}
double area(const point& a,const point& b,const point& c){return cross(b-a,c-a)/2;}

struct Line{
    point a,b;
    constexpr Line(const point& a=point(0,0),const point& b=point(0,0)):a(a),b(b){}
    constexpr Line(const double& A,const double& B,const double& C){ // Ax + By + C = 0
        if(eq(A,0)&&eq(B,0)) assert(-1);
        else if(eq(B,0)){
            a = point(-C/A,0);
            b = point(-C/A,1);
        }else{
            a = point(0,-C/B);
            b = point(1,-(A+C)/B);
        }
    }
};
bool is_orthoonal(const Line& l,const Line& r){return eq(dot(l.b-l.a,r.b-r.a),0);}
bool is_pararell(const Line& l,const Line& r){return eq(cross(l.b-l.a,r.b-r.a),0);}
bool on_Line(const Line& l,const point& p){return eq(dot(p-l.a,l.b-l.a),0);}
double dis_Lp(const Line& l,const point& p){return abs(cross(l.b-l.a,p-l.a))/abs(l.b-l.a);}
point Line_intersect(const Line& l,const Line& r){
    point lv=l.b-l.a,rv=r.b-r.a;
    return l.a + lv * abs(cross(r.b-l.a,rv) / cross(lv,rv)); 
}
Line vertical_bisector(const point& p,const point& q){return Line((p+q)/2,rotation((p+q)/2,q,point(0,1)));}
using Segment = Line;
bool on_Segment(const Segment& s,const point& p){
    return ( abs(s.a-p)+abs(s.b-p)<abs(s.a-s.b)+EPS );
}
bool is_segment_intersect(const Segment& s,const Segment& t){
    return  (cross(s.b-s.a,t.a-s.a) * cross(s.b-s.a,t.b-s.a) < EPS) &&
            (cross(t.b-t.a,s.a-t.a) * cross(t.b-t.a,s.b-t.a) < EPS);
}
point segment_intersect(const Segment& s,const Segment& t){
    double d1 = dis_Lp(s,t.a);
    double d2 = dis_Lp(s,t.b); 
    double p = d1/(d1+d2);
    return t.a + (t.b-t.a)*p;
}

bool on_circle(point cen,point p){
    return ( (cen.x-p.x)*(cen.x-p.x)+(cen.y-p.y)*(cen.y-p.y)<=1.0+EPS );
}

int main(){
    int n;
    while(cin >> n){
        if(n==0) break;
        vector<point> points(n);
        rep(i,n) cin >> points[i].x >> points[i].y;
        vector<point> centers;
        vector<Line> lines;
        rep(i,n)rep(j,i){
            point p = (points[i]+points[j])/2;
            centers.push_back(p);
            if(1-norm(p-points[i])<=0) continue;
            double dd = sqrt(1-norm(p-points[i])) / abs(p-points[i]);
            point a = rotation(p,points[i],point(0,1)*dd);
            point b = rotation(p,points[i],point(0,-1)*dd);
            centers.push_back(a);
            centers.push_back(b);
        }

        int mx = 1;
        for(point c:centers){
            int res = 0;
            rep(i,n){
                if(on_circle(c,points[i])) ++res;
            }
            mx = max(res,mx);
            //cout << n << " " << res << endl;
        }
        cout << mx << endl;
    }
    return 0;
}

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

#define loop(i,a,b) for(int i=a; i<int(b); i++)
#define rep(i,n) loop(i,0,n)

using namespace std;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef complex<double> vec;

double const pi = acos(-1);
double const eps = 1e-6;

int main(){
    int n;
    while (cin >> n && n){
        vector<vec> v(n);
        rep(i, n){
            double x, y; cin >> x >> y;
            vec p(x, y);
            v[i] = p;
        }

        if (n == 1) {
            puts("1");
            continue;
        }
        int ans = 1;
        rep(i, n)rep(j, i){
            double d = abs(v[i] - v[j]);
            if (d > 2) continue;
            ans = max(ans, 2);
            double h = sqrt(1 - d*d / 4);
            vec c1 = polar(h, arg(v[i] - v[j]) + pi / 2) + (v[i] + v[j]) / 2.;
            vec c2 = polar(h, arg(v[i] - v[j]) - pi / 2) + (v[i] + v[j]) / 2.;
            int a1 = 0, a2 = 0;
            rep(k, n){
                if (abs(c1 - v[k]) <= 1) a1++;
                if (abs(c2 - v[k]) <= 1) a2++;
            }
            ans = max(ans, a1);
            ans = max(ans, a2);
        }
        cout << ans << endl;
    }
}