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 <math.h>
#include <algorithm>
#include <iostream>
#include <map>
#include <sstream>
#include <string>
#include <vector>
using namespace std;


int circle();


typedef pair<double, double> pt;


#define x first
#define y second


pt points[305];
int n;


double dist(pt a, pt b) {
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}


pt get_circle(pt a, pt b) {
    pt mid_point = make_pair((a.x + b.x) / 2, (a.y + b.y) / 2);
    double angle = atan2(a.x - b.x, b.y - a.y);
    double d = sqrt(1 - (dist(a, mid_point) * dist(b, mid_point)));
    return make_pair(mid_point.x + d * cos(angle), mid_point.y + d * sin(angle));
}


int get_points_in_circle(pt c) {
    int num_points = 0;
    for (int i = 0; i < n; ++i) {
        if (dist(points[i], c) < 1.0 + 1e-8) {
            ++num_points;
        }
    }
    return num_points;
}


int circle() {
    int num_points = 0;

    for (int i = 0; i < n - 1; ++i) {
        for (int j = i + 1; j < n; ++j) {
            if (dist(points[i], points[j]) > 2) continue;
            pt c1 = get_circle(points[i], points[j]);
            pt c2 = get_circle(points[j], points[i]);
            num_points = max(num_points, max(get_points_in_circle(c1), get_points_in_circle(c2)));
        }
    }

    return num_points;
}


int main(int argc, char **argv) {
    ios_base::sync_with_stdio(false); cin.tie(NULL);

    while (cin >> n && n > 0) {
        for (int i = 0; i < n; ++i) {
            cin >> points[i].x >> points[i].y;
        }

        if (n == 1) {
            cout << "1" << endl;
        } else {
            cout << circle() << endl;
        }
    }
    return 0;
}

C++: #include <iostream>
#include <sstream>
#include <string>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <map>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>

using namespace std;

#define FOR(i,k,n) for(int i=(k); i<(int)n; ++i)
#define REP(i,n) FOR(i,0,n)
#define FORIT(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();++i)
const double EPS = 1e-8;
typedef pair<double,double> P;

inline double dist(P p,P q){
  return ( (p.first-q.first)*(p.first-q.first) + (p.second-q.second)*(p.second-q.second));
}

int main(){
  int N;
  double sign[2] = {+1, -1};
  while(cin>>N,N){
    vector<P> p(N);
    REP(i,N)cin>>p[i].first>>p[i].second;
    int ans = 1;
    REP(i,N){
      FOR(j,i+1,N){
        double dst = dist(p[i],p[j]);
        if(dst > 4.0) continue;
        P v((p[j].first-p[i].first)/dst, (p[j].second-p[i].second)/dst);
        P m((p[j].first+p[i].first)/2.0, (p[j].second+p[i].second)/2.0);
        double lth = sqrt(1.0-dst*dst/4.0);
        REP(k,2){
          double cx = m.first + sign[k] * v.first * lth;
          double cy = m.second - sign[k] * v.second * lth;
          P center(cx,cy);
          int count = 0;
          REP(l,N){
            if(l == i || l == j || dist(center,p[l]) < 1){
              count++;
            }
          }
          ans = max(ans, count);
        }
      }
    }
    cout<<ans<<endl;
  }
  return 0;
}

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

typedef double Num;
typedef complex<Num> Pt2;

struct Circle2 {
  Pt2 p; Num r;
  Circle2(Pt2 p, Num r) : p(p), r(r) {}
};

const double eps = 1e-10;
const double pi  = 3.14159265358979323846;

Pt2  ortho(Pt2 v)        { return v*Pt2(0,1); }
Pt2  unit (Pt2 v)        { return v/abs(v);   }

vector<Circle2> make_circle(Pt2 a, Pt2 b, Num r) {
	Num d = abs(b-a)/2.0;
	d = sqrt(r*r-d*d); // assert(r*r-d*d >= 0)
	Pt2 p = (b+a)/2.0;
	Pt2 v = unit(ortho(b-a)) * d;

	vector<Circle2> ret;
	ret.push_back(Circle2(p+v,r));
	ret.push_back(Circle2(p-v,r));
	return ret;
}

bool into_circle(Circle2 p, Pt2 a) {
	return abs(a-p.p)*abs(a-p.p) < p.r*p.r+eps;
}

int main()
{
	int N;
	while(cin>>N, N)
	{
		Pt2 ps[300];
		for(int i=0; i<N; i++) {
			Num x,y;
			cin>>x>>y;
			ps[i]=Pt2(x,y);
		}

		int ans=1;
		for(int i=  0; i<N; i++)
		for(int j=i+1; j<N; j++) {
			if(abs(ps[i]-ps[j]) > 2.0) continue;

			vector<Circle2> cs = make_circle(ps[i], ps[j], 1.0);
			for(int k=0; k<2; k++){
				int cnt=0;
				for(int l=0; l<N; l++) {
					if(into_circle(cs[k], ps[l])) cnt++;
				}
				ans=max(ans,cnt);
			}
		}
		cout << ans << endl;
	}
}

C++: #include <bits/stdc++.h>
#define PB push_back
#define MP make_pair
#define REP(i,n) for (int i=0;i<(n);i++)
#define FOR(i,a,b) for(int i=(a);i<(b);i++)
#define ALL(a) (a).begin(),(a).end()
int INF=1e9;
int MOD=1000000007;
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
double EPS = 1e-5;
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 * (P p){
		return P(x*p.x,y*p.y);
	}
};
double dist(P a,P b){
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int cnt(P pos,vector<P> dat){
	int res=0;
	REP(i,dat.size()){
		if(dist(pos,dat[i])<=1.0+EPS)res++;
	}
	return res;
}
int main(){
	int n;
	double x,y;
	while(scanf("%d",&n)&&n){
		vector<P> dat;
		REP(i,n){
			cin>>x>>y;
			dat.PB(P(x,y));
		}
		int ans=1;
		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				P a=dat[i],b=dat[j];
				P mid=P((a.x+b.x)/2.0,(a.y+b.y)/2.0);
				double dx=a.y-b.y,dy=b.x-a.x;
				double d=dist(a,b);
				if(d>2)continue;
				double dd=sqrt(dx*dx+dy*dy);
				double ddd=sqrt(1-d*d/4);
				dx=dx*ddd/dd;dy=dy*ddd/dd;
				ans=max(ans,cnt(P(mid.x+dx,mid.y+dy),dat));
				ans=max(ans,cnt(P(mid.x-dx,mid.y-dy),dat));
			}
		}
		printf("%d\n",ans);
	}
	return 0;
}