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<math.h>

using namespace std;

#define REP(i, n) for(int i=0; i<n; i++)

double p[300][2];
int n;

double dis(double x1, double y1, double x2, double y2){
    return sqrt(pow(x1-x2, 2) + pow(y1-y2, 2));
}

int circle(double x, double y, int i, int j){
    int count = 2;
    REP(k, n){
        if(k!=i && k!=j){
            if(dis(x, y, p[k][0], p[k][1]) < 1){
                count += 1;
            }
        }
    }
    return count;
}

int main(){
    while(1){
        cin >> n;

        if(n==0) break;

        REP(i, n){
            cin >> p[i][0];
            cin >> p[i][1];
        }

        int output = 0;
        REP(i, n){
            for(int j=i+1; j<n; j++){
                double cx = (p[i][0]+p[j][0]) / 2;
                double cy = (p[i][1]+p[j][1]) / 2;
                double d = dis(p[i][0], p[i][1], p[j][0], p[j][1]);

                if(d>2){
                    continue;
                } else if (d==2){
                    output = max(output, circle(cx, cy , i, j));
                } else {
                    double a1 = acos(d/2);
                    double a2 = atan2((p[j][1]-p[i][1]), p[j][0]-p[i][0]);

                    double x1 = p[i][0] + cos(a2+a1);
                    double y1 = p[i][1] + sin(a2+a1);

                    double x2 = p[i][0] + cos(a2-a1);
                    double y2 = p[i][1] + sin(a2-a1);

                    output = max(output, circle(x1, y1, i, j));
                    output = max(output, circle(x2, y2, i, j));
                }
            }
        }

        cout << output << endl;
    }
}

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

typedef complex<double> P;

const double EPS = 1e-8;
const int N = 300;
const double R = 1;

struct Circle{
  P c;
  double r;
  Circle(P c=P(0,0), double r=0):c(c),r(r){}
};

vector<P> getIntersectCC(Circle c1, Circle c2){
  vector<P> res;
  double r1 = c1.r, r2 = c2.r, d;
  P p1 = c1.c, p2 = c2.c;
  d = abs(p1-p2);

  if(d < EPS && abs(r1-r2) < EPS) return res; //????????¨?????????????????´???????????????????????????????????????
  if(r1 + r2< d - EPS || d + EPS < abs(r1-r2)) return res; //???????????????

  double a = (r1*r1 - r2*r2 + d*d) / (2*d);
  double h = sqrt(max(r1*r1 - a*a, 0.0));
  P tmp1 = p1 + a / d * (p2-p1);
  P tmp2 = h / d * (p2-p1);

  if(abs(tmp2) < EPS) res.push_back(tmp1);
  else {
    res.push_back(P(real(tmp1) - imag(tmp2), imag(tmp1) + real(tmp2)));
    res.push_back(P(real(tmp1) + imag(tmp2), imag(tmp1) - real(tmp2)));
  }
  return res;
}


int n;
vector<P> v;

int solve(){
  int ans = 0;

  for(int i=0;i<n;i++){
    for(int j=i+1;j<n;j++){
      vector<P> cc = getIntersectCC(Circle(v[i], R), Circle(v[j], R)); // v[i]??¨v[j]??????????????¢?????¨??????R??¨????????? (?????´?????????????????????????????????????????????)
      if(cc.size() == 0) continue;
      for(int k=0;k<cc.size();k++){
        int res = 0;
        for(int l=0;l<n;l++) if(abs(v[l]-cc[k]) <= R + EPS) res++;
        ans = max(ans, res);
      }
    }
  }

  return ans;
}

int main(){
  while(cin >> n && n){
    v.clear();
    for(int i=0;i<n;i++){
      double x, y;
      cin >> x >> y;
      v.push_back(P(x,y));
    }
    cout << solve() << endl;
  }
}

C++: #include <iostream>
#include <algorithm>
#include <cassert>
#include <cctype>
#include <complex>
#include <cstdio>
#include <map>
#include <math.h>
#include <queue>
#include <set>
#include <stack>
#include <string>
#include <vector>
using namespace std;

#define EPS 1e-9
#define rep(i,n) for(int i=0;i<n;i++)

int N;
double x[300],y[300];

int main(){
	while(cin>>N&&N){
		int ans=1;
		rep(i,N)cin>>x[i]>>y[i];
		rep(i,N)for(int j=i+1;j<N;j++){
			double x1=x[i],y1=y[i],x2=x[j],y2=y[j];
			if((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)<4){
				double xm=(x1+x2)/2,ym=(y1+y2)/2,dx=x2-x1,dy=y2-y1,n=sqrt(dx*dx+dy*dy),d=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)),dd=sqrt(1-d*d/4);
				double cx=xm-dy/n*dd,cy=ym+dx/n*dd;
				int res=0;
				rep(k,N)if((x[k]-cx)*(x[k]-cx)+(y[k]-cy)*(y[k]-cy)<=1+EPS)res++;
				ans=max(ans,res);
				cx=xm+dy/n*dd;cy=ym-dx/n*dd;
				res=0;
				rep(k,N)if((x[k]-cx)*(x[k]-cx)+(y[k]-cy)*(y[k]-cy)<=1+EPS)res++;
				ans=max(ans,res);
			}
		}
		cout<<ans<<endl;
	}
}

C++: #include <iostream>
#include <string>
#include <algorithm>
#include <vector>
#include <utility>
#include <cstring>
#include <cmath>
#include <queue>
#include <set>
#include <map>
using namespace std;

#define rep(i,n) for(int i=0; i<n; i++)
typedef long long ll;
typedef pair<int,int> pii;

const double EPS = 1e-6;

int main() {
	int n;
	while(cin >> n, n) {
		double x[n], y[n];
		rep(i,n) {
			cin >> x[i] >> y[i];
		}
		int maxi = 1;
		rep(i, n) {
			for(int j=i+1; j<n; j++) {
				if( ((x[i]-x[j])*(x[i]-x[j])+(y[i]-y[j])*(y[i]-y[j])) > 4.0 ) continue;
				double mx = (x[i]+x[j]) / 2.0;
				double my = (y[i]+y[j]) / 2.0;
				double px0, py0, px1, py1;
				if( fabs(y[i]-y[j]) < EPS ) {
					px0 = px1 = mx;
					py0 = my + sqrt(1.0 - (x[i]-mx)*(x[i]-mx));
					py1 = my - sqrt(1.0 - (x[i]-mx)*(x[i]-mx));
				} else if( fabs(x[i]-x[j]) < EPS ) {
					py0 = py1 = my;
					px0 = mx + sqrt(1.0 - (y[i]-my)*(y[i]-my));
					px1 = mx - sqrt(1.0 - (y[i]-my)*(y[i]-my));
				} else {
					double m = - (x[i]-x[j]) / (y[i]-y[j]);
					double l = sqrt(1.0 - ((x[i]-mx)*(x[i]-mx)+(y[i]-my)*(y[i]-my)));
					px0 = mx + sqrt(l*l/(1.0+m*m));
					px1 = mx - sqrt(l*l/(1.0+m*m));
					py0 = my + m * (px0-mx);
					py1 = my + m * (px1-mx);
				}
				int cnt0 = 0, cnt1 = 0;
				rep(k, n) {
					double d0 = (px0-x[k])*(px0-x[k]) + (py0-y[k])*(py0-y[k]);
					if( d0 <= 1.0+EPS ) {
						cnt0++;
					}
				}
				rep(k, n) {
					double d1 = (px1-x[k])*(px1-x[k]) + (py1-y[k])*(py1-y[k]);
					if( d1 <= 1.0+EPS ) {
						cnt1++;
					}
				}
				maxi = max(maxi, cnt0);
				maxi = max(maxi, cnt1);
			}
		}
		cout << maxi << endl;
	}
	return 0;
}