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++: #define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include <complex>
#include <cstdio>
#include <cstring>
#include <ctime>
#include <string>
#include <vector>
#include <algorithm>
#include <cmath>
#include <sstream>
#include <utility>
using namespace std;
typedef long long ll;
typedef vector<int> vi;
typedef vector<vi> vvi;
typedef pair<int,int> pii;
#define all(c) (c).begin(), (c).end()
#define loop(i,a,b) for(ll i=a; i<ll(b); i++)
#define rep(i,b) loop(i,0,b)
#define each(e,c) for(auto&e:c)
#define pb push_back
#define eb emplace_back
#define mp make_pair
#define mt make_tuple
#define lb lower_bound
#define ub upper_bound
#ifdef DEBUG
#define dump(...) (cerr<<#__VA_ARGS__<<" = "<<(DUMP(),__VA_ARGS__).str()<<" ["<<__LINE__<<"]"<<endl)
struct DUMP:ostringstream{template<class T>DUMP &operator,(const T&t){if(this->tellp())*this<<", ";*this<<t;return *this;}};
#else
#define dump(...)
#endif
template<class T> ostream& operator<<(ostream& os, vector<T> const& v){
    rep(i,v.size()) os << v[i] << (i+1==v.size()?"":" ");
    return os;
}

typedef double R;
typedef complex<R> P;
R const eps = 1e-8;
R const pi = acos(-1);

int main(){
    int n;
    while(cin>>n && n){
        vector<P> ps(n);
        rep(i,n){
            R x,y; cin >> x >> y;
            ps[i] = P(x,y);
        }
        int ans = 1;
        rep(i,n)rep(j,i){
            P & p = ps[i];
            P & q = ps[j];
            R d = abs(p-q)/2;
            if(d+eps > 1.0) continue;
            P m = (p+q)/R(2);
            R l = sqrt(1-d*d);
            P c = m + (p-m)/d*l*P(0,1);
            int cnt = 0;
            dump(abs(c-p));
            rep(k,n) if(abs(c-ps[k]) < 1.0+eps) cnt++;
            ans = max(ans,cnt);
        }
        cout << ans << endl;
    }
}

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

using namespace std;

int N;
pair<double,double> dots[300];

int main(){

	while(cin >> N && N != 0){
		for(int i = 0; i < N; i++){
			pair<double,double> p;
			cin >> dots[i].second >> dots[i].first;
		}
		int maxCnt = 1;
		// ~ðìé½ßÌQ_ðßé
		for(int i = 0; i < N; i++){
			for(int j = i+1; j < N; j++){
				pair<double,double> dot1,dot2;
				dot1 = dots[i];
				dot2 = dots[j];

				double v = sqrt((dot1.first - dot2.first)*(dot1.first - dot2.first)
					+(dot1.second - dot2.second)*(dot1.second - dot2.second));
				// ñ_Ô£ª2æè¨¨«¯êÎÂ­êÈ¢
				if(v > 2){
					continue;
				}

				double x = sqrt(1 - (v/2)*(v/2));

				// ñgÌPÊ@üxNg
				pair<double,double> hose[2];
				hose[0].first = -1.0/(dots[i].first - dots[j].first);
				hose[0].second = 1.0/(dots[i].second - dots[j].second);
				hose[1].first = 1.0/(dots[i].first - dots[j].first);
				hose[1].second = -1.0/(dots[i].second - dots[j].second);

				double dd1 = sqrt((hose[0].second)*(hose[0].second)+(hose[0].first)*(hose[0].first));
				double dd2 = sqrt((hose[1].second)*(hose[1].second)+(hose[1].first)*(hose[1].first));
				hose[0].first /= dd1; hose[0].second /= dd1;
				hose[1].first /= dd2; hose[1].second /= dd2;

				// x{·é
				hose[0].first *= x;hose[1].first *= x;
				hose[0].second *= x;hose[1].second *= x;

				// mð«·
				hose[0].first += (dot1.first + dot2.first)/2;
				hose[0].second += (dot1.second + dot2.second)/2;
				hose[1].first += (dot1.first + dot2.first)/2;
				hose[1].second += (dot1.second + dot2.second)/2;

				// hose[0]Æhose[1]Í~ÌSÀW
				// e_ÉÖµÄ~ÌÉ¶Ý·é©Ç¤©ðßé
				for(int k = 0; k < 2; k++){
					int cnt = 0;
					for(int l = 0; l < N; l++){
						double dist = sqrt((dots[l].first -hose[k].first)*(dots[l].first -hose[k].first)
							+ (dots[l].second -hose[k].second)*(dots[l].second -hose[k].second));
						if(dist <= 1){
							cnt++;
						}
					}
					maxCnt = max(maxCnt,cnt);
				}
			}
		}
		cout << maxCnt << endl;
	}

	return 0;
}

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 = 0;
        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 <cmath>
#include <iostream>
#include <algorithm>
using namespace std;

struct points{
    double x;
    double y;
};

points p[300];

double get_square(points a, points b){
    return pow(a.x - b.x, 2.0) + pow(a.y - b.y, 2.0);
}

double get_dist(points a, points b){
    return sqrt(get_square(a, b));
}

int find_max_cover(int n){
    int count, max_cover = 1;
    int sign[2] = {1, -1};
    double d, e;
    points m, q, v;
    for (int i = 0; i < n - 1; i++) {
        for (int j = i+1; j < n; j++) {
            d = get_dist(p[i], p[j]);
            if (d > 2.0){
                continue;
            }
            m.x = (p[i].x + p[j].x) / 2.0;
            m.y = (p[i].y + p[j].y) / 2.0;
            v.x = (p[i].x - p[j].x) / d;
            v.y = (p[i].y - p[j].y) / d;
            e = sqrt(1.0 - d * d / 4.0);
            for (int s = 0; s < 2; s++) {
                q.x = m.x + sign[s] * e * v.y;
                q.y = m.y - sign[s] * e * v.x;
                count = 0;
                for (int k = 0; k < n; k++) {
                    if (k == i || k == j || get_square(q, p[k]) <= 1.0) {
                        count++;
                    }
                }
                max_cover = max(max_cover, count);
            }
        }
    }
    return max_cover;
}

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;
        }
        cout << find_max_cover(n) << endl;
    }
    return 0;
}