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 <complex>
#include <cmath>
using namespace std;

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

#define EPS (1e-10)

typedef complex<double> P;

int main(void){
	int n;

	while(cin>>n && n){
		P p[n];
		rep(i,n) cin>>p[i].real()>>p[i].imag();

		int ans = 0;
		rep(i,n-1)REP(j,i+1,n){
			P v = p[i] - p[j];
			if(abs(v) > 2) continue;

			P q = P((p[i].real()+p[j].real())/2, (p[i].imag()+p[j].imag())/2);
			double dist = sqrt(1 - pow(abs(p[i]-p[j])/2, 2));
			P n1 = (v * P(0,+1)) / abs(v);
			P n2 = (v * P(0,-1)) / abs(v);
			P o1 = q + n1 * dist;
			P o2 = q + n2 * dist;

			int sum = 0;
			rep(k,n){
				if(abs(o1 - p[k]) < 1 + EPS) sum++;
			}

			ans = max(ans,sum);
			sum = 0;

			rep(k,n){
				if(abs(o2 - p[k]) < 1 + EPS) sum++;
			}
			ans = max(ans,sum);
		}

		cout<<(ans==0?1:ans)<<endl;
	}

	return 0;
}

C++: #include <iostream> 
#include <string> 
#include <vector> 
#include <cmath> 
#include <algorithm> 
#include <cstdlib> 
#include <ctime> 
#include <cstdio> 
#include <functional> 
#include <set> 
#include <sstream> 


using namespace std; 

int main(){

	int n;
	const double eps=1e-8;
	while(true){

		cin>>n;

		if(n==0) return 0;

		int ans=0;
		vector< double > x(n),y(n);

		for(int i=0;i<n;i++) cin>>x[i]>>y[i];


		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				const double d=hypot(x[i]-x[j],y[i]-y[j]);
				if(d>2.0-eps) continue;

				const double mx=(x[i]+x[j])/2.0;
				const double my=(y[i]+y[j])/2.0;
				const double dvx=(x[i]-x[j])/d;
				const double dvy=(y[i]-y[j])/d;

				const double cx1=mx+dvy*sqrt(1.0-d*d/4.0);
				const double cy1=my-dvx*sqrt(1.0-d*d/4.0);

				const double cx2=mx-dvy*sqrt(1.0-d*d/4.0);
				const double cy2=my+dvx*sqrt(1.0-d*d/4.0);

				int res=2;
				for(int k=0;k<n;k++){
					if(k==i || k==j) continue;
//					if(hypot(x[k]-cx1,y[k]-cy1)<1.0+eps) res++;
					if((x[k]-cx1)*(x[k]-cx1)+(y[k]-cy1)*(y[k]-cy1)<1.0+eps) res++;
				}
				
				ans=max(ans,res);

				res=2;
				for(int k=0;k<n;k++){
					if(k==i || k==j) continue;
//					if(hypot(x[k]-cx2,y[k]-cy2)<1.0+eps) res++;
					if((x[k]-cx2)*(x[k]-cx2)+(y[k]-cy2)*(y[k]-cy2)<1.0+eps) res++;
				}

				ans=max(ans,res);

			}
		}

		cout<<ans<<endl;

	}

	return 0;
}

C++: #include <bits/stdc++.h>

using namespace std;

#define X first
#define Y second

using P=pair<double, double>;

const double EPS = 1e-14;

double dist(P a, P b) {
	double dx = a.X - b.X;
	double dy = a.Y - b.Y;
	return sqrt(dx*dx+dy*dy);
}

int main() {
	int n;
	while(cin >> n, n) {
		vector<P> p(n), pl;
		for(int i = 0; i < n; i++) cin >> p[i].X >> p[i].Y;
		for(int i = 0; i < p.size(); i++) {
			P p1 = p[i];
			for(int j = 0; j < p.size(); j++) {
				if(i==j) continue;
				P p2 = p[j];
				if(dist(p1, p2)+EPS > 2) continue;

				double dxy = dist(p1, p2);
				double l = sqrt(1-(dxy/2));
				double cx = (p1.X+p2.X)/2;
				double cy = (p1.Y+p2.Y)/2;
				double dx = (p2.X - p1.X);
				double dy = (p2.Y - p1.Y);

				double dlx = l/sqrt(dx*dx+dy*dy) * -dy;
				double dly = l/sqrt(dx*dx+dy*dy) * dx;

				pl.push_back(P(cx+dlx, cy+dly));
				pl.push_back(P(cx-dlx, cy-dly));
			}
		}

		int res = 0;
		for(P p1:pl) {
			int tmp = 0;
			for(P p2:p) {
				if(dist(p1, p2)+EPS <= 1) tmp++;
			}
			res = max(res, tmp);
		}

		cout << res << endl;
	}

	return 0;
}

C++: #include <iostream>
#include <vector>
#include <complex>
#include <algorithm>
#define rep(i,end) for(int i=0;i<end;i++)
using namespace std;
typedef complex<double> EN;

int search(const EN a,vector<EN> dat){
  int count=0;
  rep(i,dat.size()){
    if(abs(a-dat[i])<=1) count++;
  }
  //  cout <<count<<" "<<a<<endl;
  return count;
}

EN tekitan(vector<EN> dat ){
  
  EN ans(0,0);
  int maxn=0;
  /*
  vector<int> chart(10000,0);
  for(int i=0;i<10000;i++)
    chart[i]=search(EN(i/1000,((int)i%100)/10.0),dat);
  for(double i=0;i<10000;i++)
    if(maxn<chart[i]){
      maxn=chart[i];
      ret=i;
    }
  */
  for(double i=0;i<=10.0;i+=0.1)
    for(double j=0;j<=10.0;j+=0.1)
      if (search(EN(i,j),dat)>maxn)
	{
	  ans=EN(i,j);
	  maxn=search(ans,dat);
	}
  return ans;
}
int pretan(vector<EN> dataset,EN nearans){
  int count=search(nearans,dataset);
  double i=0.1;
  int up,down,left,right;
  for(i=0.1;i>0.00005;i/=2){
    //    cout<<nearans<<endl;
    up=search(nearans+EN(0,i),dataset);
    down=search(nearans+EN(0,-i),dataset);
    left=search(nearans+EN(-i,0),dataset);
    right=search(nearans+EN(i,0),dataset);
    if(up==max(up,max(down,max(left,max(right,count)))))
      {count=up;nearans+=EN(0,i);continue;}
    if(down==max(up,max(down,max(left,max(right,count)))))
      {count=down;nearans+=EN(0,-i);continue;}
    if(left==max(up,max(down,max(left,max(right,count)))))
      {count=left;nearans+=EN(-i,0);continue;}
    if(right==max(up,max(down,max(left,max(right,count)))))
      {count=right;nearans+=EN(i,0);continue;}
  
  }
  
  
  return count;
}
int main(){
  int m;
  double a,b,tmp;
 
  vector<complex<double> > dataset;
  EN nearans;
  cin >>m;
  while(m){

    dataset.resize(m);
    rep(i,m){
      cin>>a>>b;
      dataset[i]=EN(a,b);
    }
    
    nearans=tekitan(dataset);

    cout <<pretan(dataset,nearans)<<endl;
    cin >>m;
  }

  return 0;
}

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

using namespace std;

const double EPS = 1e-4;

struct Coor {
	double x, y;
};

double pw(double x) { return x * x; };
int max(int n, int m) { return n > m ? n : m; };

int main()
{
	double n;
	vector<Coor> v(300);

	while (cin >> n) {
		if (n == 0) break;

		

		for (int i = 0; i < n; i++) {
			cin >> v.at(i).x >> v.at(i).y;
		}

		double x1, y1, x2, y2;
		int _max = 1;

		Coor mid;

		for (int i = 0; i < n-1; i++) {
			x1 = v.at(i).x; y1 = v.at(i).y;

			for (int j = i + 1; j < n; j++) {
				x2 = v.at(j).x; y2 = v.at(j).y;
				mid.x = (x1 + x2) / 2;
				mid.y = (y1 + y2) / 2;

				double d = sqrt(pw(x1 - x2) + pw(y1 - y2));
				double dx = sqrt(1 - pw(d / 2)) * (y1 - y2) / d;
				double dy = sqrt(1 - pw(d / 2)) * (x2 - x1) / d;

				int ct1 = 0, ct2 = 0;

				for (int k = 0; k < n; k++) {
					double tmp = pw(v.at(k).x - (mid.x + dx)) + pw(v.at(k).y - (mid.y + dy));
					if (sqrt(tmp) <= 1.0 + EPS) ct1++;

					tmp = pw(v.at(k).x - (mid.x - dx)) + pw(v.at(k).y - (mid.y - dy));
					if (sqrt(tmp) <= 1.0 + EPS) ct2++;
				}

				_max = max(_max, ct1);
				_max = max(_max, ct2);
			}
		}

		cout << _max << endl;
	}
}