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<vector>
#include<map>
#include<set>
#include<algorithm>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<string>
#include<queue>
#include<complex>
#include<numeric>
#include<bitset>
#define INF 1001001001
#define EPS 0.000000001
using namespace std;
typedef vector<int> vint;
typedef vector<vint> vvint;
typedef pair<int,int> pint;

typedef complex<double> pt;

int main(){
	pt p[310];
	int i,j,ii;
	int n;
	
	while(1){
		cin >> n;
		if(n==0) break;
		for(i=0;i<n;i++){
			double x,y;
			scanf("%lf %lf",&x,&y);
			p[i] = pt(x,y);
		}
		
		int mx=1;
		for(i=0;i<n;i++){
			for(j=i+1;j<n;j++){
				if(norm(p[i]-p[j]) > 4.0) continue;
				double l = abs(p[i]-p[j]);
				double k = sqrt(1.0-(l*l/4.0));
				pt m = (p[i]+p[j])/2.0;
				pt v(-imag(p[i]-p[j]),real(p[i]-p[j]));
				double r = abs(v);
				v= v*k/r;
				int cnt;
				
				p[n] = m+v;
				cnt=0;
				for(ii=0;ii<n;ii++){
					if(norm(p[ii]-p[n])<1.00001) cnt++;
				}
				mx = max(mx,cnt);
				
				p[n]=m-v;
				cnt=0;
				for(ii=0;ii<n;ii++){
					if(norm(p[ii]-p[n])<1.00001) cnt++;
				}
				mx = max(mx,cnt);
			}
		}
		cout << mx << endl;
	}
}

C++: #include<bits/stdc++.h>
#define eps 1e-9
using namespace std;
typedef complex<double> point;

pair<point, point> norm(point p){
  return pair<point, point>(p*point(0, 1)/abs(p),
                            p*point(0, -1)/abs(p));
}

vector<point> calcC(point a, point b){
  point p = (b - a)/2.0;
  pair<point, point> n = norm(p);
  double unit = sqrt(1 - abs(p)*abs(p));
  n.first  =  n.first*unit;
  n.second = n.second*unit;
  vector<point> res(2);
  res[0] = a + p + n.first;
  res[1] = a + p + n.second;
  return res;
}
int main(int argc, char *argv[]){
  int n;
  while(cin >> n, n){;
    point p[n];
    for (int i = 0; i < n; i++) {
      double real, imag;
      cin >> real >> imag;
      p[i] = complex<double>(real, imag);
    }
    int ans = 1, cnt = 0;
    for (int i = 0; i < n; i++) {
      for (int j = i + 1; j < n; j++) {
        point a = p[i], b = p[j];
        if(abs(b - a) > 2 - eps)continue;
        vector<point> center = calcC(a, b);
        for (int c = 0; c < 2; c++) {
          cnt = 0;
          for (int k = 0; k < n; k++) {
            if(abs(center[c] - p[k]) <= 1 + eps)cnt++;
          }
          ans = max(ans, cnt);
        }
      }
    }
    std::cout << ans << std::endl;
  }
  return 0;
}

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

#define rep(i,N) for(int i=0;i<(int)(N);i++)
#define REP(i,k,N) for(int i=k;i<(int)(N);i++)

#define vi vector<int>
#define pb push_back

#define pdd pair<double,double>
#define fi first
#define se second
const double  EPS = 1e-8;

double dist(double x,double y){
    return x*x+y*y;
};

int main(){
//    ifstream in("d_in.txt");
//    cin.rdbuf(in.rdbuf());
    int n;
    while(cin >> n,n){
        int ma = 0;
        pdd p[330];
        rep(i,n) cin >> p[i].fi >> p[i].se;
        rep(i,n-1){
            REP(j,i+1,n){
                int cnt1 = 0,cnt2=0;
                double x1,y1,x2,y2,x3,y3;
                tie(x1,y1) = p[i];
                tie(x2,y2) = p[j];
                double cx1,cy1,cx2,cy2;
                if(dist(x1-x2,y1-y2)  > 4+EPS  ) continue;

                //はりまコーナー
//                double ax,ay,bx1,by1,bx2,by2,mx,my,dx,dy;
//                ax=x1-x2;ay=y1-y2;
//                mx=(x1+x2)/2;my=(y1+y2)/2;
//                by1 = (-0.5)*pow(ax,4)+pow(ax,2)+(-0.5)*pow(ax,2)*pow(ay,2)
//                        /dist(ax,ay);
//                by1 = sqrt(abs(by1));
//                by2 = -1*by1;
//                cy1= my+by1;
//                cy2= my+by2;
//                bx1= -1*(ay*by1)/ax;
//                bx2= -1*(ay*by2)/ax;
//                cx1=bx1+mx;
//                cx2=bx2+mx;


                //あじのりコーナー
                double mx=(x1+x2)/2,my=(y1+y2)/2;

                double ld = 1-dist((x1-x2)/2,(y1-y2)/2);
                ld = sqrt(ld);
                cx1 = mx - sin(atan((y2-y1)/(x2-x1)))*ld;
                cy1 = my + cos(atan((y2-y1)/(x2-x1)))*ld;
                cx2 = mx + sin(atan((y2-y1)/(x2-x1)))*ld;
                cy2 = my - cos(atan((y2-y1)/(x2-x1)))*ld;

                for(int k = 0;k<n;k++) {
                    tie(x3, y3) = p[k];
                    if (EPS + 1> dist((cx1 - x3), (cy1 - y3))) cnt1++;
                    if (EPS + 1> dist((cx2 - x3), (cy2 - y3))) cnt2++;
                }
                ma = max({cnt1,cnt2,ma});
            }
        }
        cout<< ma << endl;
    }
    return 0;
}

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

using namespace std;

#define rep(i, n) for (int i = 0; i < (int)n; i++)
#define vec vector
#define pdd pair<double, double>
#define mp make_pair

double dis(pdd a, pdd b)
{
	pdd x = mp(a.first - b.first, a.second - b.second);
	return sqrt(x.first * x.first + x.second * x.second);
}

pdd find(pdd a, pdd b)
{
	pdd v = mp(a.first - b.first, a.second - b.second);
	double d = sqrt(v.first * v.first + v.second * v.second);
	// double d = dis(a, b);
	if (d > 2.0) return mp(1000, 1000);
	// pdd o = mp(a.first + 0.5 * v.first, a.second + 0.5 * v.second);
	pdd o = mp((a.first + b.first) / 2.0, (a.second + b.second) / 2.0);
	pdd g = mp(v.second, - v.first);
	double s = sqrt(1.0 - d * d / 4.0);
	// double g_size = sqrt(g.first * g.first + g.second * g.second);
	double g_size = d;
	pdd ret = mp(o.first + g.first / g_size * s, o.second + g.second / g_size * s);
	return ret;
}

int main(void)
{
	while(1) {
		int n;
		cin >> n;
		if (n == 0) break;
		vec<pdd> a(n);
		rep(i, n) cin >> a[i].first >> a[i].second;
		if (n == 1) {
			cout << 1 << endl;
			continue;
		}

		int MAX = 0;
		rep(i, n) rep(j, n) {
			if (i == j) continue;
			pdd c = find(a[i], a[j]);
			// cout << i << " " << j << " " << c.first << " " << c.second << endl;

			int p = 0;
			if (c.first == 1000) p = 1;
			else {
					rep(k, n) {
					// cout << dis(a[k], c) << " ";
					if (dis(a[k], c) <= 1.0001) p++;
				}
			}
			// cout << p << endl;
			MAX = max(p, MAX);
		}

		cout << MAX << endl;
	}
	return 0;
}