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 <algorithm>
#include <cassert>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <deque>
#include <functional>
#include <list>
#include <map>
#include <queue>
#include <set>
#include <sstream>
#include <string>
#include <utility>
#include <vector>

#define REP(i,s,n) for(int i=(int)(s);i<(int)(n);i++)
#define DEBUGP(val) cerr << #val << "=" << val << "\n"

using namespace std;
typedef long long int ll;
typedef vector<int> VI;
typedef vector<ll> VL;
typedef pair<int, int> PI;
typedef pair<double, double> PD;

const double EPS = 1e-11;
double sq(double x) { return x * x; }

void solve(double ax, double ay, double bx, double by, int &ma,
		 vector<double> &x, vector<double> &y) {
  double sqdist = sq(ax - bx) + sq(ay - by);
  double dist = sqrt(sqdist);
  double ndist = sqrt(1 - sqdist / 4);
  double ex = (bx - ax) / dist, ey = (by - ay) / dist;
  vector<PD> ans(2);
  int n = x.size();
  REP(i, 0, 2) {
    double coef = i == 0 ? 1 : -1;
    double cx = (ax + bx) / 2 + ey * coef * ndist;
    double cy = (ay + by) / 2 - ex * coef * ndist;
    int cnt = 0;
    REP(i, 0, n) {
      double sqdist = sq(x[i] - cx) + sq(y[i] - cy);
      if (sqdist <= 1 + EPS) {
	cnt += 1;
      }
    }
    ma = max(ma, cnt);
  }
}

int main(void) {
  int n;
  while(scanf("%d", &n) && n) {
    vector<double> x(n), y(n);
    REP(i, 0, n) scanf("%lf%lf", &x[i], &y[i]);
    int ma = 0;
    REP(i, 0, n) {
      REP(j, 0, i) {
        solve(x[i], y[i], x[j], y[j], ma, x, y);
      }
    }
    printf("%d\n", ma);
  }
}

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

#define rep(i, n) REP(i, 0, n)
#define REP(i, s, e) for(int i = (int)(s); i < (int)(e); ++i)

#define X() real()
#define Y() imag()
#define x(p) (p).X()
#define y(p) (p).Y()
#define SZ(P) (int)(P.size())
#define curr(P, i) P[(i)%SZ(P)]
#define next(P, i) P[(i+1)%SZ(P)]
#define prev(P, i) P[(i+SZ(P)-1)%SZ(P)]?´

using D = double;
using P = complex<D>;
namespace std {
  bool operator < (const P& a, const P& b) {
    return x(a) != x(b) ? x(a) < x(b) : y(a) < y(b);
  }
}
D dot(const P& a, const P& b) {
  return x(conj(a)*b);
}
D cross(const P& a, const P& b) {
  return y(conj(a)*b);
}
int ccw(P a, P b, P c) {
  b -= a;
  c -= a;
  if(cross(b, c) > 0) return +1;
  if(cross(b, c) < 0) return -1;
  if(dot(b, c) < 0) return +2;
  if(norm(b) < norm(c)) return -2;
  return 0;
}
const D EPS = 1e-8;
const D PI = acos(-1);

struct C {
  P p;
  D r;
};

P makeP(D arg, D r=1.0) {
  return r * P{cos(arg), sin(arg)};
}
P rotP(const P& p, D arg) {
  return p * makeP(arg);
} 

vector<P> getNorm(const P& p) {
  return {rotP(p, PI/2.0)/abs(p), rotP(p, -PI/2.0)/abs(p)};
}

vector<C> makeC(const P& a, const P& b, D r = 1.0) {
  P m = (a+b)/2.0;
  D d = sqrt(r*r-abs(m-a)*abs(m-a));
  auto ps = getNorm(b-a);
  rep(i, 2) ps[i] = d*ps[i]+m;
  vector<C> ret;
  rep(i, 2) ret.push_back(C{ps[i], r});
  return ret;
}
bool include(const C& c, const P& p) {
  return abs(c.p-p) < c.r+EPS;
}

int main() {
  int n;
  while(cin >> n && n) {
    vector<P> ps(n);
    rep(i, n) {
      D x, y;
      cin >> x >> y;
      ps[i] = P{x, y};
    }
    int ans = 1;
    rep(i, n) rep(j, i) {
      auto cs = makeC(ps[i], ps[j]);
      for(auto& c : cs) {
	int tmp = 0;	
	rep(k, n) if(include(c, ps[k])) tmp++;
	ans = max(ans, tmp);
      }
    }
    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;

const double EPS=1e-10;

#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+EPS);
}

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+EPS) 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));
        t*=-1.0;
        p = p2+v/2.0 + t*h;
        ans=max(ans,f(p));
    }

    cout<<ans<<endl;
}

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

C++: #include <cstdio>
#include <complex>
#include <iostream>
#include <cmath>
#include <vector>
using namespace std;

typedef complex<double> coor;

#define EPS 1e-5

double in_cir(coor o, coor x)
{
  return abs(x - o);
}

int search(vector<coor> P)
{
  int mx = 0;
  /*for (int i = 0; i < P.size(); i++) {
    printf("%lf, %lf\n", P[i].real(), P[i].imag());
    }*/
  for (int i = 0; i < P.size(); i++) {
    for (int j = i + 1; j < P.size(); j++) {
      for (int k = 0; k < 2; k++) {
	if (abs(P[i] - P[j]) > 2) {continue;}
	coor d = coor(P[i].imag() - P[j].imag(),
		      -(P[i].real() - P[j].real()));
	if (k == 1) {d = -d;}
	d /= abs(d);
	d *= sqrt(1 - norm(P[i] - P[j]) / 4);
	//printf("abs(d) = %lf\n", abs(d));
	d += (P[i] + P[j]) / 2.0;
	//printf("%lf, %lf\n", d.real(), d.imag());
	int ctr = 2;
	for (int l = 0; l < P.size(); l++) {
	  //printf("%lf\n", in_cir(d, P[l]));
	  if (l == i || l == j) {continue;}
	  if (in_cir(d, P[l]) < 1) {
	    ctr++;
	  }
	}
	//printf("_%d\n", ctr);
	mx = max(ctr, mx);
	
      }
    }
  }
  return mx;
}

int main()
{
  int N;
  while (scanf("%d ", &N) == 1) {
    if (N == 0) {
      break;
    }
    vector<coor> P(N);
    for (int i = 0; i < N; i++) {
      double x, y;
      scanf("%lf %lf ", &x, &y);
      P[i] = coor(x, y);
    }
    printf("%d\n", search(P));
  }
  return 0;
}