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

using namespace std;

typedef pair<double, double> P;

int N;
P pos[300];

double calcabs(P a, P b){
  return sqrt(pow(a.first - b.first, 2) + pow(a.second - b.second, 2));
}

int main(){
  while(1){
    cin >> N;
    if(N == 0)break;
    for(int i = 0; i < N; i++)cin >> pos[i].first >> pos[i].second;
    
    int ans = 1;

    for(int i = 0; i < N; i++){
    for(int j = i + 1; j < N; j++){
	if(calcabs(pos[i], pos[j]) > 2)continue;
	
	P center[2];
	//calc the center of the circle which we put on the selected 2 points.
	
	center[0].first = (pos[i].first + pos[j].first) / 2.0 + 
	  sqrt(1.0 - pow(calcabs(pos[i], pos[j]) , 2)/ 4.0) * 
	  1.0 / calcabs(pos[i], pos[j]) * (-1 * (pos[i].second - pos[j].second));
	
	center[0].second = (pos[i].second + pos[j].second) / 2.0 + 
	  sqrt(1.0 - pow(calcabs(pos[i], pos[j]), 2) / 4.0) * 
	  1.0 / calcabs(pos[i], pos[j]) *  (pos[i].first - pos[j].first);
	
	center[1].first = (pos[i].first + pos[j].first) / 2.0 - 
	  sqrt(1.0 - pow(calcabs(pos[i], pos[j]), 2) / 4.0) * 
	  1.0 / calcabs(pos[i], pos[j]) * (-1 * (pos[i].second - pos[j].second));
	
	center[1].second = (pos[i].second + pos[j].second) / 2.0 - 
	  sqrt(1.0 - pow(calcabs(pos[i], pos[j]), 2) / 4.0) * 
	  1.0 / calcabs(pos[i], pos[j]) *  (pos[i].first - pos[j].first);
	
	int cnt = 0;
	for(int k = 0; k < N; k++){
	  if(calcabs(center[0], pos[k]) <= 1.0001)cnt++;
	}
	ans = max(ans, cnt);
	
	cnt = 0;
	for(int k = 0; k < N; k++){
	  if(calcabs(center[1], pos[k]) <= 1.0001)cnt++;
	}
	ans = max(ans, cnt);
	
    }
    }
    cout << ans << endl;
  }
  return 0;
}

C++: #include <stdio.h>
#include <algorithm>
#include <iostream>
#include <stack>
#include <vector>
#include <math.h>
#include <queue>

using namespace std;

struct Point{
    double x,y;
};

void func(int N){

    Point point[N];
    double x_min = 100,y_min = 100,x_max = -1,y_max = -1;

    for(int i = 0; i < N; i++){
        scanf("%lf %lf",&point[i].x,&point[i].y);
        x_min = min(point[i].x,x_min);
        y_min = min(point[i].y,y_min);
        x_max = max(point[i].x,x_max);
        y_max = max(point[i].y,y_max);
    }

    if(x_min > 1)x_min -= 0.05;
    if(y_min > 1)y_min -= 0.05;
    if(x_max < 9)x_max += 0.05;
    if(y_max < 9)y_max += 0.05;

    int count,max_count = 0;

    for(double y = y_min; y <= y_max; y += 0.011){
        for(double x = x_min; x <= x_max; x += 0.011){
            count = 0;
            for(int i = 0; i < N; i++){
                if(sqrt((x-point[i].x)*(x-point[i].x)+(y-point[i].y)*(y-point[i].y)) <= 1.0)count++;
            }
            max_count = max(max_count,count);
        }
    }

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


int main(){

    int N;

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

        func(N);
    }


    return 0;
}

C++: #include <bits/stdc++.h>

using namespace std;

const double r = 1.0;
const double PI = acos(-1.0);

struct Point {
    double x, y;
    Point() {}
    Point(double x, double y) : x(x), y(y) {}
};
inline double Distance2(const Point &a, const Point &b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

int main()
{
    cin.tie(0);
    ios::sync_with_stdio(false);

    int n;
    std::random_device rnd;
    std::mt19937 mt(rnd());
    std::uniform_real_distribution<> rand100(0.0, 10.0);

    while (cin >> n, n) {
        // Input
        vector<Point> p(n);
        for (int i = 0; i < n; ++i)
            cin >> p[i].x >> p[i].y;

        // Solve
        int ans = 1;
        auto start = std::chrono::high_resolution_clock::now();
        auto end = std::chrono::high_resolution_clock::now();
        auto take_time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

        while (true) {
            end = std::chrono::high_resolution_clock::now();
            take_time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
            if (320 < take_time.count())
                break;

            for (int i = 0; i < 10; ++i) {
                Point c(rand100(mt), rand100(mt));
                int num = 0;
                for (int j = 0; j < n; ++j) {
                    double d = Distance2(p[j], c);
                    if (d <= 1.0)
                        ++num;
                }
                ans = max(ans, num);
            }
        }

        // Output
        cout << ans << endl;
    }

    return 0;
}

C++: #include <iostream>
#include <complex>
#include <cmath>
#include <algorithm>
using namespace std;
typedef complex<double> C;
#define EPS (1e-6)

int n;
C c[301];

void solve() {
  for(int i=0;i<n;i++)
  cin>>c[i].real()>>c[i].imag();
  
  int res = 0;
  for (int i = 0; i < n; i += 1) {
    for (int j = i+1; j < n; j += 1) {
      C v = c[j] - c[i];
      double absv = abs(v);
      if (absv > 2.0 + EPS) continue;
      // _
      C m = c[i] + (v / 2.0);
      // ~ÜÅÌ£
      double dist = sqrt(1.0 - absv * absv / 4.0);
      // PÊ@üxNg
      C unv  = (v * C(0, 1)) / absv;
      // S
      C nv = unv * dist;
      for (int k = 0; k < 2; k += 1) {          
        int cnt = 0;
        nv *= -1;
        C tc = m + nv;
        for (int l = 0; l < n; l += 1) {
          if (abs(tc - c[l]) < 1.0 + EPS) {
            cnt++;
          }
        }
        res = max(res, cnt);
      }
    }
  }
  cout<<res<<endl;
  //cout<<"+++++"<<endl;
}

int main (int argc, char const* argv[]) {
  while (cin>>n, n) {
    solve();
  }
  return 0;
}