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 <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <complex>
#include <string>
#include <sstream>
#include <algorithm>
#include <numeric>
#include <vector>
#include <queue>
#include <stack>
#include <functional>
#include <iostream>
#include <map>
#include <set>
using namespace std;
typedef long long ll;
typedef vector<int> vi;
typedef vector<ll> vll;
#define pu push
#define pb push_back
#define mp make_pair
#define eps 1e-7
#define INF 1000000000
#define sz(x) ((int)(x).size())
#define sq(x) ((x)*(x))
#define fi first
#define sec second
#define SORT(x) sort((x).begin(),(x).end())
#define all(x) (x).begin(),(x).end()
#define rep(i,n) for(int (i)=0;(i)<(int)(n);(i)++)
#define repn(i,a,n) for(int (i)=(a);(i)<(int)(n);(i)++)
#define EQ(a,b) (abs((a)-(b))<eps)
double add(double a,double b){
	if(abs(a+b)<eps*(abs(a)+abs(b)))return 0;
	return a+b;
}
struct P{
	double x,y;
	P() {}
	P(double x,double y) : x(x),y(y){}
	P operator + (P p){
		return P(add(x,p.x),add(y,p.y));
	}
	P operator - (P p){
		return P(add(x,-p.x),add(y,-p.y));
	}
	P operator * (double d){
		return P(x*d,y*d);
	}
	double dot(P p){
		return add(x*p.x,y*p.y);
	}
	double det(P p){
		return add(x*p.y,-y*p.x);
	}
	double norm(){
		return sqrt(x*x+y*y);
	}
	double dist2(P p){
		return sq(x-p.x)+sq(y-p.y);
	}
	P vert(){
		return P(y,-x);
	}
	void dump(){
		printf("%.12f %.12f\n",x,y);
		return;
	}
};
int N;
P pts[305];
int solve(){
	scanf("%d",&N);
	if(N==0)return 1;
	for(int i=0;i<N;i++){
		scanf("%lf %lf",&pts[i].x,&pts[i].y);
	}
	int ans = 0;
	for(int i=0;i<N;i++){
		for(int j=i+1;j<N;j++){
			P m = (pts[i]+pts[j])*0.5;
			P a = (pts[i]-pts[j]);
			double na = a.norm();
			if(na>2.0)continue;
			a = (a.vert())*(sqrt(1-sq(na/2.0))/na);
			P b = m+a;
			P c = m-a;
			/*cout << i << ' ' << j << endl;
			m.dump();
			a.dump();
			b.dump();
			c.dump();
			cout << endl;*/
			int cnt = 0;
			for(int k=0;k<N;k++){
				//if(1.0<pts[k].dist2(b)&&pts[k].dist2(b)<=1.0+eps)printf("%d:%.12f\n",k,pts[k].dist2(b));
				if(pts[k].dist2(b)<=1.0+eps)cnt++;
			}
			ans = max(ans,cnt);
			cnt = 0;
			for(int k=0;k<N;k++){
				//if(1.0<pts[k].dist2(c)&&pts[k].dist2(c)<=1.0+eps)printf("%d:%.12f\n",k,pts[k].dist2(c));
				if(pts[k].dist2(c)<=1.0+eps)cnt++;
			}
			ans = max(ans,cnt);
		}
	}
	printf("%d\n",ans);
	return 0;
}
int main(){
	while(!solve()){}
	return 0;
}

C++: #include <iostream>
#include <complex>
#include <vector>
#include <cmath>
#include <utility>
#include <algorithm>
using namespace std;


typedef complex<double> Point;
typedef pair<Point, Point> Line;
#define EPS 1e-9
#define INF 1000000000

double dis(Point a, Point b){ return abs(a - b); }
double dot(Point a, Point b){ return a.real() * b.real() + a.imag() * b.imag(); }
double cross(Point a, Point b){ return a.real() * b.imag() - a.imag() * b.real(); }
int ccw(Point a, Point b, Point c){
  if(cross(b - a, c - a) > EPS) return 1;
  if(cross(b - a, c - a) < -EPS) return -1;
  if(dot(b - a, c - a) < -EPS) return 2;
  if(dis(a, b) + EPS < dis(c, a)) return -2;
  return 0;
}
int ccw(Line a, Point b){ return ccw(a.first, a.second, b); }

bool SegmentCross(Line a, Line b){
  return ccw(a, b.first) * ccw(a, b.second) == -1
    && ccw(b, a.first) * ccw(b, a.second) == - 1;
}

Point MiddlePoint(Point a, Point b){
  double x = (a.real() + b.real()) / 2;
  double y = (a.imag() + b.imag()) / 2;
  return Point(x, y);
}

Line Equidistant(Point a, Point b, double r = -1){
  Point m = MiddlePoint(a, b);
  double m_dis = dis(m, a);
  double slope = (b.real() - a.real()) / (a.imag() - b.imag());
  
  if(r <= EPS){ r = m_dis * 2; }
  
  if(r <= m_dis + EPS){ return make_pair(Point(0, 0), Point(0, 0)); }

  double target_dis = r * r - m_dis * m_dis;
  double x = sqrt(target_dis / (1 + slope * slope));
  double y = slope * x;
  
  Point s(m.real() + x, m.imag() + y);
  Point t(m.real() - x, m.imag() - y);
  
  return make_pair(s, t);
}


int main(void){
  //double a, b, c, d;
  //cin >> a >> b >> c >> d;
  //Equidistant(Point(a, b), Point(c, d));

 
  while(true){
    int n;
    cin >> n;
    if(n == 0){ break; }

    double x, y;
    vector<Point> P;
    for(int i = 0; i < n; i++){
      cin >> x >> y;
      P.push_back(Point(x, y));
    }

    int ans = 0;
    for(int i = 0; i < n; i++){
      for(int j = i + 1; j < n; j++){
        Line tmpLine = Equidistant(P[i], P[j], 1 - EPS);
        if(tmpLine.first == tmpLine.second){ continue; }

        Point center[] = {tmpLine.first, tmpLine.second};
        for(int l = 0; l < 2; l++){
          //printf("(%f, %f)(%f, %f) -> (%f, %f)\n", P[i].real(), P[i].imag(), P[j].real(), P[j].imag(),
          //       center[l].real(), center[l].imag());
          int count = 0;
          for(int k = 0; k < n; k++){
            //printf("  (%f, %f) %f\n", P[k].real(), P[k].imag(), dis(center[l], P[k]));
            if(dis(center[l], P[k]) <= 1.0 + EPS){
              ++count;
            }
          }
          
          ans = max(ans, count);
        }
      }
    }

    // ??????
    cout << ans << endl;
     
  }

  return 0; 
}

C++: #include <bits/stdc++.h>
using namespace std;
#define EPS (1e-4)
typedef complex<double> P;

int main() {
  int n;
  while (cin >> n, n) {
    vector<P> a(n);
    for (auto &e : a) {
      double x, y;
      cin >> x >> y;
      e.real(x);
      e.imag(y);
    }
    int m = 0;
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++) {
        if (abs(a[j] - a[i]) < 2.0 + EPS) {
          P u = (a[j] - a[i]) / 2.0;
          P v = P(-u.imag(), u.real()) / abs(u) * sqrt(1.0 - pow(abs(u), 2));
          vector<P> o = {a[i] + u + v, a[i] + u - v};
          for (auto e : o) {
            int t = 2;
            for (int k = 0; k < n; k++) {
              if (k != i && k != j) {
                t += abs(a[k] - e) < 1.0 + EPS;
              }
            }
            m = max(m, t);
          }
        }
      }
    }
    cout << m << endl;
  }
}