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<bits/stdc++.h>
#define _ ios_base::sync_with_stdio(0);cin.tie(0);
#define REP(i,n) for(int i=0;i<(int)(n);i++)
#define REPP(i,j,n) for(int i=(j);i<(int)(n);i++)

using namespace std;
typedef long double ld;
typedef complex<ld> Point;
const ld eps = 1e-5, pi = acos(-1.0);
ld dot(Point a, Point b){return real(conj(a)*b);}
ld cross(Point a, Point b){ return imag(conj(a)*b);}

Point p[310];
int n;

void solve(){
  int result=1;
  REP(i,n){
    REPP(j,i+1,n){
      Point ab=p[j]-p[i];
      if(norm(ab)>4)continue;
      Point q=(p[i]+p[j])/2.L;
      Point aq=q-p[j];
      Point normal=Point{imag(aq),-real(aq)}*sqrt(1-norm(aq))/abs(aq);
      Point c=q+normal;
      Point d=q-normal;
      int cntc=0,cntd=0;
      REP(k,n){
	if(norm(p[k]-c)<=1.L)cntc++;
	if(norm(p[k]-d)<=1.L)cntd++;
      }
      result=max({result,2,cntc,cntd});
    }
  }
  cout<<result<<endl;
}

int main(){ _;
  ld r,i;
  while(cin>>n,n!=0){
    REP(j,n){
      cin>>r>>i;
      p[j]=Point{r,i};
    }
    solve();
  }
}

C++: #include <iostream>
#include <complex>
#include <cmath>
#define EPS 1e-9
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.0/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;j<n;j++)
      {
        /*if(abs(points[i]-points[j])<2.0)
        {
          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])));
        }*/
        double d=abs(points[j]-points[i])/2;
        xy_t v=(points[j]-points[i])/(2*d);
        if(d<1)
        {
            maxn=max(maxn,count_in(points[i]+d*v+sqrt(1-d*d)*v*xy_t(0,+1)));
            maxn=max(maxn,count_in(points[i]+d*v+sqrt(1-d*d)*v*xy_t(0,-1)));
        }
      }
    }
    cout<<maxn<<endl;
  }

}

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;

constexpr long double EPS = 10e-6;
constexpr long double radius = 1.0 + EPS;
constexpr long double diameter = 2.0 + EPS;

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

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, tmp;
        Point mid, midv, center;
        long double c;
        rep(n){
            repp(j, n){
                if(j == i) continue;
                mid = points[j] - points[i];
                if(mid.norm() > diameter*diameter) continue;
                mid = mid * 0.5;
                midv = Point(mid.y, -mid.x);
                c = sqrt((1.0 - mid.norm())/midv.norm());
                center = points[i] + mid + midv * c;
                tmp = 0;
                for(Point p: points){
                    Point v = p-center;
                    if(v.norm() <= radius) tmp++;
                }
                ans = max(tmp, ans);
            }
        }
        printf("%d\n", ans);
        debug(ans);
    }
}

C++: #include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <complex>
#include <string>
#include <sstream>
#include <algorithm>
#include <numeric>
#include <vector>
#include <queue>
#include <stack>
#include <functional>
#include <iostream>
#include <map>
#include <set>
using namespace std;
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vll;
#define pu push
#define pb push_back
#define mp make_pair
#define eps 1e-7
#define INF 1000000000
#define sz(x) ((int)(x).size())
#define sq(x) ((x)*(x))
#define fi first
#define sec second
#define SORT(x) sort((x).begin(),(x).end())
#define all(x) (x).begin(),(x).end()
#define rep(i,n) for(int (i)=0;(i)<(int)(n);(i)++)
#define repn(i,a,n) for(int (i)=(a);(i)<(int)(n);(i)++)
#define EQ(a,b) (abs((a)-(b))<eps)
double add(double a,double b){
	if(abs(a+b)<eps*(abs(a)+abs(b)))return 0;
	return a+b;
}
struct P{
	double x,y;
	P() {}
	P(double x,double y) : x(x),y(y){}
	P operator + (P p){
		return P(add(x,p.x),add(y,p.y));
	}
	P operator - (P p){
		return P(add(x,-p.x),add(y,-p.y));
	}
	P operator * (double d){
		return P(x*d,y*d);
	}
	double dot(P p){
		return add(x*p.x,y*p.y);
	}
	double det(P p){
		return add(x*p.y,-y*p.x);
	}
	double norm(){
		return sqrt(x*x+y*y);
	}
	double dist2(P p){
		return sq(x-p.x)+sq(y-p.y);
	}
	P vert(){
		return P(y,-x);
	}
	void dump(){
		printf("%.12f %.12f\n",x,y);
		return;
	}
};
int N;
P pts[305];
int solve(){
	scanf("%d",&N);
	if(N==0)return 1;
	for(int i=0;i<N;i++){
		scanf("%lf %lf",&pts[i].x,&pts[i].y);
	}
	int ans = 1;
	for(int i=0;i<N;i++){
		for(int j=i+1;j<N;j++){
			P m = (pts[i]+pts[j])*0.5;
			P a = (pts[i]-pts[j]);
			double na = a.norm();
			if(na>2.0)continue;
			a = (a.vert())*(sqrt(1-sq(na/2.0))/na);
			P b = m+a;
			P c = m-a;
			/*cout << i << ' ' << j << endl;
			m.dump();
			a.dump();
			b.dump();
			c.dump();
			cout << endl;*/
			int cnt = 0;
			for(int k=0;k<N;k++){
				//if(1.0<pts[k].dist2(b)&&pts[k].dist2(b)<=1.0+eps)printf("%d:%.12f\n",k,pts[k].dist2(b));
				if(pts[k].dist2(b)<=1.0+eps)cnt++;
			}
			ans = max(ans,cnt);
			cnt = 0;
			for(int k=0;k<N;k++){
				//if(1.0<pts[k].dist2(c)&&pts[k].dist2(c)<=1.0+eps)printf("%d:%.12f\n",k,pts[k].dist2(c));
				if(pts[k].dist2(c)<=1.0+eps)cnt++;
			}
			ans = max(ans,cnt);
		}
	}
	printf("%d\n",ans);
	return 0;
}
int main(){
	while(!solve()){}
	return 0;
}