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 <bits/stdc++.h>
#define PB push_back
#define MP make_pair
#define REP(i,n) for (int i=0;i<(n);i++)
#define FOR(i,a,b) for(int i=(a);i<(b);i++)
#define ALL(a) (a).begin(),(a).end()
int INF=1e9;
int MOD=1000000007;
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
double EPS = 1e-10;
double add(double a,double b){
	if(abs(a+b)<EPS*(abs(a)+abs(b)))return 0;
	return a+b;
}
struct P{
	double x,y;
	P(){}
	P(double x,double y):x(x),y(y){}
	P operator + (P p){
		return P(add(x,p.x),add(y,p.y));
	}
	P operator - (P p){
		return P(add(x,p.x),add(y,-p.y));
	}
	P operator * (P p){
		return P(x*p.x,y*p.y);
	}
};
double dist(P a,P b){
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int cnt(P pos,vector<P> dat){
	int res=0;
	REP(i,dat.size()){
		if(dist(pos,dat[i])<=1.0)res++;
	}
	return res;
}
int main(){
	int n;
	double x,y;
	while(scanf("%d",&n)&&n){
		vector<P> dat;
		REP(i,n){
			cin>>x>>y;
			dat.PB(P(x,y));
		}
		int ans=0;
		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				P a=dat[i],b=dat[j];
				P mid=P((a.x+b.x)/2,(a.y+b.y)/2);
				if(dist(a,b)>2){
					ans=max(ans,1);
					continue;
				}
				double dx=a.y-b.y,dy=b.x-a.x;
				double dd=sqrt(dx*dx+dy*dy);
				double ddd=sqrt(1-dist(a,b)/2*dist(a,b)/2);
				dx=dx/dd*ddd;dy=dy/dd*ddd;
				ans=max(ans,cnt(P(mid.x+dx,mid.y+dy),dat));
				ans=max(ans,cnt(P(mid.x-dx,mid.y-dy),dat));
			}
		}
		printf("%d\n",ans);
	}
	return 0;
}

C++: #include <iostream>
#include <complex>
#include <cmath>
using namespace std;

#define EPS 1.0e-8
typedef complex<double> P;

int main(){
  int n;
  P p[301];

  while( cin>>n && n ){
    for( int i=0;i<n;i++ )
      cin >> p[i].real() >> p[i].imag();

    int ans=0;
    for( int i=0;i<n;i++ ){
      for( int j=0;j<n;j++ ){
	if( i==j ) 	  continue;
	if( abs( p[i]-p[j] )+EPS > 2.0 ) continue;

	P m  = (p[i]+p[j]) / 2.0;         // p[i]Æp[j]Ì_
	P nn = (p[i]-p[j]) * P(0.0,1.0);  // @üxNg
	nn = nn / abs(nn);

	double x = abs( p[i]-p[j] )/2.0;
	P c = m + nn*sqrt( 1.0-x*x );     // ~ÌS

	int num=0;
	for( int k=0;k<n;k++ )
	  if( abs( c-p[k] ) < 1.0 + EPS  )
	    num++;
	ans = max( ans,num );
      }
    }
    cout << ans << endl;
  }

  return 0;
}

C++: #include <iostream>
#include <cstdio>
#include <iomanip>
#include <vector>
#include <map>
#include <set>
#include <queue>
#include <bitset>
#include <stack>
#include <utility>
#include <numeric>
#include <algorithm>
#include <functional>
#include <cctype>
#include <complex>
#include <string>
#include <sstream>

using namespace std;

#define all(c) c.begin(),c.end()
#define rall(c) c.rbegin(),c.rend()
#define mp(a,b) make_pair((a),(b))
#define eq ==

typedef long long ll;
typedef complex<double> point;
typedef pair<int,int> pii;

// →↑←↓
const int dx[] = {1,0,-1,0};
const int dy[] = {0,-1,0,1};


const double EPS = 1e-9;

int main(){
    while(true){
        int N;
        cin >> N;
        if(N == 0) break;
        vector<point> P(N);
        for(int i=0;i<N;i++){
            cin >> P[i].real() >> P[i].imag();
        }

        int ret = 1;

        for(int i=0;i<N;i++){
            for(int j=0;j<N;j++){
                if(i == j) continue;
                point p1 = P[i];
                point p2 = P[j];

                if(abs(p1-p2) <= 2.0+EPS){
                    point d = p1 - p2;
                    point m = (p1+p2)/(point(2,0));

                    point nor = d * point(0,1);
                    nor /= abs(nor);
                    nor *= point(sqrt(1-(abs(d)/2)*(abs(d)/2)),0);

                    point p = m + nor;
                    int cnt = 0;
                    for(int k=0;k<N;k++){
                        if(abs(P[k]-p) <= 1+EPS) cnt++;
                    }
                    ret = max(ret,cnt);
                }
            }
        }
        cout << ret << endl;
    }
    return 0;
}

C++: //???????????????????????????, ?????????????????????????????°????????????????????????????????????????????????2???????????¨???????????£?????\????????¨?????§????????????
//?????£???, 2??????????????????1?????????????????????????????????????????°??????????????????????£????, O(N^2)????????????????????????, ??¨?????§O(N^3)????????§?§£?????????
#include <iostream>
#include <complex>
#include <vector>
using namespace std;
typedef complex<double> Point;

int n;
double x[300], y[300];
vector<Point> points;	//????£???????

int main() {
	while (cin >> n) {
		if (n == 0) break;
		for (int i = 0; i < n; i++) cin >> x[i] >> y[i];
		
		for (int i = 0; i < n; i++) {
			for (int j = i + 1; j < n; j++) {
				Point p1 = Point(x[i], y[i]);
				Point p2 = Point(x[j], y[j]);
				if (norm(p1 - p2) < 4) {
					double h = sqrt(1.0 - norm(p1 - p2) / 4.0);
					Point t = p2 - p1; t /= abs(t);
					Point u = t * Point(0, 1);
					points.push_back(p1 + t * abs(p1 - p2) / 2.0 + u * h);
					points.push_back(p1 + t * abs(p1 - p2) / 2.0 - u * h);
				}
			}
		}
		
		double eps = 1e-10;
		int ans = 1;
		for (int i = 0; i < points.size(); i++) {
			double px = points[i].real();
			double py = points[i].imag();
			int cnt = 0;
			for (int j = 0; j < n; j++) {
				if ((x[j] - px) * (x[j] - px) + (y[j] - py) * (y[j] - py) < 1 + eps) {
					cnt++;
				}
			}
			if (ans < cnt) ans = cnt;
		}
		cout << ans << endl;
	}
	return 0;
}