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>
using namespace std;

typedef complex<double> C;

int solve(int N, vector<double>& x, vector<double>& y) {
    vector<C> point(N);

    for (int j = 0; j < N; ++j) {
        point[j] = C(x[j], y[j]);
    }

    int ret = 0;

    for (int j = 0; j < N; ++j) {
        for (int k = 0; k < N; ++k) {
            C mid = (point[j] + point[k]) * 0.5;

            C normal = (point[j] - point[k]) * C(0, 1);
            normal /= abs(point[j] - point[k]);

            double a = 1 - norm(mid - point[j]);

            C centre1 = mid + a * normal;
            int count = 0;
            for (int l = 0; l < N; ++l) {
                if (abs(centre1 - point[l]) < 1) {
                    ++count;
                }
            }
            ret = max(ret, count);

            C centre2 = mid - a * normal;
            count = 0;
            for (int l = 0; l < N; ++l) {
                if (abs(centre2 - point[l]) < 1.0001) {
                    ++count;
                }
            }
            ret = max(ret, count);

        }
    }

    return ret;
}

int main () {
    int N;
    vector<double> x, y;

    while (true) {
        cin >> N;
        if (N == 0) break;

        x.resize(N);
        y.resize(N);

        for (int j = 0; j < N; ++j) {
            cin >> x[j] >> y[j];
        }

        cout << solve(N, x, y) << endl;
    }

    return 0;
}

C++: #include<iostream>
#include<vector>
#include<complex>
#include<algorithm>
#include<cmath>
#define dot complex<double>
using namespace std;
int main(){
  int n;
  while(cin >>n,n){
    vector<dot > num;
    for(int i=0; i<n; i++){
      double a,b;
      cin >>a>>b;
      num.push_back(dot(a,b));
    }
    int ans = 1;
    for(int i=0; i<n; i++){
      for(int j=i+1; j<n; j++){
	dot v = num[j]-num[i],o;
	if(abs(v) >= 2) continue;
	v/=2;
	int f = 0;
	double r = sqrt(1-abs(v)*abs(v))/abs(v);
	o = num[i]+v+dot(v.imag(),-v.real())*r;
	for(int k=0; k<n; k++){
	  dot vv = num[k]-o;
	  if(abs(vv)<1) f++;
	}
	ans = max(ans,f);
	f = 0;
	o = num[i]+v-dot(v.imag(),-v.real())*r;
	for(int k=0; k<n; k++){
	  dot vv = num[k]-o;
	  if(abs(vv)<=1.00000001) f++;
	}
	ans = max(ans,f);
      }
    }
    cout <<ans<<endl;
  }
  return 0;
}

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

using namespace std;

// ??????????????§?????????????????§???3????????????????????????????±?????????????????????¨?????????

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 reps(i,f,n) for(int i=f;i<int(n);i++)
#define rep(i,n) reps(i,0,n)

typedef complex<double> P;

vector<P> point;

const double EPS = 0.0000000001;


void init(){
  point.clear();
}


bool input(){
  int n;
  cin>>n;
  if(n==0)return false;

  rep(i,n){
    double a ,b;
    cin>>a>>b;
    point.push_back(P(a,b));
  }
  return true;
}

vector<P> make_cand(){
  vector<P> ret;
  rep(i,point.size()){
    reps(j,i+1,point.size()){
      P dist = point[j]-point[i];
      double d = abs(dist)/2;

      if(d>1)continue;
      double h = sqrt(1-d*d);
      P e = dist/(d*2);

      P e1 = e*P(0,1);
      P e2 = e*P(0,-1);

      P c1 = point[i]+e*d+h*e1;
      P c2 = point[i]+e*d+h*e2;

      //cout<<"d="<<d<<" h="<<h<<endl;
      //cout<<point[i]<<" "<<point[j]<<" "<<c1<<" "<<c2<<endl;
      ret.push_back(c1);
      ret.push_back(c2);
    }
  }
  rep(i,point.size()){
    ret.push_back(point[i]);
  }
  return ret;
}

int count_in(P p){
  int sum = 0;
  rep(i,point.size()){
    P d  = point[i]-p;
    if(abs(d)<=1+EPS)sum++;
  }
  return sum;
}
int solve(){
  vector<P> cand = make_cand();
  int ans = 0;
  rep(i,cand.size()){
    ans = max(ans, count_in(cand[i]));
  }
  return ans;
}

int main(){
  while(init(),input()){
    printf("%d\n",solve());
  }
}

C++: #include <bits/stdc++.h>

using namespace std;

const double r = 1.0;
const double PI = acos(-1.0);

struct Point {
    double x, y;
    Point() {}
    Point(double x, double y) : x(x), y(y) {}
};
inline double Distance2(const Point &a, const Point &b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

int main()
{
    cin.tie(0);
    ios::sync_with_stdio(false);

    int n;
    std::random_device rnd;
    std::mt19937 mt(rnd());
    std::uniform_real_distribution<> rand100(0.0, 10.0);

    while (cin >> n, n) {
        // Input
        vector<Point> p(n);
        for (int i = 0; i < n; ++i)
            cin >> p[i].x >> p[i].y;

        // Solve
        int ans = 1;
        auto start = std::chrono::high_resolution_clock::now();
        auto end = std::chrono::high_resolution_clock::now();
        auto take_time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

        while (true) {
            end = std::chrono::high_resolution_clock::now();
            take_time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
            if (280 < take_time.count())
                break;

            for (int i = 0; i < 10; ++i) {
                Point c(rand100(mt), rand100(mt));
                int num = 0;
                for (int j = 0; j < n; ++j) {
                    double d = Distance2(p[j], c);
                    if (d <= 1.0)
                        ++num;
                }
                ans = max(ans, num);
            }
        }

        // Output
        cout << ans << endl;
    }

    return 0;
}