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 <string>
#include <vector>
#include <sstream>
#include <iostream>
#include <algorithm>
#include <map>
#include <list>
#include <set>
#include <numeric>
#include <queue>
#include <stack>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <cstring>
#include <climits>
#include <cfloat>
#include <ctime>
#include <complex>
#include <cassert>
#include <tuple>
using namespace std;
typedef long long LL;
typedef pair<int,int> P;

int N;
double x[300];
double y[300];
inline double dist2(double a, double b)
{
	return a*a+b*b;
}
void exec()
{
	cin >> N;
	if(N==0){
		exit(0);
	}
	for(int i=0;i<N;i++){
		cin >> x[i] >> y[i];
	}
	int ret=1;
	for(int i=0;i<N;i++){
		for(int j=i+1;j<N;j++){
			double dx=(x[i]-x[j])/2;
			double dy=(y[i]-y[j])/2;
			double r2=dx*dx+dy*dy;
			double R=sqrt(1/r2-1);
			double px=x[j]+dx-R*dy;
			double py=y[j]+dy+R*dx;
			double qx=x[j]+dx+R*dy;
			double qy=y[j]+dy-R*dx;
			int sp=0,sq=0;
			for(int k=0;k<N;k++){
				if(k==i||k==j){
					sp++;sq++;
					continue;
				}
				if(dist2(x[k]-px,y[k]-py)<=1){
					sp++;
				}
				if(dist2(x[k]-qx,y[k]-qy)<=1){
					sq++;
				}
			}
			ret=max(ret,sp);
			ret=max(ret,sq);
		}
	}
	cout << ret << endl;
}
int main() {
	while(1)
		exec();
}

C++: #include <bits/stdc++.h>
#define EPS (1e-10)
using namespace  std;

typedef complex<double> V;

int n;
vector<V> points;

pair<V, V> circles(const V& a, const V& b) {
  V v = a - b;
  double dist = abs(v);
  V n1 = (v * V(0, 1)) / dist;
  V n2 = (v * V(0, -1)) / dist;

  dist /= 2;
  double x = sqrt(1 - (dist * dist));

  return make_pair((b + v / 2.) + x * n1, (b + v / 2.) + x * n2);
}

int countInCircle(const V& center) {
  int count = 0;
  for (auto p : points) {
    if (abs(p - center) < 1 + EPS) count++;
  }
  return count;
}

int solve() {
  int count = 1;
  for (int i = 0; i < n - 1; ++i) {
    for (int j = i + 1; j < n; ++j) {
      if (abs(points[i] - points[j]) < 2 + EPS) {
	auto c = circles(points[i] , points[j]);
	count = max(count, countInCircle(c.first));
	count = max(count, countInCircle(c.second));
      }
    }
  }
  return count;
}

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

    cout << solve() << endl;
  }
  return 0;
}

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

#define fi first
#define se second
#define repl(i,a,b) for(int i=(int)(a);i<(int)(b);i++)
#define repr(i,n) for(int i=(int)(n-1);i>=0;i--)
#define rep(i,n) repl(i,0,n)
#define each(itr,v) for(auto itr:v)
#define pb(s) push_back(s)
#define all(x) (x).begin(),(x).end()
#define dbg(x) cout << #x" = " << x << endl
#define print(x) cout << x << endl
#define maxch(x,y) x=max(x,y)
#define minch(x,y) x=min(x,y)
#define uni(x) x.erase(unique(all(x)),x.end())
#define exist(x,y) (find(all(x),y)!=x.end())
#define bcnt(x) bitset<32>(x).count()

typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> P;
typedef pair<double, double> PD;
typedef pair<P, int> PPI;
typedef pair<int, P> PIP;
typedef pair<ll, ll> PL;
typedef pair<P, ll> PPL;
typedef set<int> S;

#define INF INT_MAX/3
#define MAX_N 1000000001

vector<PD> cli(double a, double b, double c, PD c1) {
  double l = a * a + b * b, k = a * c1.fi + b * c1.se + c, d = l - k * k;
  vector<PD> v;
  if (d > 0) {
    double ds = sqrt(d), apl = a / l, bpl = b / l;
    double xc = c1.fi - apl * k, xd = bpl * ds, yc = c1.se - bpl * k, yd = apl * ds;
    v.pb(PD(xc - xd, yc + yd)), v.pb(PD(xc + xd, yc - yd));
  } else if (d == 0) {
    v.pb(PD(c1.fi - a * k / l, c1.se - b * k / l));
  }
  return v;
}

vector<PD> cci(PD c1, PD c2) {
  double a = c1.fi - c2.fi, b = c1.se - c2.se;
  return cli(2.0 * a, 2.0 * b, 0 - a * (c1.fi + c2.fi) - b * (c1.se + c2.se), c1);
}

bool dist(PD p1, PD p2) {
  return (pow(p1.fi - p2.fi, 2.0) + pow(p1.se - p2.se, 2.0)) <= 1.0 + 1e-7;
}

int main(){
  cin.sync_with_stdio(false);
  int n;
  while(cin >> n, n) {
    vector<PD> v(n);
    rep(i, n) cin >> v[i].fi >> v[i].se;

    int maxi = 1;
    rep(i, n - 1)repl(j, i + 1, n) {
      vector<PD> p = cci(v[i], v[j]);
      rep(k, p.size()) {
        int cnt = 0;
        rep(l, n) {
          if (dist(p[k], v[l])) cnt++;
        }
        maxch(maxi, cnt);
      }
    }

    cout << maxi << 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-7;

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 = 1;
        pdd p[330];
        rep(i,n) cin >> p[i].fi >> p[i].se;
        rep(i,n-1){
            REP(j,i+1,n){

                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 = sqrt(1-dist((x1-x2),(y1-y2))/4);
                double theta = atan((y2-y1)/(x2-x1));
                cx1 = mx - sin(theta)*ld;
                cy1 = my + cos(theta)*ld;
                cx2 = mx + sin(theta)*ld;
                cy2 = my - cos(theta)*ld;

                int cnt1 = 0,cnt2=0;
                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;
}