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 <cmath>
#define REP(i, a, n) for(int i = (a); i < (n); i++)
using namespace std;

int main(void) {
  for(int N; cin >> N, N > 0;) {
    double x[300], y[300];
    REP(i, 0, N) cin >> x[i] >> y[i];

    int answer = 1;
    REP(i, 0, N) {
      REP(j, i + 1, N) {
        double tan = (y[i] - y[j]) / (x[i] - x[j]);
        double cot = (x[i] - x[j]) / (y[i] - y[j]);
        double dsq = (x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j]);
        if(dsq > 4.0) continue;
        double l = sqrt(dsq - 1.);
        double sx = (x[i] + x[j]) / 2.;
        double sy = (y[i] + y[j]) / 2.;

        double cx1 = sx - tan * l;
        double cy1 = sy + cot * l;
        int cnt1 = 0;
        REP(k, 0, N) if(k == i || k == j || (x[k] - cx1) * (x[k] - cx1) + (y[k] - cy1) * (y[k] - cy1) <= 1.0) cnt1++;
        answer = max(answer, cnt1);

        double cx2 = sx + tan * l;
        double cy2 = sy - cot * l;
        int cnt2 = 0;
        REP(k, 0, N) if(k == i || k == j || (x[k] - cx2) * (x[k] - cx2) + (y[k] - cy2) * (y[k] - cy2) <= 1.0) cnt2++;
        answer = max(answer, cnt2);
      }
    }

    cout << answer << endl;
  }

  return 0;
}

C++: #include <vector>
#include <map>
#include <set>
#include <stack>
#include <queue>
#include <cstring>
#include <string>
#include <sstream>
#include <algorithm>
#include <iomanip>
#include <iostream>

#define VARIABLE(x) cerr << #x << "=" << x << endl
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define REP(i,m,n) for (int i=m;i<(int)(n);i++)

const int INF = 10000000;
int dx[4]={1, 0, -1, 0}, dy[4]={0, 1, 0, -1};
using namespace std;
typedef long long  ll;
typedef pair<double, double> P;

/** Problem1132 : Circle and Point **/
int main()
{
	int N, ans;
	double X, Y;
	vector<P> points;
	
	while (cin>>N, N) {
		ans=0;
		
		points.clear();
		rep(i, N) {
			cin >> X >> Y;
			points.push_back(P(X, Y));
		}

		for (double x=0.0; x<=10.0; x+=0.1) {
			for (double y=0.0; y<=10.0; y+=0.1) {
				int num=0;
				vector<P>::iterator it = points.begin();
				for (; it!=points.end(); it++) {
					double d = (it->first-x)*(it->first-x) + (it->second-y)*(it->second-y);
					if (d<1) {
						num++;
					}
				}
				if (ans<num)
					ans = num;
			}
		}
		cout << ans << endl;
	}
	return 0;
}

C++: #include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
//#include<cctype>
#include<climits>
#include<iostream>
#include<string>
#include<vector>
#include<map>
//#include<list>
#include<queue>
#include<deque>
#include<algorithm>
//#include<numeric>
#include<utility>
//#include<memory>
#include<functional>
#include<cassert>
#include<set>
#include<stack>
#include<random>

const int dx[] = {1, 0, -1, 0};
const int dy[] = {0, 1, 0, -1};
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef vector<int> vi;
typedef vector<ll> vll;
typedef pair<int, int> pii;
typedef double Real;

Real eps = 1e-9;

Real add(Real a, Real b) {
    if (abs(a+b) < eps * (abs(a)+abs(b))) return 0;
    return a+b;
}

bool equal(Real a, Real b) {
    return add(a, -b) == 0;
}

struct P {
    Real x, y;
    P() {}
    P(Real x, Real y) : x(x), y(y) {}
    P operator+(P p) const {return P(add(x, p.x), add(y, p.y));}
    P operator-(P p) const {return P(add(x, -p.x), add(y, -p.y));}
    P operator*(Real d) const {return P(x*d, y*d);}
    Real dot(P p) const {return add(x*p.x, y*p.y);} // ??????
    Real det(P p) const {return add(x*p.y, -y*p.x);} // ??????
    Real dist(P p) const {return sqrt((x-p.x)*(x-p.x) + (y-p.y)*(y-p.y));} // ?????¢
    void normalize() {Real d = sqrt(x*x+y*y); x /= d; y /= d;} // ??£??????
    bool operator<(const P& rhs) const {
        if (x != rhs.x) return x < rhs.x;
        return y < rhs.y;
    }
    bool operator==(const P& rhs) const {
        return equal(x, rhs.x) && equal(y, rhs.y);
    }
};
const int MAXN = 333;
P pnt[MAXN];
int N;

int calc(P center) {
	int ret = 0;
	for (int i = 0; i < N; i++) {
		P vec = pnt[i]-center;
		if (vec.dot(vec) < 1+eps) {
			ret++;
		}
	}
	return ret;
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    while (cin >> N) {
    	if (N==0) break;
    	for (int i = 0; i < N; i++) 
    		cin >> pnt[i].x >> pnt[i].y;
    	int ans = min(2, N);
    	for (int i = 0; i < N; i++) for (int j = i+1; j < N; j++) {
    		if (pnt[i].dist(pnt[j]) > 2) continue;
    		// diff ?????????
    		P vec = pnt[j]-pnt[i];
    		// vec ????????´????????????????????????
    		P n = P(vec.y, -vec.x);
    		n.normalize();
    		vec = vec*0.5;
    		double len = sqrt(1-vec.dot(vec));
    		P cand1 = pnt[i] + vec + n*len;
    		P cand2 = pnt[i] + vec - n*len;
    		ans = max(ans, calc(cand1));
    		ans = max(ans, calc(cand2));
    	}
    	cout << ans << endl;
    }
    return 0;
}