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=0.0, ld y=0.0) : 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);
    }
};

struct Circle {
    ld x, y;
    ld r;
    Circle(ld x=0.0, ld y=0.0, ld r=1.0) : x(x), y(y), r(r) {}
    const static ld eps = 1e-7;
    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 {
        ld dx = x - o.x;
        ld dy = y - o.y;
        ld a = (dx*dx + dy*dy + r*r - o.r*o.r) / 2.0;
        ld d = sqrt((dx*dx+dy*dy)*r*r-a*a);
        Point p, q;
        p.x = (a*dx + dy * d) / (dx*dx+dy*dy) + o.x;
        q.x = (a*dx - dy * d) / (dx*dx+dy*dy) + o.x;
        p.y = (a*dy - dx * d) / (dx*dx+dy*dy) + o.y;
        q.y = (a*dy + dx * d) / (dx*dx+dy*dy) + o.y;
        //cout << p.x << ',' << p.y << endl;
        //cout << q.x << ',' << q.y << endl;
        return make_pair(p, q);
    }
    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;
}

int N;
int main() {
    //cout << Circle(0, 0, 5).intersectionPoint(Circle(8, 0, 5)).first.to_s() << endl;
    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);
            }
        }

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

C++: #include <iostream>
#include <complex>
#include <cmath>
#define EPS 10e-6
using namespace std;
typedef complex<double> xy_t;
int n;
xy_t points[301];
int count_in(xy_t center)
{
  int ans=0;
  for(int i=0;i<n;i++)
  {
    if (abs(center-points[i])<1.0+EPS)
      ans++;
  }
  return ans;
}
xy_t ppdcl(xy_t a)
{
  return xy_t(a.imag(),-a.real());
}
xy_t delta_v(xy_t C1,xy_t C2)
{
  xy_t half=(C1-C2)*0.5;
  double l=abs(half);
  return ppdcl(half)*sqrt(1-l*l)*(1/l);
}
int main()
{
  while(1)
  {
    cin>>n;
    if (n==0) break;
    double x,y;
    for(int i=0;i<n;i++)
    {
      cin>>x>>y;
      points[i]=xy_t(x,y);
    }
    int maxn=-1;
    for(int i=0;i<n;i++)
    {
      for(int j=i+1;j<n;j++)
      {
        if(abs(points[i]-points[j])<2.0+EPS)
        {
          xy_t mid_p=(points[i]+points[j])*0.5;
          maxn=max(maxn,count_in(mid_p+delta_v(points[i],points[j])));
          maxn=max(maxn,count_in(mid_p-delta_v(points[i],points[j])));
        }
      }
    }
    cout<<maxn<<endl;
  }

}

C++: #include <iostream>
#include <queue>
#include <vector>
#include <string>
#include <map>
#include <algorithm>
#include <complex>

#define REP(i,n) for(int i=0;i<(int)(n);i++)
#define EQ(a,b) (abs((a)-(b))<eps)

using namespace std;

typedef long double ld;
typedef complex<ld> P;
typedef vector<P> VP;
const ld eps = 1e-8, pi = acos(-1.0);

ld dot (P a, P b) { return real(conj(a) * b); }
ld cross (P a, P b) { return imag(conj(a) * b); }

namespace std {
  bool operator<(const P &lhs, const P &rhs) {
    return lhs.real() == rhs.real() ? lhs.imag() < rhs.imag() : lhs.real() < rhs.real();
  }
}

struct sq {
  P p;
  ld sz;
  sq(P p_, ld sz_) : p(p_), sz(sz_) {}
};

bool operator<(const sq& l, const sq& r) { return l.p < r.p; }

pair<int, int> f(const VP& vp, sq s) {
  ld margin = s.sz / sqrt(2.0);
  P c = s.p + P(s.sz / 2.0, s.sz / 2.0);
  int mx = 0;
  int mn = 0;
  for(P p : vp) {
    if(abs(p-c) <= 1 + margin) {
      ++mx;
      if(abs(p-c) <= 1 - margin) {
        ++mn;
      }
    }
  }
  return make_pair(mx, mn);
}

int calc(const VP& vp) {
  std::priority_queue<pair<int, sq>> q;
  int mx_of_mn = 0;
  sq start(P(), 10.0);
  q.push(make_pair(f(vp,start).first, start));
  while (!q.empty()){
    auto p = q.top(); q.pop();
    if (mx_of_mn >= p.first) return mx_of_mn;
    REP(i,2)REP(j,2){
      P c = p.second.p + P(i * p.second.sz / 2, j * p.second.sz / 2);
      sq ns(c, p.second.sz / 2);
      auto rng = f(vp, ns);
      mx_of_mn = max(mx_of_mn, rng.second);
      q.push(make_pair(rng.first, ns));
    }
  }
  return mx_of_mn;
}

int main() {
  while(1){
    int n;
    cin >> n;
    if(!n) break;
    VP vp;
    REP(i,n){
      ld x,y;
      cin >> x >> y;
      vp.emplace_back(x,y);
    }
    cout << calc(vp) << endl;
  }
  return 0;
}