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++: #define _USE_MATH_DEFINES
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<climits>
#include<string>
#include<vector>
#include<list>
#include<map>
#include<set>
#include<cmath>
#include<queue>
#include<cstring>
#include<stack>
using namespace std;
#define EPS 1e-8
#define INF 1000000
 
struct Point{
    double x,y;
    Point(){}
    Point(double _x,double _y){
        x=_x; y=_y;
    }
    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 double d)const{
        return Point(x*d,y*d);
    }
    bool operator <(const Point &p)const{
        if(x==p.x) return y<p.y;
        return x<p.x;
    }
    double norm(){
        return x*x+y*y;
    }
	bool input(){
		if(cin>>x>>y) return true;
		return false;
	}
};
 
struct Line{
    Point a,b;
    Line(){}
    Line(Point _a,Point _b){
        a=_a; b=_b;
    }
	bool input(){
		if(a.input() && b.input()) return true;
		return false;
	}
};
 
struct Circle{
    Point c;
    double r;
    Circle(){}
    Circle(Point _c,double _r){
        c=_c; r=_r;
    }
};

typedef vector<Point> Polygon;

double dot(Point p,Point q){
    return p.x*q.x+p.y*q.y;
}
 
double cross(Point p,Point q){
    return p.x*q.y-q.x*p.y;
}
 
int ccw(Point a,Point b,Point c){ 
	Point v1 = b-a;
	Point v2 = c-a;
    if(cross(v1,v2)>EPS) return +1; 
    if(cross(v1,v2)<-EPS) return -1; 
	if(dot(v1,v2)<-EPS) return +2; 
	if(v1.norm()<v2.norm()) return -2; 
    return 0; 
}
 
bool intersect_ss(Line l,Line m){
	return ccw(l.a,l.b,m.a)*ccw(l.a,l.b,m.b)<=0 && ccw(m.a,m.b,l.a)*ccw(m.a,m.b,l.b)<=0;
}
 
bool crosspoint_ll(Line l,Line m,Point &p){
    Point a1,a2,b1,b2;
    a1 = l.b-l.a; a2 = m.b-m.a;
    b1 = m.a-l.a; b2 = l.a-m.b;
    double s1,s2;
    s1 = cross(a1,b1)/2; s2 = cross(a1,b2)/2;
    if(s1+s2<EPS) return false; 
    p = Point(m.a.x+a2.x*s1/(s1+s2),m.a.y+a2.y*s1/(s1+s2));
    return true;
}
 
int crosspoint_ss(Line l,Line m,Point &p){
    if(intersect_ss(l,m)==false) return 0; 
    if(crosspoint_ll(l,m,p)==true) return 1;
	return -1; 
}
 
int crosspoint_cc(Circle c1,Circle c2,Point &p1,Point &p2){
    double d,a,t;
    d = sqrt((c2.c-c1.c).norm());
	if(abs(c1.c.x-c2.c.x)<EPS && abs(c1.c.y-c2.c.y)<EPS && abs(c1.r-c2.r)<EPS)
		return -1; 
    if(d<abs(c1.r-c2.r) || c1.r+c2.r<d) return 0; 
    a = acos((c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d));
    t = atan2(c2.c.y-c1.c.y,c2.c.x-c1.c.x);
    p1 = Point(c1.c.x+c1.r*cos(t+a),c1.c.y+c1.r*sin(t+a));
    p2 = Point(c1.c.x+c1.r*cos(t-a),c1.c.y+c1.r*sin(t-a));
    if(abs(p1.x-p2.x)<EPS && abs(p1.y-p2.y)<EPS) return 1; 
    return 2; 
}

int contains(Polygon g,Point p){
	Line l = Line(p,Point(INF,p.y));
	int cnt = 0, n = g.size();
	for(int i=0;i<n;i++){
		Point a = g[i]-p;
		Point b = g[(i+1)%n]-p;
		if(ccw(g[i],g[(i+1)%n],p)==0) return 1; 
		if(a.y>b.y) swap(a,b);
		if(a.y<=EPS && EPS<b.y && cross(a,b)>EPS) cnt++;
	}
	if((cnt&1)==1) return 2; 
	return 0; 
}

Polygon andrewScan(Polygon s){
	if(s.size()<=2) return s;
	sort(s.begin(),s.end());
	Polygon g;
	for(int i=0;i<s.size();i++){
		for(int n=g.size(); n>=2 && ccw(g[n-2],g[n-1],s[i])!=-1; n--){
			g.pop_back();
		}
		g.push_back(s[i]);
	}
	int upper_n = g.size();
	for(int i=s.size()-2;i>=0;i--){
		for(int n=g.size(); n>upper_n && ccw(g[n-2],g[n-1],s[i])!=-1; n--){
			g.pop_back();
		}
		g.push_back(s[i]);
	}
	reverse(g.begin(),g.end());
	g.pop_back();
	return g;
}

int cntInsidePoint(Circle c,Polygon s){
	int cnt=0;
	for(int i=0;i<s.size();i++){
		if((s[i]-c.c).norm()<=1+EPS) cnt++;
	}
	return cnt;
}

int main(){
	int N;
	while(cin>>N,N){
		Polygon S;
		S.resize(N);
		for(int i=0;i<N;i++) S[i].input();
		int ans = 0;
		for(int i=0;i<N;i++){
			for(int j=i+1;j<N;j++){
				Point p1,p2;
				if(crosspoint_cc(Circle(S[i],1),Circle(S[j],1),p1,p2)>0){
					ans = max(ans,max(cntInsidePoint(Circle(p1,1),S),cntInsidePoint(Circle(p2,1),S)));
				}
			}
		}
		if(ans==0) puts("1");
		else printf("%d\n",ans);
	}
}

C++: // http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1132&lang=jp

#include <iostream>
#include <complex>
#include <cmath>
#include <iomanip>
using namespace std;

typedef complex<double> xy_t;

const xy_t rot[2] = { xy_t(0, 1), xy_t(0, -1) };
const double eps = 1e-6; // 最初から誤差の見積もりがあるのでそれに甘えた。

int main() {
  int N;
  while (scanf("%d", &N) != EOF && N) {
    if (N == 1) {
      printf("1\n");
      continue;
    }
    xy_t P[310];
    for (int i=0; i<N; i++) {
      double X, Y;
      scanf("%lf%lf", &X, &Y);
      P[i] = xy_t(X, Y);
    }
    int ans = 0;
    for (int i=0; i<N; i++) {
      for (int j=i+1; j<N; j++) {
	if ( abs(P[i] - P[j]) > 2 - eps ) continue;
	for (int k=0; k<2; k++) {
	  xy_t c = (P[i] + P[j])/2.0; // 中点
	  xy_t n = (P[i] - P[j]);
	  xy_t m = n/abs(n) * rot[k] * sqrt(1 - abs(n/2.0) * abs(n/2.0));
	  xy_t center = c + m;
	  int count = 0;
	  for (int k=0; k<N; k++) {
	    if (abs(center - P[k]) < 1 + eps) count++;
	    // cerr << fixed << setprecision(10) << abs(center - P[k]) << endl;
	  }
	  ans = max(ans, count);
	}
      }
    }
    printf("%d\n", ans);
  }
}

C++: #include <bits/stdc++.h>
using namespace std;

typedef complex<double> P;

const double EPS = 1e-8;
const int N = 300;
const double R = 1;

struct Circle{
  P c;
  double r;
  Circle(P c=P(0,0), double r=0):c(c),r(r){}
};

vector<P> getIntersectCC(Circle c1, Circle c2){
  vector<P> res;
  double r1 = c1.r, r2 = c2.r, d;
  P p1 = c1.c, p2 = c2.c;
  d = abs(p1-p2);

  if(d < EPS && abs(r1-r2) < EPS) return res; //????????¨?????????????????´???????????????????????????????????????
  if(r1 + r2< d - EPS || d + EPS < abs(r1-r2)) return res; //???????????????

  double a = (r1*r1 - r2*r2 + d*d) / (2*d);
  double h = sqrt(max(r1*r1 - a*a, 0.0));
  P tmp1 = p1 + a / d * (p2-p1);
  P tmp2 = h / d * (p2-p1);

  if(abs(tmp2) < EPS) res.push_back(tmp1);
  else {
    res.push_back(P(real(tmp1) - imag(tmp2), imag(tmp1) + real(tmp2)));
    res.push_back(P(real(tmp1) + imag(tmp2), imag(tmp1) - real(tmp2)));
  }
  return res;
}


int n;
vector<P> v;

int solve(){
  int ans = 1;

  for(int i=0;i<n;i++){
    for(int j=i+1;j<n;j++){
      vector<P> cc = getIntersectCC(Circle(v[i], R), Circle(v[j], R)); // v[i]??¨v[j]??????????????¢?????¨??????R??¨????????? (?????´?????????????????????????????????????????????)
      if(cc.size() == 0) continue;
      for(int k=0;k<cc.size();k++){
        int res = 0;
        for(int l=0;l<n;l++) if(abs(v[l]-cc[k]) <= R + EPS) res++;
        ans = max(ans, res);
      }
    }
  }

  return ans;
}

int main(){
  while(cin >> n && n){
    v.clear();
    for(int i=0;i<n;i++){
      double x, y;
      cin >> x >> y;
      v.push_back(P(x,y));
    }
    cout << solve() << endl;
  }
}