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 <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <climits>
#include <cfloat>
#include <ctime>
#include <map>
#include <utility>
#include <set>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>
#include <sstream>
#include <complex>
#include <stack>
#include <queue>
#include <numeric>
#include <list>

using namespace std;

#ifdef _MSC_VER
#define __typeof__ decltype
template <class T> int __builtin_popcount(T n) { return n ? 1 + __builtin_popcount(n & (n - 1)) : 0; }
#endif

#define foreach(it, c) for (__typeof__((c).begin()) it=(c).begin(); it != (c).end(); ++it)
#define all(c) (c).begin(), (c).end()
#define rall(c) (c).rbegin(), (c).rend()
#define popcount __builtin_popcount
#define CLEAR(arr, val) memset(arr, val, sizeof(arr))

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

template <class T> void max_swap(T& a, const T& b) { a = max(a, b); }
template <class T> void min_swap(T& a, const T& b) { a = min(a, b); }

typedef long long ll;
typedef pair<int, int> pint;

const double PI = acos(-1.0);
const int dx[] = { 0, 1, 0, -1 };
const int dy[] = { 1, 0, -1, 0 };



typedef pair<double, double> point;
const double EPS = 1e-6;
int main()
{
	int n;
	while (scanf("%d", &n), n)
	{
		point p[333];
		rep (i, n)
			scanf("%lf%lf", &p[i].first, &p[i].second);
		sort(p, p + n);

		queue<point> q;
		rep (i, n)
			q.push(p[i]);

		list<point> r;
		int res = 0;
		for (int i = -10000; i <= 100000; ++i)
		{
			double cx = i * 1e-4;
			while (!r.empty() && r.front().first < cx - 1 + EPS)
				r.pop_front();
			while (!q.empty() && q.front().first < cx + 1 + EPS)
				r.push_back(q.front()), q.pop();

			vector<double> cy;
			foreach (it, r)
			{
				double x = it->first, y = it->second;
				double d = sqrt(1 - (x - cx) * (x - cx));
				cy.push_back(y - d);
				cy.push_back(y + d);
			}

			for (int i = 0; i < cy.size(); ++i)
			{
				int c = 0;
				foreach (it, r)
				{
					double dx = cx - it->first;
					double dy = cy[i] - it->second;
					if (dx * dx + dy * dy < 1 + EPS)
						++c;
				}
				max_swap(res, c);
			}
		}
		printf("%d\n", res);
	}
}

C++: #include <stdio.h>
#include <algorithm>
#include <iostream>
#include <stack>
#include <vector>
#include <math.h>
#include <queue>

using namespace std;

struct Point{
    double x,y;
};

void func(int N){

    Point point[N];
    int x_table[1001],y_table[1001],x_start,x_end,y_start,y_end;

    for(int i = 0; i <= 1000; i++){
    	x_table[i] = y_table[i] = 0;
    }

    for(int i = 0; i < N; i++){
        scanf("%lf %lf",&point[i].x,&point[i].y);
        x_start = point[i].x*100-50;
        x_start = max(49,x_start);
        x_end = point[i].x*100+50;
        x_end = min(x_end,901);
        y_start = point[i].y*100-50;
        y_start = max(49,y_start);
        y_end = point[i].y*100+50;
        y_end = min(y_end,901);

        for(int k = x_start; k <= x_end; k++){
        	x_table[k]++;
        }

        for(int k = y_start; k <= y_end;k++){
        	y_table[k]++;
        }
    }


    int count,max_count = 0;
    double calc_x,calc_y;

    for(int x = 49; x <= 951; x++){
    	for(int y = 49; y <= 951; y++){
    		if(x_table[x] > 0 && y_table[y] > 0){
    			calc_x = x,calc_y = y;
    			calc_x /= 100,calc_y/=100;
    			count = 0;
    			for(int i = 0; i < N; i++){
    				if(sqrt((calc_x-point[i].x)*(calc_x-point[i].x)+(calc_y-point[i].y)*(calc_y-point[i].y)) <= 1){
    					count++;
    				}
    			}
    			max_count = max(max_count,count);
    		}
    	}
    }
    printf("%d\n",max_count);
}


int main(){

    int N;

    while(true){
        scanf("%d",&N);
        if(N == 0)break;

        func(N);
    }


    return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;
typedef complex<double> P;

const double eps = 1e-5;
const int MAXN = 333;
int N;
P ps[MAXN];

int calc(P c) {
  int res = 0;
  for(int i = 0; i < N; ++i) {
    double d = norm(c - ps[i]);
    if(d < 1.0 + eps) ++res;
  }
  return res;
}

int main() {
  for(; cin >> N && N; ) {
    for(int i = 0; i < N; ++i) {
      cin >> ps[i].real() >> ps[i].imag();
    }
    int res = 0;
    for(int i = 0; i < N; ++i) {
      for(int j = i+1; j < N; ++j) {
        P base = ps[j] - ps[i];
        P m = ps[i] + base / 2.0;
        double d = abs(base);
        double a = sqrt(1.0 - d * d / 4.0);
        P base2(-base.imag(), base.real());
        base2 = base2 / abs(base2) * a;
        res = max(res, max(calc(m + base2), calc(m - base2)));
      }
    }
    cout << res << endl;
  }
  return 0;
}