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 <stdio.h>
#include <algorithm>
#include <iostream>
#include <stack>
#include <vector>
#include <math.h>
#include <queue>

using namespace std;

struct Point{
    double x,y;
};

void func(int N){

    Point point[N];
    double x_min = 100,y_min = 100,x_max = -1,y_max = -1;

    for(int i = 0; i < N; i++){
        scanf("%lf %lf",&point[i].x,&point[i].y);
        x_min = min(point[i].x,x_min);
        y_min = min(point[i].y,y_min);
        x_max = max(point[i].x,x_max);
        y_max = max(point[i].y,y_max);
    }

    if(x_min > 1)x_min -= 0.05;
    if(y_min > 1)y_min -= 0.05;
    if(x_max < 9)x_max += 0.05;
    if(y_max < 9)y_max += 0.05;

    int count,max_count = 0;

    for(double y = y_min; y <= y_max; y += 0.01){
        for(double x = x_min; x <= x_max; x += 0.01){
            count = 0;
            for(int i = 0; i < N; i++){
                if(sqrt((x-point[i].x)*(x-point[i].x)+(y-point[i].y)*(y-point[i].y)) <= 1.0)count++;
            }
            max_count = max(max_count,count);
        }
    }

    printf("%d\n",max_count);
}


int main(){

    int N;

    while(true){
        scanf("%d",&N);
        if(N == 0)break;

        func(N);
    }


    return 0;
}

C++: #include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <map>
#include <utility>
#include <set>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>
#include <sstream>
#include <complex>
#include <stack>
#include <queue>
#include <cstring>
#include <sstream>
#include <cassert>
#include <list>
using namespace std;
static const double EPS = 1e-5;
typedef long long ll;
typedef pair<int,int> PI;
typedef vector<int> vi;
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#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()
#define mp(a,b) make_pair(a,b)
#define pb(a) push_back(a)
#define F first
#define S second


int dx[]={0,1,0,-1};
int dy[]={1,0,-1,0};


main(){
  int n;
  while(cin>>n,n){
    double x[n],y[n];

    rep(i,n)cin>>x[i]>>y[i];

    int ans=0;

    rep(i,n){
      rep(j,i){
	double sx,sy;
	double tx=x[i]-x[j],ty=y[i]-y[j];
	double len=sqrt(tx*tx+ty*ty);
	if(len>2+EPS)continue;
	len=sqrt(1-len*len/4);
	double theta=atan2(ty,tx)+M_PI/2;
	int tans=0;
	sx=(x[i]+x[j])/2+cos(theta)*len;
	sy=(y[i]+y[j])/2+sin(theta)*len;
	rep(k,n){
	  tx=sx-x[k];
	  ty=sy-y[k];
	  if(tx*tx+ty*ty<=1+EPS)++tans;
	}
	ans=max(ans,tans);
	tans=0;
	sx=(x[i]+x[j])/2-cos(theta)*len;
	sy=(y[i]+y[j])/2-sin(theta)*len;
	rep(k,n){
	  tx=sx-x[k];
	  ty=sy-y[k];
	  if(tx*tx+ty*ty<=1+EPS)++tans;
	}
	ans=max(ans,tans);
      }
    }
    cout<<ans<<endl;
  }
}

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

using namespace std;

typedef complex<double> Point;

bool eqv(double a, double b) {
    return abs(a-b) < 1e-5;
}

pair<double,double> key_point(double y, const Point &p) {
    double len = sqrt(1.0 - pow(y-p.imag(), 2));
    return make_pair(p.real()+len, p.real()-len);
}

bool inCircle(const Point &p, const Point &center, double radius) {
    double dist2 = norm(p-center);
    double rad2 = radius*radius;
    return dist2 < radius*radius || eqv(dist2, rad2);
}

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]);
    const double hlen = abs(v)/2;

    if(hlen*hlen >= 1) return 0;
    const double plen = sqrt(1 - hlen*hlen);
    Point perp = v*Point(0, 1);
    perp /= abs(perp);
    const Point mid = points[b] + v/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));
            }
        }
        cout << ans << endl;
    }
    return 0;
}

C++: #include "bits/stdc++.h"
using namespace std;
typedef long long LL;
struct Point {
	double x;
	double 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);
	}
	double abs() {
		return sqrt(norm());
	}
	double norm() {
		return x*x + y*y;
	}
};
typedef Point Vector;
struct Circle {
	Point c;
	double r;
	Circle(Point _c = Point(), double _r = 0.0) :c(_c), r(_r) {}
};
double arg(Vector p) {
	return atan2(p.y, p.x);
}
Vector polar(double a, double r) {
	return Point(cos(r)*a, sin(r)*a);
}
pair<Point, Point> getCrossPoints(Circle c1, Circle c2) {
	double d = (c1.c - c2.c).abs();
	double a = acos((c1.r*c1.r + d*d - c2.r*c2.r) / (2 * c1.r*d));
	double t = arg(c2.c - c1.c);
	return make_pair(c1.c + polar(c1.r, t + a), c1.c + polar(c1.r, t - a));
}
double X[300], Y[300];
int main() {
	int N;
	while (cin >> N, N) {
		for (int i = 0; i < N; i++) {
			cin >> X[i] >> Y[i];
		}
		int ans = 0;
		for (int i = 0; i < N; i++) {
			for (int j = i + 1; j < N; j++) {
				Circle c1(Point(X[i], Y[i]), 1.0);
				Circle c2(Point(X[j], Y[j]), 1.0);
				if ((c1.c - c2.c).abs() >= 2.0001) continue;
				pair<Point, Point> c = getCrossPoints(c1, c2);
				int cnt1 = 0, cnt2 = 0;
				for (int k = 0; k < N; k++) {
					Point p(X[k], Y[k]);
					if ((p - c.first).abs() <= 1.0001) cnt1++;
					if ((p - c.second).abs() <= 1.0001) cnt2++;
				}
				ans = max(ans, cnt1);
				ans = max(ans, cnt2);
			}
		}
		cout << ans << endl;
	}
}