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<vector>
#include<map>
using namespace std;
#define INF 1e+10
#define EPS 1e-10
#define EQ(a,b) (abs(a-b)<EPS)

//誤差を考慮して足し算
inline double add(double a,double b){
	if(abs(a+b) < EPS*(abs(a)+abs(b)))return 0;
	return a+b;
}
struct P{//2次元ベクトル
	double x,y;
	P(){}
	P(double x,double y):x(x),y(y){}
	P operator + (P p){
		return P(add(x,p.x),add(y,p.y));
	}
	P operator - (P p){
		return P(add(x,-p.x),add(y,-p.y));
	}
	P operator * (double d){
		return P(x*d,y*d);
	}
	bool operator ==(P p){
		return x==p.x && y==p.y;
	}
	bool equal(P p){
		return (x-p.x)*(x-p.x)+(y-p.y)*(y-p.y) < EPS*EPS;
	}	
	double dot(P p){//内積
		return add(x*p.x,y*p.y);
	}
	double det(P p){//外積
		return add(x*p.y,-y*p.x);
	}
};
//線分p1-p2上に点qがあるか判定
bool on_seg(P p1,P p2,P q){
	return (p1-q).det(p2-q)==0 && (p1-q).dot(p2-q)<=0;
}
//直線p1-p2と直線q1-q2の交点
P intersection(P p1,P p2,P q1,P q2){
	return p1+(p2-p1)*( (q2-q1).det(q1-p1)/(q2-q1).det(p2-p1) );
}
//線分p1-p2と線分q1-q2が交差するか
bool intersecting(P p1,P p2,P q1,P q2){
	if(on_seg(p1,p2,q1)||on_seg(p1,p2,q2)||on_seg(q1,q2,p1)||on_seg(q1,q2,p2))return true;
	if((p1-p2).det(q1-q2)==0 )return false;//平行
	P a=intersection(p1,p2,q1,q2);
	return on_seg(p1,p2,a) && on_seg(q1,q2,a);
}
P rotate(P p,double rad){
	double s=sin(rad),c=cos(rad);
	return P(c*p.x-s*p.y, s*p.x+c*p.y);
}
double dist(P p,P q){
	p = p-q;
	return sqrt(p.x*p.x+p.y*p.y);
}

//円同士の交点を計算
pair<P,P> circle_intersection(double x1,double y1,double r1,double x2,double y2,double r2){
	x2 -= x1;
	y2 -= y1;
	double rad = atan2(y2,x2);
	P tmp = rotate(P(x2,y2),-rad);
	x2 = tmp.x;
	y2 = tmp.y;
	double x = (r1*r1-r2*r2+x2*x2)/(2*x2);
	if(abs(r1) < abs(x))return make_pair(P(0,0),P(0,0));//交点なし
	double y = sqrt(r1*r1-x*x);
	return make_pair( P(x1,y1)+rotate(P(x,y),rad) ,  P(x1,y1)+rotate(P(x,-y),rad) );
}
//2点を通る円の中心座標
pair<P,P> circle_on_2P(P p,P q,double r){
	if(2*r < dist(p,q))return make_pair(P(0,0),P(0,0));//解なし
	if(q.x < p.x)swap(p,q);
	q = q - p;
	double rad = atan2(q.y,q.x);
	P tmp = rotate(q,-rad);
	double x = tmp.x/2;
	double y = sqrt(4*r*r-tmp.x*tmp.x)/2;
	return make_pair( p+rotate(P(x,y),rad) ,  p+rotate(P(x,-y),rad) );
}

int n;
P c[300];
int calc(P p,P q){
	if (2 < dist(p,q))return 1;
	pair<P,P> pa = circle_on_2P(p,q,1.0);
	int res=0,t=0;
	for(int i=0;i<n;i++){
		if(dist(pa.first,c[i])<=1)t++;
	}
	res = t;
	t=0;
	for(int i=0;i<n;i++){
		if(dist(pa.second,c[i])<=1)t++;
	}
	res = max(res,t);
	return res;
}
int main(){
	while(cin>>n&&n){
		for(int i=0;i<n;i++){
			double x,y;
			cin>>x>>y;
			c[i] = P(x,y);
		}
		int ans=1;
		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				ans = max(ans,calc(c[i],c[j]));
			}
		}
		cout<<ans<<endl;
	}
	return 0;
}

C++: #define _USE_MATH_DEFINES
#include <algorithm>
#include <cstdio>
#include <functional>
#include <iostream>
#include <cfloat>
#include <climits>
#include <cstring>
#include <cmath>
#include <map>
#include <queue>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <time.h>
#include <vector>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> i_i;
typedef pair<ll, int> ll_i;
typedef pair<double, int> d_i;
typedef pair<ll, ll> ll_ll;
typedef pair<double, double> d_d;
struct edge { int u, v; ll w; };

ll MOD = 1000000007;
ll _MOD = 1000000009;
double EPS = 1e-5;

int main() {
	for (;;) {
		int N; cin >> N;
		if (N == 0) break;
		vector<double> x(N), y(N);
		for (int i = 0; i < N; i++)
			cin >> x[i] >> y[i];
		int maxi = 1;
		for (int i = 0; i < N; i++)
			for (int j = 0; j < N; j++) {
				if (i == j) continue;
				double dx = x[j] - x[i], dy = y[j] - y[i];
				double d = sqrt(dx * dx + dy * dy);
				if (d > 2) continue;
				double h = sqrt(1 - d * d / 4);
				double ox = x[i] + dx / 2 - dy / d * h;
				double oy = y[i] + dy / 2 + dx / d * h;
				int cnt = 0;
				for (int k = 0; k < N; k++) {
					double dx = x[k] - ox, dy = y[k] - oy;
					double d2 = dx * dx + dy * dy;
					if (d2 <= 1 + EPS) cnt++;
				}
				maxi = max(maxi, cnt);
			}
		cout << maxi << endl;
	}
}

C++: //30
#include<iostream>
#include<complex>
#include<cmath>

using namespace std;

typedef complex<double> P;

int main(){
  for(int n;cin>>n,n;){
    P p[300];
    for(int i=0;i<n;i++){
      double x,y;
      cin>>x>>y;
      p[i]=P(x,y);
    }
    int ans=1;
    for(int i=0;i<n;i++){
      for(int j=0;j<n;j++){
	if(i==j||abs(p[i]-p[j])>2)continue;
	P c=p[i]+polar(1.,arg(p[j]-p[i])+acos(abs(p[j]-p[i])/2));
	int cc=0;
	for(int k=0;k<n;k++){
	  cc+=abs(p[k]-c)<1+1e-9;
	}
	if(ans<cc){
	  ans=cc;
	}
      }
    }
    cout<<ans<<endl;
  }
  return 0;
}

C++: #include <bits/stdc++.h>
#define REP(i,n) for(int i=0; i<(int)(n); ++i)

using namespace std;

const double EPS = 1e-8;

typedef complex<double> Point;

int calc(Point p, vector<Point>& ps){
    int ans = 0;
    for(int i = 0; i < ps.size(); i++){
        if(abs(ps[i] - p) < 1.0 + EPS) {
            ans++;
        }
    }
    return ans;
}
vector<Point> touching_circle(Point a, Point b, double r){
    vector<Point> res;

    double d = abs(b - a);
    if(d > 2 * r + EPS) return res;

    Point mid = 0.5 * (a + b);
    Point dir = polar(sqrt(r * r - d * d / 4), arg(b - a) + M_PI / 2);
    res.push_back(mid + dir);
    res.push_back(mid - dir);
    return res;
}
int main(){
    int N;
    while(cin >> N && N){
        vector<Point> ps(N);
        REP(i, N){
            double x, y;
            cin >> x >> y;
            ps[i] = Point(x, y);
        }
        int ans = 0;
        REP(i, N) REP(j, N) if(abs(ps[i] - ps[j]) > EPS) {
            vector<Point> tp = touching_circle(ps[i], ps[j], 1.0);
            for(Point c : tp){
            assert(abs(1.0 - abs(c - ps[i])) < EPS);
            assert(abs(1.0 - abs(c - ps[j])) < EPS);
                ans = max(ans, calc(c, ps));
            }
        }
        cout << ans << endl;
    }
    return 0;
}