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

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 <sstream>
#include <vector>
#include <set>
#include <map>
#include <string>
#include <algorithm>
#include <cmath>
#include <climits>
#include <cstdio>
#include <cstdlib>
#include <complex>
using namespace std;

typedef long long LL;
typedef vector<int> vi;
typedef vector<string> vs;
typedef vector<vector<int> > vii;
typedef vector<vector<string> > vss;
typedef pair<int, int> pii;
#define INF 0x7fffffff
//#define INFLL 0x7fffffffffffffff
#define INFLL LLONG_MAX
#define SIZE(v) int((v).size())
#define ALL(v) (v).begin(), (v).end()
#define FOR(i,x) for(int i=0;i<x;i++)
#define FORic(i,y,x) for(int i=y;i<=x;i++)
#define FORdc(i,y,x) for(int i=y;i>=x;i--)
#define FILL(a,v) memset(a, v, sizeof(a))
#define PI 3.14159265358979323846

typedef complex<double> P; /* _ */
typedef pair<P, double> C; /* ~ */
struct L { P p, q; L(P p, P q) : p(p), q(q) {} }; /* üª */

const double EPS = 1e-13;

// complexÉå¬ÖWªKvÈê(intersCSÅgp)
bool operator < (const P &a, const P &b) {
	return make_pair(real(a),imag(a)) < make_pair(real(b),imag(b));
}


// ~Ìð_ðÔ·(æÉð·»èK{)
pair<P,P> crossCC(C c1, C c2){
	double r1 = c1.second;
	double r2 = c2.second;
	double r3 = abs(c1.first-c2.first);
	double rc = (r3*r3+r1*r1-r2*r2)/(2*r3);
	double rs = sqrt(r1*r1-rc*rc);
	P dif = (c2.first-c1.first)/r3;
	P p1 = c1.first+dif*P(rc,rs);
	P p2 = c1.first+dif*P(rc,-rs);
	return make_pair(p1,p2);
}

int main(){
	int n;
	while(cin>>n){
		if(!n) break;
		vector<P> all;
		FOR(i, n){
			double x, y;
			cin>>x>>y;
			all.push_back(P(x,y));
		}
		vector<P> couho;
		FOR(i, n-1){
			FORic(j, i+1, n-1){
				if(abs(all[j]-all[i])<=2+EPS){
					pair<P,P> temp = crossCC(C(all[i], 1), C(all[j], 1));
					couho.push_back(temp.first);
					couho.push_back(temp.second);
				}
			}
		}
		if(SIZE(couho)==0){
			cout<<1<<endl;
			continue;
		}
		int ans = 0;
		FOR(i, SIZE(couho)){
			int count = 0;
			FOR(j, n){
				if(abs(couho[i]-all[j])<=1+EPS){
					count++;
				}
			}
			ans = max(ans, count);
		}
		cout<<ans<<endl;
	}
	return 0;
}

C++: #include <iostream>
#include <string>
#include <cmath>
#include <set>
#include <utility>
#include <vector>
#include <map>
#include <queue>
#include <algorithm>

using namespace std;

constexpr long long mod = 1000000007;

constexpr double eps = 0.00001;

int n;
double p[300][2];

inline double dist2(double x1, double y1, double x2, double y2) {
    return (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2);
}

int cnt_c(double x, double y) {
    int ret = 0;
    for (int i = 0; i < n; i++) {
        if (dist2(x, y, p[i][0], p[i][1]) < (1.0 + eps) * (1.0 + eps)) ret++;
    }
    return ret;
}

int main() {
    while (1) {
        scanf("%d", &n);
        if (n == 0) return 0;
        int r = 1;
        for (int i = 0; i < n; i++) {
            scanf("%lf %lf", &p[i][0], &p[i][1]);
        }
        for (int i = 0; i < n; i++)
            for (int j = i + 1; j < n; j++) {
                if (dist2(p[i][0], p[i][1], p[j][0], p[j][1]) > (2 + eps) * (2 + eps)) continue;
                double ok, ng, dr;
                ok = dr = atan2(p[j][1] - p[i][1], p[j][0] - p[i][0]);
                ng = ok + M_PI;
                while (abs(ok - ng) > eps / 1000) {
                    double m = (ok + ng) / 2;
                    if (dist2(p[i][0] + cos(m), p[i][1] + sin(m), p[j][0], p[j][1]) <= 1)
                        ok = m;
                    else ng = m;
                }
                double is = (ok + ng) / 2;
                r = max(r, cnt_c(p[i][0] + cos(is), p[i][1] + sin(is)));
                r = max(r, cnt_c(p[i][0] + cos(dr * 2 - is), p[i][1] + sin(dr * 2 - is)));
            }
        cout << r << endl;
    }
}