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 <complex>
#include <iostream>
#include <cmath>
#include <vector>
using namespace std;

typedef complex<double> coor;

#define EPS 1e-5

double in_cir(coor o, coor x)
{
  return abs(x - o);
}

int search(vector<coor> P)
{
  int mx = 1;
  /*for (int i = 0; i < P.size(); i++) {
    printf("%lf, %lf\n", P[i].real(), P[i].imag());
    }*/
  for (int i = 0; i < P.size(); i++) {
    for (int j = i + 1; j < P.size(); j++) {
      for (int k = 0; k < 2; k++) {
	if (abs(P[i] - P[j]) - 2.0 > -EPS) {continue;}
	coor d = coor(P[i].imag() - P[j].imag(),
		      -(P[i].real() - P[j].real()));
	if (k == 1) {d = -d;}
	d /= abs(d);
	d *= sqrt(1 - norm(P[i] - P[j]) / 4);
	//printf("abs(d) = %lf\n", abs(d));
	d += (P[i] + P[j]) / 2.0;
	//printf("%lf, %lf\n", d.real(), d.imag());
	int ctr = 2;
	for (int l = 0; l < P.size(); l++) {
	  //printf("%lf\n", in_cir(d, P[l]));
	  if (l == i || l == j) {continue;}
	  if (in_cir(d, P[l]) - 1.0 < -EPS) {
	    ctr++;
	  }
	}
	//printf("_%d\n", ctr);
	mx = max(ctr, mx);
	
      }
    }
  }
  return mx;
}

int main()
{
  int N;
  while (scanf("%d ", &N) == 1) {
    if (N == 0) {
      break;
    }
    vector<coor> P(N);
    for (int i = 0; i < N; i++) {
      double x, y;
      scanf("%lf %lf ", &x, &y);
      P[i] = coor(x, y);
    }
    printf("%d\n", search(P));
  }
  return 0;
}

C++: #include <stdio.h>
#include <string.h>
#include <algorithm>
#include <iostream>
#include <math.h>
#include <assert.h>
#include <vector>

using namespace std;
typedef long long ll;
typedef unsigned int uint;
typedef unsigned long long ull;
static const double EPS = 1e-9;
static const double PI = acos(-1.0);

#define REP(i, n) for (int i = 0; i < (int)(n); i++)
#define FOR(i, s, n) for (int i = (s); i < (int)(n); i++)
#define FOREQ(i, s, n) for (int i = (s); i <= (int)(n); i++)
#define FORIT(it, c) for (__typeof((c).begin())it = (c).begin(); it != (c).end(); it++)
#define MEMSET(v, h) memset((v), h, sizeof(v))

#include <complex>
typedef complex<double> Point;
// O(n^2logn)
int MaxPointCircle(const vector<Point> points, double r) {
  const int n = points.size();
  int ret = 1;
  for (int i = 0; i < n; i++) {
    vector<pair<double, int> > radians;
    for (int j = 0; j < n; j++) {
      double d = abs(points[i] - points[j]) / 2.0;
      if (i == j || d * 2.0 - r * 2.0 > EPS) { continue; }
      Point center = (points[i] + points[j]) / 2.0;
      Point vect = (points[j] - points[i]) * Point(0, -1);
      vect /= abs(vect);
      double l = sqrt(r * r - d * d);
      Point c1 = center + vect * l;
      Point c2 = center - vect * l;
      double arg1 = arg(c1 - points[i]) - EPS;
      double arg2 = arg(c2 - points[i]) + EPS;
      radians.push_back(make_pair(arg1, (j + 1)));
      radians.push_back(make_pair(arg2, -(j + 1)));
    }
    if ((int)radians.size() / 2 < ret - 1) { continue; }
    sort(radians.begin(), radians.end());
    vector<bool> used(points.size(), false);
    int lsum = 1;
    REP(iter, 2) {
      REP(j, radians.size()) {
        int index = abs(radians[j].second) - 1;
        bool start = radians[j].second > 0 ? true : false;
        if (used[index] && !start) {
          used[index] = false;
          lsum--;
        } else if (!used[index] && start) {
          used[index] = true;
          lsum++;
        }
        ret = max(ret, lsum);
      }
    }
  }
  return ret;
}

int main() {
  int n;
  while (scanf("%d", &n), n) {
    vector<Point> points(n);
    REP(i, n) {
      double x, y;
      scanf("%lf %lf", &x, &y);
      points[i] = Point(x, y);
    }
    printf("%d\n", MaxPointCircle(points, 1.0));
  }
}

C++: #include <bits/stdc++.h>
#define EPS 1e-9
using namespace std;

typedef complex<double> P;

P rot (P x, P p, double ang) {
    ang = ang / 180.0 * M_PI;
    double p1 = p.real() - x.real();
    double p2 = p.imag() - x.imag();
    return P( p1 * cos(ang) - p2 * sin(ang) + x.real()
              , p1 * sin(ang) + p2 * cos(ang) + x.imag() );
}

int main () {
    int n;
    while (cin >> n, n) {
        vector<P> v(n);
        for (auto &i : v) {
            double x, y;
            cin >> x >> y;
            i = P(x,y);
        }

        int res = 1;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                if (i != j && abs(v[i] - v[j]) <= 2. + EPS) {
                    int count = 0;
                    P m = P(abs(v[i].real() + v[j].real()) / 2., abs(v[i].imag() + v[j].imag()) / 2.);
                    P p = m + (m - v[i]) / abs(m - v[i]) * sqrt(1. - pow(abs(m - v[i]), 2.));
                    p = rot(m, p, 90);
                    for (int k = 0; k < n; k++) {
                        if (abs(p - v[k]) <= 1. + EPS) count++;
                    }

                    res = max(res, count);
                }
            }
        }

        cout << res << endl;
    }
    return 0;
}

C++: #include <algorithm>
#include <cmath>
#include <complex>
#include <iostream>
#include <vector>
using namespace std;

// 問題名："Circle and Points"
// https://onlinejudge.u-aizu.ac.jp/problems/1132
// 所要時間… 考察：10分程度 実装：20分程度

/*
***方針***
最大数囲む半径1の円の囲み方が有るときに、囲む点を変えずに少なくとも2点を同時に円周上に持ってくることができる。したがって、周上に2点が存在するような円のみを考えれば良い。周上に置く2点を決めた時、考えられる半径1の円の囲み方は最大2通り存在する（2点間の距離が2より大きいと存在しない）。この2点の選び方はNC2=O(N^2)通りなので、それぞれについて円で囲まれる点を一つづつ数えると、計算量はO(N^3)となりN<=300で間に合いそうである。

***学んだこと***
complex型によって、平面上の点の操作が容易となる（複素平面的扱いにより変数の数を減らすことが可能）。
*/

int main() {
  int n;
  while (cin >> n, n > 0) {
    // 入力
    vector<complex<double>> p(n);
    double xt, yt;
    for (int i = 0; i < n; i++) {
      cin >> xt >> yt;
      p[i].real(xt);
      p[i].imag(yt);
    }

    int ret =
        1;  // 1<=Nなので答えは必ず1以上（0とすると、下のループに入らない場合にバグる）
    for (int i = 0; i < n; i++) {  // 円周に来る2点の組でループ
      for (int j = i + 1; j < n; j++) {
        if (abs(p[i] - p[j]) > 2)
          continue;  // 2点間の距離が2より大きかったら無視

        // 回転操作と実数倍で、2つの点p[i],p[j]から円の中心を二つ求める（c+pmとc-pm）
        complex<double> pm = (p[i] + p[j]) / 2.0;
        complex<double> c = p[i] - pm;
        double d = abs(c);
        c *= sqrt(1 - d * d) / d * 1.0i;
        int cntp = 0, cntm = 0;

        // できる2つの円について、含まれる点の個数を数える
        // ちょうど周上に来る2点を確実に入れるために、ボーダーを少し広めに設定（制約上、外の点を含むことはない）
        for (int k = 0; k < n; k++) {
          if (abs(p[k] - (c + pm)) < 1.00005) cntp++;
          if (abs(p[k] - (-c + pm)) < 1.00005) cntm++;
        }
        ret = max({ret, cntp, cntm});  // 大きいものを採用
      }
    }
    cout << ret << endl;
  }
}