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<cmath>
#include<algorithm>
#include<iostream>
#include<vector>
#include<climits>
#include<cfloat>
#define curr(P, i) P[(i) % P.size()]
#define next(P, i) P[(i+1) % P.size()]
#define prev(P, i) P[(i+P.size()-1) % P.size()]


using namespace std;

typedef double Real;

Real EPS = 1e-8;

int sgn(Real a, Real b=0){return a<b-EPS?-1:a>b+EPS?1:0;}
Real sqr(Real a){return sqrt(max(a,(Real)0));}

struct Point{  
  Real add(Real a, Real b){
    if(abs(a+b) < EPS*(abs(a)+abs(b)))return 0;
    return a+b;
  }

  Real x, y;
  Point(){}
  Point(Real x,Real y) : x(x) , y(y){}

  Point operator + (Point p){return Point(add(x,p.x), add(y,p.y));}
  Point operator - (Point p){return Point(add(x,-p.x), add(y,-p.y));}
  Point operator * (Real d){return Point(x*d,y*d);}
  Point operator / (Real d){return Point(x/d,y/d);}
  bool operator == (Point p){return !sgn(dist(p));}
  bool operator < (Point p){return (p.x!=x)?p.x<x:p.y<y;}
  Real norm(){return sqr(x*x+y*y);}
  Real dist(Point a){return (*this-a).norm();}
  Real dot(Point a){return x*a.x+y*a.y;}
  Real cross(Point a){return x*a.y-y*a.x;}
};

struct Circle{
  Point p;
  Real r;
  Circle(){}
  Circle(Point p, Real r) : p(p) , r(r){}
  
  vector<Point> intersectionPoints(Circle b){
    vector<Point> res;

    Point tmp=p-b.p;
    if(tmp.norm()>pow(r+b.r,2))return res;

    Real L = p.dist(b.p);
    Real C = atan2(b.p.y-p.y,b.p.x-p.x);
    Real a = acos((L*L+r*r-b.r*b.r)/(2*L*r));

    res.push_back(Point(p.x+cos(C+a), p.y+sin(C+a)));
    res.push_back(Point(p.x+cos(C-a), p.y+sin(C-a)));
    return res;
  }
};

int main(void){
  int n;
  Point a;
  vector<Point>P;
  vector<Circle>C;

  while(cin >> n,n){
    P.clear();
    C.clear();
    for(int i=0;i<n;i++){
      cin >> a.x >> a.y;
      P.push_back(a);
    }

    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
	if(sgn(P[i].dist(P[j]),2.0)>=0)continue;
	vector<Point>res=Circle(P[i],1.0).intersectionPoints(Circle(P[j],1.0));
	for(int k=0;k<res.size();k++)
	  C.push_back(Circle(res[k],1));
      }
    }

    int ans=1;
    for(int i=0;i<C.size();i++){
      int cnt=0;
      for(int j=0;j<P.size();j++)
	if(sgn(C[i].p.dist(P[j]),1.0)<=0)cnt++;

      ans=max(ans,cnt);
    }
    cout << ans << endl;
  }

  return 0;
}

C++: #include <bits/stdc++.h>
#define N 301
using namespace std;
typedef complex <double> P;
const double EPS = 1e-8;
P p[N];
int n;

int cnt(P c){
  int res=0;
  for(int i=0;i<n;i++)
    if(1-abs(c-p[i])>-EPS) res++;
  return res;
}

P get_C(P a,P b){
  double x=abs(b-a);
  if(x>2||a==b) return a;
  P c=P(x/2,sqrt(1-x*x/4))*(b-a)/abs(b-a)+a;

  return c;  
}


int main(){
  while(1){
  cin>>n;
  if(!n)break;
  double a,b;
  for(int i=0;i<n;i++)cin>>a>>b,p[i]=P(a,b);
  
  int ans=0;
  for(int i=0;i<n;i++)
    for(int j=0;j<n;j++){
      P c=get_C(p[i],p[j]);
      ans=max(ans,cnt(c));
    }
  cout <<ans<<endl;  
  }
  return 0;
}

C++: #include <iostream> 
#include <string> 
#include <vector> 
#include <cmath> 
#include <algorithm> 
#include <cstdlib> 
#include <ctime> 
#include <cstdio> 
#include <functional> 
#include <set> 
#include <sstream> 

using namespace std; 

double d2(double x1,double x2,double y1,double y2){
	return (x1-x2)*(x1-x2)+(y1-y2)*(y1-y2);
}

int main(){

	int n;
	const double eps=1e-8;
	while(true){

		cin>>n;

		if(n==0) return 0;

		int res=1;
		vector<double> x(n),y(n);

		for(int i=0;i<n;i++) cin>>x[i]>>y[i];

		for(int i=0;i<n;i++){
			for(int j=i+1;j<n;j++){
				
				int res2=0;
				const double X=(x[i]+x[j])/2.0;
				const double Y=(y[i]+y[j])/2.0;

				for(int k=0;k<n;k++){
//					printf("P1(%f, %f) P2(%f, %f) -> d2=%f\n",x[k],y[k],X,Y,sqrt(d2(x[k],X,y[k],Y)));
					if(d2(x[k],X,y[k],Y)<1.0+eps) res2++;
				}
//				cout<<endl;

				res=max(res,res2);

			}
		}

		cout<<res<<endl;

	}

	return 0;
}

C++: #include <bits/stdc++.h>

using namespace std;

typedef complex<double> P;

#define EPS 1e-12

int main() {
	int N;
	while (cin >> N, N) {
		vector<P> v;
		for (int i = 0; i < N; ++i) {
			double x, y; cin >> x >> y;
			v.push_back( P(x, y) );
		}

		vector<P> c;
		for (int i = 0; i < N; ++i) {
			for (int j = i+1; j < N; ++j) {
				P r = (v[j] - v[i])/2.0;
				double l = abs(r);
				if (l > 1 + EPS) {
					continue;
				}
				else if (l >= 1 - EPS) {
					c.push_back( v[i] + r );
				}
				else {
					P t = r / l * sqrt(1 - l*l);
					c.push_back(v[i] + r + t * P(0, 1));
					c.push_back(v[i] + r + t * P(0, -1));
				}
			}
		}

		int ans = 1;
		for (int i = 0; i < c.size(); ++i) {
			int cnt = 0;
			for (int j = 0; j < v.size(); ++j) {
				P p = c[i] - v[j];
				if ( p.real()*p.real()+p.imag()*p.imag() <= 1 + EPS ) {
					++cnt;
				}
			}
			ans = max(ans, cnt);
		}
		cout << ans << endl;
	}
}

C++: #include <bits/stdc++.h>
using namespace std;
using ld = long double;
using P = complex<ld>;
const ld eps = 1e-6;

int main()
{
	int N;
	while (cin >> N, N) {
		vector<ld> x(N), y(N);
		vector<P> p(N);
		for (int i = 0; i < N; i++) {
			cin >> x[i] >> y[i];
			p[i] = P(x[i], y[i]);
		}
		vector<P> ko;
		for (int i = 0; i < N; i++) {
			for (int j = i + 1; j < N; j++) {
				if (norm(p[i] - p[j]) < 4.0 + eps) {
					auto d = p[i] - p[j];
					ko.push_back(p[j] + (d * (ld)0.5) + d * P(0, 1) / abs(d) * sqrtl(1 - norm(d * (ld)0.5)));
					ko.push_back(p[j] + (d * (ld)0.5) - d * P(0, 1) / abs(d) * sqrtl(1 - norm(d * (ld)0.5)));
				}
			}
		}
		int res = 0;
		for (auto c : ko) {
			int cnt = 0;
			for (int i = 0; i < N; i++) {
				if (norm(c - p[i]) < 1.0001) {
					cnt++;
				}
			}
			res = max(res, cnt);
		}
		cout << res << endl;
	}
	return 0;
}