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<iostream>
#include<vector>
#include<complex>
#define EPS 0.0001
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 = 0;
		int n;
		scanf("%d", &n);
		if(n==0)break;
		if(n==1){
			printf("1\n");
			continue;
		}
		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 = 0; j < points.size(); ++j){
				if( i != j && dist( points[i], points[j] ) < 2+EPS ){
					int count_p = 0;
					int count_q = 0;
					elem d = sqrt(1-pow(dist(points[i], points[j])/2,2 ) );
					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;
					for(unsigned int k = 0; k < points.size(); ++k){
						if( dist(p, points[k]) < 1.0+EPS ){ count_p++; }
						if( dist(q, points[k]) < 1.0+EPS ){ count_q++; }
					}
					ans = max(ans, max(count_q, count_p) );
				}
			}
		}
		printf("%d\n", ans);
	}
	return 0;
}

C++: /**
 *  Circle and Points (PKU 1981)
 *  by Kenji Inoue, 2007-01-05
 **/

#include <iostream>
#include <complex>

using namespace std;

typedef complex<double> P;

#define EPS (1e-6)

P p[300];
P tmp;

int main() {
  int n;
  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=i+1; j<n; j++) {
        P v = p[i] - p[j];
        double absv = abs(v);
        if (absv > 2.0 + EPS) { continue; }

        P m = p[j] + (v / 2.0);
        P unv = (v * P(0, 1)) / absv; // unit normal vector
        P nv = sqrt(1.0 - absv*absv/4.0) * unv;
        for (int s=0; s<2; s++) {
          nv *= -1;
          P c = m + nv; // center point of the circle

          int cnt = 0;
          for (int k=0; k<n; k++) {
            // if (abs(p[k] - c) < 1.0 + EPS) { cnt++; } // too slow
            tmp = p[k] - c;
            if ( tmp.real()*tmp.real() + tmp.imag()*tmp.imag() < 1.0 + EPS ) { cnt++; }
          }
          if (cnt > ans) { ans = cnt; }
        }
      }
    }

    cout << ans << endl;
  }

  return 0;
}

C++: #include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <cstring>
#include <map>
#include <utility>
#include <set>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>
#include <sstream>
#include <complex>
#include <stack>
#include <queue>

#define FOR(i,b,n) for(int i=b;i<n;i++)
#define RFOR(i,b,n) for(int i=n-1;i>=b;i--)
#define CLR(mat) memset(mat, 0, sizeof(mat))
#define NCLR(mat) memset(mat, -1, sizeof(mat))

#define EPS 0.0001
#define EQ(a,b) (abs((a)-(b)) < EPS)
#define EQV(a,b) ( EQ((a).real(), (b).real()) && EQ((a).imag(), (b).imag()) )
#define PI (acos(-1.0))
#define POSARG(a) (arg((a)) > 0.0? arg((a)) : 2.0*PI + arg((a)))

using namespace std;

typedef complex<double> P;

int n;

int solve(vector<P> vc)
{	
	int res = 1;
	
	FOR(j, 0, n)
		FOR(i, j+1, n)
		{
			int tmpp = 2;
			int tmpn = 2;
			
			if(abs(vc[j]-vc[i]) > 2.0+EPS)
				continue;
			
			double cs = abs(vc[j]-vc[i])/2.0, sn = sqrt(1.0-cs*cs);		//線分IJと円の中心の距離
			
			P mid = (vc[j]+vc[i])/2.0;//点IとJの中点
			
			P nml1 = (vc[j]-vc[i])*P(0.0, sn), nml2 = (vc[j]-vc[i])*P(0.0, -sn);
														//線分IJの法線ベクトル(長さsn)
			
			P c1 = mid + nml1, c2 = mid + nml2;//円の中心2つ
						
						
			FOR(k, 0, n)
			{	
				if( k != j && k != i )
				{
					if( abs(vc[k]-c1) < 1.0+EPS)
					{
						tmpp++;
					}
					
					if( abs(vc[k]-c2) < 1.0+EPS)
					{
						tmpn++;
					}
				}
			}
			
			res = max(res, max(tmpp, tmpn));
		}
		
	return res;
}

int main()
{
	while(cin >> n, n)
	{
		vector<P> vc = vector<P>();
		
		FOR(i, 0, n)
		{
			double x, y;
			cin >> x >> y;
			vc.push_back(P(x,y));
		}

		cout << solve(vc) << endl;
	}
		
	
	return 0;
}

C++: #include <cstdio>
#include <complex>
#include <vector>
using namespace std;
typedef complex<double> P;
const double PI=3.141592653589793238462643383279;

bool inC(P p, P c, double r){return abs(p-c)<r+1e-9;}

P getE(P p){return p/abs(p);}

double sq(double x){ return  x*x;}

int main(){
    int n,i,j,k,ans,cnt;
    double x, y;
    while(scanf("%d",&n)){
        if(n==0) break;
        ans=1;
        vector<P> ps(n);
        for(i=0; i<n; i++){
            scanf("%lf%lf",&x,&y);
            ps[i] = P(x,y);
        }
        for(i=0;i<n;i++)
            for(j=i+1;j<n;j++)if(abs(ps[i]-ps[j])<2){
                double d=sqrt(1-sq(abs(ps[j]-ps[i])*0.5));
                P c1=0.5*(ps[i]+ps[j])+getE(getE(ps[j]-ps[i]))*polar(d,PI/2);
                P c2=0.5*(ps[i]+ps[j])+getE(getE(ps[j]-ps[i]))*polar(d,-PI/2);
                cnt=0;
                for(k=0;k<n;k++)if(inC(ps[k],c1,1.0))cnt++;
                if(cnt>ans)ans=cnt;
                cnt=0;
                for(k=0;k<n;k++)if(inC(ps[k],c2,1.0))cnt++;
                if(cnt>ans)ans=cnt;
            }
       // if(ans==0) ans++;
        printf("%d\n",ans);
    }
    return 0;
}