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 <cstring>
#include <cstdlib>
#include <cmath>
#include <complex>
#include <string>
#include <sstream>
#include <algorithm>
#include <numeric>
#include <vector>
#include <queue>
#include <stack>
#include <functional>
#include <iostream>
#include <map>
#include <set>
using namespace std;
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vll;
#define pu push
#define pb push_back
#define mp make_pair
#define eps 1e-9
#define INF 1000000000
#define sz(x) ((int)(x).size())
#define sq(x) ((x)*(x))
#define fi first
#define sec second
#define SORT(x) sort((x).begin(),(x).end())
#define all(x) (x).begin(),(x).end()
#define rep(i,n) for(int (i)=0;(i)<(int)(n);(i)++)
#define repn(i,a,n) for(int (i)=(a);(i)<(int)(n);(i)++)
#define EQ(a,b) (abs((a)-(b))<eps)
double add(double a,double b){
	if(abs(a+b)<eps*(abs(a)+abs(b)))return 0;
	return a+b;
}
struct P{
	double x,y;
	P() {}
	P(double x,double y) : x(x),y(y){}
	P operator + (P p){
		return P(add(x,p.x),add(y,p.y));
	}
	P operator - (P p){
		return P(add(x,-p.x),add(y,-p.y));
	}
	P operator * (double d){
		return P(x*d,y*d);
	}
	double dot(P p){
		return add(x*p.x,y*p.y);
	}
	double det(P p){
		return add(x*p.y,-y*p.x);
	}
	double norm(){
		return sqrt(x*x+y*y);
	}
	double dist2(P p){
		return sq(x-p.x)+sq(y-p.y);
	}
	P vert(){
		return P(y,-x);
	}
	void dump(){
		printf("%.12f %.12f\n",x,y);
		return;
	}
};
int N;
P pts[305];
int solve(){
	scanf("%d",&N);
	if(N==0)return 1;
	for(int i=0;i<N;i++){
		scanf("%lf %lf",&pts[i].x,&pts[i].y);
	}
	int ans = 0;
	for(int i=0;i<N;i++){
		for(int j=i+1;j<N;j++){
			P m = (pts[i]+pts[j])*0.5;
			P a = (pts[i]-pts[j]);
			double na = a.norm();
			if(na>2.0)continue;
			a = (a.vert())*(1/na);
			a = a*(sqrt(1-sq(na/2.0)));
			P b = m+a;
			P c = m-a;
			/*cout << i << ' ' << j << endl;
			m.dump();
			a.dump();
			b.dump();
			c.dump();
			cout << endl;*/
			int cnt = 0;
			for(int k=0;k<N;k++){
				//printf("%d:%.12f\n",k,pts[k].dist2(b));
				if(pts[k].dist2(b)<=1.0+eps)cnt++;
			}
			ans = max(ans,cnt);
			cnt = 0;
			for(int k=0;k<N;k++){
				//printf("%d:%.12f\n",k,pts[k].dist2(c));
				if(pts[k].dist2(c)<=1.0+eps)cnt++;
			}
			ans = max(ans,cnt);
		}
	}
	printf("%d\n",ans);
	return 0;
}
int main(){
	while(!solve()){}
	return 0;
}

C++: #include <iostream>
#include <algorithm>
#include <cmath>
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);
  int cnt = 0;
  int newCnt;
  double p1p2 = norm(p1,p2);
  if(p1p2 > 4) return 1;
  double l = sqrt(4.0/p1p2-1)/2.0;

  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(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;
  }
}

C++: #include<bits/stdc++.h>
#define range(i,a,b) for(int i = (a); i < (b); i++)
#define rep(i,b) for(int i = 0; i < (b); i++)
#define all(a) (a).begin(), (a).end()
#define show(x)  cerr << #x << " = " << (x) << endl;
//const int INF = 1e8;
using namespace std;

typedef complex<double> P;
int n;
P in[305];
map<vector<bool>, int> memo;

bool check(vector<bool> s, int p){
	rep(i,n){
		if(s[i]){
			if(abs(in[i] - in[p]) > 2){
				return false;
			}
		}
	}
	return true;
}

int dfs(vector<bool> s){
	if(memo.count(s)) return memo[s];

	int res = 0;
	rep(i,n){
		if(not s[i]){
			if(check(s, i)){
				s[i] = true;
				res = max(res, dfs(s) + 1);
				s[i] = false;
			}
		}
	}

	return memo[s] = res;
}

int main(){
	while(cin >> n,n){
		memo.clear();
		rep(i,n){
			double a,b;
			cin >> a >> b;
			in[i] = P(a,b);
		}

		cout << dfs(vector<bool>(n,0)) << endl;
	}
}

C++: #include <cstdio>
#include <complex>
#include <iostream>
#define rep(i,n) for(int i=0;i<n;i++)
#define re real()
#define im imag()
using namespace std;
int main(){
	while(true){
		int n,mx=1;
		scanf("%d",&n);
		if(n==0) break;
		complex<double> p[300];
		rep(i,n) scanf("%lf%lf",&p[i].re,&p[i].im);
		rep(i,n) rep(j,i){
			double d=abs(p[i]-p[j]);
			if(d>2) continue;
			complex<double> c=(p[i]+p[j])/2.0,dis=(p[i]-p[j])/d*sqrt(1-d*d/4);
			swap(dis.re,dis.im);
			dis.im=-dis.im;
			complex<double> cen=c+dis;
			int cnt=0;
			rep(k,n) if(abs(cen-p[k])<=1.0000000001) cnt++;
			mx=max(mx,cnt);
			cen=c-dis;
			cnt=0;
			rep(k,n) if(abs(cen-p[k])<=1.0000000001) cnt++;
			mx=max(mx,cnt);
		}
		printf("%d\n",mx);
	}
	return 0;
}

C++: #include <iostream>
#include <math.h>

using namespace std;

#define EPS (1e-6)

double x[300], y[300];

int sign[] = {+1, -1};

int main(void) {
  int n;
  while (cin >> n, n) {
    for (int i=0; i<n; i++) {
      cin >> x[i] >> y[i];
    }

    int ans = 1;
    for (int i=0; i<n; i++) {
      for (int j=i+1; j<n; j++) {
        double d = hypot(x[i]-x[j], y[i]-y[j]);
        if ( d > 2.0 + EPS ) { continue; }

        double mx = (x[i] + x[j]) / 2;
        double my = (y[i] + y[j]) / 2;
        double vx = (x[j] - x[i]) / d;
        double vy = (y[j] - y[i]) / d;
        double e = sqrt(1.0 - d*d/4.0);

        for (int s=0; s<2; s++) {
          double px = mx + sign[s]*e*vy;
          double py = my - sign[s]*e*vx;

          int cnt = 2;
          for (int k=0; k<n; k++) {
            if (k==i || k==j) { continue; }
            if ( (px-x[k])*(px-x[k]) + (py-y[k])*(py-y[k]) < 1.0 + EPS ) { cnt++; }
          }
          if (cnt > ans) { ans = cnt; }
        }
      }
    }

    cout << ans << endl;
  }
  return 0;
}