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 <string>
#include <sstream>
#include <algorithm>
#include <vector>
#include <utility>
#include <functional>
#include <stack>
#include <queue>
#include <map>
#include <set>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <deque>
#include <ctime>
  
using namespace std;
  
#define rep(i,n) REP(i,0,n)
#define REP(i,s,e) for(int i=(s); i<(int)(e); i++)
#define pb push_back
#define mp make_pair
#define all(r) r.begin(),r.end()
#define rall(r) r.rbegin(),r.rend()
#define fi first
#define se second
#define println(X) cout<<X<<endl;
#define DBG(X) cout<<#X<<" : "<<X<<endl;
  
typedef long long ll;
typedef vector<int> vi;
typedef vector<vi> vii;
typedef vector<ll> vl;
typedef vector<vl> vll;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
 
const int INF = 1e9;
const ll LINF = 1e18;
const ll MOD = 1e9 + 7;
double EPS = 1e-10;
const double PI = acos(-1);

struct Vector2{
	double x, y;
	Vector2(){};
	Vector2(double x, double y) : x(x), y(y){};

	double cross(const Vector2& p){
		return x*p.y - y*p.x;
	}

	double dot(const Vector2& p) const{
		return x*p.x + y*p.y;
	}

	const Vector2 medPoint(const Vector2& v){
		return {(x + v.x) / 2.0, (y + v.y) / 2.0};
	}

	const bool operator == (const Vector2& v) const{
		return abs(x - v.x) < EPS && abs(y - v.y) < EPS;
	}
	const bool operator != (const Vector2& v) const{
		return !(*this==(v));
	}

	const Vector2 operator - (const Vector2& p) const{
		return {x-p.x, y-p.y};
	}

	const Vector2& operator += (const Vector2& p) {
		x += p.x;
		y += p.y;
		return *this;
	}

	double dist2(const Vector2& p) const{
		return ((*this) - p).dot((*this) - p);
		//return (x-p.x)*(x-p.x) + (y-p.y)*(y-p.y);
	}

	Vector2 spin(double arg){
		double si = sin(arg);
		double co = cos(arg);
		return {x*co - y*si, x*si + y*co};
	}
};

struct Circle{
	Vector2 p;
	double r;
	Circle(){};
	Circle(const Vector2 &p, double r) : p(p), r(r){};
	Circle(double x, double y, double r) : p({x, y}), r(r){};

	bool isPointIn(const Vector2& v) const{
		return p.dist2(v) < r * r + EPS;
	}
};

vector<Circle> getCircles(const Vector2 &p1, const Vector2 &p2){
	if(p1.dist2(p2) > 4.0 + EPS) return {};
	vector<Circle> ret;
	Vector2 p = p1 - p2;
	// double arg = atan2(p.y, p.x);
	// p = p.spin(-arg);

	// ret.pb({1.0, p.y/2.0, 1.0});
	// ret.pb({-1.0, p.y/2.0, 1.0});

	double arg = atan2(p.y, p.x);
	p = p.spin(-arg);

	double x = p.x/2.0, y = sqrt(1.0-x * x);
	ret.pb({{x, y}, 1.0});
	ret.pb({{x, -y}, 1.0});

	rep(i, 2){
		ret[i].p = ret[i].p.spin(arg);
		ret[i].p += p2;
	}
	return ret;
}

int main(){
	int n;
	while(cin>>n && n){
		vector<Vector2> v;
		rep(i, n){
			double x, y;
			cin>>x>>y;
			v.pb({x, y});
		}
		int ans = 1;
		for(int i = 1; i < n; i++){
			for(int j = 0; j < i; j++){
				vector<Circle> c = getCircles(v[j], v[i]);
				for(int k = 0; k < c.size(); k++){
					int tmp = 0;
					for(int l = 0; l < n; l++){
						if(c[k].isPointIn(v[l])) tmp++;
					}
					ans = max(ans, tmp);
				}
			}
		}
		cout<<ans<<endl;
	}
}

C++: #include <iostream> 
#include <string> 
#include <vector> 
#include <cmath> 
#include <algorithm> 
#include <cstdlib> 
#include <ctime> 
#include <cstdio> 
#include <functional> 
#include <set> 
#include <sstream> 

using namespace std; 

double d2(double x1,double x2,double y1,double y2){
	return (x1-x2)*(x1-x2)+(y1-y2)*(y1-y2);
}

int main(){

	int n;
	const double eps=1e-8;
	while(true){

		cin>>n;

		if(n==0) return 0;

		int res=1;
		vector<double> x(n),y(n);

		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++){
				
				int res2=0;
				const double X=(x[i]+x[j])/2.0;
				const double Y=(y[i]+y[j])/2.0;

				for(int k=0;k<n;k++){
					printf("P1(%f, %f) P2(%f, %f) -> d2=%f\n",x[k],y[k],X,Y,sqrt(d2(x[k],X,y[k],Y)));
					if(d2(x[k],X,y[k],Y)<1.0+eps) res2++;
				}
				cout<<endl;

				res=max(res,res2);

			}
		}

		cout<<res<<endl;

	}

	return 0;
}

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

#define EPS (1e-10)

using namespace std;

int count( vector<complex<double>> &vpos );

int main(int argc, char const *argv[]) {
    int n, result;
    double x, y;
    vector<complex<double>> vpos;
    vector<int> results;

    while(1){
        cin>>n;
        if( n == 0 ) break;
        for(int i = 0; i < n; i++){
            cin>>x>>y;
            vpos.push_back( complex<double> (x,y) );
        }
        results.push_back( count( vpos ) );
        vpos.clear();
    }

    for (auto it = results.begin(); it != results.end(); it++) {
        cout<<*it<<endl;
    }

    return 0;
}

int count( vector<complex<double>> &vpos ){
    int count = 0;
    int tmp_count1, tmp_count2;
    int size = vpos.size();
    bool flag = false;
    complex<double> m, n1, n2, center1, center2;

    if( size == 1 ) return 1;

    for(int i = 0; i < size; i++){
        for(int j = i + 1; j < size ; j++){
            if( vpos[i].real() == vpos[j].real() ) continue;
            tmp_count1 = 0;
            tmp_count2 = 0;

            m = ( vpos[i] + vpos[j] ) / 2.0;
            n1 = ( vpos[i] -  vpos[j] ) * complex<double> ( 0.0, 1.0 );
            n1 = n1 / abs(n1);
            n2 = ( vpos[i] -  vpos[j] ) * complex<double> ( 0.0, -1.0 );
            n2 = n2 / abs(n2);

            double x = abs( vpos[i] -  vpos[j] ) / 2.0;
            center1 =  m + n1 * sqrt( 1.0 - x * x );
            center2 =  m + n2 * sqrt( 1.0 - x * x );

            if( center1.real() > 0.0 ) flag = true;
            if( center2.real() > 0.0 ) flag = true;

            for(int k = 0; k < size; k++){
                double dist1 = abs( center1 - vpos[k] );
                double dist2 = abs( center2 - vpos[k] );
                if( dist1 <= 1.0 + EPS ){
                    tmp_count1++;
                }
                if( dist2 <= 1.0 + EPS ){
                    tmp_count2++;
                }
            }
            if( count < max(tmp_count1, tmp_count2) ){
                count = max(tmp_count1, tmp_count2);
            }
        }
    }

    if(!flag) return 1;
    return count;
}