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 (1e8)
#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 rotate(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 n;

pair<Circle,Circle> makeCirclebySegment(Point a,Point b,int r){
  Point m=a+(b-a)/2.0;
  double d=sqrt(r*r-abs(m-b)*abs(m-b));
  Vector v=(rotate(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()
{
  double x,y;

  while(1){
    cin>>n;
    if(n==0)break;
    vp.clear();

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

C++: #include<bits/stdc++.h>
#define eps 1e-9
using namespace std;
typedef complex<double> point;

pair<point, point> norm(point p){
  double unit = 1 / sqrt(p.real()*p.real() + p.imag()*p.imag());
  return pair<point, point>(p*point(0, 1)/unit,
                            p*point(0, -1)/unit);
}

vector<point> calcC(point a, point b){
  point p = (b - a)/2.0;
  pair<point, point> n = norm(p);
  double unit = sqrt(1 - abs(p)*abs(p));
  n.first  *= unit;
  n.second *= unit;
  vector<point> res(2);
  res[0] = a + p + n.first;
  res[1] = a + p + n.second;
  return res;
}
int main(int argc, char *argv[]){
  int n;
  while(cin >> n, n){;
    point p[n];
    for (int i = 0; i < n; i++) {
      double real, imag;
      cin >> real >> imag;
      p[i] = complex<double>(real, imag);
    }
    int ans = 0, cnt = 0;
    for (int i = 0; i < n; i++) {
      for (int j = i + 1; j < n; j++) {
        point a = p[i], b = p[j];
        if(abs(b - a) > 2 - eps)continue;
        vector<point> center = calcC(a, b);
        for (int c = 0; c < 2; c++) {
          cnt = 0;
          // std::cout << center[c].real() << " " <<
          //   center[c].imag()<< std::endl;
          for (int k = 0; k < n; k++) {
            if(abs(center[c] - p[k]) <= 1 + eps)cnt++;
          }
          ans = max(ans, cnt);
        }
      }
    }
    std::cout << ans << std::endl;
  }
  return 0;
}

C++: #include <stdio.h>
#include <string.h>
#include <algorithm>
#include <iostream>
#include <math.h>
#include <assert.h>
#include <vector>
#include <complex>

using namespace std;
typedef long long ll;
typedef complex<double> Point;
static const double EPS = 1e-9;
static const double PI = acos(-1.0);

#define REP(i, n) for (int i = 0; i < (int)(n); i++)
#define FOR(i, s, n) for (int i = (s); i < (int)(n); i++)
#define FOREQ(i, s, n) for (int i = (s); i <= (int)(n); i++)
#define FORIT(it, c) for (__typeof((c).begin())it = (c).begin(); it != (c).end(); it++)
#define DEC(i, s) for (int i = (s); i >= 0; i--)

#define SIZE(v) (int)((v).size())
#define MEMSET(v, h) memset((v), h, sizeof(v))
#define FIND(m, w) ((m).find(w) != (m).end())

int n;
int ans;
double R;
Point points[3000];
pair<double, int> radian[30000];
bool used[3000];

int calc(int s) {
  int ret = 1;
  int m = 0;
  REP(i, n) {
    double d = abs(points[s] - points[i]) / 2.0;
    if (i == s || d * 2.0 > R * 2) { continue; }
    Point center = (points[s] + points[i]) / 2.0;
    Point vect = (points[i] - points[s]) * Point(0, -1);
    vect /= abs(vect);
    double l = sqrt(R * R - d * d);
    Point c1 = center + vect * l;
    Point c2 = center - vect * l;
    double arg1 = arg(c1 - points[s]) - EPS;
    double arg2 = arg(c2 - points[s]) + EPS;
    radian[m++] = pair<double, int>(arg1 + 0 * PI, (i + 1));
    radian[m++] = pair<double, int>(arg2 + 0 * PI, -(i + 1));
  }
  if (m / 2 < ans - 1) { return 0; }
  sort(radian, radian + m);
  MEMSET(used, false);
  int lsum = 1;
  REP(iter, 2) {
    REP(i, m) {
      int index = abs(radian[i].second) - 1;
      bool start = radian[i].second > 0 ? true : false;
      if (used[index] && !start) {
        used[index] = false;
        lsum--;
      } else if (!used[index] && start) {
        used[index] = true;
        lsum++;
      }
      ret = max(ret, lsum);
    }
  }
  return ret;
}

int main() {
  R = 1.0;
  while (scanf("%d %lf", &n), n) {
    REP(i, n) {
      double x, y;
      scanf("%lf %lf", &x, &y);
      points[i] = Point(x, y);
    }
    ans = 0;
    REP(i, n) {
      ans = max(ans, calc(i));
    }
    printf("It is possible to cover %d points.\n", ans);
  }
}

C++: #include<iostream>
#include<complex>
#include<vector>
using namespace std;
typedef complex<double> P;
int n;
vector<P> p;
int ret(){
	int cnt=0;
	for(int i=0;i<n;i++){
		for(int j=i+1;j<n;j++){
			if(abs(p[i]-p[j])>2.0)
				continue;
			P c=(p[i]+p[j])*0.5;
			double d=sqrt(1.0-abs(p[i]-c)*abs(p[i]-c));
			P l1=(p[i]-p[j])*P(0,1);
			P l2=(p[i]-p[j])*P(0,-1);
			l1=l1*d/abs(l1)+c;l2=l2*d/abs(l2)+c;
			int x=2,y=2;
			for(int k=0;k<n;k++){
				if(k==i||k==j)
					continue;
				if(abs(l1-p[k])<=1.0){
					x++;
				}
				if(abs(l2-p[k])<=1.0){
					y++;
				}
				cnt=max(cnt,max(x,y));
			}
		}
	}
	return cnt;
	
}
int main(){
	while(cin>>n,n){
		p.resize(n);
		for(int i=0;i<n;i++){
			double x,y;
			cin>>x>>y;
			p[i]=P(x,y);
		}
		cout<<ret()<<endl;
	}
	return 0;
}