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 <cstdio>
#include <cmath>
#include <cstring>
#include <cassert>
#include <iostream>
#include <string>
#include <vector>
#include <queue>
#include <stack>
#include <algorithm>
#include <map>

// +, -, *, /, sqrt, >, >=, <, <=, ==
template<typename T>
class Point2 {
public:
  Point2() {}
  Point2(T x_, T y_) : x(x_), y(y_) {}
  Point2 operator * (T b) { return Point2<T>(x*b, y*b); };
  Point2 operator / (T b) { return Point2<T>(x/b, y/b); };
  T x, y;
};

class Mean_ {
public:
  template<typename T>
  Point2<T> operator() (Point2<T> a, Point2<T> b) {
    return Point2<T>((a.x + b.x) / 2.0, (a.y + b.y) / 2.0);
  }
};
Mean_ Mean;

class Distance_ {
public:
  template<typename T>
  T operator () (Point2<T> a, Point2<T> b) {
    T dx = a.x - b.x;
    T dy = a.y - b.y;
    return sqrt(dx * dx + dy * dy);
  }
};
Distance_ Distance;

class Heron_ {
public:
  template<typename T>
  T operator () (T a, T b, T c) {
    T s = (a + b + c) / 2.0;
    return sqrt(s*(s-a)*(s-b)*(s-c));
  }
};
Heron_ Heron;

class CenterCoord_ {
public:
  template<typename T>
  bool check(Point2<T> a, Point2<T> b, T r) {
    T dist = Distance(a, b);
    if( dist >= r * 2.0 ) return false;
    return true;
  }
  template<typename T>
  std::pair<Point2<T>, Point2<T> > operator () (Point2<T> a, Point2<T> b, T r) {
    T dist = Distance(a, b);
    assert( dist < 2.0 );
    auto t1 = Mean(a, b);
    auto t2 = Point2<T>((a.x - b.x) / dist, (a.y - b.y) / dist);
    T h = 2 * Heron(r, r, dist) / dist;
    auto t3 = t2 * h;
    auto c1 = Point2<T>(t1.x + t3.y, t1.y - t3.x);
    auto c2 = Point2<T>(t1.x - t3.y, t1.y + t3.x);
    return std::pair<Point2<T>, Point2<T> >(c1, c2);
  }
};
CenterCoord_ CenterCoord;

class InCircle_ {
public:
  template<typename T>
  bool operator () (Point2<T> center, T r, Point2<T> target) {
    return Distance(center, target) <= r;
  }
};
InCircle_ InCircle;

int n;
Point2<double> p[512];

void init() {
  for(int i = 0; i < (int)512; ++i) {
    p[i].x = p[i].y = 0.0;
  }
} 


class solve_ {
public:
  int operator() (void) {
    int res = 1;
    for(int i = 0; i < (int)n; ++i) {
      for(int j = 0; j < (int)i; ++j) {
        res = std::max(res, CountPoints(i, j));
      }
    }
    return res;
  }
  int CountPoints(int i, int j) {
    if( not CenterCoord.check(p[i], p[j], 1.0) ) return 0;
    auto c12 = CenterCoord(p[i], p[j], 1.0);
    auto c1 = c12.first;
    auto c2 = c12.second;
    int count1 = 0;
    int count2 = 0;
    for(int k = 0; k < (int)n; ++k) {
      if( k == i or k == j ) continue;
      if( InCircle(c1, 1.0, p[k]) ) {
        count1 += 1;
      }
      if( InCircle(c2, 1.0, p[k]) ) {
        count2 += 1;
      }
    }
    int res = std::max(count1, count2) + 2;
    return res;
  }
};
solve_ solve;

int main() {
  for(;;) {
    init();
    scanf("%d", &n);
    if( n == 0 ) break;
    for(int i = 0; i < (int)n; ++i) {
      scanf("%lf %lf", &p[i].x, &p[i].y);
    }
    printf("%d\n", solve());
  }
  return 0;
}

C++: #include <bits/stdc++.h>
using namespace std;
#define rep(i,n) for(int i=0;i<n;++i)
struct Point{double x,y;};
double dist(Point a,Point b) {return sqrt(pow(a.x-b.x,2)+pow(a.y-b.y,2));}

int main(void){
    int N;
    while(cin>>N,N){
        vector<Point> p(N);
        int maxcnt=1;
        rep(i,N)cin>>p[i].x>>p[i].y;
        rep(i,N){
            for(int j=i+1;j<N;j++){
                double d = dist(p[i],p[j]);
                if(d>2.0)continue;
                int sign[]={-1,1};
                rep(s,2){
                    Point C;
                    C.x = p[i].x + cos(atan2(p[j].y-p[i].y,p[j].x-p[i].x) + sign[s]*acos(d/2.0));
                    C.y = p[i].y + sin(atan2(p[j].y-p[i].y,p[j].x-p[i].x) + sign[s]*acos(d/2.0));
                    int cnt=0;
                    rep(k,N) cnt += (pow(C.x-p[k].x,2) + pow(C.y-p[k].y,2) <= 1.0);
                    if(maxcnt<cnt)maxcnt=cnt;
                }
            }
        }
        cout<<maxcnt<<endl;
    }

    return 0;
}

C++: #include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <complex>
#include <utility>

using namespace std;

typedef complex<double> Point;

bool inCircle(const Point &p, const Point &center, double radius) {
    double dist2 = norm(p-center);
    double rad2 = radius*radius;
    return dist2 < rad2+1e-9;
}

int check(const Point &center, const vector<Point> &points) {
    const int N = points.size();
    int cnt = 0;
    for(int i = 0; i < N; ++i) {
        if(inCircle(points[i], center, 1)) ++cnt;
    }
    return cnt;
}

int test(int a, int b, const vector<Point> &points) {
    const Point v(points[a]-points[b]);
    if(abs(v) >= 2) return 1;
    const double hlen = abs(v)/2;

    const double plen = sqrt(1 - hlen*hlen);
    Point perp = v*Point(0, 1);
    perp /= abs(perp);
    const Point mid = (points[a]+points[b]) / 2.0;
    return max(check(mid + perp*plen, points), check(mid - perp*plen, points));
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(0);

    while(true) {
        int N;
        cin >> N;
        if(!N) break;

        vector<Point> points(N);
        for(int i = 0; i < N; ++i) {
            cin >> points[i].real() >> points[i].imag();
        }
        int ans = 0;
        for(int i = 0; i < N; ++i) {
            for(int j = 0; j < i; ++j) {
                ans = max(ans, test(i, j, points));
            }
        }
        if(N == 1) ans = 1;
        cout << ans << endl;
    }
    return 0;
}