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<bits/stdc++.h>
#define f first
#define s second
#define mp make_pair
#define inf 1000000000
#define eps (1e-11)
#define equals(a,b) (fabs((a)-(b))<eps)
using namespace std;

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 k){ return Point(x*k,y*k);}
  Point operator/(double k){ return Point(x/k,y/k);}
  bool operator<(Point p)const{ return (x!=p.x ? x<p.x : y<p.y);}
  bool operator==(Point p)const{ return fabs(x-p.x)<eps && fabs(y-p.y)<eps;}

  double abs(){ return sqrt(norm());}
  double norm(){ return (x*x+y*y);}
};
typedef Point Vector;
typedef vector<Point> Polygon;

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

class Circle{
public:
  Point c;
  double r;
  Circle(Point c=Point(),double r=0.0):c(c),r(r){}
};

double norm(Vector a){ return (a.x*a.x+a.y*a.y);}
double abs(Vector a){ return sqrt(norm(a));}
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 median_Point(Point a,Point b){
  Vector v=a-b;
  return b+(v/2);
}

Point PointRotation(Point base,Point a,double r){
  Point b=a-base;
  a.x=b.x*cos((r/180)*M_PI)-b.y*sin((r/180)*M_PI);
  a.y=b.x*sin((r/180)*M_PI)+b.y*cos((r/180)*M_PI);
  a=a+base;
  return a;
}

vector<Point> vp;
int ans;
int n;

pair<Circle,Circle> makeCircle(Point a,Point b){
  Point m=median_Point(a,b);
  double d=sqrt(1-abs(m-b)*abs(m-b));
  Vector v=(PointRotation(m,a,90)-m)*(d/abs(m-a));
  Point A=m+v,B=m-v;
  return mp(Circle(A,1),Circle(B,1));
}

int PointsInCircle(Circle C)
{
  int num=0;
  for(int i=0;i<n;i++){
    if(abs(vp[i]-C.c)<=1+eps)num++;
  }
  return num;
}

int main()
{
  while(1){
    cin>>n;
    if(n==0)break;
    vp.clear();
    vp.resize(n);
    for(int i=0;i<n;i++)cin>>vp[i].x>>vp[i].y;
    ans=1;
    for(int i=0;i<n;i++){
      for(int j=0;j<n;j++){
	if(i==j)continue;
	if(2<abs(vp[i]-vp[j]))continue;
	pair<Circle,Circle> pc=makeCircle(vp[i],vp[j]);
	ans=max(ans,PointsInCircle(pc.f));
	ans=max(ans,PointsInCircle(pc.s));
      }
    }
    cout<<ans<<endl;
  }
  return 0;
}

C++: // 2004à\ID Circle and Points

#include <iostream>
#include <queue>
#include <cmath>

using namespace std;

const double sqrt2 = sqrt(2.0); 
const double eps = 1e-6;
char dx[] = {-1,1,-1,1};
char dy[] = {1,1,-1,-1};

class Square{
public:
	int bound; double x, y, size;
	Square(int bound, double x, double y, double size) : bound(bound), x(x), y(y), size(size) {}
	bool operator < (const Square &s) const { return bound < s.bound; }
};

inline double sqr(double x) { return x*x; }

int main(){
	int N;
	double x[300], y[300];
	while(cin >> N){
		if(!N) break;
		for(int i=0;i<N;i++) cin >> x[i] >> y[i];
		int ans = 1;
		priority_queue<Square> qu;
		qu.push(Square(N, 5.0, 5.0, 10.0));
		while(!qu.empty()){
			Square s = qu.top(); qu.pop();
			if(s.bound <= ans) break;
			for(int i=0;i<4;i++){
				double nx = s.x + dx[i]*s.size/2.0;
				double ny = s.y + dy[i]*s.size/2.0;
				int c = 0, upper = 0;
				for(int j=0;j<N;j++){
					double dist = sqr(x[j]-nx)+sqr(y[j]-ny);
					if(dist < 1.0+eps) c++;
					if(dist < sqr(1.0+sqrt2*s.size/2.0)+eps) upper++;
				}
				ans = max(ans, c);
				if(upper > ans) qu.push(Square(upper, nx, ny, s.size/2.0));
			}
		}
		cout << ans << endl;
	}
}

C++: #include <iostream>
#include <string>
#include <sstream>
#include <cmath>
#include <vector>

using namespace std;

typedef long double ld;

struct Point {
    ld x, y;
    Point(ld x=0, ld y=0) : x(x), y(y) {}
    string to_s() {
        ostringstream os;
        os << x << ',' << y;
        return os.str();
    }
    Point operator+(const Point& o) const {
        return Point(x+o.x, y+o.y);
    }
    Point operator/(const int d) const {
        return Point(x/d, y/d);
    }
};

struct Circle {
    ld x, y;
    ld r;
    Circle(ld x=0, ld y=0, ld r=1) : x(x), y(y), r(r) {}
    const static ld eps = 1e-11;
    bool isIntersect(const Circle& o) const {
        ld dx = o.x - x;
        ld dy = o.y - y;
        return (r + o.r) * (r + o.r) > (dx*dx + dy*dy);
    }
    pair<Point, Point> intersectionPoint(const Circle& o) const {
        ld dx = x - o.x;
        ld dy = y - o.y;
        ld a = (dx*dx + dy*dy + r*r - o.r*o.r) / 2.0;
        ld d = sqrt((dx*dx+dy*dy)*r*r-a*a);
        Point p, q;
        p.x = (a*dx + dy * d) / (dx*dx+dy*dy) + o.x;
        q.x = (a*dx - dy * d) / (dx*dx+dy*dy) + o.x;
        p.y = (a*dy - dx * d) / (dx*dx+dy*dy) + o.y;
        q.y = (a*dy + dx * d) / (dx*dx+dy*dy) + o.y;
        //cout << p.x << ',' << p.y << endl;
        //cout << q.x << ',' << q.y << endl;
        return make_pair(p, q);
    }
    bool contains(ld px, ld py) const {
        return r + eps >= sqrt((px - x) * (px - x) + (py - y) * (py - y));
    }
    string to_s() {
        ostringstream os;
        os << x << ',' << y << ':' << r << endl;
        return os.str();
    }
};

int countContainingPoints(const Circle& c, const vector<Circle>& ps) {
    int cnt = 0;
    for (int i = 0; i < ps.size(); i++) 
        if (c.contains(ps[i].x, ps[i].y)) 
            cnt++;
    return cnt;
}

Point midPoint(pair<Point, Point>& pp) {
    return (pp.first + pp.second) / 2;
}

int N;
int main() {
    //cout << Circle(0, 0, 5).intersectionPoint(Circle(8, 0, 5)).first.to_s() << endl;
    vector<Circle> cs;
    while (cin >> N, N) {
        cs.clear(); cs.resize(N);
        for (int i = 0; i < N; i++) cin >> cs[i].x >> cs[i].y;

        vector<Point> ips; // Intersection PointS
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                if (i == j) continue;
                if (!cs[i].isIntersect(cs[j])) continue;
                pair<Point, Point> ip = cs[i].intersectionPoint(cs[j]);
                ips.push_back(ip.first);
                ips.push_back(ip.second);
                //ips.push_back(midPoint(ip));
            }
        }

        int maxPoint = 0;
        for (int i = 0; i < ips.size(); i++) {
            maxPoint = max(maxPoint, countContainingPoints(Circle(ips[i].x, ips[i].y, 1), cs));
        }
        cout << maxPoint << endl;
    }
    return 0;
}