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>

using namespace std;

#define EPS (1e-6)


void calc(double x1,double y1,double x2,double y2,double *c_x1,double *c_y1,double *c_x2,double *c_y2){
  double midx,midy,hox,hoy,dx,dy,dis;
  midx = (x1+x2) / 2.0 , midy = (y1+y2) / 2.0;
  dx = x1-x2 , dy = y1-y2;
  dis = sqrt(1 - pow(x1-midx,2.0) - pow(y1-midy,2.0));
  hox = sqrt(pow(dy,2.0) / (pow(dy,2.0) + pow(dx,2.0)));
  hoy = hox * (-1.0) * (dx / dy);
  *c_x1 = midx + hox*dis,*c_y1 = midy + hoy*dis;
  *c_x2 = midx - hox*dis,*c_y2 = midy - hoy*dis;
}

int main(){
  int N;

  while(cin >> N){
    double x[300],y[300],c_x1,c_y1,c_x2,c_y2;
    int count,m;
    m = 0;
    if(!(N)) break;
    for(int i=0;i<N;i++) cin >> x[i] >> y[i];

    for(int i=0;i<N;i++){
      for(int j=0;j<N;j++){
        if(i == j) continue;
        if(sqrt((x[i] - x[j])*(x[i] - x[j]) + (y[i] - y[j])*(y[i] - y[j])) > 2.0 + EPS  ) continue;
        count = 2;
        calc(x[i],y[i],x[j],y[j],&c_x1,&c_y1,&c_x2,&c_y2);
        for(int k=0;k<N;k++){
          if(k == i || k == j) continue;
          if((x[k] - c_x1)*(x[k] - c_x1) + (y[k] - c_y1)*(y[k] - c_y1) < 1.0 + EPS  ) count++;
        }
        if(m < count) m = count;
        count = 2;
        for(int k=0;k<N;k++){
          if(k == i || k == j) continue;
          if((x[k] - c_x2)*(x[k] - c_x2) + (y[k] - c_y2)*(y[k] - c_y2) < 1.0 + EPS  ) count++;
        }
        if(m < count) m = count;
      }
    }

    printf("%d\n",m);
  }

  return 0;
}

C++: #include <iostream>
#include <vector>
#include <string>
#include <complex>
#include <algorithm>
using namespace std;

typedef complex<double> p;
const double eps = 1e-9;

int n;
p ps[300];

int main() {
	while (cin >> n , n != 0) {
		for (int i = 0; i < n; i++) {
			double x, y;
			cin >> x >> y;
			ps[i] = p(x, y);
		}

		if (n == 1) {
			cout << 1 << endl;
			continue;
		}

		int ans = 0;
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				p v = ps[i] - ps[j], ve = v / abs(v);
				double l = sqrt(1 - norm(v) * 0.25);
				p vn1 = ve * p(0, 1) * l, vn2 = ve * p(0, -1) * l;
				p c1 = ps[j] + v * 0.5 + vn1, c2 = ps[j] + v * 0.5 + vn2;
				int count1 = 2, count2 = 2;
				for (int k = 0; k < n; k++) {
					if (k == i || k == j) continue;
					if (norm(c1 - ps[k]) < 1.0 + eps) count1++;
					if (norm(c2 - ps[k]) < 1.0 + eps) count2++;
				}
				ans = max(ans, count1);
				ans = max(ans, count2);
			}
		}

		cout << ans << endl;
	}
}

C++: #include <iostream>
#include <stdio.h>
#include <algorithm>
#include <vector>
#include <queue>
#include <numeric>
#include <string>
#include <string.h>
#include <map>
#include <set>
#include <functional>
#include <complex>
#include <stack>
#include <tuple>
#include <array>
using namespace std;

#define FOR(i,a,b) for (int i=(a);i<(b);i++)
#define RFOR(i,a,b) for (int i=(b)-1;i>=(a);i--)
#define REP(i,n) for (int i=0;i<(n);i++)
#define RREP(i,n) for (int i=(n)-1;i>=0;i--)
#define ALL(u) begin(u),end(u)
#define PB push_back
#define LE(n,m) ((n) < (m) + EPS)
#define GE(n,m) ((n) + EPS > (m))
#define EQ(n,m) (abs((n)-(m)) < EPS)

typedef long long int ll;
const int INF = (1<<30) - 1;
const double EPS = 1e-9;
const int MOD = 1000000007;
typedef double D;      // ??§?¨?????????????double???long double?????????
typedef complex<D> P;  // Point
typedef pair<P, P> L;  // Line
typedef vector<P> VP;
#define X real()
#define Y imag()


// ???a??¨???b?????????????????????r????????????????????????
VP circlesPointsRadius(P a, P b, D r) {
    VP cs;
    P abH = (b-a)*0.5;
    D d = abs(abH);
    if (d == 0 || d > r) return cs;	// ???????????? !LE(d,r) ??¨?????????1??????????????´????????????
    D dN = sqrt(r*r - d*d);					// ???????????? max(r*r - d*d, 0) ??¨??????
    P n = abH * P(0,1) * (dN / d);
    cs.push_back(a + abH + n);
    if (dN > 0) cs.push_back(a + abH - n);
    return cs;
}

int N;

//g++ -std=c++0x -msse4.2 -O3
//#include <bits/stdc++.h>
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);
    //cout.precision(16);
    //cout.setf(ios::fixed);
    
    while(cin >> N, N){
        vector<P> points(N);
        for(int i=0; i<N; i++){
            D x, y; cin >> x >> y;
            points[i] = P(x, y);
        }
        
        int ans = 0;
        for(int i=0; i<N; i++)
            for(int j=i+1; j<N; j++){
                VP vs = circlesPointsRadius(points[i], points[j], 1);
                for(auto v: vs){
                    int tans = 0;
                    for(auto p: points)
                        if(LE(abs(p-v), 1))
                            tans++;
                    ans = max(ans, tans);
                }
            }
        cout << ans << endl;
    }
    
    return 0;
}

C++: #include <cstdio>
#include <cmath>


int main() {
	double x[300];
	double y[300];

	while (true) {
		int N;

		scanf("%d", &N);
		if (N == 0) break;

		for (int i = 0; i < N; i++) {
			scanf("%lf%lf", &x[i], &y[i]);
		}

		int max = 1;
		for (int i = 0; i < N; i++) {
			for (int j = 0; j < N; j++) {
				double dx = x[i] - x[j];
				double dy = y[i] - y[j];
				double dist = hypot(dx, dy);
				if (dist > 2.0) continue;

				double mx = (x[i] + x[j]) / 2;
				double my = (y[i] + y[j]) / 2;
				double l = sqrt(1.0 - dist * dist / 4);
				double cx = mx + (-dy) * l / dist;
				double cy = my + dx * l / dist;
				int count = 0;
				for (int k = 0; k < N; k++) {
					if (k == i || k == j ||
						(cx - x[k]) * (cx - x[k]) + (cy - y[k]) * (cy - y[k]) <= 1.0) count++;
				}

				if (max < count) max = count;
			}
		}

		printf("%d\n", max);
	}
	return 0;
}