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;

#define EPS 1e-8
typedef complex<double> Point;

int main(){
  int N;
  while(cin >> N,N){
    Point p[N];
    for(int i = 0 ; i < N ; i++){
      cin >> p[i].real() >> p[i].imag();
    }
    int res = 1;
    for(int i = 0 ; i < N-1 ; i++){
      for(int j = i+1 ; j < N ; j++){
        Point np = p[j]-p[i];
        double d = abs(np);
        if(d + EPS > 2.0){ continue; }
        np /= d;
        np *= Point(0,1);
        Point cp = (p[i]+p[j])/2.0;
        Point v = sqrt(1.0-d*d/4.0)*np;
        int cnt = 2;
        Point center = cp + v;
        for(int k = 0 ; k < N ; k++){
          if(i == k || j == k) continue;
          Point diff = center-p[k];
          if(norm(diff) < 1.0 + EPS){
            cnt++;
          }
        }
        res = max(res,cnt);
      }
    }
    cout << res << endl;
  }
  return 0;
}

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

using namespace std;

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

double x[400], y[400];

double EPS = 1e-8;

int main() {
  int n;
  for (;;) {
    cin >> n;
    if (n == 0) break;
    int res = 1;
    rep (i, n) cin >> x[i] >> y[i];
    rep (i, n) rep (j, n) if (i < j && sqrt(pow(x[i] - x[j], 2) + pow(y[i] - y[j], 2)) <= 2 + EPS) {
      double mx = (x[i] + x[j]) / 2, my = (y[i] + y[j]) / 2;
      double nx = -(my - y[i]), ny = (mx - x[i]);
      double d = sqrt(max(1 - (pow(x[i] - mx, 2) + pow(y[i] - my, 2)), 0.));
      double dx = nx / sqrt(nx * nx + ny * ny) * d;
      double dy = ny / sqrt(nx * nx + ny * ny) * d;
      double xx = mx + dx, yy = my + dy;
      //cout << i <<  " " << j << endl;
      //cout << d << endl;
      //cout << dx << " " << dy << endl;
      //cout << xx << " " << yy << endl;
      int count = 0;
      rep (k, n) {
	//cout << pow(x[k] - xx, 2) + pow(y[k] - yy, 2) << endl;
	if (pow(x[k] - xx, 2) + pow(y[k] - yy, 2) <= 1 + EPS * 10) ++count;
      }
      res = max(res, count);
      xx = mx - dx, yy = mx - dy;
      count = 0;
      rep (k, n) {
	//cout << pow(x[k] - xx, 2) + pow(y[k] - yy, 2) << endl;
	if (pow(x[k] - xx, 2) + pow(y[k] - yy, 2) <= 1 + EPS * 10) ++count;
      }
      res = max(res, count);
    }
    cout << res << endl;
  }
}