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 <bits/stdc++.h>
using namespace std;
using ll=long long;
using vi=vector<int>;
using vvi=vector<vi>;
using pii=pair<int,int>;
#define rep(i,n) for(int i=0;i<n;i++)
#define range(i,a,n) for(int i=a;i<n;i++)
#define all(a) a.begin(),a.end()
#define rall(a) a.rbegin(),a.rend()
#define INF 1e9
#define EPS 1e-9
#define MOD (1e9+7)
void put(string d){}template<class H,class...T>void put(string d,H&h,T&...t){cout<<h;if(sizeof...(t))cout<<d;put(d,t...);}
template<class T>void puti(T&a,string d=" "){bool f=1;for(auto&_:a)cout<<(exchange(f,0)?"":d)<<_;cout<<endl;}
template<class T>void putii(T&a,string d=" "){for(auto&_:a)puti(_,d);}

using P=complex<double>;

int main(){
    int n;
    while(cin>>n,n){
        int ans=1;
        vector<P> c(n);
        rep(i,n) { double p,q; cin>>p>>q; c[i]=P(p,q); }
        rep(i,n)range(j,i+1,n){
            P v=c[i]-c[j];
            double d=abs(v);
            if(d>2.0) continue;
            P med=(c[i]+c[j])/2.0;
            P nml=v/d*P(0,1)*sqrt(1-d*d/4);
            P r1=med+nml;
            P r2=med-nml;
            int s1=0,s2=0;
            rep(k,n){
                s1+=(abs(r1-c[k])<=1.0+1e-6);
                s2+=(abs(r2-c[k])<=1.0+1e-6);
            }
            ans=max(ans,max(s1,s2));
        }
        cout<<ans<<endl;
    }
	return 0;
}

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

const double EPS = 1.0e-8;

int n;
double xs[311];
double ys[311];

void solve() {
    int maxval = 1;
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
            double dx = xs[i] - xs[j];
            double dy = ys[i] - ys[j];
            if (dx * dx + dy * dy >= 4.0 - EPS) continue;

            double mx = 0.5 * (xs[i] + xs[j]);
            double my = 0.5 * (ys[i] + ys[j]);
            double nx = -(ys[i] - ys[j]);
            double ny = xs[i] - xs[j];
            double nl = sqrt(nx * nx + ny * ny);
            double cl = sqrt(1.0 - 0.25 * nl * nl);
            nx = cl * nx / nl;
            ny = cl * ny / nl;


            int cnt;
            double cx, cy;

            cx = mx + nx;
            cy = my + ny;
            cnt = 0;
            for (int k = 0; k < n; k++) {
                dx = cx - xs[k];
                dy = cy - ys[k];
                if (dx * dx + dy * dy <= 1.0 + EPS) cnt++;
            }
            maxval = max(maxval, cnt);

            cx = mx - nx;
            cy = my - ny;
            cnt = 0;
            for (int k = 0; k < n; k++) {
                dx = cx - xs[k];
                dy = cy - ys[k];
                if (dx * dx + dy * dy <= 1.0 + EPS) cnt++;
            }
            maxval = max(maxval, cnt);
        }
    }
    printf("%d\n", maxval);
}

int main(int argc, char** argv) {
    while (cin >> n) {
        if (n == 0) break;
        for (int i = 0; i < n; i++) {
            cin >> xs[i] >> ys[i];
        }
        solve();
    }
}

C++: #include <bits/stdc++.h>
using namespace std;
using ll = long long;
using pii = pair<int, int>;

const double EPS = 1e-9;
using Point = pair<double, double>;

vector<Point> getCenter(Point a, Point b) {
  double xm = (a.first + b.first) / 2;
  double ym = (a.second + b.second) / 2;
  double dx = b.first - a.first;
  double dy = b.second - a.second;
  double l = hypot(dx, dy) / 2;
  double t = sqrt((1 - l*l) / (dx*dx + dy*dy));
  vector<Point> res;
  res.emplace_back(xm - (-dy)*t, ym - dx*t);
  res.emplace_back(xm + (-dy)*t, ym + dx*t);
  return res;
}

int N;
double X[310], Y[310];
Point p[310];

bool contain(Point a, Point b) {
  return hypot(a.first-b.first, a.second-b.second) - EPS < 1;
}

int main() {
  cin.tie(0);
  ios_base::sync_with_stdio(false);
  cout << fixed << setprecision(10);
  
  while (cin >> N, N) {
    for (int i = 0; i < N; i++) {
      cin >> X[i] >> Y[i];
      p[i] = Point(X[i], Y[i]);
    }
    int ans = 1;
    for (int i = 0; i < N; i++) {
      for (int j = i+1; j < N; j++) {
        auto c = getCenter(p[i], p[j]);
        for (int ci = 0; ci < c.size(); ci++) {
          int sum = 0;
          for (int k = 0; k < N; k++) {
            if (contain(c[ci], p[k])) sum++;
          }
          ans = max(ans, sum);
        }
      }
    }
    cout << ans << endl;
  }

  return 0;
}

C++: #include<bits/stdc++.h>
#define _ ios_base::sync_with_stdio(0);cin.tie(0);
#define REP(i,n) for(int i=0;i<(int)(n);i++)
#define REPP(i,j,n) for(int i=(j);i<(int)(n);i++)

using namespace std;
typedef long double ld;
typedef complex<ld> Point;
const ld eps = 1e-5, pi = acos(-1.0);
ld dot(Point a, Point b){return real(conj(a)*b);}
ld cross(Point a, Point b){ return imag(conj(a)*b);}

Point p[310];
int n;

void solve(){
  int result=1;
  REP(i,n){
    REPP(j,i+1,n){
      Point ab=p[j]-p[i];
      if(norm(ab)>4)continue;
      Point q=(p[i]+p[j])/2.L;
      Point aq=q-p[j];
      Point normal=Point{imag(aq),-real(aq)}*sqrt(1-norm(aq))/abs(aq);
      Point c=q+normal;
      Point d=q-normal;
      int cntc=0,cntd=0;
      REP(k,n){
	if(norm(p[k]-c)<=1.001L)cntc++;
	if(norm(p[k]-d)<=1.001L)cntd++;
      }
      result=max({result,2,cntc,cntd});
    }
  }
  cout<<result<<endl;
}

int main(){ _;
  ld r,i;
  while(cin>>n,n!=0){
    REP(j,n){
      cin>>r>>i;
      p[j]=Point{r,i};
    }
    solve();
  }
}