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 <cmath>

using namespace std;

#define EPS (1e-6)


void calc(double x1,double y1,double x2,double y2,double *c_x1,double *c_y1,double *c_x2,double *c_y2){
  double midx,midy,hox,hoy,dx,dy,dis;
  midx = (x1+x2) / 2.0 , midy = (y1+y2) / 2.0;
  dx = x1-x2 , dy = y1-y2;
  dis = sqrt(1 - pow(x1-midx,2.0) - pow(y1-midy,2.0));
  hox = sqrt(pow(dy,2.0) / (pow(dy,2.0) + pow(dx,2.0)));
  hoy = hox * (-1.0) * (dx / dy);
  *c_x1 = midx + hox*dis,*c_y1 = midy + hoy*dis;
  *c_x2 = midx - hox*dis,*c_y2 = midy - hoy*dis;
}

int main(){
  int N;

  while(cin >> N){
    double x[300],y[300],c_x1,c_y1,c_x2,c_y2;
    int count,m;
    m = 0;
    if(!(N)) break;
    for(int i=0;i<N;i++) cin >> x[i] >> y[i];

    for(int i=0;i<N;i++){
      for(int j=0;j<N;j++){
        if(i == j) continue;
        if(sqrt((x[i] - x[j])*(x[i] - x[j]) + (y[i] - y[j])*(y[i] - y[j])) > 2.0 + EPS  ) continue;
        count = 2;
        calc(x[i],y[i],x[j],y[j],&c_x1,&c_y1,&c_x2,&c_y2);
        for(int k=0;k<N;k++){
          if(k == i || k == j) continue;
          if((x[k] - c_x1)*(x[k] - c_x1) + (y[k] - c_y1)*(y[k] - c_y1) < 1.0  ) count++;
        }
        if(m < count) m = count;
        count = 2;
        for(int k=0;k<N;k++){
          if(k == i || k == j) continue;
          if((x[k] - c_x2)*(x[k] - c_x2) + (y[k] - c_y2)*(y[k] - c_y2) < 1.0) count++;
        }
        if(m < count) m = count;
      }
    }

    printf("%d\n",m);
  }

  return 0;
}

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 = 0;
    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 "bits/stdc++.h"
using namespace std;
typedef long long LL;
struct Point {
	double x;
	double y;
	Point(double _x = 0, double _y = 0) :x(_x), y(_y) {}
	Point operator+(Point p) {
		return Point(x + p.x, y + p.y);
	}
	Point operator-(Point p)  {
		return Point(x - p.x, y - p.y);
	}
	double abs() {
		return sqrt(norm());
	}
	double norm() {
		return x*x + y*y;
	}
};
typedef Point Vector;
struct Circle {
	Point c;
	double r;
	Circle(Point _c = Point(), double _r = 0.0) :c(_c), r(_r) {}
};
double arg(Vector p) {
	return atan2(p.y, p.x);
}
Vector polar(double a, double r) {
	return Point(cos(r)*a, sin(r)*a);
}
pair<Point, Point> getCrossPoints(Circle c1, Circle c2) {
	double d = (c1.c - c2.c).abs();
	double a = acos((c1.r*c1.r + d*d - c2.r*c2.r) / (2 * c1.r*d));
	double t = arg(c2.c - c1.c);
	return make_pair(c1.c + polar(c1.r, t + a), c1.c + polar(c1.r, t - a));
}
double X[300], Y[300];
int main() {
	int N;
	while (cin >> N, N) {
		for (int i = 0; i < N; i++) {
			cin >> X[i] >> Y[i];
		}
		int ans = 0;
		for (int i = 0; i < N; i++) {
			for (int j = i + 1; j < N; j++) {
				Circle c1(Point(X[i], Y[i]), 1.0);
				Circle c2(Point(X[j], Y[j]), 1.0);
				if ((c1.c - c2.c).abs() > 2.0) continue;
				pair<Point, Point> c = getCrossPoints(c1, c2);
				int cnt1 = 0, cnt2 = 0;
				for (int k = 0; k < N; k++) {
					Point p(X[k], Y[k]);
					if ((p - c.first).abs() <= 1.0) cnt1++;
					if ((p - c.second).abs() <= 1.0) cnt2++;
				}
				ans = max(ans, cnt1);
				ans = max(ans, cnt2);
			}
		}
		cout << ans << endl;
	}
}

C++: #include <bits/stdc++.h>
#define rep(i,n) for(int i = 0; i < n; i++)
#define INF 100000000
#define EPS 1e-10
#define MOD 1000000007
using namespace std;
typedef pair<double,double> P;

int n;
int ans;
vector<double> x, y;
vector<P> v;

void solve(){
	x.clear(); y.clear();
	x.resize(n); y.resize(n);
	ans = 1;
	rep(i,n) cin >> x[i] >> y[i];
	for(int i = 0; i <= 1000; i++) for(int j = 0; j <= 1000; j++){
		double I = i/100.0, J = j/100.0;
		int cnt = 0;
		rep(k,n){
			if((x[k]-I)*(x[k]-I)+(y[k]-J)*(y[k]-J) < 1) cnt++;
		}
		ans = max(ans,cnt);
	}
	cout << ans << endl;
}

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

C++: #include <iostream>
#include <vector>
#include <string>
#include <complex>
#include <algorithm>
using namespace std;

typedef complex<double> p;

int n;
p ps[300];

int main() {
	while (cin >> n , n != 0) {
		for (int i = 0; i < n; i++) {
			double x, y;
			cin >> x >> y;
			ps[i] = p(x, y);
		}

		if (n == 1) {
			cout << 1 << endl;
			continue;
		}

		int ans = 0;
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				p v = ps[i] - ps[j], ve = v / abs(v);
				double l = sqrt(1 - norm(v) * 0.25);
				p vn1 = ve * p(0, 1) * l, vn2 = ve * p(0, -1) * l;
				p c1 = ps[j] + v * 0.5 + vn1, c2 = ps[j] + v * 0.5 + vn2;
				int count1 = 2, count2 = 2;
				for (int k = 0; k < n; k++) {
					if (k == i || k == j) continue;
					if (norm(c1 - ps[k]) <= 1.0) count1++;
					if (norm(c2 - ps[k]) <= 1.0) count2++;
				}
				ans = max(ans, count1);
				ans = max(ans, count2);
			}
		}

		cout << ans << endl;
	}
}