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 <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <stack>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <string>
#include <cstring>
#include <ctime>
#include <queue>
#include <complex>

 
using namespace std;
 
#define INF 1000000000
#define EPS 1e-9
#define PI acos(-1)
 
typedef long long ll;
typedef complex<double> vec;

#define MAX_N 300
int N;
double x[MAX_N], y[MAX_N];
vector<int> ans_list;

int par[MAX_N];

int main(){

  while(true){
    cin >> N;
    if(N == 0)
      break;
    for(int i = 0; i < N; i++){
      cin >> x[i] >> y[i];
    }

    int ans = 0;

    for(int i = 0; i < N; i++){
      for(int j = i+1; j < N; j++){
	if(sqrt((x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j])) > 2.0001)
	  continue;
	int tmp = 0;
	vec p1 = vec(x[i], y[i]), p2 = vec(x[j], y[j]);
	vec mid = (p2-p1);
	mid /= 2;
	vec r1 = (mid * vec(0, 1)) / abs(mid);
	vec r2 = (mid * vec(0, -1)) / abs(mid);

	vec R = p1 + mid + r1*sqrt((1-abs(mid)*abs(mid)));
	double rx = R.real();
	double ry = R.imag();
	tmp = 2;
	for(int k = 0; k < N; k++){
	  if(k == i || k == j)
	    continue;
	  if((rx - x[k])*(rx-x[k]) + (ry - y[k])*(ry - y[k]) < 1)
	    tmp++;
	}
	ans = max(ans, tmp);
	R = p1 + mid + r2*sqrt((1-abs(mid)*abs(mid)));
	rx = R.real();
	ry = R.imag();
	tmp = 0;
	for(int k = 0; k < N; k++){
	  if(k == i || k == j)
	    continue;
	  if((rx - x[k])*(rx-x[k]) + (ry - y[k])*(ry - y[k]) < 1)
	    tmp++;
	}
	ans = max(ans, tmp);
      }
    }

    ans_list.push_back(ans);
    
  }
  
  for(int i = 0; i < ans_list.size(); i++){
    cout << ans_list[i] << endl;
  }

  return 0;

}

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

#define EPS (1e-9)
#define prev(P,i) P[(i+P.size()-1)%P.size()]
#define curr(P,i) P[i]
#define next(P,i) P[(i+1)%P.size()]
typedef complex < double > P; // P:点
struct L:public vector< P >{ // L:線分
  L(){resize(2);}
  L(const P& a,const P& b){push_back(a);push_back(b);}
};
struct C{ //C:円
  P p; double r;
  C():p(0),r(0){}
  C(const P& p, double r):p(p),r(r){};
};

double cross(const P& a,const P& b){
  return imag(conj(a) * b);
}
double dot(const P& a,const P& b){
  return real(conj(a) * b);
}
P projection(const L& l,const P& p){
  double t = dot(p-l[0],l[0]-l[1]) / norm(l[0]-l[1]);
  return l[0] + t * (l[0] - l[1]);
}
int ccw(P a,P b,P c){
  b -= a; c -= a;
  if(cross(b,c) > 0) return 1;
  if(cross(b,c) < 0) return -1;
  if(dot(b,c) < 0) return 2; // c--a--b
  if(norm(b) < norm(c)) return -2; // a--b--c
  return 0; // a--c--b
}
bool interectCC(const C& c,const C& d){
  return (norm(c.p - d.p) - (c.r + d.r)*(c.r + d.r) < EPS) &&
    (norm(c.p - d.p) - (c.r - d.r)*(c.r - d.r) > -EPS);
}
pair< P , P > crosspointCC(const C& c,const C& q){
  P a = conj(q.p - c.p);
  P b = q.r  * q.r - c.r * c.r - (q.p - c.p) * a;
  P C = c.r * c.r * (q.p - c.p);
  P d = b * b - 4.0 * a * C;
  P z1 = (-b * sqrt(d)) / (2.0 * a) + c.p;
  P z2 = (-b - sqrt(d)) / (2.0 * a) + c.p;
  return pair< P , P >( z1, z2);
}
int main(){
  int n;
  while(cin >> n , n){
    vector < P > pts(n);
    int ans = 1;
    for(int i = 0 ; i < n ; i++ ){
    cin >> pts[i].real() >> pts[i].imag();
    }
    for(int i = 0 ; i < n ; i++ ){
      C c1 = C( pts[i], 1);
      for(int j = 0 ; j < i ; j++ ){
        C c2 = C( pts[j], 1);
        if(interectCC(c1,c2)){
          pair< P , P > nnr = crosspointCC( c1, c2);
          int ret = 0, rec = 0;
          for(int k = 0 ; k < n ; k++ ){
            if(abs(nnr.first - pts[k]) <= 1 + EPS) ret++;
            if(abs(nnr.second - pts[k]) <= 1 + EPS) rec++;
          }
          ans = max( ans, max( ret, rec));
        }
      }
    }
    cout << ans << endl;
  }
}

C++: #include<iostream>
#include<vector>
#include<complex>
#include<algorithm>
#include<cmath>
using namespace std;
#define REP(i, j) for(int i = 0; i < j; i++)
#define FOR(i, j, k) for(int i = j; i < k; i++)
#define P complex<double>
#define X real()
#define Y imag()
const double EPS = 1e-9;

int n;
vector<P> vv;

int count_points(P a, P b, double mx, double my){
    P v = (a - b);
    double d_v = abs(v);
    if(d_v >= 2.0 + EPS) return 0;
    P m = v / 2.0 + b;
    double half_d_v = d_v / 2;
    double d = sqrt(1.0 - half_d_v * half_d_v);
    P u = (v * P(mx, my)) / abs(v) * d;
    int cnt = 0;
    REP(i, n){
        P c = m + u;
        if(abs(c - vv[i]) <= 1.0 + EPS) cnt++;
    }
    return cnt;
}

int main(){
    while(cin >>n && n){
        vv = vector<P>(n);
        REP(i, n){
          double a, b; cin >>a >>b;
          vv[i] = P(a, b);
        }

        int ans = 1;
        REP(a, n){
            FOR(b, a + 1, n){
                ans = max(ans, count_points(vv[a], vv[b], 0, 1));
                ans = max(ans, count_points(vv[a], vv[b], 0, -1));
            }
        }
        cout <<ans <<endl;
    }
    return 0;
}