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>
#define pb push_back
#define mp make_pair
#define INF 120000000
#define EPS 1e-6

using namespace std;

typedef complex<double> P;
typedef long long int LL;
vector<P> in;

void read(int &n) {
    in.resize(n);
    for (int i = 0; i < n; i++) {
        double x, y;
        cin >> x >> y;
        in[i] = P(x, y);
    }
}

int solve(int &n) {
    int mx = 1;

    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
            P c[2];
            double v, w;
            v = abs(in[i] - in[j]);
            if (v > 2.0) continue;
            w = sqrt(1.0 - v * v / 4.0);
            // cout << w << endl;
            P un1 = ((in[i] - in[j]) * P(0, 1)) / v;
            P un2 = ((in[i] - in[j]) * P(0, -1)) / v;
            // cout << real(un2) << " " << imag(un2) << endl;
            c[0] = un1 * w + (in[i] + in[j]) / 2.0;
            c[1] = un2 * w + (in[i] + in[j]) / 2.0;
            // cout << real((in[i] + in[j]) / 2.0) << " " << imag((in[i] + in[j]) / 2.0) << endl;
            // cout << real(c[0]) << " " << imag(c[0]) << endl;

            for (int k = 0; k < 2; k++) {
                // cout << real(c[k]) << " " << imag(c[k]) << endl;
                int cnt = 0;
                for (int l = 0; l < n; l++) {
                    // cout << i << " " << j << " " << l << " " << abs(c[k] - in[l]) << endl;
                    if (abs(c[k] - in[l]) <= 1.0 + EPS) {
                        cnt++;
                    }
                }
                mx = max(mx, cnt);
            }
        }
    }

    return mx;
}

int main() {
    int n;
    while (cin >> n, n) {
        read(n);
        cout << solve(n) << endl;
        in.clear();
    }
    return 0;
}

C++: #include<bits/stdc++.h>
#define N 300
using namespace std;
typedef complex<double> P;
const double EPS = 1e-10;

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;
}

double getdistancepp(P a,P b){
  return sqrt((a.real()-b.real())*(a.real()-b.real())+(a.imag()-b.imag())*(a.imag()-b.imag()));
}

int main(){
  int n,ans=0;
  double a,b;
  P p[N];
  while(1){
    cin>>n;
    if(!n)break;
    for(int i=0;i<n;i++){
      cin>>a>>b;
      p[i]=P(a,b);
    }
    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
	Circle c1,c2;
	c1.c=p[i];
	c2.c=p[j];
	c1.r=c2.r=1;
	vector<P> r=getIntersectCC(c1,c2);
	if(!r.size())continue;
	for(int k=0;k<r.size();k++){
	  int cnt=0;
	  for(int l=0;l<n;l++)
	    if(getdistancepp(r[k],p[l])<=1+EPS)cnt++;	    
	  ans=max(ans,cnt);
	}
      }
    }
    cout<<ans<<endl;
  }
  return 0;
}

C++: #include<iostream>
#include<vector>
#include<complex>
#define EPS 0.0001
using namespace std;
typedef double elem;
typedef complex<elem> point,vec;
inline elem dist(const point &a, const point &b){return abs(a-b);}
inline vec normr(const vec &v){return v*vec(0,-1);}
inline vec norml(const vec &v){return v*vec(0,1);}
inline vec unormr(const vec &v){return normr(v) / abs(v);}
inline vec unorml(const vec &v){return norml(v) / abs(v);}
int main(){
	while(true){
		int ans = 0;
		int n;
		scanf("%d", &n);
		if(n==0)break;
		vector< point > points;
		for(int i=0;i<n;++i){
			elem x, y;
			scanf("%Lf%Lf", &x, &y);
			points.push_back( point(x,y) );
		}
		for(unsigned int i = 0; i < points.size(); ++i){
			for(unsigned int j = 0; j < points.size(); ++j){
				if( i != j && dist( points[i], points[j] ) < 2 ){
					int count_p = 0;
					int count_q = 0;
					elem d = sqrt(1-pow(dist(points[i], points[j])/2,2 ) );
					point p,q;
					vec ab = points[j] - points[i];
					point t = 0.5 * ab + points[i];  
					p = d * unormr(ab) + t;
					q = d * unorml(ab) + t;
					for(unsigned int k = 0; k < points.size(); ++k){
						if( dist(p, points[k]) <= 1.0+EPS ){ count_p++; }
						if( dist(q, points[k]) <= 1.0+EPS ){ count_q++; }
					}
					ans = max(ans, max(count_q, count_p) );
				}
			}
		}
		printf("%d\n", ans);
	}
	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;
  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 = sqrt( 1-dis(p1,m) );
	  // ツ陳心ツ点 cツづーツ仰づ淞づゥ
	  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;

	  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<bits/stdc++.h>
using namespace std;
typedef complex<double> point;

point cal(point a,point b,double c){
  double w=abs(b-a);
  double h=sqrt(c*c-w*w/4);
  point cen=point(w/2,h);
  return cen*(b-a)/w+a;
}

int main(){
  int n;
  while(cin>>n,n){
    int ans=1;
    double x,y;
    point p[300];
    for(int i=0;i<n;i++)cin>>x>>y,p[i]=point(x,y);

    for(int i=0;i<n;i++)
      for(int j=i+1;j<n;j++){
	if(abs(p[i]-p[j])>2)continue;
	point center=cal(p[i],p[j],1);
	int can=0;
	for(int k=0;k<n;k++)
	  if(abs(p[k]-center)<=1+1e-8)can++;
	ans=max(ans,can);
      }
    
    cout<<ans<<endl;
  }
  return 0;
}