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 <vector>
#include <cmath>
#include <algorithm>
#include <queue>
#include <stack>
#include <set>

using namespace std;


int main(){
	int n;
	double x1,y1,xc,yc,dis,vx,vy,ex,ey;
	while(1){
		cin >> n;
		if(n==0)break;
		vector<double> x(n);
		vector<double> y(n);
		for(int i=0;i<n;i++){
			cin >> x[i] >> y[i];
		}
		if(n==1){
			cout << 1 << endl;
			continue;
		}
		int ans =0;
		for(int i=0;i<n;i++){
			for(int j=0;j<n;j++){
				if(i==j)break;
				xc = (x[i]+x[j])/2;
				yc = (y[i]+y[j])/2;
				dis = sqrt(1.0-(xc-x[i])*(xc-x[i])-(yc-y[i])*(yc-y[i]));
				vx = x[i]-x[j];
				vy = y[i]-y[j];
				ex = dis * vy / sqrt(vx*vx+vy*vy);
				ey = -dis* vx / sqrt(vx*vx+vy*vy);
				x1 = xc + ex;
				y1 = yc + ey;
				int count=0;
				for(int k=0;k<n;k++){
					if(k==i||k==j)count +=1;
					else{
						if((x1-x[k])*(x1-x[k])+(y1-y[k])*(y1-y[k])<=1.0001){
							count +=1;
						}
					}
				}
				ans = max(ans,count);
				x1 = xc - ex;
				y1 = yc - ey;
				count = 0;
				for(int k=0;k<n;k++){
					if(k==i||k==j)count +=1;
					else{
						if((x1-x[k])*(x1-x[k])+(y1-y[k])*(y1-y[k])<=1.0001){
							count +=1;
						}
					}
				}
				ans = max(ans,count);
			}
		}
		cout << ans << endl;
	}
	return 0;
}

C++: #include <cstdio>
#include <cstring>
#include <string>
#include <iostream>
#include <cmath>
#include <vector>
#include <map>
#include <set>
#include <queue>
#include <deque>
#include <algorithm>
#define EPS 1e-10
 
using namespace std;
typedef long long int ll;
 
double add(double a, double b){
	if(abs(a+b)<EPS*(abs(a)+abs(b))) return 0;
	return a+b;
}
 
struct P{
	double x, y;
	P() {}
	P(double x, double y): x(x), y(y){
	}
	P operator+ (P p){
		return P(add(x, p.x), add(y, p.y));
	}
	P operator- (P p){
		return P(add(x, -p.x), add(y, -p.y));
	}
	P operator* (double d){
		return P(x*d, y*d);
	}
	double dot(P p){
		return add(x*p.x, y*p.y);
	}
	double det(P p){
		return add(x*p.y, -y*p.x);
	}
};

int main()
{
	while(1){
      int n;
      scanf("%d", &n);
      if(n==0) return 0;
      P p[300];
      for(int i=0; i<n; i++){
        double x, y;
        scanf("%lf %lf", &x, &y);
        p[i]=P(x, y);
      }
      
      int ans=1;
      for(int i=0; i<n; i++){
        for(int j=i+1; j<n; j++){
          double d=sqrt((p[i]-p[j]).dot(p[i]-p[j]));
          if(d>2.0) continue;
          P q;
          q.x=-(p[i]-p[j]).y;
          q.y=(p[i]-p[j]).x;
          double t=sqrt((1-((p[i]-p[j]).dot(p[i]-p[j]))/4.0)/(q.dot(q)));
          P r1=(p[i]+p[j])*0.5+q*t;
          P r2=(p[i]+p[j])*0.5-q*t;
          int c=2;
          for(int k=0; k<n; k++){
            if(k==i || k==j) continue;
            if(sqrt((p[k]-r1).dot(p[k]-r1))<1.0) c++;
          }
          if(ans<c) ans=c;
          c=2;
          for(int k=0; k<n; k++){
            if(k==i || k==j) continue;
            if(sqrt((p[k]-r2).dot(p[k]-r2))<1.0) c++;
          }
          if(ans<c) ans=c;
        }
      }
      printf("%d\n", ans);
    }
	return 0;
}

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

double EPS=1e-5;

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

		int res=1;
		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				if(abs(p[i]-p[j])>=2) continue;
				complex<double> mid=(p[i]+p[j])/complex<double>(2.0,0.0);
				auto a=(p[i]-p[j]);
				auto v=a*complex<double>(0,1);
				v=v/abs(v);
				double d=abs(mid-p[i]);
				double t=sqrt(1-d*d);
				complex<double> c[2]={mid+t*v,mid-t*v};
				for(int k=0;k<2;k++){
					int kai=0;
					for(int i=0;i<n;i++){
						if(abs(c[k]-p[i])<1+EPS) kai++;
					}
					res=max(res,kai);
				}
			}
		}
		cout<<res<<endl;
	}

}

C++: #include <bits/stdc++.h>
#define EPS 1e-6
using namespace std;
using p = complex<double>;

int n;
vector<p> v;
p rot[2] = {p(0, 1), p(0, -1)};

int solve();
int count(p now);

int main() {
  while(1) {
    cin >> n;
    if(n == 0) break;
    v.resize(n);
    for(int i = 0; i < n; ++i) {
      double x, y;
      cin >> x >> y;
      v[i] = p(x, y);
    }
    cout << solve() << endl;
  }
  return 0;
}

int solve() {
  int ans = 1;
  for(int i = 0; i < n; ++i)
    for(int j = i + 1; j < n; ++j)
      if(abs(v[i] - v[j]) <= 2.0 + EPS)
        for(int t = 0; t < 2; ++t) {
          p rotatep = (v[j] - v[i]) * rot[t];
          double len = 0;
          len = sqrt(4.0 - pow(abs(v[i] - v[j]), 2));
          rotatep *= len / abs(rotatep);
          rotatep += v[i];
          p nowp = (rotatep - v[j]);
          nowp /= 2;
          ans = max(ans, count(v[j] + nowp));
        }
  return ans;
}

int count(p now) {
  int ans = 0;
  for(int i = 0; i < n; ++i)
    if(abs(now - v[i]) <= 1.0 + EPS) ++ans;
  return ans;
}

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 (1000 < 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;
}