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 <string>
#include <sstream>
#include <cmath>
#include <vector>

using namespace std;

typedef long double ld;

struct Point {
    ld x, y;
    Point(ld x, ld y) : x(x), y(y) {}
    string to_s() {
        ostringstream os;
        os << x << ',' << y;
        return os.str();
    }
    Point operator+(const Point& o) const {
        return Point(x+o.x, y+o.y);
    }
    Point operator/(const int d) const {
        return Point(x/d, y/d);
    }
};

struct Circle {
    ld x, y;
    ld r;
    Circle(ld x=0, ld y=0, ld r=1) : x(x), y(y), r(r) {}
    const static ld eps = 1e-10;
    bool isIntersect(const Circle& o) const {
        ld dx = o.x - x;
        ld dy = o.y - y;
        return (r + o.r) * (r + o.r) > (dx*dx + dy*dy);
    }
    pair<Point, Point> intersectionPoint(const Circle& o) const {
        ld dx = o.x - x;
        ld dy = o.y - y;
        ld d = sqrt(dx*dx + dy*dy);

        ld slope = dy / dx;
        ld theta = atan(slope);
        ld theta_p = acos((r*r - o.r*o.r + d*d) / (2*r*d));

        return make_pair(
                Point(x + r*cos(theta - theta_p), y + r*sin(theta - theta_p)),
                Point(x + r*cos(theta + theta_p), y + r*sin(theta + theta_p)));
    }
    bool contains(ld px, ld py) const {
        return r * r + eps > (px - x) * (px - x) + (py - y) * (py - y);
    }
    string to_s() {
        ostringstream os;
        os << x << ',' << y << ':' << r << endl;
        return os.str();
    }
};

int countContainingPoints(const Circle& c, const vector<Circle>& ps) {
    int cnt = 0;
    for (int i = 0; i < ps.size(); i++) 
        if (c.contains(ps[i].x, ps[i].y)) 
            cnt++;
    return cnt;
}

Point midPoint(pair<Point, Point>& pp) {
    return (pp.first + pp.second) / 2;
}

int N;
int main() {
    vector<Circle> cs;
    while (cin >> N, N) {
        cs.clear(); cs.resize(N);
        for (int i = 0; i < N; i++) cin >> cs[i].x >> cs[i].y;

        vector<Point> ips; // Intersection PointS
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                if (i == j) continue;
                if (!cs[i].isIntersect(cs[j])) continue;
                pair<Point, Point> ip = cs[i].intersectionPoint(cs[j]);
                ips.push_back(ip.first);
                ips.push_back(ip.second);
                //ips.push_back(midPoint(ip));
            }
        }

        int maxPoint = 0;
        for (int i = 0; i < ips.size(); i++) {
            maxPoint = max(maxPoint, countContainingPoints(Circle(ips[i].x, ips[i].y, 1), cs));
        }
        cout << maxPoint << endl;
    }
    return 0;
}

C++: #include <bits/stdc++.h>
#define N 301
using namespace std;
typedef complex <double> P;
P p[N];
int n;

int cnt(P c){
  int res=0;
  for(int i=0;i<n;i++)if(abs(c-p[i])<1) res++;
  return res;
}

P get_C(P a,P b){
  double x=abs(b-a);
  if(x>2||a==b) return a;
  return P(x/2,sqrt(1-x*x/4))*(b-a)/abs(b-a)+a;  
}


int main(){
  while(1){
  cin>>n;
  if(!n)break;
  double a,b;
  for(int i=0;i<n;i++)cin>>a>>b,p[i]=P(a,b);
  
  int ans=0;
  for(int i=0;i<n;i++)
    for(int j=0;j<n;j++)
      ans=max(ans,cnt(get_C(p[i],p[j])));
  cout <<ans<<endl;  
  }
  return 0;
}

C++: #include <algorithm>
#include <functional>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <string>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <vector>
#include <list>
#include <stack>
#include <queue>
#include <map>
#include <set>
#include <bitset>
#include <climits>
#include <complex>

#define all(c) (c).begin(), (c).end()
#define rep(i,n) for(int i=0;i<(n);i++)
#define pb(e) push_back(e)
#define mp(a, b) make_pair(a, b)
#define fr first
#define sc second

const int INF=100000000;
int dx[4]={1,0,-1,0};
int dy[4]={0,1,0,-1};
using namespace std;
typedef pair<int ,int > P;
typedef long long ll;

#define equals(a, b) (fabs((a)-(b))<EPS)
#define X real()
#define Y imag()

typedef complex<double> Point;
typedef complex<double> Vector;


struct Circle {
    Point c;
    double r;
    Circle(Point c=Point(), double r=0.0) :
        c(c),r(r){}
};

struct Segment {
    Point p1, p2;
};

bool include(Circle c,Point p) {
    return !(norm(c.c-p) > c.r*c.r);
}

int n;
double xs[301],ys[301];

Point getMidpoint(Point p1, Point p2, Point p3) {
    Point p=p1+p2+p3;
    p/=3;

    return p;
}

int f(Point p) {
    int ret=0;
    Circle c=Circle(p,1);
    rep(i,n) {
        if(include(c,Point{xs[i],ys[i]})) ret++;
    }

    return ret;
}

void solve() {
    rep(i,n) cin>>xs[i]>>ys[i];
    int ans=1;
    rep(i,n) rep(j,n) if(i<j) {
        Point p1=Point{xs[i],ys[i]};
        Point p2=Point{xs[j],ys[j]};
        Vector v=p1-p2;
        if(abs(v) > 2.0) continue;
        
        Vector t=(v/abs(v))*Vector(0,1);
        double h=sqrt(1-abs(v)*abs(v)/4.0);

        Point p = p2+v/2.0 + t*h;


        ans=max(ans,f(p));
        p*=-1.0;
        ans=max(ans,f(p));
    }

    cout<<ans<<endl;
}

int main() {
    while(cin>>n) {
        if(!n) break;
        solve();
    }
    return 0;
}

C++: #include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <map>
#include <utility>
#include <set>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>
#include <sstream>
#include <complex>
#include <stack>
#include <queue>
#include <cstring>
#include <sstream>
#include <cassert>
#include <list>
using namespace std;
static const double EPS = 1e-5;
typedef long long ll;
typedef pair<int,int> PI;
typedef vector<int> vi;
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#define rep(i,n) for(int i=0;i<(int)n;++i)
#define FOR(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();++i)
#define ALL(c) (c).begin(), (c).end()
#define mp(a,b) make_pair(a,b)
#define pb(a) push_back(a)
#define F first
#define S second


int dx[]={0,1,0,-1};
int dy[]={1,0,-1,0};


main(){
  int n;
  while(cin>>n,n){
    double x[n],y[n];

    rep(i,n)cin>>x[i]>>y[i];

    int ans=1;

    rep(i,n){
      rep(j,i){
	double sx,sy;
	double tx=x[i]-x[j],ty=y[i]-y[j];
	double len=sqrt(tx*tx+ty*ty);
	if(len>2)continue;
	len=sqrt(1-len*len/4);
	double theta=atan2(ty,tx)+M_PI/2;
	int tans=0;
	sx=(x[i]+x[j])/2+cos(theta)*len;
	sy=(y[i]+y[j])/2+sin(theta)*len;
	rep(k,n){
	  tx=sx-x[k];
	  ty=sy-y[k];
	  if(tx*tx+ty*ty<=1+EPS)++tans;
	}
	ans=max(ans,tans);
	tans=0;
	sx=(x[i]+x[j])/2-cos(theta)*len;
	sy=(y[i]+y[j])/2-sin(theta)*len;
	rep(k,n){
	  tx=sx-x[k];
	  ty=sy-y[k];
	  if(tx*tx+ty*ty<=1+EPS)++tans;
	}
	ans=max(ans,tans);
      }
    }
    cout<<ans<<endl;
  }
}