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 <bits/stdc++.h>
using namespace std;

typedef long long ll;
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define each(itr,c) for(__typeof(c.begin()) itr=c.begin(); itr!=c.end(); ++itr)
#define all(x) (x).begin(),(x).end()
#define mp make_pair
#define pb push_back
#define fi first
#define se second

const double eps=1e-8;
typedef pair<double,double> p;

double dist(p a, p b)
{
    return sqrt((a.fi-b.fi)*(a.fi-b.fi)+(a.se-b.se)*(a.se-b.se));
}

double size(p a)
{
    return sqrt(a.fi*a.fi+a.se*a.se);
}

int main()
{
    int n;
    while(cin >>n,n)
    {
        vector<p> v(n);
        rep(i,n) scanf(" %lf %lf", &v[i].fi, &v[i].se);

        int ans=1;
        //select 2 points
        rep(i,n)rep(j,i)
        {
            double d=dist(v[i],v[j]);
            if(d>2.0) continue;

            double r=sqrt(1.0 - (d*d/4.0));

            double X=v[j].fi-v[i].fi, Y=v[j].se-v[i].se;

            p v1=p(X,Y);
            double v1size=size(v1);
            v1.fi/=v1size; v1.se/=v1size;
            v1.fi*=d/2; v1.se*=d/2;

            p v2;
            double v2size;
            p c;

            //2 candidate for center of circle
            //first
            int ct=0;
            v2=p(-Y,X);
            v2size=size(v2);
            v2.fi/=v2size; v2.se/=v2size;
            v2.fi*=r; v2.se*=r;

            c.fi=v[i].fi+v1.fi+v2.fi; c.se=v[i].se+v1.se+v2.se;
            assert(dist(c,v[i]));
            rep(k,n)
            {
                if(k==i || k==j) ++ct;
                else
                {
                    if(dist(c,v[k])<1.0) ++ct;
                }
            }
            ans=max(ans,ct);

            //second
            ct=0;
            v2=p(Y,-X);
            v2size=size(v2);
            v2.fi/=v2size; v2.se/=v2size;
            v2.fi*=r; v2.se*=r;

            c.fi=v[i].fi+v1.fi+v2.fi; c.se=v[i].se+v1.se+v2.se;
            assert(dist(c,v[i]));
            rep(k,n)
            {
                if(k==i || k==j) ++ct;
                else
                {
                    if(dist(c,v[k])<1.0) ++ct;
                }
            }
            ans=max(ans,ct);
        }

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

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-10;
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)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=0;
		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,(a.y+b.y)/2);
				double dx=a.y-b.y,dy=b.x-a.x;
				double dd=sqrt(dx*dx+dy*dy);
				double ddd=sqrt(1-dist(a,b)/2*dist(a,b)/2);
				dx=dx/dd*ddd;dy=dy/dd*ddd;
				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;
}

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

#define rep(i, n)       rep2(i, 0, n)
#define rep2(i, m, n)   for(int i = m; i < (int)(n); ++i)
#define X real()
#define Y imag()

typedef complex<double> Pt;

const double EPS = 1e-6;

int main()
{
    int n;
    while (cin >> n, n) {
        double x, y;
        vector<Pt> p;
        rep (i, n) {
            cin >> x >> y;
            p.push_back(Pt(x, y));
        }
        
        int ans = 1;
        rep (i, n) rep2 (j, i + 1, n) {
            double d = abs(p[i] - p[j]);
            if (d > 2.0 + EPS) continue;
            
            Pt m = (p[i] + p[j]) / 2.0;
            Pt v = Pt(p[i].Y - p[j].Y, p[j].X - p[i].X);
            v /= abs(v);
            v *= sqrt(1.0 - abs(p[i] - p[j]) / 2.0);
            
            int t1 = 0, t2 = 0;
            Pt c1 = m + v, c2 = m - v;
            rep (k, n) {
                if (abs(p[k] - c1) < 1.0 + EPS) ++t1;
                if (abs(p[k] - c2) < 1.0 + EPS) ++t2;
            }
            ans = max(max(t1, t2), ans);
        }
        
        cout << ans << endl;        
    }
}