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

using namespace std;
typedef complex<double> Complex;

static const double E = 0.000001;

inline bool collide_circle_circle(const Complex p1, const double r1, const Complex p2, const double r2) {
  assert(r1 >= 0.0 && r2 >= 0.0);
  Complex dif = p1 - p2;
  return dif.real() * dif.real() + dif.imag() * dif.imag() <= (r1 + r2) * (r1 + r2) + E;
}

// include circle inside
inline bool collide_circle_point(const Complex p1, const double r1, const Complex p2) {
  return collide_circle_circle(p1, r1, p2, 0.0);
}

inline void intersection_circle_circle(const Complex p1, const double r1, const Complex p2, const double r2, int *size, Complex *ret1, Complex *ret2) {
  assert(p1 != p2);
  assert(ret1 != NULL && ret2 != NULL);
  if (!collide_circle_circle(p1, r1, p2, r2)) {
    if (size != NULL) { *size = 0; }
    *ret1 = *ret2 = Complex(NAN, NAN);
    return;
  }
  Complex center = (p1 * r2 + p2 * r1) / (r1 + r2);
  Complex vect = (center - p1) * Complex(0.0, 1.0);
  if (abs(vect) < E) {
    *ret1 = *ret2 = center;
    if (size != NULL) { *size = 1; }
    return;
  }
  vect *= sqrt(r1 * r1 - norm(vect)) / abs(vect);
  *ret1 = center + vect;
  *ret2 = center + vect * Complex(-1.0, 0.0);
  if (size != NULL) { *size = 2; }
}

Complex point[301];

int main() {
  int n;
  while (scanf("%d", &n), n != 0) {
    for (int i = 0; i < n; i++) {
      double a, b;
      scanf("%lf %lf", &a, &b);
      point[i] = Complex(a, b);
    }
    int max_ans = 1;
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < i; j++) {
        if (!collide_circle_circle(point[i], 1.0, point[j], 1.0)) { continue; }
        int size;
        Complex m[2];
        intersection_circle_circle(point[i], 1.0, point[j], 1.0, &size, &m[0], &m[1]);
        for (int s = 0; s < size; s++) {
          int ans = 0;
          for (int k = 0; k < n; k++) {
            if (collide_circle_point(m[s], 1.0, point[k])) {
              ans++; }
          }
          max_ans = max(max_ans, ans);
        }
      }
    }
    printf("%d\n", max_ans);
  }
}

C++: #include<iostream>
#include<vector>
#include<complex>
#define EPS 1.0e-10
using namespace std;
typedef double elem;
typedef complex<elem> point,vec;
inline elem dist(const point &a, const point &b){return abs(a-b);}
inline vec normr(const vec &v){return v*vec(0,-1);}
inline vec norml(const vec &v){return v*vec(0,1);}
inline vec unormr(const vec &v){return normr(v) / abs(v);}
inline vec unorml(const vec &v){return norml(v) / abs(v);}
int main(){
	while(true){
		int ans = 1;
		int n;
		scanf("%d", &n);
		if(n==0)break;
		vector< point > points;
		for(int i=0;i<n;++i){
			elem x, y;
			scanf("%lf%lf", &x, &y);
			points.push_back( point(x,y) );
		}
		for(unsigned int i = 0; i < points.size(); ++i){
			for(unsigned int j = i + 1; j < points.size(); ++j){
				if( dist( points[i], points[j] ) <= 2 ){
					int count_p = 0;
					int count_q = 0;
					elem d = (elem)sqrt( 1.0 - pow( dist(points[i], points[j]) / 2.0, 2.0 ) );
					point p,q;
					vec ab = points[j] - points[i];
					point t = 0.5 * ab + points[i];  
					p = d * unormr(ab) + t;
					q = d * unorml(ab) + t;
					//cout << i << "," << j << endl;
					for(unsigned int k = 0; k < points.size(); ++k){
						//cout << dist(p,points[k]) << " " << dist(q,points[k]) << endl;

						if( dist(p, points[k]) <= 1+EPS ){ count_p++; }
						if( dist(q, points[k]) <= 1+EPS ){ count_q++; }
					}
					ans = max(ans, max(count_q, count_p) );
				}
			}
		}
		printf("%d\n", ans);
	}
	return 0;
}