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 <bits/stdc++.h>
struct UnionFind {
    int sz;
    std::vector<int> par;
    UnionFind(int sz) {
        this->sz = sz;
        par.clear();
        par.resize(sz, -1);
    }
    int root(int x) {
        return par[x] < 0 ? x : par[x] = root(par[x]);
    }
    void unite(int a, int b) {
        int roota = root(a);
        int rootb = root(b);
        if(roota == rootb) return;
        if(par[rootb] < par[roota]) std::swap(roota, rootb);
        par[roota] += par[rootb];
        par[rootb] = roota;
    }
    bool same(int a, int b) {
        return (root(a) == root(b));
    }
    int size(int a) {
        return -par[root(a)];
    }
};
using namespace std;
using ll = long long;
using ld = long double;
using P = pair<int, int>;
using vi = vector<int>;
using vvi = vector<vector<int>>;
using vll = vector<ll>;
using vvll = vector<vector<ll>>;
const int MOD = 1000000007;
const int INF = 1000000000;
const ll LINF = 1ll<<50;
template<typename T>
void printv(const vector<T>& s) {
  for(int i=0;i<(int)(s.size());++i) {
    cout << s[i];
    if(i == (int)(s.size())-1) cout << endl;
    else cout << " ";
  }
}
using ld = long double;
using Point = std::complex<ld>;
const ld eps = 1e-9, pi = acos(-1.0);
namespace std
{
bool operator<(const Point &lhs, const Point &rhs)
{
 if (lhs.real() < rhs.real() - eps)
  return true;
 if (lhs.real() > rhs.real() + eps)
  return false;
 return lhs.imag() < rhs.imag();
}
}
Point input_point()
{
 ld x, y;
 std::cin >> x >> y;
 return Point(x, y);
}
bool eq(ld a, ld b)
{
 return (abs(a - b) < eps);
}
ld dot(Point a, Point b)
{
 return real(conj(a) * b);
}
ld cross(Point a, Point b)
{
 return imag(conj(a) * b);
}
int ccw(Point a, Point b, Point c)
{
 b -= a;
 c -= a;
 if (cross(b, c) > eps)
  return 1;
 if (cross(b, c) < -eps)
  return -1;
 if (dot(b, c) < 0)
  return 2;
 if (norm(b) < norm(c))
  return -2;
 return 0;
}
int main() {
  cin.tie(0);
  cout << fixed << setprecision(10);
  while(1) {
    int n; cin >> n;
    if(n == 0) break;
    vector<Point> v(n);
    for(int i=0;i<n;++i) {
      double x, y; cin >> x >> y;
      v[i] = Point(x, y);
    }
    int ma = 1;
    for(int i=0;i<n;++i) {
      for(int j=i+1;j<n;++j) {
        Point mid = (v[i] + v[j]) / (ld)2.0;
        Point diff = (v[i] - v[j]) / (ld)2.0;
        ld l2 = abs(diff);
        if(l2 + eps > 1) continue;
        ld tang = 1.0 / (l2 * l2) - 1;
        Point d(-imag(diff), real(diff));
        Point test1 = mid + d * sqrt(tang);
        Point test2 = mid - d * sqrt(tang);
        int cnt1 = 0, cnt2 = 0;
        for(int k=0;k<n;++k) {
          if(abs(v[k] - test1) < 1 + eps) {
            cnt1++;
          }
          if(abs(v[k] - test2) < 1 + eps) {
            cnt2++;
          }
        }
        ma = max(ma, max(cnt1, cnt2));
      }
    }
    cout << ma << endl;
  }
}

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) 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 = 1;
        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){
                ans = max(ans, calc(c, ps));
            }
        }
        cout << ans << endl;
    }
    return 0;
}

C++: #include <iostream>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <cmath>

using namespace std;

double pow2(double d){
	return d * d;
}

struct Vector2{
	double x;
	double y;
};

Vector2 P[500];

int main(){
	int n;
	while(true){
		cin >> n;
		if(n == 0){
			break;
		}
		for(int i = 0; i < n; i++){
			cin >> P[i].x >> P[i].y;
		}
		double k, L_2, a, b;
		
		int max = 0;
		
		for(int num1 = 0; num1 < n; num1++){
			for(int num2 = num1 + 1; num2 < n; num2++){
				L_2 = pow2(P[num1].x - P[num2].x) + pow2(P[num1].y - P[num2].y);
				if(L_2 >= 4){
					continue;
				}
				k = sqrt((4 - L_2) / L_2);
				
				a = (k * (P[num1].y - P[num2].y) + P[num1].x + P[num2].x) / 2;
				b = (k * (P[num2].x - P[num1].x) + P[num1].y + P[num2].y) / 2;
				
				int count;
				
				count = 0;
				
				for(int i = 0; i < n; i++){
					if(i == num1 || i == num2){
						count++;
						continue;
					}
					if(pow2(P[i].x - a) + pow2(P[i].y - b) < 1){
						count++;
					}
				}
				
				if(max < count){
					max = count;
				}
				
				a = (k * (P[num2].y - P[num1].y) + P[num1].x + P[num2].x) / 2;
				b = (k * (P[num1].x - P[num2].x) + P[num1].y + P[num2].y) / 2;
				
				count = 0;
				
				for(int i = 0; i < n; i++){
					if(i == num1 || i == num2){
						count++;
						continue;
					}
					if(pow2(P[i].x - a) + pow2(P[i].y - b) < 1){
						count++;
					}
				}
				
				if(max < count){
					max = count;
				}
			}
		}
		printf("%d\n", max);
	}
	return 0;
}

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

#define rep(i, n) for (ll i = 0; i < (ll)(n); i++)
#define rept(i,k,n) for(int i = (k); i < (int)(n); i++)
#define Sort(v) sort(v.begin(),v.end())
#define Reverse(v) reverse(v.begin(),v.end())
#define X first
#define Y second

typedef long long ll;
typedef pair<double,double> Pdd;
typedef pair<int,int> Pii;

int n;
vector<vector<int> > tree;

int main(){
    while(1){
        int n;
        cin>>n;
        if(!n) return 0;
        vector<Pdd> xy(n);
        rep(i,n) cin>>xy[i].X>>xy[i].Y;
        int ans =1;
        rep(i,n) rep(j,n){
            if(i==j) continue;
            double dx=xy[i].X-xy[j].X, dy=xy[i].Y-xy[j].Y;
            double d = sqrt(dx*dx+dy*dy);
            if(d>2.000001)continue;
            double dr = sqrt(1.0-d*d/4);
            double xx=(xy[i].X+xy[j].X)/2+-dy/d*dr;
            double yy=(xy[i].Y+xy[j].Y)/2+dx/d*dr;
            int cnt=0;
            rep(k,n){
                if(pow(xy[k].X-xx,2.0)+pow((xy[k].Y-yy),2.0)<1.000001)cnt++;
            }
            ans = max(ans,cnt);
        }
        cout<<ans<<endl;
    }
    return 0;
}