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 <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for(int i=a;i<b;i++)

typedef complex<double> P;
#define PI acos(-1)
#define EPS 1e-10
#define INF 1e9

int N;
P p[300];

int solve()
{
	int ans = 1;

	rep(i, 0, N) rep(j, i + 1, N)
	{
		double dd = abs(p[i] - p[j]);

		if (2.0 < dd - EPS) continue;

		P vec = p[i] - p[j];
		P m = p[i] + p[j]; m /= 2.0;
		P grad = P(-imag(vec), real(vec));
		grad /= abs(grad);

		P c = m + grad * sqrt(1 - abs(vec) * abs(vec) / 4.0);

		int cnt = 0;
		rep(k, 0, N)
		{
			double d = abs(p[k] - c);
			if (d <= 1.0 + EPS) cnt++;
		}
		ans = max(ans, cnt);

		c = m - grad * sqrt(1 - abs(vec) * abs(vec) / 4.0);

		cnt = 0;
		rep(k, 0, N)
		{
			double d = abs(p[k] - c);
			if (d <= 1.0 + EPS) cnt++;
		}
		ans = max(ans, cnt);
	}

	return ans;
}

int main()
{
	cin.tie(0);
	ios::sync_with_stdio(false);

	while (1)
	{
		cin >> N;
		if (N == 0) return 0;

		rep(i, 0, N)
		{
			double x, y; cin >> x >> y;
			p[i] = P(x, y);
		}

		cout << solve() << endl;
	}
}

C++: #include <cstdio>
#include <cstring>
#include <string>
#include <cmath>
#include <cassert>
#include <iostream>
#include <algorithm>
#include <stack>
#include <queue>
#include <vector>
#include <set>
#include <map>
#include <bitset>
#include <functional>
#include <numeric>
using namespace std;

#define repl(i,a,b) for(int i=(int)(a);i<(int)(b);i++)
#define rep(i,n) repl(i,0,n)
#define mp(a,b) make_pair((a),(b))
#define pb(a) push_back((a))
#define all(x) (x).begin(),(x).end()
#define dbg(x) cout<<#x"="<<((x))<<endl
#define fi first
#define se second

#define INF 2147483600
#define long long long
#define EPS 1e-8

// ???????°???°???????????????
inline double add(double a, double b){
  if(abs(a+b) < EPS*(abs(a) + abs(b))) return 0;
  return a+b;
}

struct Point{
  double x,y;
  Point() {}
  Point(double nx, double ny) : x(nx), y(ny) {}
  inline Point operator + (const Point & p){ return Point(add(x, p.x), add(y, p.y)); }
  inline Point operator - (const Point & p){ return Point(add(x,-p.x), add(y,-p.y)); }
  inline Point operator * (double d){ return Point(x*d, y*d); }
  inline double dot(const Point & p){ return add(x * p.x, y*p.y); }  //??????
  inline double det(const Point & p){ return add(x * p.y, -y*p.x); } //??????
  inline double dist(const Point & p){ return sqrt((x-p.x)*(x-p.x) + (y-p.y)*(y-p.y)); }
  inline bool operator < (const Point & p) const {
    if(x != p.x) return x < p.x;
    else return y < p.y;
  }
  inline bool operator == (const Point & p) const {
    return (add(x, -p.x)==0) && (add(y, -p.y)==0);
  }
  friend ostream& operator<<(ostream& os, const Point& p) {
    os << "[" << p.x << "," << p.y << "]";
    return os;
  }
};

pair<Point, Point> calcCenter(Point &p1, Point &p2){
  if(add(p1.x, -p2.x)==0){
    double y = (p1.y + p2.y)/2;
    double dx = sqrt(1 - (p1.y - y)*(p1.y - y));
    return mp(Point(p1.x+dx, y), Point(p1.x-dx,y));
  }
  // 2???????????´???????????? Ax + By + C = 0 , A != 0
  double A = 2*(p1.x - p2.x);
  double B = 2*(p1.y - p2.y);
  double C = - p1.x*p1.x + p2.x*p2.x - p1.y*p1.y + p2.y*p2.y;
  // x = py+q
  double p = -B/A; double q = -C/A;

  // p1, p2 ??¨????????¢???1?????? ay^2 + 2by + c = 0
  double a = 1 + p*p;
  double b = p*q - p1.x*p - p1.y;
  double c = q*q - 2*p1.x*q + p1.x*p1.x + p1.y*p1.y - 1;

  double rt = sqrt(b*b-a*c);
  double y1 = (-b+rt)/a;
  double y2 = (-b-rt)/a;

  return mp(Point(p*y1+q, y1), Point(p*y2+q, y2));
}

Point vec[300];

int main(){
  int n;
  while(scanf("%d", &n), n){
    rep(i,n){
      scanf("%lf %lf", &vec[i].x, &vec[i].y);
    }

    int res=1;
    rep(i,n) repl(j,i+1,n) if(vec[i].dist(vec[j]) < 2){
      pair<Point, Point> pr = calcCenter(vec[i], vec[j]);
      int cnt=0;
      rep(k,n) if(pr.fi.dist(vec[k]) < 1) cnt++;
      res = max(res, cnt);
      cnt=0;
      rep(k,n) if(pr.se.dist(vec[k]) < 1) cnt++;
      res = max(res, cnt);
    }

    printf("%d\n", res);
  }
  return 0;
}

C++: #include <cmath>
#include <cstdio>
#include <algorithm>

#define MAX_N 300
const double PI=acos(-1.0);

struct Point{
    double x,y;
    Point(){
        x=y=0;
    }
    bool operator<(const Point& p)const{
        return this->x < p.x;
    }
};

struct comp{
    bool operator()(const Point &P1, const Point &P2){
        return P1.x < P2.x;
    }
    bool operator()(const Point &P, double d){
        return P.x < d;
    }
    bool operator()(double d, const Point &P){
        return d < P.x;
    }
};

Point points[MAX_N];
int N;

double distance(Point *a, Point *b){
    return std::sqrt(pow(a->x-b->x,2)+pow(a->y-b->y,2));
}

double distance2(Point *a, Point *b){
    return pow(a->x-b->x,2)+pow(a->y-b->y,2);
}

int solve(){
    Point C, *P1, *P2, *P3;
    int max=1;
    for(int i=0;i<N;i++){
        P1=&points[i];
        Point *LB,*UB;
        LB=std::lower_bound(points,points+N,P1->x-2.0001,comp());
        UB=std::upper_bound(points,points+N,P1->x+2.0001,comp());
        for(int j=(LB-points);j<(UB-points);j++){
            if(i==j) continue;
            P2=&points[j];
            if(distance(P1,P2)>2.0) continue;

            double alpha,beta;
            alpha=acos(distance(P1,P2)/2.0);
            if((P2->x-P1->x)>=0){
                if((P2->y - P1->y)>=0){
                    beta=atan((P2->y - P1->y) / (P2->x - P1->x));
                }else{
                    beta=atan((P2->y - P1->y) / (P2->x - P1->x))+2.0*PI;
                }
            }else{
                beta=atan((P2->y - P1->y) / (P2->x - P1->x))+PI;
            }
            C.x=P1->x+cos(alpha+beta);
            C.y=P1->y+sin(alpha+beta);

            int enclosed=0;
            Point *lb,*ub;
            lb=std::lower_bound(points,points+N,C.x-1.0001,comp());
            ub=std::upper_bound(points,points+N,C.x+1.0001,comp());
            for(int k=(lb-points);k<(ub-points);k++){
                if(k==i||k==j){
                    enclosed++;
                    continue;
                }
                P3=&points[k];
                if(distance2(P3,&C)<1.0) enclosed++;
            }
            if(enclosed>max) max=enclosed;
            if(max==N) return max;
        }
    }
    return max;
}

int main(void){
    while(scanf("%d",&N), N){
        for(int i=0;i<N;i++){
            scanf("%lf %lf", &points[i].x, &points[i].y);
        }
        std::sort(points,points+N);
        printf("%d\n",solve());
    }
    return 0;
}