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

// * XYツ催?標
#define X real()
#define Y imag()

// * ツ点ツづ個表ツ個サ
typedef complex<double> P;

// * ツ仰鳴容ツつキツづゥツ古ォツ債キツε?
const double EPS = 1e-10;
// * Infinity
const double INF = 1e12;
// * ツ円ツ偲シツ猟ヲ
const double PI = acos(-1.0);
// * 
const P INF_P(INF,INF);

// * complex<double> ツづ個渉?渉伉づーツ津ィツ義ツつキツづゥ
namespace std {
	bool operator < (const P& a, const P& b) {
		return (real(a) != real(b))? real(a) < real(b) : imag(a) < imag(b);
	}
}

// * ツ点ツ催?標ツづ個デツバツッツグツ出ツ療?
void print(P p) {
	cout << "point" << p << ";" << endl;
}

// * 2ツ湘ヲツつキツづゥ
double sq(double x) {
	return x * x;
}

// * 2ツ点ツ甘板づ個仰猟猟」
double dist(P p1, P p2) {
	return abs(p1 - p2);
}

// * 2ツづつづ個スツカツδ可ーツつェツ督卍つオツつ「ツつゥツづ?つ、ツつゥ
bool equal(double a, double b) {
	return ( fabs( a - b ) < EPS );
}
// * 2ツづつづ個ベツクツトツδ仰つェツ督卍つオツつ「ツつゥツづ?つ、ツつゥ
bool eq_v(P a, P b) {
	return ( equal( a.X, b.X ) && equal( a.Y, b.Y ) );
}

// * ツ禿?静?(dot product) : aツ・b = |a||b|cosツδヲ
double dot(P a, P b) {
	return real( conj(a) * b );
}

// * ツ外ツ静?(cross product) : aツ×b = |a||b|sinツδヲ
double cross(P a, P b) {
	return imag( conj(a) * b ) ;
}

// ツベツクツトツδ?a ツづ個単ツ暗環法ツ静シツベツクツトツδ仰づーツ仰?づ淞づゥ
vector<P> normal_unit(P a) {
	vector<P> vp;
	vp.push_back( (a * P(0,  1)) / abs(a) );
	vp.push_back( (a * P(0, -1)) / abs(a) ) ;
	return vp;
}

// ツ個エツ点ツづーツ篠イツづ可点 p ツづーツ角ツ度 a (ツδ可ジツアツδ?ツつセツつッツ嘉アツ転ツつウツつケツつスツ点ツづーツ陛板つキ
P rot(P p, double a){
	double x = p.X * cos(a) - p.Y * sin(a);
	double y = p.X * sin(a) + p.Y * cos(a);
	return P(x,y);
}

int main(){
	int n;
	while( cin >> n , n ){
		vector<P> v;
		int ans=1;
		for(int i=0 ; i < n ; i++ ){
			double x, y;
			scanf("%lf %lf", &x, &y);
			v.push_back( P(x,y) );
		}
		for(int i=0 ; i < v.size() ; i++ ){
			for(int j=i+1 ; j < v.size() ; j++ ){
				if( i == j ) continue;
				P p1 = v[i];
				P p2 = v[j];
				double d = dist(p1,p2);
				if( d >= 2.0 ) continue;
				
				double angle = atan2( p2.Y-p1.Y , p2.X-p1.X );
				P mid = rot( p2-p1 , -angle ) / 2.0;
				double h = sqrt(1.0 - sq(d/2.0));
				P p3 = P( mid.X , h );
				P p4 = P( mid.X , -h );
				int cnt1=0, cnt2=0;
				for(int k=0 ; k < v.size() ; k++ ){
					P p = rot( v[k]-p1 , -angle );
					if( dist(p3,p) <= 1.0 ) cnt1++;
					if( dist(p4,p) <= 1.0 ) cnt2++;
				}
				ans = max( ans , max(cnt1,cnt2) );
			}
		}
		printf("%d\n", ans);
	}
}

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

using namespace std;

inline int toInt(string s) {int v; istringstream sin(s);sin>>v;return v;}
template<class T> inline string toString(T x) {ostringstream sout;sout<<x;return sout.str();}
template<class T> inline T sqr(T x) {return x*x;}

typedef vector<int> vi;
typedef vector<vi> vvi;
typedef vector<string> vs;
typedef pair<int, int> pii;
typedef long long ll;

#define all(a)  (a).begin(),(a).end()
#define rall(a) (a).rbegin(), (a).rend()
#define pb push_back
#define mp make_pair
#define each(i,c) for(typeof((c).begin()) i=(c).begin(); i!=(c).end(); ++i)
#define exist(s,e) ((s).find(e)!=(s).end())
#define range(i,a,b) for(int i=(a);i<(b);++i)
#define rep(i,n)  range(i,0,n)
#define clr(a,b) memset((a), (b) ,sizeof(a))
#define dump(x)  cerr << #x << " = " << (x) << endl;
#define debug(x) cerr << #x << " = " << (x) << " (L" << __LINE__ << ")" << " " << __FILE__ << endl;

const double eps = 1e-10;
const double pi  = acos(-1.0);
const ll INF =1LL << 62;
const int inf =1 << 29;

// ??\??? ????????¨???
template<class T> bool operator==(T a, T b){return abs(a-b)< eps;}
template<class T> bool operator< (T a, T b){return a < b-eps;}
template<class T> bool operator<=(T a, T b){return a < b+eps;}
template<class T> int sig(T r) {return (r==0||r==-0) ? 0 : r > 0 ? 1 : -1;}
#define X real()
#define Y imag()

typedef bool B;
typedef long double D;
typedef complex<D> P;

P p[310];

int main(void){
	int n;
	while(cin >> n,n){
		D ix,iy;
		rep(i,n){
			cin >> ix >> iy;
			p[i].real(ix),p[i].imag(iy);
		}
		int ans=1;
		rep(i,n)rep(j,n){
			if(i<=j) continue;
			D d=norm(p[j]-p[i]);
			if(d>4.0) continue;
			P v=(p[j]-p[i])/P(2.0,0.0);
			D s=sqrt((1.0-norm(v))/norm(v));
			P cnt1=p[i]+v+s*P(v.imag(),-v.real());
			int num=2;
			rep(k,n){
				if(k==i||k==j) continue;
				if(norm(p[k]-cnt1)<1.0) num++;
			}
			ans=max(ans,num);
			P cnt2=p[i]+v+s*P(-v.imag(),v.real());
			int num2=2;
			rep(k,n){
				if(k==i||k==j) continue;
				if(norm(p[k]-cnt2)<1.0) num2++;
			}
			ans=max(ans,num2);
		}
		cout << ans << endl;
	}
	return 0;
}

C++: #include <bits/stdc++.h>
typedef long long LL;
#define SORT(c) sort((c).begin(),(c).end())

#define FOR(i,a,b) for(int i=(a);i<(b);++i)
#define REP(i,n)  FOR(i,0,n)

using namespace std;
pair<double,double> smul(pair<double,double> a,double k){
  return make_pair(a.first*k,a.second*k);
}
pair<double,double> sub2(pair<double,double> a,pair<double,double> b){
  return make_pair(a.first-b.first,a.second-b.second);
}
pair<double,double> mean2(pair<double,double> a,pair<double,double> b){
  return make_pair((a.first+b.first)/2,(a.second+b.second)/2);
}
double dist(pair<double,double> a,pair<double,double> b){
  return sqrt(pow(a.first-b.first,2)+pow(a.second-b.second,2));
}
int main(void)
{
  for(;;){
    int n;
    cin >> n;
    if(!n) return 0;
    vector<pair<double,double> > ps;
    ps.resize(n);
    REP(i,n) cin >> ps[i].first >> ps[i].second;
    vector<pair<double,double> > co;
    co.clear();
    REP(j,n) REP(i,j){
      if(dist(ps[i],ps[j])>2) continue;
      pair<double,double> m=mean2(ps[i],ps[j]);
      double k=sqrt(1-pow(dist(ps[i],ps[j])/2.0,2.0));
      pair<double,double> tmp=sub2(ps[i],ps[j]);
      tmp=smul(make_pair(tmp.second,-tmp.first),k/dist(ps[i],ps[j]));
      co.push_back(sub2(m,tmp));
      tmp=smul(tmp,-1.0);
      co.push_back(sub2(m,tmp));
    }
    int answer=1;
    for(auto x:co){
      int tmpans=0;
      REP(i,n) if(dist(x,ps[i])<1.000001) ++tmpans;
      answer=max(answer,tmpans);
    }
    cout << answer << endl;

  }
}