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<math.h>
#include<vector>
using namespace std;

// 点
class Point {
  public:
    float x, y;
    Point(double ix = 0.0, double iy = 0.0) {
      set(ix, iy);
    }
    void set(double ix = 0.0, double iy = 0.0) {
      x = ix;
      y = iy;
    }
    double length() {
      return sqrt(x * x + y * y);
    }
    double distanceTo(Point p) {
      double dx = p.x - x, dy = p.y - y;
      return sqrt(dx * dx + dy * dy);
    }
};

int main() {
  vector<int> results;
  while (true) {
    // 入力
    int n;
    cin >> n;
    if (n == 0) break;
    vector<Point> p(n);
    for (int i=0; i<n; i++) {
      double x, y;
      cin >> x >> y;
      p[i].set(x, y);
    }
    
    // とりあえず格子状に調べてみる（）
    const float X = 1000;
    Point pc;
    int max = 0;
    for (int xi=0; xi<=X; xi++) {
      for (int yi=0; yi<=X; yi++) {
	float x = xi * 10 / X, y = yi * 10 / X; // 0 ~ 10に補正
	// 距離が1以下のものをカウント
	pc.set(x, y);
	int count = 0;
	for (int i=0; i<n; i++) {
	  if (pc.distanceTo(p[i]) <= 1) count++;
	}
	if (count > max) max = count;
      }
    }
    results.push_back(max);
  }
  
  int size = results.size();
  for (int i=0; i<size; i++) {
    cout << results[i] << endl;
  }
  return 0;
}

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

#define rep(i,N) for(int i=0;i<(int)(N);i++)
#define REP(i,k,N) for(int i=k;i<(int)(N);i++)

#define vi vector<int>
#define pb push_back

#define pdd pair<double,double>
#define fi first
#define se second
const double  EPS = 1e-9;

double dist(double x,double y){
    return x*x+y*y;
};

int main(){
    int n;
    while(cin >> n,n){
        int ma = 0;
        pdd p[330];
        rep(i,n) cin >> p[i].fi >> p[i].se;
        rep(i,n-1){
            REP(j,i+1,n){
                int cnt1 = 0,cnt2=0;
                double x1,y1,x2,y2,x3,y3;
                tie(x1,y1) = p[i];
                tie(x2,y2) = p[j];
                double cx1,cy1,cx2,cy2;
                if(dist(x1-x2,y1-y2)>4) continue;

                //はりまコーナー
//                double ax,ay,bx1,by1,bx2,by2,mx,my,dx,dy;
//                ax=x1-x2;ay=y1-y2;
//                mx=(x1+x2)/2;my=(y1+y2)/2;
//                by1 = (-0.5)*pow(ax,4)+pow(ax,2)+(-0.5)*pow(ax,2)*pow(ay,2)
//                        /dist(ax,ay);
//                by1 = sqrt(abs(by1));
//                by2 = -1*by1;
//                cy1= my+by1;
//                cy2= my+by2;
//                bx1= -1*(ay*by1)/ax;
//                bx2= -1*(ay*by2)/ax;
//                cx1=bx1+mx;
//                cx2=bx2+mx;


                //あじのりコーナー
                double mx=(x1+x2)/2,my=(y1+y2)/2;

                double ld = 1-dist((x1-x2)/2,(y1-y2)/2);
                ld = sqrt(ld);
                cx1 = mx - sin(atan((y2-y1)/(x2-x1)))*ld;
                cy1 = my + cos(atan((y2-y1)/(x2-x1)))*ld;
                cx2 = mx + sin(atan((y2-y1)/(x2-x1)))*ld;
                cy2 = my - cos(atan((y2-y1)/(x2-x1)))*ld;

                for(int k = 0;k<n;k++) {
                    tie(x3, y3) = p[k];
                    if (EPS > dist((cx1 - x3), (cy1 - y3))-1) cnt1++;
                    if (EPS > dist((cx2 - x3), (cy2 - y3))-1) cnt2++;
                }
                ma = max({cnt1,cnt2,ma});
            }
        }
        cout<< ma << endl;
    }
    return 0;
}