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<complex>
#include<vector>
#include<algorithm>
#include<cmath>
#include<map>
#include<list>

#define EPS (1e-10)
#define EQ(a,b) (abs((a)-(b)) < EPS)
#define fs first
#define sc second
#define pb push_back
#define sz size()
#define FOR(i,a,b) for(int i=(int)(a);i<(int)(b);i++)

using namespace std;

typedef complex<double> P;
typedef pair<P,P> L;
typedef vector<P> Poly;
typedef pair<P,double> C;

const double PI = acos(-1);

//for vector
P unit(P p){return p / abs(p);}

pair<P,P> norm(P p){return make_pair(p*P(0,1),p*P(0,-1));}

double dot(P x,P y){return real(conj(x)*y);}

double cross(P x,P y){return imag(conj(x)*y);}

//rotate a point counter-clockwise on the origin
P rotate(P v,double s){
  return P(real(v)*cos(s) - imag(v)*sin(s), real(v)*sin(s) + imag(v)*cos(s) );
}


//return seta A
double arg(double a,double b,double c){return acos( (b*b+c*c-a*a)/(2*b*c) );}


//for circle
double area_cir(C c){return PI * c.sc * c.sc;}

bool in_cir(C c,P x){return (abs(x-c.fs) +EPS < c.sc);}
bool on_cir(C c,P x){return EQ(abs(x-c.fs),c.sc);}

//circle position relation
int cpr(C a,C b){
  double d = abs(a.fs-b.fs);
  if(a.sc+b.sc + EPS < d)return -1;      //no cross point (outside)
  if(b.sc+d + EPS < a.sc)return 1;       //no cross point (inside,B in A)
  if(a.sc+d + EPS < b.sc)return 2;       //no cross point (inside,A in B)
  //------------above verified----------------// 
  if(abs(a.sc+b.sc - d) < EPS)return -3; //one cross point (outside)
  if(abs(b.sc+d - a.sc) < EPS)return 3;  //one cross point (inside,B in A)
  if(abs(a.sc+d - b.sc) < EPS)return 4;  //one cross point (inside,A in B)
  return 0;                              //two cross point
}
  
//not verified
vector<P> cp_cir_to_cir(C a, C b){
  vector<P> v;
  int pos = cpr(a,b);
  if(pos==0){
    double s = arg(b.sc,abs(b.fs-a.fs),a.sc);
    P x = a.sc * unit(b.fs - a.fs);
    v.pb(a.fs + rotate(x,s));
    v.pb(a.fs + rotate(x,-s));
  }else if(abs(pos) == 3){
    v.pb(a.fs + a.sc * unit(b.fs-a.fs));
  }else if(pos == 4){
    v.pb(a.fs + a.sc * unit(b.fs-a.fs));
  }
  return v;
}

int n;
C c[500];

int main(){
  while(scanf("%d",&n),n){
    for(int i=0;i<n;i++){
      scanf("%lf%lf",&c[i].fs.real(),&c[i].fs.imag());
      c[i].sc = 1.0;
    }

    int pos;
    vector<P> s,cp;

    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
	pos = cpr(c[i],c[j]);
	if(pos==1)s.pb(c[i].fs);
	else if(pos==2)s.pb(c[j].fs);
	else{
	  cp = cp_cir_to_cir(c[i],c[j]);
	  for(int x=0;x<cp.sz;x++)s.pb(cp[x]);
	}
      }
    }

    int ans = 1;
    for(int i=0;i<s.sz;i++){
      int num = 0;
      for(int j=0;j<n;j++){
	if(abs(s[i]-c[j].fs) < 1.0 + EPS)num++;
      }
      if(ans<num)ans = num;
    }
    printf("%d\n",ans);
  }
}

C++: #include <iostream>
#include <complex>
#include <cmath>

using namespace std;

typedef complex<double> xy_t;

const xy_t INF_VECT = xy_t(100, 100);
const xy_t ZERO_VECT = xy_t(0, 0);

xy_t center(xy_t p, xy_t q);

int main()
{
    int num;
    while(1){
        cin >> num;
        if(num == 0) break;
        
        xy_t point[num];
        for(int i = 0; i < num; i++){
            double x, y;
            cin >> x >> y;
            point[i] = xy_t(x, y);
        }

        int max = 0;
        for(int i = 0; i < num-1; i++){
            for(int j = i+1; j < num; j++){
                xy_t vec = center(point[i], point[j]);
                if(vec != ZERO_VECT && vec != INF_VECT){
                    xy_t mid = (point[i] + point[j]) * xy_t(0.5, 0);
                    xy_t cent = mid + vec;
                    int cnt = 0;
                    for(int k =0; k < num; k++){
                        xy_t r = point[k] - cent;
                        if(abs(r) <= 1) cnt++;
                    }
                    if(cnt > max) max = cnt;

                    cent = mid - vec;
                    cnt = 0;
                    for(int k =0; k < num; k++){
                        xy_t r = point[k] - cent;
                        if(abs(r) <= 1) cnt++;
                    }
                    if(cnt > max) max = cnt;
                }
            }
        }
        
        cout << max << endl;
        
    }
    return 0;
}


xy_t center(xy_t p, xy_t q)
{    
    xy_t r = q-p;
    double d = abs(r);

    if(d > 2){
        return INF_VECT;
    }
    else if(d == 2){
        return ZERO_VECT;
    }
    else{
        double l = sqrt(4 - d*d) / 2;
        xy_t r_ = r * xy_t(0, 1);
        r_ *= xy_t(l/d, 0);
        return r_;
   }
}

C++: #include<iostream>
#include<vector>
#include<complex>
#include<algorithm>
#include<map>
using namespace std;
#define EPS (1e-9)
typedef complex < double > P; // P:点
struct C{ //C:円
  P p; double r;
  C():p(0),r(0){}
  C(const P& p, double r):p(p),r(r){};
};
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;
  }
}