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 <cmath>

using namespace std;
const double EPS = 1e-8;
typedef pair<double, double> pt;
#define x first
#define y second

pt operator-(pt a, pt b) {
  return pt(a.x - b.x, a.y - b.y);
}

bool zero(double x) {
  return fabs(x) <= EPS;
}

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

double dist(pt p, pt q) {
  return sqrt(sq(p.x - q.x) + sq(p.y - q.y));
}

pt midpoint(pt p1, pt p2){
    double midx = (p1.x+p2.x)/2;
    double midy = (p1.y+p2.y)/2;
    return make_pair(midx, midy);
}

int get_max_points(pt centre, vector<pt> points){

    int count = 0;
    for (int i = 0; i < points.size(); i++){
        pt current = points[i];
        double num = sq(current.x - centre.x) + sq(current.y - centre.y);
        // cout << num << endl;
        // cout << fabs(num - 1) << endl;
        if (sqrt(num) < 1+EPS) {
            count++;
        }
    }
    return count;
}

int find_max_points(vector<pt> points){
    // for every pair of two points, create a cricle and check num points in circle

    int max_points = 1;
    for (int i = 0; i < points.size(); ++i){
        pt p1 = points[i];
        for (int j = i+1; j < points.size(); ++j){
            pt p2 = points[j];

            pt mid = midpoint(p1, p2);
            double q = dist(p1, p2);

            double x1 = mid.x + sqrt(1-sq(q/2))*(p1.y-p2.y)/q;
            double y1 = mid.y + sqrt(1-sq(q/2))*(p2.x-p1.x)/q;
            pt c1(x1, y1);
            int m = get_max_points(c1, points);
            max_points = max(m, max_points);
            // cout << mid.x << " " << mid.y << endl;
            // int m = get_max_points(mid, points);
            double x2 = mid.x - sqrt(1-sq(q/2))*(p1.y-p2.y)/q;
            double y2 = mid.y - sqrt(1-sq(q/2))*(p2.x-p1.x)/q;
            pt c2(x1, x2);
            m = get_max_points(c2, points);

            max_points = max(max_points, m);
        }
    }
    // cout << max_points 
    return max_points;
}

int main(){

    int n;

    while (1){
        cin >> n;
        if (n == 0) break;
        vector<pt> points;
        for (int i = 0; i < n; ++i){
            double x, y;
            cin >> x >> y;
            points.push_back(make_pair(x, y));
        }
        cout << find_max_points(points) << endl;
    }


    return 0;
}

C++: #include "bits/stdc++.h"
using namespace std;

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

const double EPS = 1e-8;

int n;
double x[303], y[303];

int main() {
        while (cin >> n && n) {
                rep(i, n) cin >> x[i] >> y[i];
                int ans = 1;
                rep(i, n) rep(j, n) {
                        if (i == j) continue;
                        double dx = x[i] - x[j], dy = y[i] - y[j];
                        double d = sqrt(dx * dx + dy * dy);
                        if (d > 2.0 + EPS) continue;
                        double r = sqrt(1.0 - (d / 2.0) * (d / 2.0));
                        double vx = -dy / d * r;
                        double vy = dx / d * r;
                        double xx = (x[i] + x[j]) / 2.0 + vx;
                        double yy = (y[i] + y[j]) / 2.0 + vy;
                        int cnt = 0;
                        rep(k, n) if ((x[k] - xx) * (x[k] - xx) + (y[k] - yy) * (y[k] - yy) < 1.0 + EPS) cnt ++;
                        ans = max(ans, cnt);
                }
                cout << ans << endl;
        }
        return 0;
}

C++: #include<bits/stdc++.h>
#define N 300
using namespace std;
typedef complex<double> P;
const double EPS = 1e-7;

struct Circle{
  P c;
  double r;
  Circle(P c=P(0,0),double r=0):c(c),r(r){}
};

vector<P> getIntersectCC(Circle c1, Circle c2){
  vector<P> res;
  double r1=c1.r,r2=c2.r,d;
  P p1=c1.c,p2=c2.c;
  d=abs(p1-p2);
  if(d<EPS&&abs(r1-r2)<EPS) return res;
  if(r1+r2<d-EPS||d+EPS<abs(r1-r2)) return res; 
  double a=(r1*r1-r2*r2+d*d)/(2*d);
  double h=sqrt(max(r1*r1-a*a,0.0));
  P tmp1=p1+a/d*(p2-p1);
  P tmp2=h/d*(p2-p1);
  if(abs(tmp2)<EPS)res.push_back(tmp1);
  else{
    res.push_back(P(real(tmp1)-imag(tmp2),imag(tmp1)+real(tmp2)));
    res.push_back(P(real(tmp1)+imag(tmp2),imag(tmp1)-real(tmp2)));
  }
  return res;
}

double getdistancepp(P a,P b){
  return sqrt((a.real()-b.real())*(a.real()-b.real())+(a.imag()-b.imag())*(a.imag()-b.imag()));
}

int main(){
  int n,ans=0;
  double a,b;
  P p[N];
  while(1){
    cin>>n;
    if(!n)break;
    for(int i=0;i<n;i++){
      cin>>a>>b;
      p[i]=P(a,b);
    }
    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
	Circle c1,c2;
	c1.c=p[i];
	c2.c=p[j];
	c1.r=c2.r=1;
	vector<P> r=getIntersectCC(c1,c2);
	if(!r.size())continue;
	for(int k=0;k<r.size();k++){
	  int cnt=0;
	  for(int l=0;l<n;l++)
	    if(getdistancepp(r[k],p[l])<=1+EPS)cnt++;	    
	  ans=max(ans,cnt);
	}
      }
    }
    cout<<ans<<endl;
  }
  return 0;
}

C++: #include <iostream>
#include <cmath>

using namespace std;

#define EPS (1e-6)


void calc(double x1,double y1,double x2,double y2,double *c_x1,double *c_y1,double *c_x2,double *c_y2){
  double midx,midy,hox,hoy,dx,dy,dis;
  midx = (x1+x2) / 2.0 , midy = (y1+y2) / 2.0;
  dx = x1-x2 , dy = y1-y2;
  dis = sqrt(1 - pow(x1-midx,2.0) - pow(y1-midy,2.0));
  hox = sqrt(pow(dy,2.0) / (pow(dy,2.0) + pow(dx,2.0)));
  hoy = hox * (-1.0) * (dx / dy);
  *c_x1 = midx + hox*dis,*c_y1 = midy + hoy*dis;
  *c_x2 = midx - hox*dis,*c_y2 = midy - hoy*dis;
}

int main(){
  int N;

  while(cin >> N){
    double x[300],y[300],c_x1,c_y1,c_x2,c_y2;
    int count,m;
    m = 1;
    if(!(N)) break;
    for(int i=0;i<N;i++) cin >> x[i] >> y[i];
    for(int i=0;i<N;i++){
      for(int j=0;j<N;j++){
        if(i == j) continue;
        if(sqrt((x[i] - x[j])*(x[i] - x[j]) + (y[i] - y[j])*(y[i] - y[j])) > 2.0 + EPS  ) continue;
        count = 2;
        calc(x[i],y[i],x[j],y[j],&c_x1,&c_y1,&c_x2,&c_y2);
        for(int k=0;k<N;k++){
          if(k == i || k == j) continue;
          if((x[k] - c_x1)*(x[k] - c_x1) + (y[k] - c_y1)*(y[k] - c_y1) < 1.0 + EPS  ) count++;
        }
        if(m < count) m = count;
        count = 2;
        for(int k=0;k<N;k++){
          if(k == i || k == j) continue;
          if((x[k] - c_x2)*(x[k] - c_x2) + (y[k] - c_y2)*(y[k] - c_y2) < 1.0 + EPS) count++;
        }
        if(m < count) m = count;
      }
    }

    printf("%d\n",m);
  }

  return 0;
}