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 <string>
#include <vector>
#include <sstream>
#include <iostream>
#include <algorithm>
#include <map>
#include <list>
#include <set>
#include <numeric>
#include <queue>
#include <stack>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <cstring>
#include <climits>
#include <cfloat>
#include <ctime>
#include <complex>
#include <cassert>
#include <tuple>
using namespace std;
typedef long long LL;
typedef pair<int,int> P;

int N;
double x[300];
double y[300];
inline double dist2(double a, double b)
{
	return a*a+b*b;
}
void exec()
{
	cin >> N;
	if(N==0){
		exit(0);
	}
	for(int i=0;i<N;i++){
		cin >> x[i] >> y[i];
	}
	int ret=1;
	for(int i=0;i<N;i++){
		for(int j=i+1;j<N;j++){
			double dx=(x[i]-x[j])/2;
			double dy=(y[i]-y[j])/2;
			double r2=dx*dx+dy*dy;
			if(r2>1){
				continue;
			}
			double R=sqrt(1/r2-1);
			double px=x[j]+dx-R*dy;
			double py=y[j]+dy+R*dx;
			double qx=x[j]+dx+R*dy;
			double qy=y[j]+dy-R*dx;
			int sp=0,sq=0;
			for(int k=0;k<N;k++){
				if(k==i||k==j){
					sp++;sq++;
					continue;
				}
				if(dist2(x[k]-px,y[k]-py)<=1){
					sp++;
				}
				if(dist2(x[k]-qx,y[k]-qy)<=1){
					sq++;
				}
			}
			ret=max(ret,sp);
			ret=max(ret,sq);
		}
	}
	cout << ret << endl;
}
int main() {
	while(1)
		exec();
}

C++: #include<iostream>
#include<complex>
#include<algorithm>
#include<cmath>
#define EPS 1e-8
using namespace std;

typedef complex<double> P;

P center(P a,P b){
  P m = (a+b)/2.0;
  P n = (a-b)*P(0,1);
  n = n/abs(n);

  double x = abs(a-b)/2.0;
  return m + n*sqrt(1.0 - x*x);
}

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

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

    int ans = 1;
    for(int i=0;i<n;i++){
      for(int j=0;j<n;j++){
	if(i!=j){
	  if(abs(p[i]-p[j])+EPS>2.0)continue;
	  P c = center(p[i],p[j]);
	  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;
  }
}

C++: #include "bits/stdc++.h"
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
const int INF = 1e9;
const ll LINF = 1e18;
template<class S,class T> ostream& operator << (ostream& out,const pair<S,T>& o){ out << "(" << o.first << "," << o.second << ")"; return out; }
template<class T> ostream& operator << (ostream& out,const vector<T> V){ for(int i = 0; i < V.size(); i++){ out << V[i]; if(i!=V.size()-1) out << " ";} return out; }
template<class T> ostream& operator << (ostream& out,const vector<vector<T> > Mat){ for(int i = 0; i < Mat.size(); i++) { if(i != 0) out << endl; out << Mat[i];} return out; }
template<class S,class T> ostream& operator << (ostream& out,const map<S,T> mp){ out << "{ "; for(auto it = mp.begin(); it != mp.end(); it++){ out << it->first << ":" << it->second; if(mp.size()-1 != distance(mp.begin(),it)) out << ", "; } out << " }"; return out; }

/*
 <url:http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1132>
 問題文============================================================
 =================================================================
 解説=============================================================
 ================================================================
 */

typedef complex<double> Point;
int n;
int solve(){
    int res = 1;
    vector<Point> ps(n);
    for(int i = 0; i < n;i++){
        double x,y; cin >> x >> y;
        ps[i] = Point(x,y);
    }
    for(int i = 0; i < n;i++){
        for(int j = i+1; j < n;j++){
            double dist = abs(ps[i]-ps[j]);
            if(dist > 2) continue;
            Point mp = (ps[i]+ps[j])/2.;
            Point dp = ps[i]-ps[j];
            double x = sqrt(1-dist*dist/4);
            Point cp = Point(-dp.imag()*x/dist+mp.real(),dp.real()*x/dist+mp.imag());
            int cnt = 0;
            for(int k = 0; k < n;k++){
                if(abs(cp-ps[k]) <= 1.0) cnt++;
            }
            res = max(res,cnt);
        }
    }
    return res;
}
int main(void) {
    cin.tie(0); ios_base::sync_with_stdio(false);
    while(cin >> n,n){
        cout << solve() << endl;
    }
    return 0;
}

C++: #include <iostream>
#include <cstdio>
#include <vector>
#include <algorithm>
#include <cmath>

using namespace std;

// ??????????????§?????????????????§???3????????????????????????????±?????????????????????¨?????????

long double pow2(long double d){
	return d * d;
}

struct Vector2{
	long double x;
	long double y;
};

Vector2 P[500];

int main(){
	int n;
	while(true){
		cin >> n;
		if(n == 0){
			break;
		}
		for(int i = 0; i < n; i++){
			cin >> P[i].x >> P[i].y;
		}
		long double k, L_2, a, b;
		
		int max = 0;
		
		for(int num1 = 0; num1 < n; num1++){
			for(int num2 = num1 + 1; num2 < n; num2++){
				L_2 = pow2(P[num1].x - P[num2].x) + pow2(P[num1].y - P[num2].y);
				if(L_2 >= 4){
					continue;
				}
				k = sqrt((4 - L_2) / L_2);
				
				a = (k * (P[num1].y - P[num2].y) + P[num1].x + P[num2].x) / 2;
				b = (k * (P[num2].x - P[num1].x) + P[num1].y + P[num2].y) / 2;
				
				int count;
				
				count = 0;
				
				for(int i = 0; i < n; i++){
					if(i == num1 || i == num2){
						count++;
						continue;
					}
					if(pow2(P[i].x - a) + pow2(P[i].y - b) <= 1){
						count++;
					}
				}
				
				if(max < count){
					max = count;
				}
				
				a = (k * (P[num2].y - P[num1].y) + P[num1].x + P[num2].x) / 2;
				b = (k * (P[num1].x - P[num2].x) + P[num1].y + P[num2].y) / 2;
				
				count = 0;
				
				for(int i = 0; i < n; i++){
					if(i == num1 || i == num2){
						count++;
						continue;
					}
					if(pow2(P[i].x - a) + pow2(P[i].y - b) <= 1){
						count++;
					}
				}
				
				if(max < count){
					max = count;
				}
			}
		}
		printf("%d\n", max);
	}
	return 0;
}

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

int main() {
  const complex<double> I(0,1);
  complex<double> d2(2,0);
  for(int tci = 0; ; tci++) {
    int n; scanf("%d", &n);
    if(!n) break;
    complex<double> *pts = new complex<double>[n];
    for(int i = 0; i < n; i++) {
      scanf("%lf%lf", &pts[i].real(), &pts[i].imag());
    }
    int maxcnt = 0;
    for(int i = 0; i < n; i++) {
      for(int j = 0; j < n; j++) {
        if(i==j) continue;
        complex<double> d = (pts[i] - pts[j]) / d2;
        double absd = abs(d);
        complex<double> k = d*I/absd * sqrt(1-absd*absd);
        complex<double> c = (pts[i] + pts[j]) / d2 + k;
        int cnt = 0;
        for(int k = 0; k < n; k++) {
          cnt += (abs(pts[k]-c) <= 1 + 1e-9);
        }
        maxcnt = max(maxcnt, cnt);
      }
    }
    delete[] pts;
    printf("%d\n", maxcnt);
  }
  return 0;
}