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-10;
    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 * r + eps > (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 <bits/stdc++.h>
#define ALL(V) (V).begin(),(V).end()
using namespace std;
using ll = int64_t;
template <typename T> using V = vector<T>;

pair<long double, long double> lr_theta[300][300];
ll default_cnt[300];
pair<long double, long double> circle[300];

long double eps = 1e-5;
const long double inf = 5e15;
auto error = make_pair(-inf, -inf);

long double length(long double x, long double y) { return sqrt(x * x + y * y); }

ll solve(ll N) {
    for(ll i = 0; i < N; i++) {
        long double x, y;
        cin >> x >> y;
        circle[i] = make_pair(x, y);
    }

    sort(circle, circle + N);

    auto fix_theta = [](long double theta) { return theta < 0 ? theta + 2 * M_PI : theta; };

    ll ans = 1;
    for(ll i = 0; i < N; i++) {
        V<pair<long double, bool>> lr_lis;
        ll cur = 1;
        long double xi, yi, xj, yj;
        tie(xi, yi) = circle[i];
        for(ll j = 0; j < N; j++) {
            if(i == j) continue;
            tie(xj, yj) = circle[j];
            long double dx = xj - xi;
            long double dy = yj - yi;
            long double cen_len = length(dx, dy);
            if(2.0 < cen_len) continue;
            long double orx = -dy, ory = dx;
            long double cur_orlen = length(orx, ory);
            long double orlen = sqrt(1 - cen_len * cen_len / 4);
            (orx *= orlen) /= cur_orlen;
            (ory *= orlen) /= cur_orlen;
            long double theta1 = fix_theta(atan2(dy / 2 + ory, dx / 2 + orx));
            long double theta2 = fix_theta(atan2(dy / 2 - ory, dx / 2 - orx));
            if(theta2 < theta1) swap(theta1, theta2);
            if((theta1 + 2 * M_PI) - theta2 <= theta2 - theta1) {
                cur++;
                lr_lis.emplace_back(theta1, false);
                lr_lis.emplace_back(theta2, true);
            } else {
                lr_lis.emplace_back(theta1, true);
                lr_lis.emplace_back(theta2, false);
            }
        }
        ans = max(ans, cur);
        sort(ALL(lr_lis));
        for(auto &&ele : lr_lis) {
            cur += ele.second ? 1 : -1;
            ans = max(ans, cur);
        }
    }
    return ans;
}

int main() {
    while(true) {
        ll N;
        cin >> N;
        if(!N) break;
        cout << solve(N) << endl;
    }
    return 0;
}