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 <stdio.h>
#include <string.h>
#include <algorithm>
#include <iostream>
#include <math.h>
#include <assert.h>
#include <vector>
#include <complex>

using namespace std;
typedef long long ll;
typedef complex<double> Point;
static const double EPS = 1e-9;
static const double PI = acos(-1.0);

#define REP(i, n) for (int i = 0; i < (int)(n); i++)
#define FOR(i, s, n) for (int i = (s); i < (int)(n); i++)
#define FOREQ(i, s, n) for (int i = (s); i <= (int)(n); i++)
#define FORIT(it, c) for (__typeof((c).begin())it = (c).begin(); it != (c).end(); it++)
#define DEC(i, s) for (int i = (s); i >= 0; i--)

#define SIZE(v) (int)((v).size())
#define MEMSET(v, h) memset((v), h, sizeof(v))
#define FIND(m, w) ((m).find(w) != (m).end())

int n;
int ans;
double R;
Point points[3000];
pair<double, int> radian[30000];
bool used[3000];

int calc(int s) {
  int ret = 1;
  int m = 0;
  REP(i, n) {
    double d = abs(points[s] - points[i]) / 2.0;
    if (i == s || d * 2.0 > R * 2) { continue; }
    Point center = (points[s] + points[i]) / 2.0;
    Point vect = (points[i] - points[s]) * Point(0, -1);
    vect /= abs(vect);
    double l = sqrt(R * R - d * d);
    Point c1 = center + vect * l;
    Point c2 = center - vect * l;
    double arg1 = arg(c1 - points[s]) - EPS;
    double arg2 = arg(c2 - points[s]) + EPS;
    radian[m++] = pair<double, int>(arg1 + 0 * PI, (i + 1));
    radian[m++] = pair<double, int>(arg2 + 0 * PI, -(i + 1));
  }
  if (m / 2 < ans - 1) { return 0; }
  sort(radian, radian + m);
  MEMSET(used, false);
  int lsum = 1;
  REP(iter, 2) {
    REP(i, m) {
      int index = abs(radian[i].second) - 1;
      bool start = radian[i].second > 0 ? true : false;
      if (used[index] && !start) {
        used[index] = false;
        lsum--;
      } else if (!used[index] && start) {
        used[index] = true;
        lsum++;
      }
      ret = max(ret, lsum);
    }
  }
  return ret;
}

int main() {
  R = 1.0;
  while (scanf("%d", &n), n) {
    REP(i, n) {
      double x, y;
      scanf("%lf %lf", &x, &y);
      points[i] = Point(x, y);
    }
    ans = 0;
    REP(i, n) {
      ans = max(ans, calc(i));
    }
    printf("It is possible to cover %d points.\n", ans);
  }
}

C++: #define _USE_MATH_DEFINES
#include <algorithm>
#include <cstdio>
#include <functional>
#include <iostream>
#include <cfloat>
#include <climits>
#include <cstring>
#include <cmath>
#include <map>
#include <queue>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <time.h>
#include <vector>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> i_i;
typedef pair<ll, int> ll_i;
typedef pair<double, int> d_i;
typedef pair<ll, ll> ll_ll;
typedef pair<double, double> d_d;
struct edge { int u, v; ll w; };

ll MOD = 1000000007;
ll _MOD = 1000000009;
double EPS = 1e-8;

int main() {
	for (;;) {
		int N; cin >> N;
		if (N == 0) break;
		vector<double> x(N), y(N);
		for (int i = 0; i < N; i++)
			cin >> x[i] >> y[i];
		int maxi = 1;
		for (int i = 0; i < N; i++)
			for (int j = 0; j < N; j++) {
				if (i == j) continue;
				double dx = x[j] - x[i], dy = y[j] - y[i];
				double d = sqrt(dx * dx + dy * dy);
				if (d > 2) continue;
				double h = sqrt(1 - d * d / 4);
				double ox = x[i] + dx / 2 - dy / d * h;
				double oy = y[i] + dy / 2 + dx / d * h;
				int cnt = 0;
				for (int k = 0; k < N; k++) {
					double dx = x[k] - ox, dy = y[k] - oy;
					double d2 = dx * dx + dy * dy;
					if (d <= 1 + EPS) cnt++;
				}
				maxi = max(maxi, cnt);
			}
		cout << maxi << endl;
	}
}

C++: #include <iostream> 
#include <string> 
#include <vector> 
#include <cmath> 
#include <algorithm> 
#include <cstdlib> 
#include <ctime> 
#include <cstdio> 
#include <functional> 
#include <set> 
#include <sstream> 


using namespace std; 

int main(){

	int n;
	const double eps=1e-6;
	while(true){

		cin>>n;

		if(n==0) return 0;

		int ans=0;
		vector< double > x(n),y(n);

		for(int i=0;i<n;i++) cin>>x[i]>>y[i];


		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				const double d=hypot(x[i]-x[j],y[i]-y[j]);
				if(d>2.0-eps) continue;

				const double mx=(x[i]+x[j])/2.0;
				const double my=(y[i]+y[j])/2.0;
				const double dvx=(x[i]-x[j])/d;
				const double dvy=(y[i]-y[j])/d;

				const double cx1=mx+dvy*sqrt(1.0-d*d/4.0);
				const double cy1=my-dvx*sqrt(1.0-d*d/4.0);

				const double cx2=mx-dvy*sqrt(1.0-d*d/4.0);
				const double cy2=my+dvx*sqrt(1.0-d*d/4.0);

				int res=2;
				for(int k=0;k<n;k++){
					if(k==i || k==j) continue;
//					if(hypot(x[k]-cx1,y[k]-cy1)<1.0+eps) res++;
					if((x[k]-cx1)*(x[k]-cx1)+(y[k]-cy1)*(y[k]-cy1)<1.0+eps) res++;
				}
				
				ans=max(ans,res);

				res=2;
				for(int k=0;k<n;k++){
					if(k==i || k==j) continue;
//					if(hypot(x[k]-cx2,y[k]-cy2)<1.0+eps) res++;
					if((x[k]-cx2)*(x[k]-cx2)+(y[k]-cy2)*(y[k]-cy2)<1.0+eps) res++;
				}

				ans=max(ans,res);

			}
		}

		cout<<ans<<endl;

	}

	return 0;
}

C++: #include<bits/stdc++.h>
#define eps 1e-9
using namespace std;
typedef complex<double> point;

double dist(point a, point b){
  return sqrt((a.real() - b.real())*(a.real() - b.real()) +
              (a.imag() - b.imag())*(a.imag() - b.imag()));
}

pair<point, point> norm(point p){
  double unit = 1 / sqrt(p.real()*p.real() + p.imag()*p.imag());
  return pair<point, point>(p*point(0, 1)/unit,
                            p*point(0, -1)/unit);
}

vector<point> calcC(point a, point b){
  point p = (b - a)/2.0;
  pair<point, point> n = norm(p);
  double unit = sqrt(1 - abs(p)*abs(p));
  n.first  *= unit;
  n.second *= unit;
  vector<point> res(2);
  res[0] = a + p + n.first;
  res[1] = a + p + n.second;
  return res;
}
int main(int argc, char *argv[]){
  int n;
  while(cin >> n, n){;
    point p[n];
    for (int i = 0; i < n; i++) {
      double real, imag;
      cin >> real >> imag;
      p[i] = complex<double>(real, imag);
    }
    int ans = 0, cnt = 0;
    for (int i = 0; i < n; i++) {
      for (int j = i + 1; j < n; j++) {
        point a = p[i], b = p[j];
        if(dist(a, b) > 2 - eps)continue;
        vector<point> center = calcC(a, b);
        for (int c = 0; c < 2; c++) {
          cnt = 0;
          // std::cout << center[c].real() << " " <<
          //   center[c].imag()<< std::endl;
          for (int k = 0; k < n; k++) {
            if(dist(center[c], p[k]) <= 1 - eps)cnt++;
          }
          ans = max(ans, cnt);
        }
      }
    }
    std::cout << ans << std::endl;
  }
  return 0;
}