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 <cstdio>
#include <cmath>
#include <iostream>
#include <vector>
using namespace std;

struct Point {
    double x, y;
    Point(double x_ = 0.0, double y_ = 0.0) : x(x_), y(y_) {}
    Point(const Point& p) : x(p.x), y(p.y) {}

    Point& operator=(const Point& p) {
        x = p.x;
        y = p.y;
        return *this;
    }

    double det(const Point& p) {
        return x * p.y - y * p.x;
    }

    Point nrm() {
        Point ret;
        double n = length();
        ret.x = x / n;
        ret.y = y / n;
        return ret;
    }

    double length() {
        return sqrt(x * x + y * y);
    }
};

Point operator+(const Point& p1, const Point& p2) {
    return Point(p1.x + p2.x, p1.y + p2.y);
}

Point operator-(const Point& p1, const Point& p2) {
    return Point(p1.x - p2.x, p1.y - p2.y);
}

Point operator*(const Point& p1, double s) {
    return Point(p1.x * s, p1.y * s);
}

Point operator/(const Point& p1, double s) {
    return Point(p1.x / s, p1.y / s);
}

Point calcCenter(const Point& p1, const Point& p2, const Point& p3) {
    double a = (p1 - p2).length();
    double b = (p2 - p3).length();
    double c = (p3 - p1).length();
    double A = a * a;
    double B = b * b;
    double C = c * c;
    double S = 0.5 * std::abs((p2 - p1).det(p3 - p1));
    Point cent = (p1 * (A * (B + C - A)) + p2 * (B * (C + A - B)) + p3 * (C * (A + B - C))) / (16.0 * S * S);
    return p1 + (cent - p1).nrm();
}

void solve(const std::vector<Point>& pts) {
    int n = (int)pts.size();
    int maxval = 1;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            for (int k = j + 1; k < n; k++) {
                Point cnt = calcCenter(pts[i], pts[j], pts[k]);
                int s = 0;
                for (int l = 0; l < n; l++) {
                    if ((cnt - pts[l]).length() <= 1.0) {
                        s++;
                    }
                }

                if (maxval < s) {
                    maxval = s;
                }
            }
        }
    }
    printf("%d\n", maxval);
}

int main(int argc, char** argv) {
    int n;
    double x, y;
    while (cin >> n) {
        if (n == 0) break;

        vector<Point> pts;
        for (int i = 0; i < n; i++) {
            cin >> x >> y;
            pts.push_back(Point(x, y));
        }
        solve(pts);
    }
}

C++: #include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vl;
typedef complex<double> P;
typedef pair<int,int> pii;
#define REP(i,n) for(ll i=0;i<n;++i)
#define REPR(i,n) for(ll i=1;i<n;++i)
#define FOR(i,a,b) for(ll i=a;i<b;++i)

#define DEBUG(x) cout<<#x<<": "<<x<<endl
#define DEBUG_VEC(v) cout<<#v<<":";REP(i,v.size())cout<<" "<<v[i];cout<<endl
#define ALL(a) (a).begin(),(a).end()

#define MOD (ll)(1e9+7)
#define ADD(a,b) a=((a)+(b))%MOD
#define FIX(a) ((a)%MOD+MOD)%MOD

int n;
P pts[353];

int check(P c){
  int ret = 0;
  REP(i,n){
    if(abs(c-pts[i])<=1.00001)ret++;
  }
  return ret;
}

int main(){
  while(true){
    scanf("%d",&n);
    if(n==0)break;
    REP(i,n){
      double x,y;
      scanf("%lf%lf",&x,&y);
      pts[i] = P(x,y);
    }
    int ans = 0;
    REP(i,n)REP(j,i){
      P a,b;
      a = pts[i];
      b = pts[j];
      double l = abs(b-a);
      if(l>2.0)continue;
      double h = sqrt(1.0-l*l/4.0);
      P c;
      c = a + (b-a)/2.0 + (b-a)*P(0,1)*h/l;
      ans = max(ans,check(c));
      c = a + (b-a)/2.0 + (b-a)*P(0,-1)*h/l;
      ans = max(ans,check(c));
    }
    printf("%d\n",ans);
  }
  return 0;
}

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

using namespace std;

int N;
pair<double,double> dots[300];
#define EPS (1e-6)


int main(){

	while(cin >> N && N != 0){
		for(int i = 0; i < N; i++){
			pair<double,double> p;
			cin >> dots[i].second >> dots[i].first;
		}
		int maxCnt = 1;
		// ~ðìé½ßÌQ_ðßé
		for(int i = 0; i < N; i++){
			for(int j = i+1; j < N; j++){
				pair<double,double> dot1,dot2;
				dot1 = dots[i];
				dot2 = dots[j];

				double v = sqrt((dot1.first - dot2.first)*(dot1.first - dot2.first)
					+(dot1.second - dot2.second)*(dot1.second - dot2.second));
				// ñ_Ô£ª2æè¨¨«¯êÎÂ­êÈ¢
				if(v > 2.0 + EPS){
					continue;
				}

				double x = sqrt(1 - (v/2)*(v/2));

				// ñgÌPÊ@üxNg
				pair<double,double> hose[2];
				hose[0].first = -1.0/(dots[i].first - dots[j].first);
				hose[0].second = 1.0/(dots[i].second - dots[j].second);
				hose[1].first = 1.0/(dots[i].first - dots[j].first);
				hose[1].second = -1.0/(dots[i].second - dots[j].second);

				double dd1 = sqrt((hose[0].second)*(hose[0].second)+(hose[0].first)*(hose[0].first));
				double dd2 = sqrt((hose[1].second)*(hose[1].second)+(hose[1].first)*(hose[1].first));
				hose[0].first /= dd1; hose[0].second /= dd1;
				hose[1].first /= dd2; hose[1].second /= dd2;

				// x{·é
				hose[0].first *= x;hose[1].first *= x;
				hose[0].second *= x;hose[1].second *= x;

				// mð«·
				hose[0].first += (dot1.first + dot2.first)/2;
				hose[0].second += (dot1.second + dot2.second)/2;
				hose[1].first += (dot1.first + dot2.first)/2;
				hose[1].second += (dot1.second + dot2.second)/2;

				// hose[0]Æhose[1]Í~ÌSÀW
				// e_ÉÖµÄ~ÌÉ¶Ý·é©Ç¤©ðßé
				for(int k = 0; k < 2; k++){
					int cnt = 0;
					for(int l = 0; l < N; l++){
						double dist = sqrt((dots[l].first -hose[k].first)*(dots[l].first -hose[k].first)
							+ (dots[l].second -hose[k].second)*(dots[l].second -hose[k].second));
						if(dist <= 1.0 + EPS){
							cnt++;
						}
					}
					maxCnt = max(maxCnt,cnt);
				}
			}
		}
		cout << maxCnt << endl;
	}

	return 0;
}

C++: #include <iostream>
#include <cmath>
#define eps 1.0e-9

using namespace std;

int N;
double x[305], y[305];

int count(double cx, double cy)
{
	int ret = 0;
	for(int i = 0; i < N; i++){
		if( (cx - x[i])*(cx - x[i]) + (cy - y[i])*(cy - y[i]) <= 1 + eps) ret++;
	}
	return ret;
}

int main(void)
{
	while(1){
		cin >> N;
		if(N == 0) break;
		for(int i = 0; i < N; i++) cin >> x[i] >> y[i];
		
		int ans = -1;
		double cx, cy;
		for(int i = 0; i < N; i++){
			for(int j = 0; j < N; j++){
				if(i >= j) continue;
				
				double px = x[j] - x[i];
				double py = y[j] - y[i];
				double pnorm = sqrt(px*px + py*py);
				px /= pnorm;
				py /= pnorm;
				
				double mx = (x[i] + x[j]) / 2.0;
				double my = (y[i] + y[j]) / 2.0;
				
				double qnorm = 1 - pnorm*pnorm/4.0;
				if(qnorm < eps) continue;
				qnorm = sqrt(qnorm);
				double qx = qnorm * py;
				double qy = qnorm * -px;
				
				ans = max(ans, count(mx + qx, my + qy));
				ans = max(ans, count(mx - qx, my - qy));
			}
		}
		cout << ans << endl;
	}
	
	return 0;
}