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>
#include <complex>

using namespace std;

int main(){
	int n;
	cin >> n;
	while(n){
		complex<long double> p[300];
		int s = 1;
		for(int i = 0; i < n; i++){
			long double a,b;
			cin >> a >> b;
			p[i] = complex<long double>(a,b);
		}
		for(int i  = 0; i < n; i++){
			for(int j = i+1; j < n; j++){
				if(abs(p[i]-p[j]) < 2){
					complex<long double> _p = polar(pow(1-norm(p[i]-p[j])/4,0.5),arg(p[i]-p[j])+M_PI_2);
					complex<long double> o[2];
					for(int k = 0; k < 2; k++){
						o[k] = complex<long double>(real(p[i]+p[j])/2,imag(p[i]+p[j])/2);
					}
					o[0] += _p;
					o[1] -= _p;
					for(int k = 0; k < 2; k++){
						int t = 0;
						for(int l = 0; l < n; l++){
							if(abs(p[l]-o[k]) <= 1.0001) ++t;
						}
						s = max(s,t);
					}
				}
			}
		}
		cout << s << endl;
		cin >> n;
	}
}

C++: #include<bits/stdc++.h>
using namespace std;
typedef pair<double,double> P;
int main(){
  int n;
  while(cin>>n,n){
    vector<P> v;
    double x,y,mx1,mx2,my1,my2,p,q,a,b,c;
    int ans=1;
    int i,j,k,l;
    for(i=0;i<n;i++){
      cin >> x >> y;
      v.push_back(P(x,y));
    }
    for(i=0;i<n;i++){
      for(j=i+1;j<n;j++){
        mx1=v[i].first;my1=v[i].second;
	mx2=v[j].first;my2=v[j].second;
	//cout << (mx1-mx2)*(mx1-mx2)+(my1-my2)*(my1-my2) << endl;
	if((mx1-mx2)*(mx1-mx2)+(my1-my2)*(my1-my2)>2.0*2.0) continue;
	p=(my1-my2)/(mx2-mx1);
	q=((mx1+mx2)-p*(my1+my2))/2;
	//cout << p << ":" << q << endl;
	a=p*p+1;b=2*p*q-2*my1-2*mx1*p;c=q*q-2*mx1*q+mx1*mx1+my1*my1-1.0;
	//cout << a << ":" << b << ":" << c << endl;
	l=0;
	//cout << b*b-4*a*c << "/";
	y=(-b+sqrt(b*b-4*a*c))/(2*a);x=p*y+q;
	//cout << x << ":" << y << endl;
	for(k=0;k<n;k++){
	  if((v[k].first-x)*(v[k].first-x)
	     +(v[k].second-y)*(v[k].second-y)<=1.0001) l++;
	}
	l=0;
	y=(-b-sqrt(b*b-4*a*c))/(2*a);x=p*y+q;
	for(k=0;k<n;k++){
	  if((v[k].first-x)*(v[k].first-x)
	     +(v[k].second-y)*(v[k].second-y)<=1.0001) l++;
	}
	ans=max(ans,l);
      }
    }
    cout << ans << endl;
  }
  return 0;
}

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

#define MAX_N 310

double x[MAX_N], y[MAX_N];
int n, maxn;

int main() {
	while (true) {
		cin >> n;
		if (n == 0) { break; }
		for (int i = 0; i < n; i++) {
			cin >> x[i] >> y[i];
		}
		maxn = 0;
		for (int i = 0; i < n; i++) {
			for (int j = i + 1; j < n; j++) {
				int sum = 0;
				double X = (x[i] + x[j]) / 2.0;
				double Y = (y[i] + y[j]) / 2.0;
				for (int k = 0; k < n; k++) {
					if (((x[k] - X)*(x[k] - X) + (y[k] - Y)*(y[k] - Y)) <= 1) {
						sum++;
					}
				}
				maxn = max(maxn, sum);
			}
		}
		for (int i = 0; i < n; i++) {
			int sum = 0;
			double X = x[i];
			double Y = y[i];
			for (int k = 0; k < n; k++) {
				if (((x[k] - X)*(x[k] - X) + (y[k] - Y)*(y[k] - Y)) <= 1) {
					sum++;
				}
			}
			maxn = max(maxn, sum);
		}
		cout << maxn << endl;
	}
	return 0;
}

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

#define EPS (1e-7)
#define X second
#define Y first

int main(){

  while(true){
    int N, ans = 1, n = 0;
    long double x, y;
    vector<long double> line;
    vector<pair<long double,int> > ip;
    vector<pair<long double, long double> > points;
    
    cin>>N;

    if(!N) break;

    
    for(int i = 0; i < N; i++){
      cin>>x>>y;
      
      points.push_back(make_pair(y,x));
      line.push_back(y+1.0);
      line.push_back(y-1.0);
    }

    for(int i = 0; i < N; i++){
      for(int j = i+1; j < N; j++){
       	if((points[i].X-points[j].X)*(points[i].X-points[j].X) + (points[i].Y-points[j].Y)*(points[i].Y-points[j].Y) > 4.0) continue;
       	long double X = points[i].X - points[j].X, Y = points[i].Y - points[j].X, a, x1, x2, y1, y2;
			
       	a = (X*X + Y*Y)/2.0;
       	y1 = (a*Y - X*sqrt(X*X+Y*Y - a*a))/(X*X+Y*Y) + points[j].Y;	
       	y2 = (a*Y + X*sqrt(X*X+Y*Y - a*a))/(X*X+Y*Y) + points[j].Y;	
       	line.push_back(y1);
       	line.push_back(y2);
       		
      }
    }

    sort(line.begin(), line.end());

    for(int i = 0; i < line.size()-1; i++){
      line[i] = (line[i+1]+line[i])/2.0;
    }
    line.pop_back();
    sort(points.begin(),points.end());

    for(int i = 0; i < line.size(); i++){
      
      for(int j = n; j < N; j++){
	long double a = fabs(line[i] - points[j].Y), b;
	if(a < 1.0){
	  b = sqrt(1.0 - a*a);
	  ip.push_back(make_pair(points[j].X-b,0));
	  ip.push_back(make_pair(points[j].X+b,1));
	}
      }

      sort(ip.begin(), ip.end());
      
      for(int j = 0, num = 0; j < ip.size(); j++,ans = max(ans,num) ) num += ip[j].second?-1:1;
      
      ip.clear();
    }

    cout<<ans<<endl;
  }

  return 0;
}