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 <algorithm>
#include <vector>
#include <utility>
#include <cstring>
#include <cmath>
#include <queue>
#include <set>
#include <map>
using namespace std;

#define rep(i,n) for(int i=0; i<n; i++)
typedef long long ll;
typedef pair<int,int> pii;

const double EPS = 1e-8;

int main() {
	int n;
	while(cin >> n, n) {
		double x[n], y[n];
		rep(i,n) {
			cin >> x[i] >> y[i];
		}
		int maxi = 0;
		rep(i, n) {
			for(int j=i+1; j<n; j++) {
				if( ((x[i]-x[j])*(x[i]-x[j])+(y[i]-y[j])*(y[i]-y[j])) > 4.0 ) continue;
				double mx = (x[i]+x[j]) / 2.0;
				double my = (y[i]+y[j]) / 2.0;
				double px0, py0, px1, py1;
				if( fabs(y[i]-y[j]) < EPS ) {
					px0 = px1 = mx;
					py0 = my + sqrt(1.0 - (x[i]-mx)*(x[i]-mx));
					py1 = my - sqrt(1.0 - (x[i]-mx)*(x[i]-mx));
				} else if( fabs(x[i]-x[j]) < EPS ) {
					py0 = py1 = my;
					px0 = mx + sqrt(1.0 - (y[i]-my)*(y[i]-my));
					px1 = mx - sqrt(1.0 - (y[i]-my)*(y[i]-my));
				} else {
					double m = - (x[i]-x[j]) / (y[i]-y[j]);
					double l = sqrt(1.0 - ((x[i]-mx)*(x[i]-mx)+(y[i]-my)*(y[i]-my)));
					px0 = mx + sqrt(l*l/(1.0+m*m));
					px1 = mx - sqrt(l*l/(1.0+m*m));
					py0 = my + m * (px0-mx);
					py1 = my + m * (px1-mx);
				}
				int cnt0 = 0, cnt1 = 0;
				rep(k, n) {
					double d0 = (px0-x[k])*(px0-x[k]) + (py0-y[k])*(py0-y[k]);
					if( d0 <= 1.0+EPS ) {
						cnt0++;
					}
				}
				rep(k, n) {
					double d1 = (px1-x[k])*(px1-x[k]) + (py1-y[k])*(py1-y[k]);
					if( d1 <= 1.0+EPS ) {
						cnt1++;
					}
				}
				maxi = max(maxi, cnt0);
				maxi = max(maxi, cnt1);
			}
		}
		cout << maxi << endl;
	}
	return 0;
}

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

using namespace std;

typedef complex<double> Point;

bool eqv(double a, double b) {
    return abs(a-b) < 1e-10;
}

bool inCircle(const Point &p, const Point &center, double radius) {
    return abs(p-center) < radius+1e-10;
}

int check(const Point &center, const vector<Point> &points) {
    const int N = points.size();
    int cnt = 0;
    for(int i = 0; i < N; ++i) {
        if(inCircle(points[i], center, 1)) ++cnt;
    }
    return cnt;
}

int test(int a, int b, const vector<Point> &points) {
    const Point v(points[a]-points[b]);
    const double hlen = abs(v)/2;

    if(hlen >= 1) return 0;
    const Point unit = v / abs(v);
    const double angle = acos(hlen);
    return max(check(points[b] + unit*polar(1.0, angle), points), check(points[b] + unit*polar(1.0, -angle), points));
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(0);

    while(true) {
        int N;
        cin >> N;
        if(!N) break;

        vector<Point> points(N);
        for(int i = 0; i < N; ++i) {
            cin >> points[i].real() >> points[i].imag();
        }
        int ans = 0;
        for(int i = 0; i < N; ++i) {
            for(int j = 0; j < i; ++j) {
                ans = max(ans, test(i, j, points));
            }
        }
        if(N == 1) ans = 1;
        cout << ans << endl;
    }
    return 0;
}

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

}