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<complex>
#include<vector>
using namespace std;
typedef complex<double> P;
int n;
vector<P> p;
int ret(){
	int cnt=0;
	for(int i=0;i<n;i++){
		for(int j=i+1;j<n;j++){
			if(abs(p[i]-p[j])>2.0)
				continue;
			P c=(p[i]+p[j])*0.5;
			double d=sqrt(1.0-abs(p[i]-c)*abs(p[i]-c));
			P l1=(p[i]-p[j])*P(0,1);
			P l2=(p[i]-p[j])*P(0,-1);
			l1=l1*d/abs(l1)+c;l2=l2*d/abs(l2)+c;
			int x=2,y=2;
			for(int k=0;k<n;k++){
				if(k==i||k==j)
					continue;
				if(abs(l1-p[k])<1.0+0.00000001){
					x++;
				}
				if(abs(l2-p[k])<1.0+0.00000001){
					y++;
				}
				
			}
			cnt=max(cnt,max(x,y));
		}
	}
	return cnt;
	
}
int main(){
	while(cin>>n,n){
		p.resize(n);
		for(int i=0;i<n;i++){
			double x,y;
			cin>>x>>y;
			p[i]=P(x,y);
		}
		cout<<ret()<<endl;
	}
	return 0;
}

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

const double EPS = 1e-8;
using P = complex<double>;
int n;
P p[300];

int main(){
  while (cin >> n, n){
    for (int i = 0; i < n; i++){
      double x, y;
      cin >> x >> y;
      p[i] = P(x, y);
    }
    vector<P> ps;
    for (int i = 0; i < n; i++){
      for (int j = 0; j < i; j++){
        if (abs(p[i] - p[j]) >= 2.0) continue;
        P v = p[j] - p[i];
        P u = v * P(0.0, 1.0);
        u *= sqrt(1.0 - abs(v) / 2.0) / abs(u);
        ps.push_back(p[i] + v / 2.0 + u);
        ps.push_back(p[i] + v / 2.0 - u);
      }
    }
    int res = 1;
    for (P c : ps){
      int s = 0;
      for (int i = 0; i < n; i++){
        s += abs(c - p[i]) <= 1.0 + EPS;
      }
      res = max(res, s);
    }
    cout << res << endl;
  }
}

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

#define rep(n) for(int i=0;i<n;i++)
#define repp(j, n) for(int j=0;j<n;j++)
#define reppp(i, m, n) for(int i=m;i<n;i++)
#define all(c) c.begin(), c.end()
#define rall(c) c.rbegin(), c.rend()
#define debug(x) cerr << #x << ": " << x << endl

using namespace std;

typedef long long ll;
typedef pair<ll, ll> Pll;
typedef pair<int, int> Pii;

constexpr long double EPS = 10e-6;
constexpr long double radius = 1.0 + EPS;
constexpr long double diameter = 2.0 + EPS;

class Point{
    public: 
    long double x,y;
    Point(long double x = 0.0, long double y = 0.0): x(x), y(y) {}
    Point operator + (Point p) const {
        return Point(x+p.x, y+p.y);
    }
    Point operator - (Point p) const {
        return Point(x-p.x, y-p.y);
    }
    Point operator * (long double a) const {
        return Point(x*a, y*a);
    }
    Point operator / (long double a) const {
        return Point(x/a, y/a);
    }
    long double norm() const {
        return x*x + y*y;
    }
    long double dot(Point p) const {
        return x*p.x + y*p.y;
    }
};

int main() {
    int n;
    while(scanf("%d", &n) && n){
        vector<Point> points(n);
        rep(n) scanf("%Lf %Lf", &points[i].x, &points[i].y);

        int ans = 1, tmp;
        Point mid, midv, center;
        long double c;
        rep(n){
            repp(j, n){
                if(j == i) continue;
                mid = points[j] - points[i];
                if(mid.norm() > diameter*diameter) continue;
                mid = mid * 0.5;
                midv = Point(mid.y, -mid.x);
                c = sqrt((1.0 - mid.norm())/midv.norm());
                center = points[i] + mid + midv * c;
                tmp = 0;
                for(Point p: points){
                    Point v = p-center;
                    if(v.norm() <= radius) tmp++;
                }
                ans = max(tmp, ans);
            }
        }
        printf("%d\n", ans);
    }
}

C++: #include <iostream>
#include <cstdio>
#include <cassert>
#include <cstring>
#include <vector>
#include <valarray>
#include <array>
#include <queue>
#include <set>
#include <unordered_set>
#include <map>
#include <unordered_map>
#include <algorithm>
#include <cmath>
#include <complex>
#include <random>

using namespace std;
typedef long long ll;
typedef unsigned long long ull;

typedef double R;
typedef complex<R> P;

const R EPS = 1e-10;
const R PI = acos((R)(-1));

R radNorP(R x) {
    return fmod(fmod(x, 2*PI)+2*PI, 2*PI);
}

const int MN = 330;

P p[MN];
typedef pair<R, int> Pi;

bool solve() {
    int n;
    scanf("%d", &n);
    if (!n) return false;

    for (int i = 0; i < n; i++) {
        R x, y;
        scanf("%lf %lf", &x, &y);
        p[i] = P(x, y);
    }

    int res = 1;
    vector<Pi> v;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            if (i == j) continue;
            if (2.0 < abs(p[j]-p[i])) continue;
            R th = acos(abs(p[j]-p[i])/2.0);
            R l = radNorP(arg(p[j]-p[i])-th-EPS);
            R r = radNorP(arg(p[j]-p[i])+th+EPS);
            if (r < l) {
                v.push_back(Pi(l-EPS, 1));
                v.push_back(Pi(2*PI+EPS, -1));
                v.push_back(Pi(0-EPS, 1));
                v.push_back(Pi(r+EPS, -1));
            } else {
                v.push_back(Pi(l-EPS, 1));
                v.push_back(Pi(r+EPS, -1));
            }

        }

        sort(v.begin(), v.end());
        int sm = 1;
        for (Pi p: v) {            
            sm += p.second;
            res = max(res, sm);
        }
        v.clear();
    }
    cout << res << endl;
    return true;
}

int main() {
    while (solve()) {}
    return 0;
}