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<vector>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cmath>
using namespace std;

using ll = long long;
#define rep(i, j) for(int i=0; i < (int)(j); i++)
#define repeat(i, j, k) for(int i = (j); i < (int)(k); i++)
#define all(v) v.begin(),v.end()
#define debug(x) cerr << #x << " : " << x << endl

const double EPS = 1e-8;

// vector
template<class T> istream& operator >> (istream &is , vector<T> &v) {
    for(T &a : v) is >> a; return is;
}
using Coordinate = double;
bool eq(Coordinate a, Coordinate b) {
    return abs(a - b) < EPS;
}
/*
 * Point
 */
struct Point {
    Coordinate x,y;
    Point():x(0),y(0){}
    Point(Coordinate x,Coordinate y):x(x),y(y){}
    pair<Coordinate,Coordinate> to_pair() const { return make_pair(x,y); }
    bool operator == (const Point &a) const { return eq(x, a.x) && eq(y, a.y); }
    bool operator != (const Point &a) const { return !(*this == a); }
    
    Point operator + (const Point& p) const { return Point(x + p.x, y + p.y); }
    Point operator - (const Point& p) const { return Point(x - p.x, y - p.y); }
    Point operator - () const { return Point(-x, -y); }
    Point& operator += (const Point& p) { x += p.x; y += p.y; return *this; }
    Point& operator -= (const Point& p) { x -= p.x; y -= p.y; return *this; }
};
template<typename T> Point operator * (T k, const Point &p) { return Point(k * p.x, k * p.y); }
template<typename T> Point operator / (const Point &p, T k) { return Point(p.x / k, p.y / k); }
double dot(const Point &a, const Point &b) { return a.x * b.x + a.y * b.y; }
double cross(const Point &a, const Point &b) { return a.x * b.y - a.y * b.x; }
double norm(const Point& a) { return sqrt( dot(a, a) ); }

istream &operator >> (istream& is, Point& p){ return is >> p.x >> p.y;}

double distance(const Point& a, const Point& b) { return norm(a - b); }
double distance2(const Point& a, const Point& b) { return dot(a - b, a - b); }

// ??£?????´???????????????
Point orthonormal(Point &p) {
    Point p2 = Point(-1 * p.y, p.x); //rotate(p, 90);
    return p2 / norm(p2);
}

Point midPoint(const Point &a, const Point &b) {
    return Point((a.x + b.x) / 2, (a.y + b.y) / 2);
}


bool solve() {
    int n; cin >> n;
    if(n == 0) return 0;
    vector<Point> P(n); cin >> P;
    int ans = 1;
    rep(i, n) rep(j, i) if(distance(P[i], P[j]) <= 2.0 + EPS){
        Point a = P[i] - P[j];
        Point mid = midPoint(P[i], P[j]);
        double len = sqrt(1 - distance2(mid, P[j]));
        Point c = len * orthonormal(a) + mid;
        int cnt = 0;
        rep(l, n) {
            cnt += distance2(P[l], c) <= 1.0 + EPS;
        }
        ans = max(ans, cnt);
    }
    cout << ans << endl;
    return 1;
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    while(solve());
    return 0;
}

C++: #include <algorithm>
#include <cassert>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <deque>
#include <functional>
#include <list>
#include <map>
#include <queue>
#include <set>
#include <sstream>
#include <string>
#include <utility>
#include <vector>

#define REP(i,s,n) for(int i=(int)(s);i<(int)(n);i++)
#define DEBUGP(val) cerr << #val << "=" << val << "\n"

using namespace std;
typedef long long int ll;
typedef vector<int> VI;
typedef vector<ll> VL;
typedef pair<int, int> PI;
typedef pair<double, double> PD;

const double EPS = 1e-11;
double sq(double x) { return x * x; }

void solve(double ax, double ay, double bx, double by, int &ma,
		 const vector<double> &x, const vector<double> &y) {
  double sqdist = sq(ax - bx) + sq(ay - by);
  if (sqdist >= 4 + EPS) return;
  double dist = sqrt(sqdist);
  double ndist = sqrt(1 - sqdist / 4);
  double ex = (bx - ax) / dist, ey = (by - ay) / dist;
  int n = x.size();
  REP(i, 0, 2) {
    double coef = i == 0 ? 1 : -1;
    double cx = (ax + bx) / 2 + ey * coef * ndist;
    double cy = (ay + by) / 2 - ex * coef * ndist;
    int cnt = 0;
    REP(j, 0, n) {
      double sqdist = sq(x[j] - cx) + sq(y[j] - cy);
      if (sqdist <= 1 + EPS) {
	cnt += 1;
      }
    }
    ma = max(ma, cnt);
  }
}

int main(void) {
  int n;
  while(scanf("%d", &n) && n) {
    vector<double> x(n), y(n);
    REP(i, 0, n) scanf("%lf%lf", &x[i], &y[i]);
    int ma = 1;
    REP(i, 0, n) {
      REP(j, 0, i) {
        solve(x[i], y[i], x[j], y[j], ma, x, y);
      }
    }
    printf("%d\n", ma);
  }
}

C++: #include <iostream>
#include <deque>
#include <string>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <string>
#include <cmath>
#include <complex>

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

using namespace std;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef complex<double> vec;

double const pi = acos(-1);
double const eps = 1e-9;

int main(){
    int n;
    while(cin >> n && n){
        vector<vec> v(n);
        rep(i, n){
            double x, y; cin >> x >> y;
            vec p(x, y);
            v[i] = p;
        }
        int ans = 0;
        rep(i, n)rep(j, i){
            double d = abs(v[i] - v[j]);
            if (d - eps > 2) continue;
            double h = 1 - d*d/4;
            vec c1 = polar(h, arg(v[i] - v[j]) + pi / 2) + (v[i] + v[j]) / 2.;
            vec c2 = polar(h, arg(v[i] - v[j]) - pi / 2) + (v[i] + v[j]) / 2.;
            int a1 = 0, a2 = 0;
            rep(k, n){
                if (abs(c1 - v[k]) <= 1 + eps) a1++;
                if (abs(c2 - v[k]) <= 1 + eps) a2++;
            }
            ans = max(ans, a1);
            ans = max(ans, a2);
        }
        cout << ans << endl;
    }
}