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 complex<double> C;

int solve(int N, vector<double>& x, vector<double>& y) {
    vector<C> point(N);

    for (int j = 0; j < N; ++j) {
        point[j] = C(x[j], y[j]);
    }

    int ret = 0;

    for (int j = 0; j < N; ++j) {
        for (int k = j+1; k < N; ++k) {
            C mid = (point[j] + point[k]) * 0.5;

            C normal = (point[j] - point[k]) * C(0, 1);
            normal /= abs(point[j] - point[k]);

            double a = 1 - norm(mid - point[j]);

            C centre1 = mid + a * normal;
            int count = 0;
            for (int l = 0; l < N; ++l) {
                if (abs(centre1 - point[l]) < 1) {
                    ++count;
                }
            }
            ret = max(ret, count);

            C centre2 = mid - a * normal;
            count = 0;
            for (int l = 0; l < N; ++l) {
                if (abs(centre2 - point[l]) < 1.0001) {
                    ++count;
                }
            }
            ret = max(ret, count);

        }
    }

    return ret;
}

int main () {
    int N;
    vector<double> x, y;

    while (true) {
        cin >> N;
        if (N == 0) break;

        x.resize(N);
        y.resize(N);

        for (int j = 0; j < N; ++j) {
            cin >> x[j] >> y[j];
        }

        cout << solve(N, x, y) << endl;
    }

    return 0;
}

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 <complex>
#include <vector>
#include <cmath>
#define EPS 1.0e-10
using namespace std;

#define X real()
#define Y imag()
typedef complex<double> P;
 
vector<P> p;
int check(P a){
	int ret=0;
	for(int i=0;i<p.size();i++) if(abs(a-p[i]) < 1+EPS) ret ++;
	return ret;
}

int adjcircle(P a, P b){
	P dir = a-b, pos=(a+b)*0.5;
	dir = P(dir.Y, -dir.X);
	dir = dir*(1/abs(dir))*sqrt(1-abs(pos-a)*abs(pos-a));
	return max(check(pos+dir), check(pos-dir));
}

main(){
	int i,j,N;
	while(cin >> N, N){
		p = vector<P>(N);
		int ans = 1;
		for(i=0;i<N;i++) cin >> p[i].X >> p[i].Y;
		for(i=1;i<N;i++)
			for(j=0;j<i;j++)
				if(abs(p[i]-p[j])<=2.0) ans = max(ans, adjcircle(p[i], p[j]));
		cout << ans << endl;
	}
}