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 <vector>
#include <stack>

using namespace std;

template<class T> bool operator<(const complex<T> &a, const complex<T> &b){
    return a.real() == b.real() ? a.imag() < b.imag() : a.real() < b.real();
}

typedef complex<double> P;
#define X real()
#define Y imag()

double dot(P a, P b){
    return (a * conj(b)).X;
}
double cross(P a, P b){
    return (a * conj(b)).Y;
}

vector<P> points;
deque<P> selected;
int max_count;

void search(int index) {
    if(index > points.size()) return;
    int selected_size = selected.size();

    bool is_all_in_circle = true;
    if(selected_size >= 2 && selected_size > max_count) {
        // calc is all overlaps
        vector<P> copied_selected;
        copy(selected.begin(), selected.end(), back_inserter(copied_selected));

        for (int i = 0; i < selected_size ; ++i) {
            for (int j = 0; j < selected_size; ++j) {
                if(i == j) continue;
                // two point check
                double d12 = abs(copied_selected[i] - copied_selected[j]);
                if(d12 > 2.0 + 1e-10) {
                    is_all_in_circle = false;
                    break;
                }

                for (int k = 0; k < selected_size; ++k) {
                    if(i == k || j == k) continue;
                    // three point check
                    P p1 = copied_selected[i];
                    P p2 = copied_selected[j];
                    P p3 = copied_selected[k];

                    P center12 = P((p1.X + p2.X) / 2, (p1.Y + p2.Y) / 2);
                    double allow_distance = 1 + sqrt(1 - ((d12 * d12) / 4));
                    if(abs(center12 - p3)  > allow_distance + 1e-10) {
                        is_all_in_circle = false;
                        break;
                    }
                }

            }
            if(!is_all_in_circle) {
                break;
            }
        }

        copied_selected.clear();
        copied_selected.shrink_to_fit();
    }

    if (is_all_in_circle) {
        // update max count
        if(selected_size > max_count) {
            max_count = selected_size;
        }

        P p = points[index];
        // remove bad situation
        if(selected.empty() || abs(p - selected.front()) <= 2.0 + 1e-10) {
            // take
            selected.push_front(p);
            search(index + 1);
            selected.pop_front();
        }

        // for next dispatch
        // not take
        search(index + 1);
    }
}

int main() {
    int n;
    while (cin >> n) {
        if(n == 0) break;
        points.clear();
        points.shrink_to_fit();
        selected.clear();
        selected.shrink_to_fit();
        max_count = 0;

        // load points
        for (char i = 0; i < n; ++i) {
            double x;
            double y;
            cin >> x >> y;
            P p = P(x, y);
            points.push_back(p);
        }

        search(0);
        cout << max_count << endl;
    }
    return 0;
}

C++: #include<iostream>
#include<cmath>
#include<sstream>
using namespace std;
#define REP(i,b,n) for(int i=b;i<n;i++)
#define rep(i,n)   REP(i,0,n)
#define eps 1e-10
const double rs =1.;

int count(double x,double y,double *dx,double *dy,int n){
  int ret =0;
  rep(i,n){
    if ( (x-dx[i])*(x-dx[i])+(y-dy[i])*(y-dy[i]) <=1.+eps)ret++;
  }
  return ret;
}



int solve(double *x,double *y,int n){
  int cnt=1;
  rep(i,n){
    REP(j,i+1,n){
      double midx=(x[i]+x[j])/2.,midy=(y[i]+y[j])/2.;
      double tx=x[i]-midx,ty=y[i]-midy;
      double tr=sqrt(1-tx*tx-ty*ty);
      if (isnan(tr))continue;
      if ( tx==0 &&ty==0)continue;
      double rate = tr/sqrt(tx*tx+ty*ty);
      cnt=max(cnt,max(count(-ty*rate+midx,tx*rate+midy,x,y,n),
		      count(ty*rate+midx,-tx*rate+midy,x,y,n)));
    }
  }
  return cnt;
}





main(){
  int n;
  while(cin>>n&&n){
    double x[n],y[n];
    rep(i,n)cin>>x[i]>>y[i];
    cout << solve(x,y,n)<<endl;
  }

}

C++: #include "bits/stdc++.h"
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
const int INF = 1e9;
const ll LINF = 1e18;
template<class S,class T> ostream& operator << (ostream& out,const pair<S,T>& o){ out << "(" << o.first << "," << o.second << ")"; return out; }
template<class T> ostream& operator << (ostream& out,const vector<T> V){ for(int i = 0; i < V.size(); i++){ out << V[i]; if(i!=V.size()-1) out << " ";} return out; }
template<class T> ostream& operator << (ostream& out,const vector<vector<T> > Mat){ for(int i = 0; i < Mat.size(); i++) { if(i != 0) out << endl; out << Mat[i];} return out; }
template<class S,class T> ostream& operator << (ostream& out,const map<S,T> mp){ out << "{ "; for(auto it = mp.begin(); it != mp.end(); it++){ out << it->first << ":" << it->second; if(mp.size()-1 != distance(mp.begin(),it)) out << ", "; } out << " }"; return out; }

/*
 <url:http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1132>
 問題文============================================================
 =================================================================
 解説=============================================================
 ================================================================
 */

typedef complex<double> Point;
int n;
int solve(){
    int res = 1;
    vector<Point> ps(n);
    for(int i = 0; i < n;i++){
        double x,y; cin >> x >> y;
        ps[i] = Point(x,y);
    }
    for(int i = 0; i < n;i++){
        for(int j = i+1; j < n;j++){
            double dist = abs(ps[i]-ps[j]);
            if(dist > 2) continue;
            Point mp = (ps[i]+ps[j])/2.;
            Point dp = ps[i]-ps[j];
            double x = sqrt(1-dist*dist/4);
            Point cp = Point(-dp.imag()*x/dist+mp.real(),dp.real()*x/dist+mp.imag());
            int cnt = 0;
            for(int k = 0; k < n;k++){
                if(abs(cp-ps[k]) < 1.0 + 1e-9) cnt++;
            }
            res = max(res,cnt);
        }
    }
    return res;
}
int main(void) {
    cin.tie(0); ios_base::sync_with_stdio(false);
    while(cin >> n,n){
        cout << solve() << endl;
    }
    return 0;
}