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 <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;

const long double EPS = 10e-6;

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) {
        return Point(x+p.x, y+p.y);
    }
    Point operator - (Point p) {
        return Point(x-p.x, y-p.y);
    }
    Point operator * (long double a) {
        return Point(x*a, y*a);
    }
    Point operator / (long double a) {
        return Point(x/a, y/a);
    }
    
    long double norm() {
        return x*x + y*y;
    }
    
    long double dot(Point p) {
        return x*p.x + y*p.y;
    }
};

Point calc_center(Point p, Point mid, Point midv){
    long double norm = mid.norm();
    midv = midv * sqrt(1.0 - norm / norm);
    return p+mid+midv;
}

int count(vector<Point> points, Point center){
    int tmp = 0;
    for(Point p: points){
        if((p-center).norm() <= 1.0 + EPS) tmp++;
    }
    return tmp;
}

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;
        rep(n){
            repp(j, n){
                if(j == i) continue;
                Point mid = (points[j] - points[i]) * 0.5;
                Point center = calc_center(points[i], mid, Point(mid.y, -mid.x));
                ans = max(count(points, center), ans);

                center = calc_center(points[i], mid, Point(-mid.y, mid.x));
                ans = max(count(points, center), ans);
            }
        }
        printf("%d\n", ans);
    }
}

C++: #include <bits/stdc++.h>

using namespace std;

using ld = long double;
using Point = complex<ld>;

struct Line{
    Point a, b;
    Line (Point a, Point b) : a(a), b(b) {}
    Line () : Line(Point(), Point()) {}
};

struct Circle{
    Point p;
    ld r;
    Circle (Point p, ld r) : p(p), r(r) {}
    Circle () : Circle(Point(),0.0) {}
};

constexpr ld eps=1e-9, pi=acos(-1.0);

namespace std {
  bool operator<(const Point &lhs, const Point &rhs) {
    if (lhs.real() < rhs.real() - eps) return true;
    if (lhs.real() > rhs.real() + eps) return false;
    return lhs.imag() < rhs.imag();
  }
}

namespace Geometry{
    bool eq(ld a, ld b){
        return abs(a-b) < eps;
    }
    //内積
    ld dot(Point a, Point b){
        return real(conj(a) * b);
    }
    //外積
    ld cross(Point a, Point b){
        return imag(conj(a) * b);
    }
    //3点の位置関係
    int ccw (Point a, Point b, Point c) {
        b -= a; c -= a;
        if (cross(b, c) > eps) return 1;//a,b,cで反時計周り
        if (cross(b, c) < -eps) return -1;//a,b,cで時計周り
        if (dot(b, c) < 0) return 2;//c,a,bで直線
        if (norm(b) < norm(c)) return -2;//a,b,cで直線
        return 0;//a,c,bで直線
    }

//====================================================
    Point inputPoint(){
        ld x, y;
        cin>>x>>y;
        return Point(x, y);
    }
    //2直線の交差判定
    bool isCrossed_ll(Line l, Line m){
        return !eq(cross(l.b-l.a, m.b-m.a), 0);
    }
    //直線と線分の交差判定
    bool isCrossed_ls(Line l, Line s){
        return isCrossed_ll(l, s) && 
            cross(l.b - l.a, s.a - l.a) * cross(l.b - l.a, s.b - l.a) < eps;
    }
    //線分と線分の交差判定
    bool isCrossed_ss(Line s, Line t){
        return ccw(s.a, s.b, t.a) * ccw(s.a, s.b, t.b) <= 0 &&
            ccw(t.a, t.b, s.a) * ccw(t.a, t.b, s.b) <= 0;
    }
    //点が直線上にあるか
    bool isON_l(Line l, Point p){
        return abs(cross(l.b-p, l.a-p)) < eps;
    }
    //点が線分上にあるか
    bool isON_s(Line s, Point p){
        return abs(s.a-p)+abs(s.b-p)-abs(s.b-s.a) < eps;
    }

    //点から直線への垂線の足
    Point foot(Line l, Point p){
        ld t = dot(p-l.a,l.a-l.b)/norm(l.a-l.b);
        return l.a+t*(l.a-l.b);
    }
    //直線と直線の交点
    Point intersection_ll(Line l, Line m){
        Point lv=l.b-l.a, mv=m.b-m.a;
        assert(cross(lv,mv)!=0);//平行
        return l.a+lv*cross(mv,m.a-l.a)/cross(mv,lv);
    }
    //線分と線分の交点
    Point intersection_ss(Line s, Line t){
        assert(isCrossed_ll(s, t));
        return intersection_ll(s,t);
    }
    //点と点の距離
    ld dist_pp(Point p, Point q){
        ld x=p.real()-q.real(), y=p.imag()-q.imag();
        return sqrt(x*x+y*y);
    }
    //点と直線の距離
    ld dist_lp(Line l, Point p) {
        return abs(p - foot(l, p));
    }
    //直線と直線の距離
    ld dist_ll(Line l, Line m){
        return isCrossed_ll(l, m)?0:dist_lp(l, m.a);
    }
    //直線と線分の距離
    ld dist_ls(Line l, Line s){
        return isCrossed_ls(l, s)?0:min(dist_lp(l, s.a), dist_lp(l, s.b));
    }
    //線分と点の距離
    ld dist_sp (Line s, Point p) {
        Point r = foot(s, p);
        return isON_s(s, r) ? abs(r - p) : min(abs(s.a - p), abs(s.b - p));
    }
    //線分と線分の距離
    ld dist_ss (Line s, Line t) {
        if (isCrossed_ss(s, t)) return 0;
        return min({dist_sp(s, t.a), dist_sp(s, t.b), dist_sp(t, s.a), dist_sp(t, s.b)});
    }
    //円と円の交点
    vector<Point>intersection_cc(Circle c1, Circle c2){
        vector<Point> res;
        ld d = abs(c1.p - c2.p);
        ld rc = (d * d + c1.r * c1.r - c2.r * c2.r) / (2 * d);
        ld dfr = c1.r * c1.r - rc * rc;
        if (abs(dfr) < eps) dfr = 0.0;
        else if (dfr < 0.0) return res;
        ld rs = sqrt(dfr);
        Point diff = (c2.p - c1.p) / d;
        res.push_back(c1.p + diff * Point(rc, rs));
        if (dfr != 0.0) res.push_back(c1.p + diff * Point(rc, -rs));
        return res;
    }
};

int main(){
    int n;
    using namespace Geometry;
    while(cin>>n,n){
        vector<Point>g(n);
        vector<Circle>Cs(n);
        for(int i=0;i<n;++i){
            g[i]=inputPoint();
            Cs[i]={g[i],1.0};
        }
        int ans=0;
        for(int i=0;i<n;++i){
            for(int j=0;j<n;++j){
                vector<Point>intersections(intersection_cc(Cs[i],Cs[j]));
                for(auto intersection:intersections){
                    int ret=0;
                    for(int k=0;k<n;++k){
                        if(dist_pp(intersection,g[k]) < 1.0+eps){
                            ++ret;
                        }
                    }
                    ans=max(ans,ret);
                }
            }
        }
        cout<<ans<<endl;
    }
}