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<cfloat>
#include<cassert>
#include<cmath>
#include<vector>

using namespace std;

#define EPS (1e-10)
#define equals(a, b) (fabs((a) - (b)) < EPS )
#define dle(a, b) (equals(a, b) || a < b )
static const double PI = acos(-1);
double sc(double a){return sqrt(1.0-a*a);}

class Point{
public:
 double x, y;
 Point ( double x = 0, double y = 0): x(x), y(y){}
 Point operator + ( Point p ){ return Point(x + p.x, y + p.y); }
 Point operator - ( Point p ){ return Point(x - p.x, y - p.y); }
 Point operator * ( double a ){ return Point(x*a, y*a); }
 Point operator / ( double a ){ return Point(x/a, y/a); }
 double abs() { return sqrt(norm());}
 double norm() { return x*x + y*y; }
 bool operator < ( const Point &p ) const {
   return x != p.x ? x < p.x : y < p.y;
 }
 bool operator == ( const Point &p ) const {
   return fabs(x-p.x) < EPS && fabs(y-p.y) < EPS;
 }
};

typedef Point Vector;

class Segment{
   public:
   Point p1, p2;
   Segment(Point s = Point(), Point t = Point()): p1(s), p2(t){}
};
typedef Segment Line;

double norm( Vector a ){ return a.x*a.x + a.y*a.y; }
double abs( Vector a ){ return sqrt(norm(a)); }
double getDistance( Vector a, Vector b ){ return abs(a - b); }
double dot( Vector a, Vector b ){ return a.x*b.x + a.y*b.y; }
double cross( Vector a, Vector b ){ return a.x*b.y - a.y*b.x; }
Point rotate( Point a, double si,double co){
 return Point(a.x*co-a.y*si,a.x*si+a.y*co);
}
Point rotate( Point a, double n){
 return Point( a.x*cos(n)-a.y*sin(n) , a.x*sin(n)+a.y*cos(n) );
}

Point getCrossPointCC(Point p1,double r1,Point p2,double r2){
 double a=r2,b=r1,c=getDistance(p1,p2);
 double cos_A = (b*b+c*c-a*a) / (2.0*b*c);
 double sin_A = sqrt(1.0-cos_A*cos_A);
 Vector v=p2-p1;
 v = v * r1 / c;
 return p1+rotate(v,sin_A,cos_A);
}

Segment getSegmentCC(Point p1,double r1,Point p2,double r2,bool flg=false){
 double b=getDistance(p1,p2),c=r1+r2;
 double a=sqrt(b*b-c*c);
 double cos_A=(b==0?1:c/b),sin_A=(b==0?0:a/b);

 if(flg)sin_A=-sin_A;
 
 Vector v1 = (p2-p1),v2 = (p1-p2);
 v1 = v1 * r1 / b;
 v1 = rotate(v1,sin_A,cos_A);
 v2 = v2 * r2 / b;
 v2 = rotate(v2,sin_A,cos_A);
 return Segment(p1+v1,p2+v2);
}

Segment getSegmentCC2(Point p1,double r1,Point p2,double r2,bool flg=false){
 if(r1>r2){
   swap(p1,p2);
   swap(r1,r2);
 }
 double c=abs(r1-r2),b=getDistance(p1,p2);
 double a=sqrt( b*b - c*c );

 double cos_A=c/b,sin_A=a/b;
 if(flg)sin_A=-sin_A;

 Vector v1 = p1-p2;
 v1 = v1 * r1 / b;
 v1 = rotate( v1, sin_A,cos_A);
 Vector v2 = v1 * r2 / r1;
 return Segment(p1+v1,p2+v2);
}

double getDistanceLP(Line s,Point p){
 return abs(cross(s.p2-s.p1,p-s.p1)/abs(s.p2-s.p1));
}

Point t[500];
int n;

int calc(Point p){
  int res=0;
  for(int i=0;i<n;i++){
    double c = getDistance(p,t[i]);
    if(c < 1.0 + EPS )res++;
  }
  return res;
}

int main(){
  while(cin>>n&&n){
    int ans=1;
    for(int i=0;i<n;i++){
      cin>>t[i].x>>t[i].y;
    }
    for(int i=0;i<n-1;i++){
      for(int j=i+1;j<n;j++){

	if( getDistance(t[i],t[j]) < 2.0 +EPS){
	  Point A = getCrossPointCC(t[i],1.0,t[j],1.0);
	  Point B = getCrossPointCC(t[j],1.0,t[i],1.0);
	  
	  ans=max(ans, calc(A));
	  ans=max(ans, calc(B));

	}
      }
    }
    cout<<ans<<endl;
  }
 return 0;
}

C++: #include <cstdio>
#include <vector>
#include <cmath>

using namespace std;

const double r = 1.0;
const double PI = acos(-1.0);

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, double b) {
    return Point(a.x * b, a.y * b);
}
inline double Distance2(const Point &a, const Point &b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
// p1??¨??????p1p2???theta?????¢?????????????????£????????????????????????????????????????????????
Point Circle(const Point &p1, const Point &p2, double theta) {
    Point p12 = p2 - p1;
    Point rp(p12.x * cos(theta) - p12.y * sin(theta), p12.x * sin(theta) + p12.y * cos(theta));
    double len = sqrt(p12.x * p12.x + p12.y * p12.y);

    return (rp * (1.0 / len)) + p1;
}

int main()
{
    int n;

    while (true) {
        // Input
        scanf("%d", &n);
        if (n == 0) break;
        vector<Point> p(n);
        for (int i = 0; i < n; ++i)
            scanf("%lf %lf", &p[i].x, &p[i].y);

        // Solve
        // 2???????????????????§£???????£???¨??????
        int ans = 1;
        for (int i = 0; i < n; ++i)
            for (int j = i + 1; j < n; ++j) {
                double d = Distance2(p[i], p[j]);

                if (d <= 4.0 * r * r) {
                    double theta = (PI - acos((2.0 - d) * 0.5)) * 0.5;

                    for (int s = -1; s <= 1; s += 2) {
                        Point c = Circle(p[i], p[j], s * theta);

                        int num = 2;
                        for (int k = 0; k < n; ++k) {
                            if (k != i && k != j && Distance2(c, p[k]) <= 1.0)
                                ++num;
                        }

                        ans = max(ans, num);
                    }
                }
            }

        // Output
        printf("%d\n", ans);
    }

    return 0;
}