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<cmath>
#include<algorithm>
#include<iostream>
#include<vector>
#include<climits>
#include<cfloat>
#define curr(P, i) P[(i) % P.size()]
#define next(P, i) P[(i+1) % P.size()]
#define prev(P, i) P[(i+P.size()-1) % P.size()]
 
 
using namespace std;
 
typedef double Real;
 
Real EPS = 1e-8;
 
int sgn(Real a, Real b=0){return a<b-EPS?-1:a>b+EPS?1:0;}
Real sqr(Real a){return sqrt(max(a,(Real)0));}
 
struct Point{  
  Real add(Real a, Real b){
    if(abs(a+b) < EPS*(abs(a)+abs(b)))return 0;
    return a+b;
  }
 
  Real x, y;
  Point(){}
  Point(Real x,Real y) : x(x) , y(y){}
 
  Point operator + (Point p){return Point(add(x,p.x), add(y,p.y));}
  Point operator - (Point p){return Point(add(x,-p.x), add(y,-p.y));}
  Point operator * (Real d){return Point(x*d,y*d);}
  Point operator / (Real d){return Point(x/d,y/d);}
  bool operator == (Point p){return !sgn(dist(p));}
  bool operator < (Point p){return (p.x!=x)?p.x<x:p.y<y;}
  Real norm(){return sqr(x*x+y*y);}
  Real dist(Point a){return (*this-a).norm();}
  Real dot(Point a){return x*a.x+y*a.y;}
  Real cross(Point a){return x*a.y-y*a.x;}
};
 
struct Circle{
  Point p;
  Real r;
  Circle(){}
  Circle(Point p, Real r) : p(p) , r(r){}
   
  vector<Point> intersectionPoints(Circle b){
    vector<Point> res;
 
    Point tmp=p-b.p;
    if(tmp.norm()>pow(r+b.r,2))return res;
 
    Real L = p.dist(b.p);
    Real C = atan2(b.p.y-p.y,b.p.x-p.x);
    Real a = acos((L*L+r*r-b.r*b.r)/(2*L*r));
 
    res.push_back(Point(p.x+cos(C+a), p.y+sin(C+a))*2);
    res.push_back(Point(p.x+cos(C-a), p.y+sin(C-a))*2);
    return res;
  }
};
 
int main(void){
  int n;
  Point a;
  vector<Point>P;
  vector<Circle>C;
 
  while(cin >> n,n){
    P.clear();
    C.clear();
    for(int i=0;i<n;i++){
      cin >> a.x >> a.y;
      P.push_back(a);
    }
 
    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
    if(sgn(P[i].dist(P[j]),2.0)>=0)continue;
    vector<Point>res=Circle(P[i],1.0).intersectionPoints(Circle(P[j],1.0));
    for(int k=0;k<res.size();k++)
      C.push_back(Circle(res[k],1));
      }
    }
 
    int ans=1;
    for(int i=0;i<C.size();i++){
      int cnt=0;
      for(int j=0;j<P.size();j++)
    if(sgn(C[i].p.dist(P[j]),1.0)<=0)cnt++;
 
      ans=max(ans,cnt);
    }
    cout << ans << endl;
  }
 
  return 0;
}

C++: #include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;

#define MAX_N 310

double x[MAX_N], y[MAX_N];
int n, maxn;

int main() {
	while (true) {
		cin >> n;
		if (n == 0) { break; }
		for (int i = 0; i < n; i++) {
			cin >> x[i] >> y[i];
		}
		maxn = 1;
		for (double i = -10; i <= 10; i += 0.02) {
			for (double j = -10; j <= 10; j += 0.02) {
				double X = i, Y = j;
				int sum = 0;
				for (int k = 0; k < n; k++) {
					if (((x[k] - X)*(x[k] - X) + (y[k] - Y)*(y[k] - Y)) <= 1) { sum++; }
				}
				maxn = max(maxn, sum);
			}
		}
		cout << maxn << 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(dsq - 1.);
        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 <algorithm>
#include <cmath>
#define COS45 0.707107
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 a) {
    return Point(x*a, y*a);
  }
  Point operator /(double a) {
    return Point(x/a, y/a);
  }
  bool operator ==(const Point &p) const {
    return x == p.x && y == p.y;
  }
  bool operator <(const Point &p) const {
    double l1 = x+y;
    double l2 = p.x+p.y;
    if(l1 != l2) return l1 < l2;
    if(y != p.y) return y < p.y;
    return x < p.x;
  }
};

typedef Point Vector;

double norm(Point p1, Point p2) {
  return (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y);
}

int n;
Point p[300];

int calc(Point p1, Point p2) {
  Point m = Point((p1.x+p2.x)/2.0, (p1.y+p2.y)/2.0);
  Point o;
  Vector moto = Vector(p2.y-p1.y, p1.x-p2.x);
  double l = sqrt(4.0/norm(p1,p2)-1)/2.0;
  int cnt = 0;
  int newCnt;
  for(int t = 0; t < 2; t++) {
    newCnt = 0;
    if(t == 1) {
      moto = moto * (-1);
    }
    o = m + moto * l;

    for(int i = 0; i < n; i++) {
      if(p[i].x > o.x + COS45 && p[i].y > o.y + COS45) break;
      if(newCnt + n-i <= cnt) break;
      if(p[i] == p1 || p[i] == p2) continue;
      if((p[i].x-o.x)*(p[i].x-o.x) + (p[i].y-o.y)*(p[i].y-o.y) <= 1.0) {
	newCnt++;
      }
    }
    if(newCnt > cnt) cnt = newCnt;
  }
  return cnt + 2;
}

main() {
  while(1) {
    int ans;
    cin >> n;
    if(n == 0) break;
    for(int i = 0; i < n; i++) {
      cin >> p[i].x >> p[i].y;
    }
    sort(p, &p[n]);
    ans = -1;
    for(int i = 0; i < n; i++) {
      for(int j = i+1; j < n; j++) {
	ans = max(ans, calc(p[i],p[j]));
      }
    }
    cout << ans << endl;
  }
}