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

class Point {
public:
  double x, y;
  Point(double x = 0, double y = 0) : x(x), y(y) {}
  Point operator +(Point p) {
    return Point(x+p.x, y+p.y);
  }
  Point operator -(Point p) {
    return Point(x-p.x, y-p.y);
  }
  Point operator *(double a) {
    return Point(x*a, y*a);
  }
  Point operator /(double a) {
    return Point(x/a, y/a);
  }
  bool operator ==(const Point &p) const {
    return x == p.x && y == p.y;
  }
  bool operator <(const Point &p) const {
    double l1 = x+y;
    double l2 = p.x+p.y;
    if(l1 != l2) return l1 < l2;
    if(y != p.y) return y < p.y;
    return x < p.x;
  }
};

typedef Point Vector;

double norm(Point p1, Point p2) {
  return (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y);
}

int n;
Point p[300];

int calc(Point p1, Point p2) {
  Point m = Point((p1.x+p2.x)/2.0, (p1.y+p2.y)/2.0);
  Point o;
  Vector moto = Vector(p2.y-p1.y, p1.x-p2.x);
  int cnt = 0;
  int newCnt;
  double p1p2 = norm(p1,p2);
  if(p1p2 > 4) return 1;
  double l = sqrt(4.0/p1p2-1)/2.0;

  for(int t = 0; t < 2; t++) {
    newCnt = 0;
    if(t == 1) {
      moto = moto * (-1);
    }
    o = m + moto * l;

    for(int i = 0; i < n; i++) {
      if(newCnt + n-i <= cnt) break;
      if(p[i] == p1 || p[i] == p2) continue;
      if((p[i].x-o.x)*(p[i].x-o.x) + (p[i].y-o.y)*(p[i].y-o.y) <= 1.0) {
	newCnt++;
      }
    }
    if(newCnt > cnt) cnt = newCnt;
  }
  return cnt + 2;
}

main() {
  while(1) {
    int ans;
    cin >> n;
    if(n == 0) break;
    for(int i = 0; i < n; i++) {
      cin >> p[i].x >> p[i].y;
    }
    sort(p, &p[n]);
    ans = 1;
    for(int i = 0; i < n; i++) {
      for(int j = i+1; j < n; j++) {
	ans = max(ans, calc(p[i],p[j]));
      }
    }
    cout << ans << endl;
  }
}

C++: #include <iostream>
#include <cstdio>
#include <algorithm>
#include <array>
#include <vector>
#include <cmath>

using namespace std;

typedef pair<double, double> Point;

double d2( Point p1, Point p2 ) {
  double dx = p1.first - p2.first;
  double dy = p1.second - p2.second;
  return dx * dx + dy * dy;
}

Point sub( Point p1, Point p2 ) {
  return Point( p1.first - p2.first, p1.second - p2.second );
}

Point add( Point p1, Point p2 ) {
  return Point( p1.first + p2.first, p1.second + p2.second );
}

Point div( Point p, double val ) {
  return Point( p.first / val, p.second / val );
}

Point middle( Point p1, Point p2 ) {
  return div( add( p1, p2 ), 2.0 );
}

double norm( Point p ) {
  return sqrt( p.first * p.first + p.second * p.second );
}

Point normalize( Point p ) {
  double n = norm( p );
  return div( p, n );
}

Point rotate( Point p, double th ) {
  double x = cos( th ) * p.first - sin( th ) * p.second;
  double y = sin( th ) * p.first + cos( th ) * p.second;
  return Point( x, y );
}								

void cross( Point p1, Point p2, Point &a1, Point &a2 ) {
  Point v = normalize( sub( p1, p2 ) );
  double th = acos( sqrt( d2( p1, p2 ) ) / 2.0 );
  a1 = add( rotate( v, th ), p2 );
  a2 = add( rotate( v, -th ), p2 );
}

void print( Point p ) {
  cerr << p.first << ", " << p.second << endl;
}

int main() {
  while( true ) {
	int n;
	cin >> n;
	if ( n == 0 ) break;
	
	vector<Point> vs( n );
	for ( int i = 0; i < n; i++ ) {
	  double x, y;
	  cin >> x >> y;
	  vs[ i ] = make_pair( x, y );
	}
	vector<Point> cands;
	for ( int i = 0; i < n; i++ ) {
	  for ( int j = i + 1; j < n; j++ ) {
		double d = d2( vs[ i ], vs[ j ] );
		if ( d == 4 ) {
		  cands.push_back( middle( vs[ i ], vs[ j ] ) );
		} else if ( d < 4 ) {
		  Point p1, p2;
		  cross( vs[ i ], vs[ j ], p1, p2 );
		  cands.push_back( p1 );
		  cands.push_back( p2 );
  		}
	  }
	}
	int maxCount = 1;
	for ( Point p : cands ) {
	  int c = 0;
	  for ( Point v : vs ) {
		if ( d2( p, v ) <= 1.0 ) {
		  c++;
		}
	  }
	  maxCount = max( maxCount, c );
	}
	cout << maxCount << endl;
  }
}

C++: #define _USE_MATH_DEFINES

#include <cstdio>
#include <cstdlib>
#include <string.h>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <queue>
#include <stack>
#include <map>

using namespace std;

typedef pair<long long int, long long int> P;
long long int INF = 1e18;
long long int MOD = 1000000007;

int main(){
	while(true){
		int N;
		cin >> N;
		if(N == 0){
			break;
		}
		
		double x[500], y[500];
		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++){
			double vx = x[j] - x[i];
			double vy = y[j] - y[i];
			double norm = sqrt(vx * vx + vy * vy);
			if(norm > 2 || i == j){
				continue;
			}
			double vnx = vx / norm;
			double vny = vy / norm;
			norm /= 2;
			double norm2 = sqrt(1 - norm * norm);
			double v90x = -vny;
			double v90y = vnx;
			double centx = x[i] + vx / 2 + v90x * norm2;
			double centy = y[i] + vy / 2 + v90y * norm2;
			int S = 2;
			for(int k = 0; k < N; k++){
				if(i == k || j == k){
					continue;
				}
				double xx = centx - x[k];
				double yy = centy - y[k];
				if(xx * xx + yy * yy <= 1){
					S++;
				}
			}
			ans = max(S, ans);
		}
		cout << max(1, ans) << endl;
	}
	return 0;
}