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<cmath>

using namespace std;

int n;
double x[300];
double y[300];

double dist(double x1, double y1, double x2, double y2) {
    double dx = abs(x1 - x2);
    double dy = abs(y1 - y2);
    return sqrt(dx * dx + dy * dy);
}

int count(double cx, double cy, int i, int j) {
    int ans = 0;
    for (int k = 0; k < n; ++k) {
        if (k == i || k == j || dist(cx, cy, x[k], y[k]) < 1) ans++;
    }
    return ans;
}

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

    int ans = 0;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            if (i == j) continue;
            double d = dist(x[i], y[i], x[j], y[j]);
            if (d > 2.) continue;
            double dx = x[j] - x[i];
            double dy = y[j] - y[i];
            double r = sqrt(1 - (d * d / 4));
            // cout << r << endl;
            double cx = x[i] + dx / 2 - dy * r / d;
            // cout << d << endl;
            double cy = y[i] + dy / 2 + dx * r / d;
            // cout << x[i] << ", " << y[i] << " - " << x[j] << ", " << y[j] << " : " << cx << ", "<< cy << endl;
            ans = max(ans, count(cx, cy, i, j));
        }
    }
    cout << ans << endl;
}

int main() {
    while (1) {
        cin >> n;
        if (n == 0) break;
        solve();
    }
}

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-6;
	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 <iostream>
#include <cmath>
#define EPS 1.0e-4
using namespace std;

typedef struct {
  double x,y;
} P;

double dis(P p ){
  return sqrt( p.x*p.x + p.y*p.y );
}
double dis(P p1,P p2 ){
  P p;
  p.x = p1.x-p2.x;
  p.y = p1.y-p2.y;
  return dis( p );
}


int main(){
  int n, mx;
  // mv : @üxNg
  P p[301], m, p1,p2, mv,dmv, c[2];
  double x;

  while( cin>>n && n ){
    for( int i=0;i<n;i++ )
      cin >> p[i].x >> p[i].y ;

    mx=0;
    for( int i=0;i<n;i++ ){
      for( int j=i+1;j<n;j++ ){
	if( dis(p[i],p[j]) < 2.0+EPS ){
	  p1=p[i]; p2=p[j];
	  m.x = (p1.x+p2.x)/2.0;
	  m.y = (p1.y+p2.y)/2.0;

	  x=dis(p1,m);
	  x=sqrt( 1-x*x );
	  // S_ cðßé
	  if( fabs( p1.y-p2.y ) < EPS ){
	    mv.x=0.0;  mv.y=1.0;
	  }else if( fabs( p1.x-p2.x ) < EPS ){
	    mv.x=1.0;  mv.y=0.0;
	  }else{
	    double a1 = ( p2.y-p1.y )/( p2.x-p1.x );
	    mv.x=1.0; mv.y=-1.0/a1;
	  }
	  double scale = dis( mv );
	  dmv.x = x*mv.x / scale;
	  dmv.y = x*mv.y / scale;

	  c[0].x=m.x+dmv.x; c[0].y=m.y+dmv.y;
	  c[1].x=m.x-dmv.x; c[2].y=m.y-dmv.y;
	  /*
	  cout << "m  " << m.x    << "," << m.y    << endl;
	  cout << "mv " << mv.x   << "," << mv.y   << endl;
	  cout << "dm " << dmv.x  << "," << dmv.y  << endl;
	  cout << "c0 " << c[0].x << "," << c[0].y << endl;
	  cout << "c1 " << c[1].x << "," << c[1].y << endl;
	  */

	  for( int cn=0;cn<2;cn++ ){
	    int sum=0;
	    for( int k=0;k<n;k++ ){
	      if( dis( c[cn],p[k] ) < 1.0+EPS )
		sum++;
	    }
	    if( mx<sum ) mx=sum;
	  }
	}
      }
    }
    cout << mx << endl;
  }

  return 0;
}

C++: #include<iostream>
#include<cmath>
#define MAX 10000
#include<complex>
using namespace std;
typedef complex<double> P;
int main(){
	int n;
	double p[MAX],q[MAX];
	while(cin>>n,n){
		P tab[MAX];
		int ans=1;
		for(int i=0;i<n;i++){
			cin>>p[i]>>q[i];
			tab[i]=P(p[i],q[i]);
		}
		for(int i=0;i<n;i++)
		for(int j=i+1;j<n;j++){
			if(   norm(tab[i]-tab[j]) >4.0){
				//cout<<88<<endl;
				continue;
			}
			P temp=(tab[i]-tab[j])*sqrt(1-(((p[i]-p[j])*(p[i]-p[j])+(q[i]-q[j])*(q[i]-q[j]))*0.25))/abs(tab[i]-tab[j]);
			P pre=(tab[i]+tab[j])*0.5+temp*P(0,1);
			int r=2;
			for(int k=0;k<n;k++){
				if(norm(tab[k]-pre)<=1.0&&(i!=k&&j!=k))
				r++;
			}
			ans=max(ans,r);
			temp=(tab[i]-tab[j])*sqrt(1-(   ( (p[i]-p[j])*(p[i]-p[j]) + (q[i]-q[j])*(q[i]-q[j]) )*0.25))/abs(tab[i]-tab[j]);
			pre=(tab[i]+tab[j])*0.5+temp*P(0,1)*P(0,1)*P(0,1);
			r=2;
			for(int k=0;k<n;k++){
				if(norm(tab[k]-pre)<=1.0&&(i!=k&&j!=k))
				r++;
			}
			ans=max(ans,r);
			
		}
		cout<<ans<<endl;
	}
	return 0;
}