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 <cmath>
#include <cstdio>
#include <algorithm>
??
#define MAX_N 300
const double PI=acos(-1.0);
??
struct Point{
????????double x,y;
????????Point(){
????????????????x=y=0;
????????}
????????bool operator<(const Point& p)const{
????????????????return this->x < p.x;
????????}
};
??
struct comp{
????????bool operator()(const Point &P1, const Point &P2){
????????????????return P1.x < P2.x;
????????}
????????bool operator()(const Point &P, double d){
????????????????return P.x < d;
????????}
????????bool operator()(double d, const Point &P){
????????????????return d < P.x;
????????}
};
??
Point points[MAX_N];
int N;
??
double distance(Point *a, Point *b){
????????return std::sqrt(pow(a->x-b->x,2)+pow(a->y-b->y,2));
}
??
double distance2(Point *a, Point *b){
????????return pow(a->x-b->x,2)+pow(a->y-b->y,2);
}
??
int solve(){
????????Point C, *P1, *P2, *P3;
????????int max=1;
????????for(int i=0;i<N;i++){
????????????????P1=&points[i];
????????????????Point *LB,*UB;
????????????????LB=std::lower_bound(points,points+N,P1->x-2.0001,comp());
????????????????UB=std::upper_bound(points,points+N,P1->x+2.0001,comp());
????????????????for(int j=(LB-points);j<(UB-points);j++){
????????????????????????if(i==j) continue;
????????????????????????P2=&points[j];
????????????????????????if(distance(P1,P2)>2.0) continue;
??
????????????????????????double alpha,beta;
????????????????????????alpha=acos(distance(P1,P2)/2.0);
????????????????????????if((P2->x-P1->x)>=0){
????????????????????????????????if((P2->y - P1->y)>=0){
????????????????????????????????????????beta=atan((P2->y - P1->y) / (P2->x - P1->x));
????????????????????????????????}else{
????????????????????????????????????????beta=atan((P2->y - P1->y) / (P2->x - P1->x))+2.0*PI;
????????????????????????????????}
????????????????????????}else{
????????????????????????????????beta=atan((P2->y - P1->y) / (P2->x - P1->x))+PI;
????????????????????????}
????????????????????????C.x=P1->x+cos(alpha+beta);
????????????????????????C.y=P1->y+sin(alpha+beta);
??
????????????????????????int enclosed=0;
????????????????????????Point *lb,*ub;
????????????????????????lb=std::lower_bound(points,points+N,C.x-1.0001,comp());
????????????????????????ub=std::upper_bound(points,points+N,C.x+1.0001,comp());
????????????????????????for(int k=(lb-points);k<(ub-points);k++){
????????????????????????????????if(k==i||k==j){
????????????????????????????????????????enclosed++;
????????????????????????????????????????continue;
????????????????????????????????}
????????????????????????????????P3=&points[k];
????????????????????????????????if(distance2(P3,&C)<1.0) enclosed++;
????????????????????????}
????????????????????????if(enclosed>max) max=enclosed;
????????????????????????if(max==N) return max;
????????????????}
????????}
????????return max;
}
??
int main(void){
????????while(scanf("%d",&N), N){
????????????????for(int i=0;i<N;i++){
????????????????????????scanf("%lf %lf", &points[i].x, &points[i].y);
????????????????}
????????????????std::sort(points,points+N);
????????????????printf("%d\n",solve());
????????}
????????return 0;
}

C++: #include<cstdio>
#include<algorithm>
#include<math.h>

using namespace std;
typedef pair<double, double> P;

int N;
double X[300], Y[300];
P d[300];

const double ESP = 0.0001;


bool near(double x, double y) {
  return abs(x - y) < ESP;
}

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

P minas(P p1, P p2) {
  return P(p1.first - p2.first, p1.second - p2.second);
}

P mul(P p, double c) {
  return P(p.first * c, p.second * c);
}

P ave(P p1, P p2) {
  return P((p1.first + p2.first) / 2, (p1.second + p2.second) / 2);
}

void printP(P p) {
  printf("P:(%f, %f)\n", p.first, p.second);
}

double square(P p) {
  return p.first * p.first + p.second * p.second;
}

double size(P p) {
  return sqrt(square(p));
}

int count(P p) {
  // printP(p);
  int counter = 0;
  for (int i = 0; i < N; i++) {
    if (square(minas(p, d[i])) <= 1.00001) counter++;
  }
  return counter;
}

int check(int i, int j) {
  if (square(minas(d[i], d[j])) > 4.0) return 0;
  int ans = 0;
  P average = ave(d[i], d[j]);
  P dir = mul(minas(d[i], d[j]), 0.5);
  // printP(average);
  double s = sqrt(1.0 - square(dir));
  // printf("!!!!!!%f\n",s);
  P dir1 = P((Y[i] - Y[j]), - (X[i] - X[j]));
  //printf("!!!!!!%f\n",size(dir1));
  double r = s / size(dir1);
  P d1 = add(average, mul(dir1, r));
  P d2 = add(average, mul(dir1, -r));
  return max(count(d1), count(d2));
}

int main() {
  while (1) {
    scanf("%d", &N);
    if (N == 0) break;
    for (int i = 0; i < N; i++) scanf("%lf%lf", &X[i], &Y[i]), d[i] = P(X[i], Y[i]);
    int ans = 0;
    for (int i = 0; i < N; i++)
      for (int j = 0; j < i; j++)
        ans = max(ans, check(i, j));
    printf("%d\n", ans);
  }
}

C++: #include <bits/stdc++.h>
using namespace std;
using ld = long double;
using P = complex<ld>;
const ld eps = 1e-6;

int main()
{
	int N;
	while (cin >> N, N) {
		vector<ld> x(N), y(N);
		vector<P> p(N);
		for (int i = 0; i < N; i++) {
			cin >> x[i] >> y[i];
			p[i] = P(x[i], y[i]);
		}
		vector<P> ko;
		for (int i = 0; i < N; i++) {
			for (int j = i + 1; j < N; j++) {
				if (norm(p[i] - p[j]) < 4.0 + eps) {
					auto d = p[i] - p[j];
					ko.push_back(p[j] + (d * (ld)0.5) + d * P(0, 1) / abs(d) * sqrtl(1 - norm(d * (ld)0.5)));
					ko.push_back(p[j] + (d * (ld)0.5) - d * P(0, 1) / abs(d) * sqrtl(1 - norm(d * (ld)0.5)));
				}
			}
		}
		int res = 0;
		for (auto c : ko) {
			int cnt = 0;
			for (int i = 0; i < N; i++) {
				if (norm(c - p[i]) < 1.0 + eps) {
					cnt++;
				}
			}
			res = max(res, cnt);
		}
		cout << res << endl;
	}
	return 0;
}