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(a-b)/2)*(abs(a-b)/2));
  Point A=m+(PointRotation(m,a,90)-m)*d,B=m+(PointRotation(m,b,90)-m)*d;
  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)<=C.r+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++: #include <iostream>
#include <cmath>
#define REP(i, a, n) for(int i = (a); i < (n); i++)
using namespace std;

int main(void) {
  for(int N; cin >> N, N > 0;) {
    double x[300], y[300];
    REP(i, 0, N) cin >> x[i] >> y[i];

    int answer = 1;
    REP(i, 0, N) {
      REP(j, i + 1, N) {
        double tan = (y[i] - y[j]) / (x[i] - x[j]);
        double cot = (x[i] - x[j]) / (y[i] - y[j]);
        double dsq = (x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j]);
        if(dsq > 4.0) continue;
        double l = sqrt(1. - dsq / 4.);
        double sx = (x[i] + x[j]) / 2.;
        double sy = (y[i] + y[j]) / 2.;

        double cx1 = sx - tan * l;
        double cy1 = sy + cot * l;
        int cnt1 = 0;
        REP(k, 0, N) if(k == i || k == j || (x[k] - cx1) * (x[k] - cx1) + (y[k] - cy1) * (y[k] - cy1) <= 1.0) cnt1++;
        answer = max(answer, cnt1);

        double cx2 = sx + tan * l;
        double cy2 = sy - cot * l;
        int cnt2 = 0;
        REP(k, 0, N) if(k == i || k == j || (x[k] - cx2) * (x[k] - cx2) + (y[k] - cy2) * (y[k] - cy2) <= 1.0) cnt2++;
        answer = max(answer, cnt2);
      }
    }

    cout << answer << endl;
  }

  return 0;
}

C++: #include<iostream>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<queue>
#include<algorithm>
#include<cstring>
#include<string>
#include<vector>
#define mp make_pair
#define pb push_back
#define rep(i,n) for(int i = 0;i < n;i++)
using namespace std;

double dist(double x1,double y1,double x2,double y2){
    return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}

double dist(double x1,double y1){
    return dist(x1,y1,0,0);
}

void calccenter(double c[][2],double u1[],double u2[]){
    double v[2]={u2[0]-u1[0],u2[1]-u1[1]};
    rep(i,2) {
        rep(j,2){
            c[i][j]=(u1[j]+u2[j])/2.0;
        }
    }
    double lenv = dist(v[0],v[1]);
    double lent = sqrt(1-pow((lenv/2.0),2));
    rep(i,2){
        double dg = i?1.0:-1.0;
        c[i][0]+=dg*(v[1]*lent/lenv);
        c[i][1]-=dg*(v[0]*lent/lenv);
    }
}

int main(){
    int N;
    while(1){
        cin>>N;
        if(!N)break;
        double a[300][2];
        rep(i,N){
            cin>>a[i][0]>>a[i][1];
        }
        int ans = 0;
        int temp = 0;
        rep(u1,N){
            for(int u2 = u1+1;u2<N;u2++){
                double c[2][2];
                calccenter(c,a[u1],a[u2]);
                rep(k,2){
                    temp = 0;
                    rep(l,N){
                        if(l==u1||l==u2){
                            temp++;
                        }else if(dist(c[k][0],c[k][1],a[l][0],a[l][1])<=1.0){
                            temp++;
                        }
                    }
                    ans = max(ans,temp);
                }
            }
        }
        cout<<ans<<endl;
    }
}

C++: #include <bits/stdc++.h>
using namespace std;
const double EPS = 0.0000001;
struct point{
  double x, y;
  point(){
  }
  point(double x, double y): 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);
  }
};
double abs(point P){
  return sqrt(pow(P.x, 2) + pow(P.y, 2));
}
double dist(point P, point Q){
  return abs(Q - P);
}
point midpoint(point P, point Q){
  return (P + Q) / 2;
}
point rotate90(point P){
  return point(P.y, - P.x);
}
pair<point, point> equal_distance(point P, point Q, double d){
  point M = midpoint(P, Q);
  double d2 = sqrt(pow(d, 2) - pow(dist(P, Q) / 2, 2));
  point N = rotate90(Q - P) / dist(P, Q);
  return make_pair(M + N * d2, M - N * d2);
}
int main(){
  while (1){
    int N;
    cin >> N;
    if (N == 0){
      break;
    }
    vector<point> P(N);
    for (int i = 0; i < N; i++){
      cin >> P[i].x >> P[i].y;
    }
    vector<point> Q;
    for (int i = 0; i < N; i++){
      for (int j = i + 1; j < N; j++){
        if (dist(P[i], P[j]) <= 2){
          pair<point, point> tmp = equal_distance(P[i], P[j], 1);
          Q.push_back(tmp.first);
          Q.push_back(tmp.second);
        }
      }
    }
    int sz = Q.size();
    int ans = 1;
    for (int i = 0; i < sz; i++){
      int cnt = 0;
      for (int j = 0; j < N; j++){
        if (dist(Q[i], P[j]) <= 1 + EPS){
          cnt++;
        }
      }
      ans = max(ans, cnt);
    }
    cout << ans << endl;
  }
}