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 <math.h>
 
 
int main() {
	int point_num, ans;
	float point[300][2];   
	for(;;){
		scanf("%d\n", &point_num);
		if(point_num == 0)
			break;
		int i;
		float x, y;
		for (i = 0; i < point_num; i++)
			scanf("%f %f\n", &point[i][0], &point[i][1]);
		ans = 0;
		
		for (x = 1; x < 10; x += 0.001) {
			for (y = 1; y < 10; y += 0.001) {
				int count = 0;
				for (i = 0; i < point_num; i++) {
					if (1.0 >=  (x - point[i][0]) * (x - point[i][0]) + (y - point[i][1]) * (y - point[i][1]))
						count++;
				if (count > ans)
					ans = count;
				}
			}
		}
	printf("%d\n", ans);
	}
	return 0;
}

C++: #include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <queue>
#include <stack>
#include <set>

using namespace std;


int main(){
	int n;
	double x1,y1,xc,yc,dis,vx,vy,ex,ey;
	while(1){
		cin >> n;
		if(n==0)break;
		vector<double> x(n);
		vector<double> y(n);
		for(int i=0;i<n;i++){
			cin >> x[i] >> y[i];
		}
		if(n==1){
			cout << 1 << endl;
			continue;
		}
		int ans =0;
		for(int i=0;i<n;i++){
			for(int j=0;j<n;j++){
				if(i==j)break;
				xc = (x[i]+x[j])/2;
				yc = (y[i]+y[j])/2;
				dis = sqrt(1.0-(xc-x[i])*(xc-x[i])-(yc-y[i])*(yc-y[i]));
				vx = x[i]-x[j];
				vy = y[i]-y[j];
				ex = dis * vy * sqrt(vx*vx+vy*vy);
				ey = -dis* vx * sqrt(vx*vx+vy*vy);
				x1 = xc + ex;
				y1 = yc + ey;
				int count=0;
				for(int k=0;k<n;k++){
					if(k==i||k==j)count +=1;
					else{
						if((x1-x[k])*(x1-x[k])+(y1-y[k])*(y1-y[k])<=1.0001){
							count +=1;
						}
					}
				}
				ans = max(ans,count);
				x1 = xc - ex;
				y1 = yc - ey;
				count = 0;
				for(int k=0;k<n;k++){
					if(k==i||k==j)count +=1;
					else{
						if((x1-x[k])*(x1-x[k])+(y1-y[k])*(y1-y[k])<=1.0001){
							count +=1;
						}
					}
				}
				ans = max(ans,count);
			}
		}
		cout << ans << endl;
	}
	return 0;
}

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

#define P complex<double>
#define EPS 1.0e-8

int main(){
  int n;
  P p[301];

  while( cin>>n && n ){
    for( int i=0;i<n;i++ )
      cin >> p[i].real() >> p[i].imag();

    int ans=0;
    for( int i=0;i<n;i++ ){
      for( int j=0;j<n;j++ ){
	if( i==j ) 	  continue;
	if( abs( p[i]-p[j] )+EPS > 2.0 ) continue;

	P m  = (p[i]+p[j]) / 2.0; // p[i]Æp[j]Ì_
	P nn = (p[i]-p[j]) * P(0,1);       // @üxNg
	nn = nn / abs(nn);

	double x = abs( p[i]-p[j] )/2.0;
	P c = m + nn*sqrt( 1.0-x*x );        // ~ÌS

	int num=0;
	for( int k=0;k<n;k++ )
	  if( abs( c-p[k] ) < 1.0 + EPS  )
	    num++;
	ans = max( ans,num );
      }
    }
    cout << ans << endl;
  }

  return 0;
}

C++: #include<iostream>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<algorithm>
#include<functional>
#include<iomanip>
#include<queue>
#include<ciso646>
#include<random>
#include<map>
#include<set>
#include<complex>
using namespace std;
typedef long long ll;
const ll MOD = 1000000007;
const ll INF = (ll)1000000007 * 1000000007;
const double EPS = 1e-9;
typedef pair<int, int> P;
typedef unsigned int ui;
#define stop char nyaa;cin>>nyaa;
#define rep(i,n) for(int i=0;i<n;i++)
#define per(i,n) for(int i=n-1;i>=0;i--)
#define Rep(i,sta,n) for(int i=sta;i<n;i++)
#define rep1(i,n) for(int i=1;i<=n;i++)
#define per1(i,n) for(int i=n;i>=1;i--)
#define Rep1(i,sta,n) for(int i=sta;i<=n;i++)
typedef long double ld;
typedef complex<ld> Point;
const ld eps = 1e-8;
const ld pi = acos(-1.0);
bool eq(ld a, ld b) {
	return (abs(a - b) < eps);
}
ld dot(Point a, Point b) {
	return real(conj(a)*b);
}
ld cross(Point a, Point b) {
	return imag(conj(a)*b);
}
class Line {
public:
	Point a, b;
};
class Circle {
public:
	Point p;
	ld r;
};
vector<Point> is_cc(Circle c1, Circle c2) {
	vector<Point> res;
	ld d = abs(c1.p - c2.p);
	ld rc = (d*d + c1.r*c1.r - c2.r*c2.r) / (2 * d);
	ld dfr = c1.r*c1.r - rc * rc;
	if (abs(dfr) < eps)dfr = 0.0;
	if (dfr < 0.0)return res;
	ld rs = sqrt(dfr);
	Point diff = (c2.p - c1.p) / d;
	res.push_back(c1.p + diff * Point(rc, rs));
	if (dfr != 0.0)res.push_back(c1.p + diff * Point(rc, -rs));
	return res;
}
ld dis(Point x, Point y) {
	return sqrt(dot(y - x, y - x));
}
int main() {
	int n;
	while (cin >> n, n) {
		Point a[300];
		rep(i, n) {
			ld x, y; cin >> x >> y; a[i] = { x,y };
		}
		int ma = 0;
		rep(i, n) {
			Rep(j, i + 1, n) {
				Circle s = { a[i],1.0 };
				Circle t = { a[j],1.0 };
				vector<Point> v = is_cc(s, t);
				if (v.size() < 2)continue;
				int cnt = 0;
				rep(k, n) {
					if (dis(a[k], v[0]) <= 1.0+eps)cnt++;
				}
				ma = max(ma, cnt); cnt = 0;
				rep(k, n) {
					if (dis(a[k], v[1]) <= 1.0+eps)cnt++;
				}
				ma = max(ma, cnt);
			}
		}
		cout << ma << endl;
	}
	return 0;
}

C++: #include <bits/stdc++.h>
#define REP(i,n) for(int i=0; i<(int)(n); ++i)

using namespace std;

#define X real()
#define Y imag()

typedef double D;
D EPS = 1e-8;

typedef complex<D> P;
struct L{
    P a, b;
    P vec(){ return b - a; }
};
typedef vector<P> Pol;

int sig(D a) { return a < -EPS ? -1 : a > EPS ? 1 : 0; }
D dot(P a, P b){ return a.X * b.X + a.Y * b.Y; }
D cross(P a, P b){ return a.X * b.Y - a.Y * b.X; }

D area(Pol pol) {
    int n = pol.size();
    D res = 0;
    REP(i, n) {
        res += cross(pol[i], pol[ (i + 1) % n ]);
    }
    return abs(res) / 2.0;
}

// 2点a, bを通る半径rの円の中心
vector<P> tc2(P a, P b, D r){
    vector<P> res;

    D d = abs(b - a);
    if(sig(d - 2 * r) > 0) return res;

    P mid = 0.5 * (a + b);
    P dir = polar(sqrt(r * r - d * d / 4), arg(b - a) + M_PI / 2);
    res.push_back(mid + dir);
    res.push_back(mid - dir);
    return res;
}

int main(){
    int N;
    while(cin >> N && N > 0) {
        vector<P> ps(N);
        REP(i, N) {
            D x, y;
            cin >> x >> y;
            ps[i] = P(x, y);
        }

        int ans = 1;
        REP(i, N) REP(j, N) if(i != j) {
            vector<P> cv = tc2(ps[i], ps[j], 1.0);
            for(auto c : cv) {
                int sum = 0;
                REP(k, N) {
                    if(sig(abs(c - ps[k]) - 1.0) <= 0) {
                        sum++;
                    }
                }
                ans = max(ans, sum);
            }
        }
        cout << ans << endl;
    }
    return 0;
}