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 <bits/stdc++.h>
using namespace std;
#define rep(i,n) for(int i=0;i<n;++i)
typedef complex<double> Point;
typedef vector<Point> VP;
#define X real()
#define Y imag()

// ???a??¨???b?????????????????????r????????????????????????
VP circlesPointsRadius(Point a, Point b, double r) {
  VP cs;
  Point abH = (b-a)*0.5;
  double d = abs(abH);
  if (d == 0 || d > r) return cs;
  double dN = sqrt(r*r - d*d);
  Point n = abH * Point(0,1) * (dN / d);
  cs.push_back(a + abH + n);
  if (dN > 0) cs.push_back(a + abH - n);
  return cs;
}

int main(void){
    int N;
    while(cin>>N,N){
        VP p(N);
        int res=1;
        rep(i,N){
            double x,y;
            cin>>x>>y;
            p[i] = {x,y};
        }
        rep(i,N)rep(j,i){
            VP cs = circlesPointsRadius(p[i],p[j],1.0);
            for(Point c:cs){
                int cnt=0;
                rep(k,N) if(abs(c-p[k]) <= 1.0001)cnt++;
                res=max(res,cnt);
            }
        }
        cout<<res<<endl;
    }

    return 0;
}

C++: #include <iostream>
#include <complex>
#include <cmath>
#define EPS 10e-6
using namespace std;
typedef complex<double> xy_t;
int n;
xy_t points[301];
int count_in(xy_t center)
{
  int ans=0;
  for(int i=0;i<n;i++)
  {
    if (abs(center-points[i])<=1+EPS)
      ans++;
  }
  return ans;
}
xy_t ppdcl(xy_t a)
{
  return xy_t(a.imag(),-a.real());
}
xy_t delta_v(xy_t C1,xy_t C2)
{
  xy_t half=(C1-C2)*0.5;
  double l=abs(half);
  return ppdcl(half)*sqrt(1-l*l)*(1/l);
}
int main()
{
  while(1)
  {
    cin>>n;
    if (n==0) break;
    double x,y;
    for(int i=0;i<n;i++)
    {
      cin>>x>>y;
      points[i]=xy_t(x,y);
    }
    int maxn=-1;
    for(int i=0;i<n;i++)
    {
      for(int j=i+1;j<n;j++)
      {
        if(abs(points[i]-points[j])<=4+EPS)
        {
          xy_t mid_p=(points[i]+points[j])*0.5;
          maxn=max(maxn,count_in(mid_p+delta_v(points[i],points[j])));
          maxn=max(maxn,count_in(mid_p-delta_v(points[i],points[j])));
        }
      }
    }
    cout<<maxn<<endl;
  }

}

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

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

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

inline 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=0;
        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;
}

C++: #include <iostream>
#include <string>
#include <set>
#include <cmath>
#include <algorithm>
#define llint long long
#define eps 1e-9
#define inf 1e18
#define PI 3.14159265358979323846264

using namespace std;

struct vec2d{
	double x, y;
	vec2d(){}
	vec2d(double x, double y){
		this->x = x, this->y = y;
	}
	double add(double a, double b){
		if(fabs(a+b) < eps * (fabs(a) + fabs(b))) return 0.0;
		return a+b;
	}
	vec2d operator+(vec2d ope){
		return vec2d(add(x, ope.x), add(y, ope.y));
	}
	vec2d operator-(vec2d ope){
		return vec2d(add(x, -ope.x), add(y, -ope.y));
	}
	vec2d operator*(double t){
		return vec2d(x*t, y*t);
	}
	vec2d operator/(double t){
		return vec2d(x/t, y/t);
	}
	double dot(vec2d ope){
		return add(x*ope.x, y*ope.y);
	}
	double cross(vec2d ope){
		return add(x*ope.y, -y*ope.x);
	}
	double norm(){
		double d2 = dot(*this);
		if(d2 > 0) return sqrt(d2);
		return 0.0;
	}
};

double distPP(vec2d p, vec2d q){
	return (p-q).norm();
}

llint n;
vec2d p[305];

llint check(vec2d c)
{
	//cout << c.x << " " << c.y << endl;
	
	llint ret = 0;
	for(int i = 1; i <= n; i++){
		if(distPP(c, p[i]) <= 1+eps) ret++;
	}
	return ret;
}

int main(void)
{
	while(1){
		cin >> n;
		if(n == 0) break;
		for(int i = 1; i <= n; i++) cin >> p[i].x >> p[i].y;
		
		llint ans = 1;
		for(int i = 1; i <= n; i++){
			for(int j = 1; j <= n; j++){
				if(i >= j) continue;
				if(distPP(p[i], p[j]) > 2) continue;
				vec2d v = p[i]-p[j];
				double d = v.norm(), x = sqrt(1 - d*d/4);
				
				vec2d n = vec2d(v.y, -v.x);
				n = n / n.norm() * x;
				
				ans = max(ans, check(p[j]+v/2+n));
				ans = max(ans, check(p[j]+v/2-n));
			}
		}
		cout << ans << endl;
	}
	
	return 0;
}