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<cstdio>
#include<cmath>
#include<algorithm>
#include<complex>
#include<vector>

using namespace std;

typedef complex<double> point;

#define dist2(p1,p2) (real(((p1)-(p2))*conj((p1)-(p2))))

#define mid(p1,p2) ((p1)+(p2))/2.0

#define PI (3.1415926535)
#define EPS (1e-6)

const point ang90 = point(0,1.0);

bool test (const point &a,const point &b){
    return real(a) < real(b);
}

int main(){
  int n;
  point p[300];
  while(scanf("%d",&n),n){
    for(int i=0;i<n;i++){
      double x,y;
      scanf("%lf%lf",&x,&y);
      p[i]=point(x,y);
    }
    sort(p,p+n,test);
    int ans=1;
    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
        if(real(p[j])-real(p[i]) > 2.0+EPS) break;
        if(dist2(p[i],p[j])>4.0) continue;
        point tmp = (p[i]-p[j]) * ang90;
        tmp=(tmp/abs(tmp))*sqrt(1.0-dist2(p[i],p[j])/4.0);
        point p1 = mid(p[i],p[j]) + tmp;
        int count=0;
        double rp=real(p1);
        bool f=false;
        for(int k=0;k<n;k++){
          if(k==i||k==j){
            count++;
          }else{
            if(f||(rp - real(p[k]) > 1.0 + EPS)) continue;
            f=true;
            if(real(p[k]) - rp > 1.0 + EPS) break;
            if(dist2(p1,p[k])<1.0+EPS){
              count++;
            }
          }
        }
        ans = max(count,ans);
        p1 = mid(p[i],p[j]) - tmp;
        rp=real(p1);
        count=0;
        f=false;
        for(int k=0;k<n;k++){
          if(k==i||k==j){
            count++;
          }else{
            if(f||(rp - real(p[k]) > 1.0 + EPS)) continue;
            if(real(p[k]) - rp > 1.0 + EPS) break;
            if(dist2(p1,p[k])<1.0+EPS){
              count++;
            }
          }
        }
        ans = max(count,ans);
      }
    }
    cout << ans << endl;
  }
  return 0;
}

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

using namespace std;

#define Sq(x) ((x)*(x))

#define Foreach(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();i++)

const double PI = M_PI;
const double EPS = 0.0001;

struct Point{
  double x,y;

  Point(){}
  Point(double x, double y) : x(x), y(y) {}
};

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

Point operator-(const Point &a, const Point &b){
  return Point(a.x - b.x, a.y - b.y);
}

Point operator*(const Point &a, const double k){
  return Point(a.x * k, a.y * k);
}

Point operator*(const Point &a,const Point &b){
  return Point(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);
}

Point operator/(const Point &p, const double 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 &p){
  return dot(p, p);
}

double abs(const Point &p){
  return sqrt(norm(p));
}

typedef Point Vector;

struct Circle : Point{
  double r;
  Circle(Point p = Point(0,0), double r = 0) : Point(p), r(r) {}
};

double distPP(const Point &a, const Point &b) {
  return sqrt(Sq(a.x-b.x) + Sq(a.y-b.y));
}

bool isPointInCircle(const Point &p, const Circle &c) {
  return abs(p - c) <= c.r + EPS;
}
 
vector<Point> crossPointCC(const Circle &c1, const Circle &c2) {
  vector<Point> res;
  
  if(abs(c1-c2) < EPS) return vector<Point>();
  
  double d = abs( c1 - c2 );
  double rc = ( d*d+c1.r * c1.r-c2.r * c2.r ) / ( 2*d );
  double rs = sqrt(c1.r*c1.r-rc*rc);
  Point diff = ( c2 - c1 ) / d;
  res.push_back( Point(c1 + diff * Point(rc,rs)) );
  res.push_back( Point(c1 + diff * Point(rc,-rs)) );
  return res;
}

Point rotate90(const Point &p) {
  return Point(-p.y, p.x);
}

/*
テ」ツ??」ツ?セテ」ツ?湘・ツ仰陛」ツ?凝」ツ?ェテ」ツ??
Point rotateVector(const Vector &p, double th) {
  return Point(p.x * cos(th) - p.y*sin(th), p.x * sin(th) + p.y * cos(th));
}
*/

int main() {
  int n;
  Point p[300];
  
  while(cin >> n && n) {
    for(int i=0; i<n; i++) {
      cin >> p[i].x >> p[i].y;
    }
    
    int ans = 1;
    for(int i=0; i<n-1; i++) {
      for(int j=i+1; j<n; j++) {
	if(distPP(p[i], p[j]) + EPS > 2.0) continue;
	
	Point cent((p[i]+p[j])/2);
	double a = distPP(p[i], cent);
	double b = sqrt(1.0 - a*a);
	
	Vector v1 = p[i] - cent;
	Vector v2 = rotate90(v1);
	Point center = cent + ((v2/abs(v2)) * b);
	Circle c1(center, 1.0);
	int cnt = 0;
	for(int k=0; k<n; k++) {
	  cnt += isPointInCircle(p[k], c1);
	}
	ans = max(ans, cnt);
	
	center = cent + ((v2/abs(v2)) * -b);
	Circle c2(center, 1.0);
	cnt = 0;
	for(int k=0; k<n; k++) {
	  cnt += isPointInCircle(p[k], c2);
	}
	ans = max(ans, cnt);
      }
    }
    cout << ans << endl;
  }
  
  return 0;
}

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

using namespace std;

int n;
double x[300];
double y[300];

double dist(double x1, double y1, double x2, double y2) {
    double dx = abs(x1 - x2);
    double dy = abs(y1 - y2);
    return sqrt(dx * dx + dy * dy);
}

int count(double cx, double cy, int i, int j) {
    int ans = 0;
    for (int k = 0; k < n; ++k) {
        if (k == i || k == j || dist(cx, cy, x[k], y[k]) < 1) ans++;
    }
    return ans;
}

void solve() {
    for (int i = 0; i < n; ++i) {
        cin >> x[i] >> y[i];
    }

    int ans = 0;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            if (i == j) continue;
            double d = dist(x[i], y[i], x[j], y[j]);
            if (d > 2.) continue;
            double dx = x[j] - x[i];
            double dy = y[j] - y[i];
            double r = sqrt(1 - (d * d / 4));
            //cout << r << endl;
            double cx = x[i] + dx / 2 - dy * r / d;
            double cy = y[i] + dy / 2 + dx * r / d;
            //cout << x[i] << ", " << y[i] << " - " << x[j] << ", " << y[j] << " : " << cx << ", "<< cy << endl;
            ans = max(ans, count(cx, cy, i, j));
        }
    }
    cout << ans << endl;
}

int main() {
    while (1) {
        cin >> n;
        if (n == 0) break;
        solve();
    }
}