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 <string>
#include <sstream>
#include <cmath>
#include <vector>

using namespace std;

typedef long double ld;

struct Point {
    ld x, y;
    Point(ld x, ld y) : x(x), y(y) {}
    string to_s() {
        ostringstream os;
        os << x << ',' << y;
        return os.str();
    }
    Point operator+(const Point& o) const {
        return Point(x+o.x, y+o.y);
    }
    Point operator/(const int d) const {
        return Point(x/d, y/d);
    }
};

struct Circle {
    ld x, y;
    ld r;
    Circle(ld x=0, ld y=0, ld r=1) : x(x), y(y), r(r) {}
    const static ld eps = 1e-5;
    bool isIntersect(const Circle& o) const {
        ld dx = o.x - x;
        ld dy = o.y - y;
        return (r + o.r) * (r + o.r) > (dx*dx + dy*dy);
    }
    pair<Point, Point> intersectionPoint(const Circle& o) const {
        //cout << x << ',' << y << " and " << o.x << ',' << o.y << endl;
        ld dx = o.x - x;
        ld dy = o.y - y;
        //cout << "dx: " << dx << " dy: " << dy << endl;
        ld d = sqrt(dx*dx + dy*dy);
        //cout << "d: " << d << endl;

        ld slope;
        ld theta;
        if (dx < eps) {
            if (y < o.y) {
                theta = M_PI / 2;
            } else {
                theta = - M_PI / 2;
            }
        } else {
            slope = dy / dx;
            theta = atan(slope);
        }
        //cout << "theta: " << theta << endl;
        ld theta_p = acos((r*r - o.r*o.r + d*d) / (2*r*d));
        //cout << "theta_p: " << theta_p << endl;

        return make_pair(
                Point(x + r*cos(theta - theta_p), y + r*sin(theta - theta_p)),
                Point(x + r*cos(theta + theta_p), y + r*sin(theta + theta_p)));
    }
    bool contains(ld px, ld py) const {
        return r + eps >= sqrt((px - x) * (px - x) + (py - y) * (py - y));
    }
    string to_s() {
        ostringstream os;
        os << x << ',' << y << ':' << r << endl;
        return os.str();
    }
};

int countContainingPoints(const Circle& c, const vector<Circle>& ps) {
    int cnt = 0;
    for (int i = 0; i < ps.size(); i++) 
        if (c.contains(ps[i].x, ps[i].y)) 
            cnt++;
    return cnt;
}

Point midPoint(pair<Point, Point>& pp) {
    return (pp.first + pp.second) / 2;
}

int N;
int main() {
    vector<Circle> cs;
    while (cin >> N, N) {
        cs.clear(); cs.resize(N);
        for (int i = 0; i < N; i++) cin >> cs[i].x >> cs[i].y;

        vector<Point> ips; // Intersection PointS
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                if (i == j) continue;
                if (!cs[i].isIntersect(cs[j])) continue;
                pair<Point, Point> ip = cs[i].intersectionPoint(cs[j]);
                ips.push_back(ip.first);
                ips.push_back(ip.second);
                //ips.push_back(midPoint(ip));
            }
        }

        int maxPoint = 0;
        for (int i = 0; i < ips.size(); i++) {
            maxPoint = max(maxPoint, countContainingPoints(Circle(ips[i].x, ips[i].y, 1), cs));
        }
        cout << maxPoint << endl;
    }
    return 0;
}

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

using namespace std;

int N;
pair<double,double> dots[300];
#define EPS (1e-6)


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.0 + EPS){
					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.0 + EPS){
							cnt++;
						}
					}
					maxCnt = max(maxCnt,cnt);
				}
			}
		}
		cout << maxCnt << endl;
	}

	return 0;
}