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++: //http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1132&lang=jp
#include <bits/stdc++.h>
using namespace std ;

#define pb(n) push_back(n)
#define fi first
#define se second
#define all(r) begin(r),end(r)
#define mp(a,b) make_pair(a,b)

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

using ll = long long;

const double EPS = 1e-6;
const double INF = 1e12;
  typedef complex<double> P;//????´???°????????¢?????????????????¨?§£???
  namespace std {
    bool operator < (const P& a, const P& b) {//x????????????
      return real(a) != real(b) ? real(a) < real(b) : imag(a) < imag(b);
    }
  }
  double cross(const P& a, const P& b) {//??????
    return imag(conj(a)*b);
  }
  double dot(const P& a, const P& b) {//??????
    return real(conj(a)*b);
  }

  struct L : public vector<P> {//??´???
    L(){};
    L(const P &a, const P &b) {
      push_back(a); push_back(b);
    }
  };

  typedef vector<P> G;

  struct C {//???
    P p; double r;
    C(){};
    C(const P &p, double r) : p(p), r(r) { }
  };
  int ccw(P a, P b, P c) {//3????????????????????§????????§??????????????????
    b -= a; c -= a;
    if (cross(b, c) > 0)   return +1;       // counter clockwise
    if (cross(b, c) < 0)   return -1;       // clockwise
    if (dot(b, c) < 0)     return +2;       // c--a--b on line
    if (norm(b) < norm(c)) return -2;       // a--b--c on line
    return 0;
  }
  // L:line,S:segment,P:point
  bool intersectLL(const L &l, const L &m) {
    return abs(cross(l[1]-l[0], m[1]-m[0])) > EPS || // non-parallel
      abs(cross(l[1]-l[0], m[0]-l[0])) < EPS;   // same line
  }
  bool intersectLS(const L &l, const L &s) {
    return cross(l[1]-l[0], s[0]-l[0])*       // s[0] is left of l
      cross(l[1]-l[0], s[1]-l[0]) < EPS; // s[1] is right of l
  }
  bool intersectLP(const L &l, const P &p) {
    return abs(cross(l[1]-p, l[0]-p)) < EPS;
  }
  bool intersectSS(const L &s, const L &t) {
    return ccw(s[0],s[1],t[0])*ccw(s[0],s[1],t[1]) <= 0 &&
      ccw(t[0],t[1],s[0])*ccw(t[0],t[1],s[1]) <= 0;
  }
  bool intersectSP(const L &s, const P &p) {
    return abs(s[0]-p)+abs(s[1]-p)-abs(s[1]-s[0]) < EPS; // triangle inequality
  }
  // L & P intersection
  P projection(const L &l, const P &p) {
    double t = dot(p-l[0], l[0]-l[1]) / norm(l[0]-l[1]);
    return l[0] + t*(l[0]-l[1]);
  }
  // L????????????P???????°????
  P reflection(const L &l, const P &p) {
    return p + 2.0 * (projection(l, p) - p);
  }
  double distanceLP(const L &l, const P &p) {
    return abs(p - projection(l, p));
  }
  double distanceLL(const L &l, const L &m) {
    return intersectLL(l, m) ? 0 : distanceLP(l, m[0]);
  }
  double distanceLS(const L &l, const L &s) {
    if (intersectLS(l, s)) return 0;
    return min(distanceLP(l, s[0]), distanceLP(l, s[1]));
  }
  double distanceSP(const L &s, const P &p) {
    const P r = projection(s, p);
    if (intersectSP(s, r)) return abs(r - p);
    return min(abs(s[0] - p), abs(s[1] - p));
  }
  double distanceSS(const L &s, const L &t) {
    if (intersectSS(s, t)) return 0;
    return min(min(distanceSP(s, t[0]), distanceSP(s, t[1])),
        min(distanceSP(t, s[0]), distanceSP(t, s[1])));
  }
  double distancePP(const P &p,const P &q){
    return dot((p - q), (p - q));
  }
  P crosspoint(const L &l, const L &m) {
    double A = cross(l[1] - l[0], m[1] - m[0]);
    double B = cross(l[1] - l[0], l[1] - m[0]);
    if (abs(A) < EPS && abs(B) < EPS) return m[0]; // same line
    if (abs(A) < EPS) assert(false); // !!!PRECONDITION NOT SATISFIED!!!
    return m[0] + B / A * (m[1] - m[0]);
  }
  P verticalvector(const P &a,const P &b){//2???????????????????????????????????????
    if(a.real() == b.real()){
      return {1.0,0};
    }
    if(a.imag() == b.imag()){
      return {0,1.00};
    }
    P p = a - b;
    P res = {-p.imag(),p.real()};
    return res / abs(p);
  }

int main(){
  cin.tie(0);
  ios::sync_with_stdio(false);
  ll n;
  double x,y;
  while(cin >> n && n){
    vector<P> ps(n);
    rep(i,n){
      cin >> x >> y;
      ps[i] = {x,y};
    }
    ll ans = 1;
    rep(i,n){
      REP(j,i+1,n){
        if (norm(ps[i] - ps[j]) > 4.0 + EPS) continue;
        P m = (ps[i]+ps[j])/2.0;
        P vv = verticalvector(ps[i],ps[j]);
        double x = sqrt(1 - norm((ps[i] - ps[j]) / 2.0));
        P c1 = vv * x + m,c2 = vv * (-x) + m;
        ll cnt1 = 2,cnt2 = 2;
        rep(k,n){
          if(i == k || j == k) continue;
          if(norm(ps[k] - c1) <= 1.0 + EPS){
            cnt1++;
          }
          if(norm(ps[k] - c2) <= 1.0 + EPS){
            cnt2++;
          }
        }
        ans = max(ans,max(cnt1,cnt2));
      }
    }
    cout << ans << endl;
  }
  return 0;
}

C++: #include <iostream>
#include <sstream>
#include <iomanip>
#include <algorithm>
#include <cmath>
#include <string>
#include <vector>
#include <list>
#include <queue>
#include <stack>
#include <set>
#include <map>
#include <bitset>
#include <numeric>
#include <cfloat>
using namespace std;

const double EPS = 1.0e-9;

class Point
{
public:
    double y, x;
    Point(){
        y = x = 0.0;
    }
    Point(double y0, double x0){
        y = y0;
        x = x0;
    }
    const Point operator+(const Point& p) const{
        return Point(y + p.y, x + p.x);
    }
    const Point operator-(const Point& p) const{
        return Point(y - p.y, x - p.x);
    }
    const Point operator*(double a) const{
        return Point(y * a, x * a);
    }
    double dist(const Point& p) const{
        return sqrt(pow(y-p.y, 2.0) + pow(x-p.x, 2.0));
    }
};

int solve(vector<Point>& vp)
{
    int n = vp.size();
    int ret = 1;

    for(int i=0; i<n; ++i){
        for(int j=i+1; j<n; ++j){
            double dist = vp[i].dist(vp[j]);
            if(dist > 2.0)
                continue;

            for(int k=0; k<2; ++k){
                Point p = vp[j] - vp[i];
                Point p1(p.x, p.y);
                if(k == 0)
                    p1.x *= -1;
                else
                    p1.y *= -1;
                p = vp[i] + p * 0.5 + p1 * (sqrt(pow(dist, -2) - 0.25));
                int num = 0;
                for(int l=0; l<n; ++l){
                    if(p.dist(vp[l]) < 1.0 + EPS)
                        ++ num;
                }
                ret = max(ret, num);
            }
        }
    }

    return ret;
}

int main()
{
    for(;;){
        int n;
        cin >> n;
        if(n == 0)
            break;

        vector<Point> p(n);
        for(int i=0; i<n; ++i)
            cin >> p[i].x >> p[i].y;
        cout << solve(p) << endl;
    }
}