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 _ 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.0001L)cntc++;
	if(norm(p[k]-d)<=1.0001L)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();
  }
}

C++: #include<cstdio>
#include<complex>

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

using namespace std;

typedef complex<double> Point;

const double EPS=1e-9;

int n;
Point p[300];

int f(const Point &c){
	int res=0;
	rep(i,n) if(abs(p[i]-c)<1+EPS) res++;
	return res;
}

int main(){
	for(;scanf("%d",&n),n;){
		rep(i,n){ double x,y; scanf("%lf%lf",&x,&y); p[i]=Point(x,y); }

		int ans=0;
		rep(j,n) rep(i,j) {
			Point v=0.5*(p[j]-p[i]);
			double d=abs(v);
			if(d<1-EPS){
				ans=max(ans,f(p[i]+v+(sqrt(1-d*d)/d)*v*Point(0,+1)));
				ans=max(ans,f(p[i]+v+(sqrt(1-d*d)/d)*v*Point(0,-1)));
			}
		}
		printf("%d\n",ans);
	}

	return 0;
}

C++: #include <stdio.h>
#include <cmath>
#include <algorithm>
#include <stack>
#include <queue>
#include <vector>
typedef long long int ll;
#define BIG_NUM 2000000000
#define MOD 1000000007
#define EPS 0.000001
using namespace std;

struct Point{
	bool operator<(const struct Point &arg) const{
			return x < arg.x;
	};
	double x,y;
};

double a,b,c,d;

double calc(double x){
	return ((c-a)*x)/(b-d)+(b+d)/2-(c-a)*(a+c)/(2.0*(b-d));
}


double dist(double x,double y){
	double ans1 = sqrt((x-a)*(x-a)+(y-b)*(y-b));
	double ans2 = sqrt((x-c)*(x-c)+(y-d)*(y-d));

	if(fabs(1.0-ans1) < fabs(1.0-ans2)){
		return ans1;
	}else{
		return ans2;
	}
}

int main(){

	int N,maximum,count;
	Point points[300];
	double left,right,m,tmp_x,tmp_y,tmp_dist;

	bool FLG;

	while(true){
		scanf("%d",&N);
		if(N == 0)break;

		maximum = 1;

		for(int loop = 0; loop < N; loop++){
			scanf("%lf %lf",&points[loop].x,&points[loop].y);
		}

		sort(points,points+N);

		for(int i = 0; i <= N-2; i++){
			for(int k = i+1; k <= N-1; k++){

				a = points[i].x;
				b = points[i].y;
				c = points[k].x;
				d = points[k].y;

				left = 0.0;
				right = (a+c)/2.0;
				m = (left+right)/2.0;

				FLG = false;

				while(left <= right){
					tmp_y = calc(m);
					tmp_dist = dist(m,tmp_y);

					if(fabs(tmp_dist-1.0) < EPS){
						FLG = true;
						tmp_x = m;
						break;
					}else if(tmp_dist > 1.0){
						left = m + EPS;
					}else{
						right = m - EPS;
					}
					m = (left+right)/2.0;
				}

				if(FLG){
					count = 0;
					for(int p = 0; p < N; p++){
						if(sqrt((tmp_x-points[p].x)*(tmp_x-points[p].x)+(tmp_y-points[p].y)*(tmp_y-points[p].y)) < 1.0){
							count++;
						}else if(fabs(sqrt((tmp_x-points[p].x)*(tmp_x-points[p].x)+(tmp_y-points[p].y)*(tmp_y-points[p].y))-1.0) <= EPS){
							count++;
						}
					}
					maximum = max(maximum,count);
				}

				left = (a+c)/2.0;
				right = 10.0;
				m = (left+right)/2.0;
				FLG = false;

				while(left <= right){
					tmp_y = calc(m);
					tmp_dist = dist(m,tmp_y);

					if(fabs(tmp_dist-1.0) < EPS){
						FLG = true;
						tmp_x = m;
						break;
					}else if(tmp_dist > 1.0){
						right = m - EPS;
					}else{
						left = m + EPS;
					}
					m = (left+right)/2.0;
				}

				if(FLG){
					count = 0;
					for(int p = 0; p < N; p++){
						if(sqrt((tmp_x-points[p].x)*(tmp_x-points[p].x)+(tmp_y-points[p].y)*(tmp_y-points[p].y)) < 1.0){
							count++;
						}else if(fabs(sqrt((tmp_x-points[p].x)*(tmp_x-points[p].x)+(tmp_y-points[p].y)*(tmp_y-points[p].y))-1.0) <= EPS){
							count++;
						}
					}
					maximum = max(maximum,count);
				}

			}
		}

		printf("%d\n",maximum);
	}

	return 0;
}

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

complex<double> p[300];

complex<double> core(int a,int b){
  complex<double> res,mid,norm;

  mid = (p[a]+p[b]) / 2.0;

  norm = (p[a] - p[b]) * complex<double>(0,1);
  norm = norm / abs(norm);
  norm = norm * sqrt(1.0 - abs(mid-p[a])*abs(mid-p[a]));

  res = norm + mid;
  return res;
}
  

int main(){
  int n;
  complex<double> o;
  int max,tmp;

  while(1){
    cin >> n;
    if(!n)break;

    for(int i=0;i<n;i++)cin >> p[i].real() >> p[i].imag();

    max = 1;
    for(int i=0;i<n;i++){
      for(int j=0;j<n;j++){
	if(i!=j && abs(p[i]-p[j]) < 2.0){
	  o = core(i,j);
	  tmp = 2;
	  for(int k=0;k<n;k++){
	    if(k!=i && k!=j){
	      if(abs(o-p[k])<1.0)tmp++;
	    }
	  }
	  if(max < tmp)max = tmp;
	}
      }
    }
    cout << max << endl;
  }
}

C++: #include <bits/stdc++.h>
#define N 301
using namespace std;
typedef complex <double> P;
const double EPS = 1e-8;
P p[N];
int n;

int cnt(P c){
  int res=0;
  for(int i=0;i<n;i++)if(abs(c-p[i])<1+EPS) res++;
  return res;
}

P get_C(P a,P b){
  double x=abs(b-a);
  if(x>2||a==b) return a;
  return P(x/2,sqrt(1-x*x/4))*(b-a)/abs(b-a)+a;  
}


int main(){
  while(1){
  cin>>n;
  if(!n)break;
  double a,b;
  for(int i=0;i<n;i++)cin>>a>>b,p[i]=P(a,b);
  
  int ans=0;
  for(int i=0;i<n;i++)
    for(int j=0;j<n;j++)
      ans=max(ans,cnt(get_C(p[i],p[j])));
  cout <<ans<<endl;  
  }
  return 0;
}