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>
#define EPS 1.0e-4
using namespace std;

typedef struct {
  double x,y;
} P;

double dis(P p ){
  return sqrt( p.x*p.x + p.y*p.y );
}
double dis(P p1,P p2 ){
  P p;
  p.x = p1.x-p2.x;
  p.y = p1.y-p2.y;
  return dis( p );
}


int main(){
  int n, mx;
  // mv : @üxNg
  P p[301], m, p1,p2, mv,dmv, c[2];
  double x;

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

    mx=0;

    for( int i=0;i<n;i++ ){
      for( int j=i+1;j<n;j++ ){
	if( dis(p[i],p[j]) < 2.0 + EPS ){
	  p1=p[i]; p2=p[j];
	  m.x = (p1.x+p2.x)/2.0;
	  m.y = (p1.y+p2.y)/2.0;

	  x = sqrt( 1-dis(p1,m) );
	  // S_ cðßé
	  if( fabs( p1.y-p2.y ) < EPS ){
	    mv.x=1.0;  mv.y=0.0;
	  }else if( fabs( p1.x-p2.x ) < EPS ){
	    mv.x=0.0;  mv.y=1.0;
	  }else{
	    double a1 = ( p2.y-p1.y )/( p2.x-p1.x );
	    mv.x=1.0; mv.y=-1.0/a1;
	  }
	  double scale = dis( mv );
	  dmv.x = x*mv.x / scale;
	  dmv.y = x*mv.y / scale;

	  c[0].x=m.x+dmv.x; c[0].y=m.y+dmv.y;
	  c[1].x=m.x-dmv.x; c[2].y=m.y-dmv.y;

	  cout << "m  " << m.x    << "," << m.y    << endl;
	  cout << "c0 " << c[0].x << "," << c[0].y << endl;
	  cout << "c1 " << c[1].x << "," << c[1].y << endl;


	  for( int cn=0;cn<2;cn++ ){
	    int sum=0;
	    for( int k=0;k<n;k++ ){
	      if( dis( c[cn],p[k] ) < 1.0 + EPS )
		sum++;
	    }
	    if( mx<sum ) mx=sum;
	  }
	}
      }
    }
    cout << mx << endl;
  }

  return 0;
}

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

using namespace std;
typedef double D;
const D EPS= (1e-10);

#define X real()
#define Y imag()
#define SZ size()
#define PB(x) push_back(b)
#define PREV(x,i)  (x[(i+x.SZ-1)%x.SZ])
#define NEXT(x,i) (x[(i+1)%x.SZ])

template<class T> using CR=const  T&;
using P=complex<D>;
using G=vector<P>;

int sgn(D a,D b=0){
    if(a>b+EPS)return 1;
    if(a<b-EPS)return -1;
    return 0;
}

D dist(P a,P b){
    return abs(a-b);
}

tuple<P,P> normal_direction(P v){
    D x=v.X,y=v.Y;
    return tuple<P,P>(P(y,-x),P(-y,x));
}

P normalize(P v){
    D x=v.X,y=v.Y;
    D d=abs(v);
    return P(x/d,y/d);
}

tuple<P,P> getC(P a,P b){
    auto v=b-a;
    auto pq=normal_direction(v);
    P p,q;tie(p,q)=pq;
    D d=sqrt(1-(dist(a,b)/2)*(dist(a,b)/2));
    return tuple<P,P>(a+0.5*v+d*normalize(p),a+0.5*v+d*normalize(q));
}

int cnt(G& g,P c){
    int n=g.size(),res=0;
    for(int i=0;i<n;i++)
	if(sgn(dist(c,g[i]),1.0)<=0){
	    res++;
	}
    return res;
}

int main(){
    for(int n;cin>>n,n;){
	G g(n);
	for(int i=0;i<n;i++){\
	    D x,y;
	    cin>>x>>y;
	    g[i]=P(x,y);
	}
	int res=1;
	for(int i=0;i<n;i++)
	    for(int j=i+1;j<n;j++)
		if(sgn(dist(g[i],g[j]),2.0)<=0){
		    auto pq=getC(g[i],g[j]);
		    P p,q;tie(p,q)=pq;
		    res=max({res,cnt(g,p),cnt(g,q)});
		}
	cout<<res<<endl;
    }
}

C++: #include <stdio.h>
#include <algorithm>
#include <iostream>
#include <stack>
#include <vector>
#include <math.h>
#include <queue>

using namespace std;

struct Point{
	double x,y;
};

void func(int N){

	Point point[N];
	double x_min = 100,y_min = 100,x_max = -1,y_max = -1;

	for(int i = 0; i < N; i++){
		scanf("%lf %lf",&point[i].x,&point[i].y);
		x_min = min(point[i].x,x_min);
		y_min = min(point[i].y,y_min);
		x_max = max(point[i].x,x_max);
		y_max = max(point[i].y,y_max);
	}

	if(x_min > 1)x_min -= 1;
	if(y_min > 1)y_min -= 1;
	if(x_max < 9)x_max += 1;
	if(y_max < 9)y_max += 1;

	int count,max_count = 0;

	for(double y = y_min; y <= y_max; y += 0.015){
		for(double x = x_min; x <= x_max; x += 0.015){
			count = 0;
			for(int i = 0; i < N; i++){
				if(sqrt((x-point[i].x)*(x-point[i].x)+(y-point[i].y)*(y-point[i].y)) <= 1.0)count++;
			}
			max_count = max(max_count,count);
		}
	}

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


int main(){

	int N;

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

		func(N);
	}


    return 0;
}

C++: #include <iostream>
#include <algorithm>
#include <cmath>
#include <vector>

using namespace std;

#define Sq(x) ((x)*(x))

#define Foreach(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();i++)

const double PI = M_PI;
const double EPS = 0.0001;

struct Point{
  double x,y;

  Point(){}
  Point(double x, double y) : x(x), y(y) {}
};

Point operator+(const Point &a, const Point &b){
  return Point(a.x + b.x, a.y + b.y);
}

Point operator-(const Point &a, const Point &b){
  return Point(a.x - b.x, a.y - b.y);
}

Point operator*(const Point &a, const double k){
  return Point(a.x * k, a.y * k);
}

Point operator*(const Point &a,const Point &b){
  return Point(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);
}

Point operator/(const Point &p, const double k){
  return Point(p.x / k, p.y / k);
}

double dot(const Point &a, const Point &b){
  return a.x * b.x + a.y * b.y;
}

double cross(const Point &a, const Point &b){
  return a.x * b.y - a.y * b.x;
}

double norm(const Point &p){
  return dot(p, p);
}

double abs(const Point &p){
  return sqrt(norm(p));
}

typedef Point Vector;

struct Circle : Point{
  double r;
  Circle(Point p = Point(0,0), double r = 0) : Point(p), r(r) {}
};

double distPP(const Point &a, const Point &b) {
  return sqrt(Sq(a.x-b.x) + Sq(a.y-b.y));
}

bool isPointInCircle(const Point &p, const Circle &c) {
  return abs(p - c) <= c.r + EPS;
}

Point rotate90(const Point &p) {
  return Point(-p.y, p.x);
}

/* テ・ツソツ?ィツヲツ?」ツ?ェテ」ツ?凝」ツ?」テ」ツ??*/
vector<Point> getCrossPointCC_R1(const Circle &c1, const Circle &c2) {
  vector<Point> ret;

  double dist = abs( c1 - c2 );
  if( abs( c1 - c2 ) < EPS || abs( c1 + c2 ) > dist ) return vector<Point>();
  
  Point cent = ( c1 + c2 ) / 2;
  
  Vector v1 = c1 - cent;
  Vector v2 = rotate90(v1);
  ret.push_back(cent + v1);

  v1 = c2 - cent;
  v2 = rotate90(v2);
  ret.push_back(cent + v2);

  return ret;
}

/*
テ」ツ??」ツ?セテ」ツ?湘・ツ仰陛」ツ?凝」ツ?ェテ」ツ??
Point rotateVector(const Vector &p, double th) {
  return Point(p.x * cos(th) - p.y*sin(th), p.x * sin(th) + p.y * cos(th));
}
*/

int main() {
  int n;
  Point p[300];
  
  while(cin >> n && n) {
    for(int i=0; i<n; i++) {
      cin >> p[i].x >> p[i].y;
    }
    
    int ans = 1;
    for(int i=0; i<n-1; i++) {
      for(int j=i+1; j<n; j++) {
	if(distPP(p[i], p[j]) + EPS > 2.0) continue;
	
	Point cent((p[i]+p[j])/2);
	double a = distPP(p[i], cent);
	double b = sqrt(1.0 - a*a);
	
	Vector v1 = p[i] - cent;
	Vector v2 = rotate90(v1);
	Point center = cent + ((v2/abs(v2)) * b);
	Circle c1(center, 1.0);
	int cnt = 0;
	for(int k=0; k<n; k++) {
	  cnt += isPointInCircle(p[k], c1);
	}
	ans = max(ans, cnt);
	
	center = cent + ((v2/abs(v2)) * -b);
	Circle c2(center, 1.0);
	cnt = 0;
	for(int k=0; k<n; k++) {
	  cnt += isPointInCircle(p[k], c2);
	}
	ans = max(ans, cnt);
      }
    }
    cout << ans << endl;
  }
  
  return 0;
}