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 <cmath>
#define EPS 1.0e-4
using namespace std;

typedef struct {
  double x,y;
} P;

double dis(P p ){
  return sqrt( p.x*p.x + p.y*p.y );
}
double dis(P p1,P p2 ){
  P p;
  p.x = p1.x-p2.x;
  p.y = p1.y-p2.y;
  return dis( p );
}


int main(){
  int n, mx;
  // mv : @üxNg
  P p[301], m, p1,p2, mv,dmv, c[2];
  double x;

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

    mx=0;

    for( int i=0;i<n;i++ ){
      for( int j=i+1;j<n;j++ ){
	if( dis(p[i],p[j]) < 2.0 ){
	  //	if( dis(p[i],p[j]) < 2.0 + EPS ){
	  p1=p[i]; p2=p[j];
	  m.x = (p1.x+p2.x)/2.0;
	  m.y = (p1.y+p2.y)/2.0;

	  x = sqrt( 1-dis(p1,m) );
	  // S_ cðßé
	  if( fabs( p1.y-p2.y ) < EPS ){
	    mv.x=0.0;  mv.y=1.0;
	  }else if( fabs( p1.x-p2.x ) < EPS ){
	    mv.x=1.0;  mv.y=0.0;
	  }else{
	    double a1 = ( p2.y-p1.y )/( p2.x-p1.x );
	    mv.x=1.0; mv.y=-1.0/a1;
	  }
	  double scale = dis( mv );
	  dmv.x = x*mv.x / scale;
	  dmv.y = x*mv.y / scale;

	  c[0].x=m.x+dmv.x; c[0].y=m.y+dmv.y;
	  c[1].x=m.x-dmv.x; c[2].y=m.y-dmv.y;
	  /*
	  cout << "m  " << m.x    << "," << m.y    << endl;
	  cout << "mv " << mv.x   << "," << mv.y   << endl;
	  cout << "dm " << dmv.x  << "," << dmv.y  << endl;
	  cout << "c0 " << c[0].x << "," << c[0].y << endl;
	  cout << "c1 " << c[1].x << "," << c[1].y << endl;
	  */

	  for( int cn=0;cn<2;cn++ ){
	    int sum=0;
	    for( int k=0;k<n;k++ ){
	      if( dis( c[cn],p[k] ) < 1.0 )
		//	      if( dis( c[cn],p[k] ) < 1.0 + EPS )
		sum++;
	    }
	    if( mx<sum ) mx=sum;
	  }
	}
      }
    }
    cout << mx << endl;
  }

  return 0;
}

C++: #include<iostream>
#include<math.h>

using namespace std;

#define REP(i, n) for(int i=0; i<n; i++)

double p[300][2];
int n;

double dis(double x1, double y1, double x2, double y2){
    return sqrt(pow(x1-x2, 2) + pow(y1-y2, 2));
}

int circle(double x, double y, int i, int j){
    int count = 2;
    REP(k, n){
        if(k!=i && k!=j){
            if(dis(x, y, p[k][0], p[k][1]) < 1){
                count += 1;
            }
        }
    }
    return count;
}

int main(){
    while(1){
        cin >> n;

        if(n==0) break;

        REP(i, n){
            cin >> p[i][0];
            cin >> p[i][1];
        }

        int output = 1;
        REP(i, n){
            for(int j=i+1; j<n; j++){
                double cx = (p[i][0]+p[j][0]) / 2;
                double cy = (p[i][1]+p[j][1]) / 2;
                double d = dis(p[i][0], p[i][1], p[j][0], p[j][1]);

                if(d>2){
                    continue;
                } else if (d==2){
                    output = max(output, circle(cx, cy , i, j));
                } else {
                    double a1 = acos(d/2);
                    double a2 = atan2((p[j][1]-p[i][1]), p[j][0]-p[i][0]);

                    double x1 = p[i][0] + cos(a2+a1);
                    double y1 = p[i][1] + sin(a2+a1);

                    double x2 = p[i][0] + cos(a2-a1);
                    double y2 = p[i][1] + sin(a2-a1);

                    output = max(output, circle(x1, y1, i, j));
                    output = max(output, circle(x2, y2, i, j));
                }
            }
        }

        cout << output << endl;
    }
}

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

typedef complex<double> p;
const double eps = 1e-9;

int n;
p ps[300];

int main() {
	while (cin >> n , n != 0) {
		for (int i = 0; i < n; i++) {
			double x, y;
			cin >> x >> y;
			ps[i] = p(x, y);
		}

		int ans = 1;
		for (int i = 0; i < n; i++) {
			for (int j = i + 1; j < n; j++) {
				if (!(norm(ps[i] - ps[j]) < 4.0)) continue;
				p v = ps[i] - ps[j], ve = v / abs(v);
				double l = sqrt(1 - norm(v) * 0.25);
				p vn1 = ve * p(0, 1) * l, vn2 = ve * p(0, -1) * l;
				p c1 = ps[j] + v * 0.5 + vn1, c2 = ps[j] + v * 0.5 + vn2;
				int count1 = 2, count2 = 2;
				for (int k = 0; k < n; k++) {
					if (k == i || k == j) continue;
					if (norm(c1 - ps[k]) < 1.0 + eps) count1++;
					if (norm(c2 - ps[k]) < 1.0 + eps) count2++;
				}
				ans = max(ans, count1);
				ans = max(ans, count2);
			}
		}

		cout << ans << endl;
	}
}

C++: #include<iostream>
#include<vector>
#include<complex>
#include<algorithm>
#include<cmath>
using namespace std;
#define REP(i, j) for(int i = 0; i < j; i++)
#define FOR(i, j, k) for(int i = j; i < k; i++)
#define P complex<double>
#define X real()
#define Y imag()
const double EPS = 1e-9;

int n;
vector<P> vv;

int count_points(P a, P b, double mx, double my){
    P v = (a - b);
    double d_v = abs(v);
    if(d_v >= 2.0 + EPS) return 0;
    P m = v / 2.0 + b;
    double half_d_v = d_v / 2;
    double d = sqrt(1.0 - half_d_v * half_d_v);
    P u = (v * P(mx, my)) / abs(v) * d;
    //debug
    //cout <<"s " <<a <<", " <<b <<endl;
    //cout <<v <<", " <<d_v <<", " <<m <<", " <<d <<", " <<u <<endl;
    int cnt = 0;
    REP(i, n){
        P c = m + u;
        //cout <<"c = " <<c <<endl;
        //cout <<"! " <<i <<", " <<vv[i] <<": " <<abs(c - vv[i]) <<endl;
        if(abs(c - vv[i]) <= 1.0 + EPS) cnt++;
    }
    return cnt;
}

int main(){
    while(cin >>n && n){
        vv = vector<P>(n);
        REP(i, n) cin >>vv[i].X >>vv[i].Y;

        int ans = 1;
        REP(a, n){
            FOR(b, a + 1, n){
                ans = max(ans, count_points(vv[a], vv[b], 0, 1));
                ans = max(ans, count_points(vv[a], vv[b], 0, -1));
            }
        }
        cout <<ans <<endl;
    }
    return 0;
}

C++: #include <iostream>
#include <string>
#include <vector>

using namespace std;

int main() {
    int n;
    double x, y;
    while(cin >> n && n) {
        vector<double> xs, ys;
        for (int i=0; i<n; ++i) {
            cin >> x >> y;
            xs.push_back(x);
            ys.push_back(y);
        }
        vector<double> cxs, cys;
        for (int i=0; i<n; ++i) {
            for (int j=i+1; j<n; ++j) {
                cxs.push_back((xs[i] + xs[j]) / 2.0);
                cys.push_back((ys[i] + ys[j]) / 2.0);
                for (int k=j+1; k<n; ++k) {
                    cxs.push_back((xs[i] + xs[j] + xs[k]) / 3.0);
                    cys.push_back((ys[i] + ys[j] + ys[k]) / 3.0);
                }
            }
        }
        int ret = 1;
        for (int i=0; i<(int)cxs.size(); ++i) {
            double cx = cxs[i], cy = cys[i];
            int cnt = 0;
            for (int j=0; j<n; ++j) {
                if ((xs[j] - cx) * (xs[j] - cx) + (ys[j] - cy) * (ys[j] - cy) < 1.0 + 1e-8) {
                    cnt ++;
                }
            }
            ret = max(ret, cnt);
        }
        cout << ret << endl;
    }
    return 0;
}