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<math.h>
#include<algorithm>
#include<vector>
#include<queue>
#include<string>
#include<set>
 
using namespace std;
#define rep(i,n) for(int i=0;i<n;i++)
#define INF 1001001001
#define LLINF 1001001001001001001
#define mp make_pair
#define LLIandI pair<long long int , int>
#define ll long long




//xy平面上にN個の点が与えられる．
//半径1の円をxy平面 上で動かして，それらの点をなるべくたくさん囲むようにする．
//このとき，最 大でいくつの点を同時に囲めるかを答えなさい．
//ここで，ある円が点を「囲む」 とは，その点が円の内部または円周上にあるときをいう．

int main(){
    while(1){
        //1 <= N <= 300, 0.0 <= X <= 10.0, および0.0 <= Y <= 10.0と仮定してよい
        int N;
        scanf("%d",&N);
        if(N==0)break;
        float X[303],Y[303];
        rep(i,N){
            scanf("%f %f",&X[i],&Y[i]);
        }
        //二つの点を通る縁を全て考える
        //二つの縁の組み合わせが300*299/2通り、
        //それぞれについて2通りの縁が考えられる。
        //3つの点が円を通過することがないという条件から、重複はない。
        int ans=0;
        rep(i,N){
            rep(j,i){
                //円の中心を計算
                float nx=(X[i]+X[j])/2.0;//二等分点
                float ny=(Y[i]+Y[j])/2.0;
                float len1=sqrt(pow(X[i]-nx,2.0)+pow(Y[i]-ny,2.0));
                if(len1>1.0)continue;
                float ex=(Y[i]-ny)/len1;
                float ey=-(X[i]-nx)/len1;
                float len=sqrt(1.0-pow(len1,2.0));
                for(float swi=-1.0;swi<=1.0;swi+=2.0){
                    float cx=nx+ex*len*(swi);
                    float cy=ny+ey*len*(swi);
                    int sans=0;
                    //単位ベクトル
                    rep(k,N){
                        if(k==i||k==j){sans++;}
                        else{
                            if(sqrt(pow(X[k]-cx,2.0)+pow(Y[k]-cy,2.0))<1.0)sans++;
                        }
                    }
                    ans=max(ans,sans);
                }
            }
        }
        printf("%d\n",max(1,ans));
    }
	return 0;
}

C++: #include <cstdio>
#include <cstring>
#include <string>
#include <cmath>
#include <cassert>
#include <iostream>
#include <algorithm>
#include <stack>
#include <queue>
#include <vector>
#include <set>
#include <map>
#include <bitset>
#include <functional>
#include <numeric>
using namespace std;

#define repl(i,a,b) for(int i=(int)(a);i<(int)(b);i++)
#define rep(i,n) repl(i,0,n)
#define mp(a,b) make_pair((a),(b))
#define pb(a) push_back((a))
#define all(x) (x).begin(),(x).end()
#define dbg(x) cout<<#x"="<<((x))<<endl
#define fi first
#define se second

#define INF 2147483600
#define long long long
#define EPS 1e-8

// ???????°???°???????????????
inline double add(double a, double b){
  if(abs(a+b) < EPS*(abs(a) + abs(b))) return 0;
  return a+b;
}

struct Point{
  double x,y;
  Point() {}
  Point(double nx, double ny) : x(nx), y(ny) {}
  inline Point operator + (const Point & p){ return Point(add(x, p.x), add(y, p.y)); }
  inline Point operator - (const Point & p){ return Point(add(x,-p.x), add(y,-p.y)); }
  inline Point operator * (double d){ return Point(x*d, y*d); }
  inline double dot(const Point & p){ return add(x * p.x, y*p.y); }  //??????
  inline double det(const Point & p){ return add(x * p.y, -y*p.x); } //??????
  inline double dist(const Point & p){ return sqrt((x-p.x)*(x-p.x) + (y-p.y)*(y-p.y)); }
  inline bool operator < (const Point & p) const {
    if(x != p.x) return x < p.x;
    else return y < p.y;
  }
  inline bool operator == (const Point & p) const {
    return (add(x, -p.x)==0) && (add(y, -p.y)==0);
  }
  friend ostream& operator<<(ostream& os, const Point& p) {
    os << "[" << p.x << "," << p.y << "]";
    return os;
  }
};

pair<Point, Point> calcCenter(Point &p1, Point &p2){
  if(add(p1.x, -p2.x)==0){
    double y = (p1.y + p2.y)/2;
    double dx = sqrt(1 - (p1.y - y)*(p1.y - y));
    return mp(Point(p1.x+dx, y), Point(p1.x-dx,y));
  }
  // 2???????????´???????????? Ax + By + C = 0 , A != 0
  double A = 2*(p1.x - p2.x);
  double B = 2*(p1.y - p2.y);
  double C = - p1.x*p1.x + p2.x*p2.x - p1.y*p1.y + p2.y*p2.y;
  // x = py+q
  double p = -B/A; double q = -C/A;

  // p1, p2 ??¨????????¢???1?????? ay^2 + 2by + c = 0
  double a = 1 + p*p;
  double b = p*q - p1.x*p - p1.y;
  double c = q*q - 2*p1.x*q + p1.x*p1.x + p1.y*p1.y - 1;

  double rt = sqrt(b*b-a*c);
  double y1 = (-b+rt)/a;
  double y2 = (-b-rt)/a;

  return mp(Point(p*y1+q, y1), Point(p*y2+q, y2));
}

int main(){
  int n;
  while(scanf("%d", &n), n){
    vector<Point> vec;
    rep(i,n){
      double x,y;
      scanf("%lf %lf", &x, &y);
      vec.pb(Point(x,y));
    }

    int res=1;
    rep(i,n) repl(j,i+1,n) if(vec[i].dist(vec[j]) <= 2+EPS){
      pair<Point, Point> pr = calcCenter(vec[i], vec[j]);
//      cout << pr.fi << pr.se << endl;
      int cnt=0;
      rep(k,n) if(pr.fi.dist(vec[k]) <= 1+EPS) cnt++;
      res = max(res, cnt);
      cnt=0;
      rep(k,n) if(pr.se.dist(vec[k]) <= 1+EPS) cnt++;
      res = max(res, cnt);
    }

    printf("%d\n", res);

  }
  return 0;
}

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

using namespace std;

#define int long long
#define MKP make_pair
#define fi first
#define se second
typedef pair<int, int> pii;

constexpr int MOD = 1e9 + 7;
constexpr int MAX_N = 205;
constexpr int MAX_W = 10005;
constexpr int INF = 1LL << 60;

/* 基本要素 */

typedef double D;      // 座標値の型。doubleかlong doubleを想定
typedef complex<D> P;  // Point
typedef pair<P, P> L;  // Line
typedef vector<P> VP;
const D EPS = 1e-9;    // 許容誤差。問題によって変える
#define X real()
#define Y imag()
#define LE(n,m) ((n) < (m) + EPS)
#define GE(n,m) ((n) + EPS > (m))
#define EQ(n,m) (abs((n)-(m)) < EPS)

// 内積　dot(a,b) = |a||b|cosθ
D dot(P a, P b) {
  return (conj(a)*b).X;
}
// 外積　cross(a,b) = |a||b|sinθ
D cross(P a, P b) {
  return (conj(a)*b).Y;
}

VP crosspointCC(P a, D ar, P b, D br) {
  VP ps;
  P ab = b-a;
  D d = abs(ab);
  D crL = (norm(ab) + ar*ar - br*br) / (2*d);
  if (EQ(d, 0) || ar < abs(crL)) return ps;

  P abN = ab * P(0, sqrt(ar*ar - crL*crL) / d);
  P cp = a + crL/d * ab;
  ps.push_back(cp + abN);
  if (!EQ(norm(abN), 0)) ps.push_back(cp - abN);
  return ps;
}

int n;
vector<P> v;

bool solve() {
    v.clear();
    cin >> n;
    if (!n) return false;

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

    int res = 0;
    for (int i = 0; i <= n; ++i) {
        for (int j = 0; j <= n; ++j) {
            auto crs = crosspointCC(v[i], 1.0, v[j], 1.0);
            for (int k = 0; k < crs.size(); ++k) {
                int cnt = 0;
                for (int l = 0; l <= n; ++l) {
                    double dist = abs(crs[k] - v[l]);
                    if (dist <= 1.0 + EPS) cnt++;
                }

                res = max(res, cnt);
            }
        }
    }

    cout << res << endl;

    return true;
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    while (solve());
}