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 EPS (1e-10)
#define equals(a, b) (fabs((a) - (b)) < EPS)
const double pi = acos(-1);

struct Point{
    double x, y;
    
    Point(double x = 0, double y = 0): x(x), y(y) {}

    Point operator + (Point p) { return Point(x + p.x, y + p.y); }
    Point operator - (Point p) { return Point(x - p.x, y - p.y); }
    Point operator * (double a) { return Point(a * x, a * y); }
    Point operator / (double a) { return Point(x / a, y / a); }

    double norm() {return x * x + y * y; }
    double abs() {return sqrt(norm()); }

    bool operator < (const Point &p) const{
        return (x != p.x ? x < p.x : y < p.y);
    }

    bool operator == (const Point &p) const{
        return (fabs(x - p.x) < EPS && fabs(y - p.y) < EPS);
    }
};

typedef Point Vector;

double norm(Vector a){
    return a.x * a.x + a.y * a.y;
}

double abs(Vector a){
    return sqrt(norm(a));
}

double arg(Vector p){
    return atan2(p.y, p.x);
}

Vector polar(double a, double r){
    return Point(cos(r) * a, sin(r) * a);
}

int main(){
    while(1){
        int n; cin>>n;
        if(n == 0) return 0;
        vector<Point> p;
        for(int i=0; i<n; i++){
            double x, y; cin>>x>>y;
            p.emplace_back(x, y);
        }

        int ans = 1;
        for(int i=0; i<n; i++){
            for(int j=i+1; j<n; j++){
                Point p1 = p[i];
                Point p2 = p[j];

                double dist = abs(p1 - p2);

                if(2.0 < dist) continue;

                Point mp = {(p1.x + p2.x) / 2, (p1.y + p2.y) / 2};

                double rad = arg(p2 - p1);
                
                double rad1 = rad + pi / 2;
                double rad2 = rad - pi / 2;

                double sideLength = sqrt(1.0 * 1.0 - abs(mp - p1) * abs(mp - p1));
                
                Point c1 = mp + polar(sideLength, rad1);
                int sum = 0;
                for(int k=0; k<n; k++){
                    dist = abs(c1 - p[k]);
                    if(1.0 + EPS < dist) continue;
                    sum++;
                }

                ans = max(ans, sum);

                Point c2 = mp + polar(sideLength, rad2);
                sum = 0;
                for(int k=0; k<n; k++){
                    dist = abs(c2 - p[k]);
                    if(1.0 + EPS < dist) continue;
                    sum++;
                }

                ans = max(ans, sum);
            }
        }

        cout << ans << endl;
    }
}

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

using namespace std;

typedef complex<double> Point;

const double EPS = 1e-8;

int n;
Point data[300];

int calc(Point a, Point b){
  int res = 0;
  double d = abs(b - a) / 2.0;
  Point ab = (b - a) / (2.0 * d);
  if(d < 1){
    Point o = ab * d;
    d = sqrt(1.0 - d * d);
    o += Point(d * ab.imag() + a.real(), -d * ab.real() + a.imag());
    for(int i=0;i<n;i++) if(abs(o - data[i]) < 1.0 + EPS) res++;
  }
  return res;
}

main(){
  while(cin >> n && n){
    for(int i=0;i<n;i++) cin >> data[i].imag() >> data[i].real();
    int ans = 1;
    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
        ans = max(ans, calc(data[i], data[j]));
        ans = max(ans, calc(data[j], data[i]));
      }
    }
    cout << ans << endl;
  }
}

C++: #include<iostream>
#include<cmath>
using namespace std;

class P
{
public:
	double x,y;
	P(double x=0.0,double y=0.0):x(x),y(y){}
	P operator+(const P p)const{return P(x+p.x,y+p.y);}
	P operator-(const P p)const{return P(x-p.x,y-p.y);}
};

double dis(P p1,P p2)
{
	return hypot(p1.x-p2.x,p1.y-p2.y);
}

P getc(P p1,P p2)
{
	double theta=acos(dis(p1,p2)/2.0);
	double alpha=atan2(p2.y-p1.y,p2.x-p1.x);
	return P(p1.x+cos(alpha-theta),p1.y+sin(alpha-theta));
}

bool in(P c,P p)
{
	return dis(c,p)<1.0+1e-9;
}

int main()
{
	int n;
	while(cin>>n,n)
	{
		P p[300];
		for(int i=0;i<n;i++)
		{
			cin>>p[i].x>>p[i].y;
		}
		int m=0;
		for(int i=0;i<n;i++)
		{
			for(int j=i+1;j<n;j++)
			{
				P a=p[i],b=p[j];
				P c=getc(a,b);
				int cnt=0;
				for(int k=0;k<n;k++)
				{
					if(k==i||k==j)cnt++;
					else if(in(c,p[k]))cnt++;
				}
				m>cnt?:m=cnt;
			}
		}
		cout<<m<<endl;
	}
}

C++: #include <iostream>
#include <iomanip>
#include <complex>
#include <vector>
#include <algorithm>
#include <cmath>
#define EPS 1e-8
using namespace std;
typedef complex<double> P;
#define X real()
#define Y imag()

bool equidistantpoint(const P& a, const P& b, double r, vector<P>& ep){
    if(abs(b-a)>r*2+EPS) return false;
    P mid = (a+b)/P(2,0);
    P len = P(sqrt(r*r-norm(b-mid)),0);
    P rot = P((b-a).Y,(a-b).X) / P(abs(b-a),0);
    ep.resize(2);
    ep[0] = len*rot + mid;
    ep[1] = -len*rot + mid;
    return true;
}

int main(){
    while(1){
        int n;
        cin >> n;
        if(n==0) break;

        vector<P> pset(n);
        for(int i=0; i<n; i++){
            double x,y;
            cin >> x >> y;
            pset[i] = P(x,y);
        }

        int mnop=0;
        for(int i=0; i<n-1; i++){
            for(int j=i+1; j<n; j++){
                vector<P> endpoints;
                bool ret = equidistantpoint(pset[i], pset[j], 1, endpoints);
                if(!ret)continue;
                for(int lr=0; lr<2; lr++){
                    int count=0;
                    for(int k=0; k<n; k++){
                        if(abs(pset[k]-endpoints[lr]) < 1+EPS){
                            count++;
                        }
                    }
                    mnop = max(mnop, count);
                }
            }
        }

        if(mnop==0) mnop=1;
        cout << mnop << endl;
    }
            
    return 0;
}