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 <iostream>
#include <algorithm>
#include <cmath>
#include <complex>

using namespace std;

typedef complex<double> Point;

const double EPS = 1e-8;

int n;
Point data[300];

int calc(Point& a, Point& b){
  int res = 0;
  double d = abs(b - a) / 2.0;
  Point ab = (b - a) / (2.0 * d);
  if(d < 1){
    Point o = ab * d;
    d = sqrt(1.0 - d * d);
    o += Point(d * ab.imag() + a.real(), -d * ab.real() + a.imag());
    for(int i=0;i<n;i++) if(abs(o - data[i]) < 1.0 + EPS) res++;
  }
  return res;
}

main(){
  while(cin >> n && n){
    for(int i=0;i<n;i++) cin >> data[i].imag() >> data[i].real();
    int ans = 1;
    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
        ans = max(ans, calc(data[i], data[j]));
        ans = max(ans, calc(data[j], data[i]));
      }
    }
    cout << ans << endl;
  }
}

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

using namespace std;

typedef double ld;

struct Point {
    ld x, y;
    Point(ld x, ld y) : 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-10;
    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 = o.x - x;
        ld dy = o.y - y;
        ld d = sqrt(dx*dx + dy*dy);

        ld slope = dy / dx;
        ld theta = atan(slope);
        ld theta_p = acos((r*r - o.r*o.r + d*d) / (2*r*d));

        return make_pair(
                Point(x + r*cos(theta - theta_p), y + r*sin(theta - theta_p)),
                Point(x + r*cos(theta + theta_p), y + r*sin(theta + theta_p)));
    }
    bool contains(ld px, ld py) const {
        return r * r + eps > (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() {
    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;
}

C++: #include <bits/stdc++.h>
using namespace std;

const double EPS = 1e-8;
using P = complex<double>;
int n;
P p[300];

int main(){
  while (cin >> n, n){
    for (int i = 0; i < n; i++){
      double x, y;
      cin >> x >> y;
      p[i] = P(x, y);
    }
    vector<P> ps;
    for (int i = 0; i < n; i++){
      for (int j = 0; j < i; j++){
        if (abs(p[i] - p[j]) >= 2.0) continue;
        P v = p[j] - p[i];
        P u = v * P(0.0, 1.0);
        u *= sqrt(1.0 - abs(v) * abs(v) / 4.0) / abs(u);
        ps.push_back(p[i] + v / 2.0 + u);
        ps.push_back(p[i] + v / 2.0 - u);
      }
    }
    int res = 1;
    for (P c : ps){
      int s = 0;
      for (int i = 0; i < n; i++){
        s += abs(c - p[i]) <= 1.0 + EPS;
      }
      res = max(res, s);
    }
    cout << res << endl;
  }
}

C++: 
#include <iostream>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <complex>
#include <cstring>
#include <cstdlib>
using namespace std;

#define REP(i,n) for(int i=0;i<(int)n;++i)
#define FOR(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();++i)
#define ALL(c) (c).begin(), (c).end()

typedef complex<double> P;
const double EPS = 1e-8;

int main() {
  int n;
  while(cin>>n,n) {
    vector<P> v;
    REP(i,n) {
      double x,y;
      cin >> x >> y;
      v.push_back(P(x,y));
    }
    int ans = 0;
    REP(i,n-1) {
      for (int j=i+1;j<n;j++) {
        P m = (v[i]+v[j])*P(0.5,0);
        //    cout << "m "<< m << endl;
        P a = v[i]-v[j];
        double dis = abs(a);
        if (1-dis*dis/4<0) continue;
        double d = sqrt(1-dis*dis/4);
        int hoge[]={1,-1};
        REP(k,2) {
          P un = (a * P(0, hoge[k])) / dis;
          P o1 = m + un*d;
          //cout << o1 << endl;
          int sum = 0;
          REP(l,n) {
            if (norm(v[l]-o1)<1+EPS) sum++;
          }
          ans = max(ans,sum);
        }
      }
    }
    cout << ans << endl;
    
  }
}