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>
#define chmin(a, b) ((a)=min((a), (b)))
#define chmax(a, b) ((a)=max((a), (b)))
#define fs first
#define sc second
#define eb emplace_back
using namespace std;
 
typedef long long ll;
typedef pair<int, int> P;
typedef tuple<int, int, int> T;
 
const ll MOD=1e9+7;
const ll INF=1e18;
 
int dx[]={1, -1, 0, 0};
int dy[]={0, 0, 1, -1};

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

double dist(int a,int b){
    return sqrt((p[a][0]-p[b][0])*(p[a][0]-p[b][0])+(p[a][1]-p[b][1])*(p[a][1]-p[b][1]));
}

double get_o(int i,double st,double en,double th,double b){
    double diff = 1.0;
    while(diff > 0.00001){
        double half[2];
        half[0] = (en+st)/2.0;
        half[1] = th*half[0] + b;
        if( (p[i][0]-half[0])*(p[i][0]-half[0]) + (p[i][1]-half[1])*(p[i][1]-half[1]) >1){
            en = half[0];
        }else{
            st = half[0];
        }
        diff = en - st;
        if(diff < 0 ) diff = -1.0 * diff;
    }
  //  cout << "st:" << st << ", en:" << en << endl;
    return st;
}

int find_ans(double x,double y,int a,int b){
    int ans = 2;
    for(int i = 0;i < n;i++){
        if(i == a || i == b)continue;
        if( (p[i][0]-x)*(p[i][0]-x) + (p[i][1]-y)*(p[i][1]-y) <= 1.0) ans++;
    }
    return ans;
}

int main(){
    while(true){
        cin >> n;
        if(n == 0)break;
        for(int i = 0;i < n;i++){
            cin >> p[i][0] >> p[i][1];
        }
        int ans = 1;
        for(int i = 0;i < n;i++){
            for(int j = i+1;j < n;j++){
                double d = dist(i,j);
                if(d > 2.0)continue;
                if(true){
            //    if(p[i][1] != p[j][1]){
                    double th = -(p[j][0]-p[i][0])/(p[j][1]-p[i][1]);
                    double half[2];
                    half[0] = (p[i][0]+p[j][0])/2.0;
                    half[1] = (p[i][1]+p[j][1])/2.0;
                    double b = half[1] - th*half[0];
                    double o[2];
                    o[0] = get_o(i,half[0],half[0]+1,th,b);
                    o[1] = th*o[0] + b;
                //    cout << "x1:" << p[i][0] << ",y1" << p[i][1] << ",x2" << p[j][0] << ",y2" << p[j][1] << ",ox" << o[0] << ",oy" << o[1] << endl;
                    ans = max(ans,find_ans(o[0],o[1],i,j));
                    o[0] = get_o(i,half[0],half[0]-1,th,b);
                    o[1] = th*o[0] + b;
                //    cout << "x1:" << p[i][0] << ",y1" << p[i][1] << ",x2" << p[j][0] << ",y2" << p[j][1] << ",ox" << o[0] << ",oy" << o[1] << endl;
                    ans = max(ans,find_ans(o[0],o[1],i,j));
                }
            }
        }
        cout << ans << endl;
    }
}

C++: #include <bits/stdc++.h>
using namespace std;
typedef complex<double> P;

const double eps = 1e-5;
const int MAXN = 333;
int N;
P ps[MAXN];

int calc(P c) {
  int res = 0;
  for(int i = 0; i < N; ++i) {
    double d = abs(c - ps[i]);
    if(d < 1+eps) ++res;
  }
  return res;
}

int main() {
  for(; cin >> N && N; ) {
    for(int i = 0; i < N; ++i) {
      cin >> ps[i].real() >> ps[i].imag();
    }
    int res = 0;
    for(int i = 0; i < N; ++i) {
      for(int j = i+1; j < N; ++j) {
        double d = abs(ps[i] - ps[j]);
        P base = ps[i] - ps[j];
        P m = ps[j] + base / 2.0;
        P base2(-base.imag(), base.real());
        double a = sqrt(1.0 - d * d / 4.0);
        base2 = base2 / abs(base2) * a;
        res = max(res, max(calc(m + base2), calc(m - base2)));
      }
    }
    cout << res << endl;
  }
  return 0;
}

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

using namespace std;

typedef long double ld;

struct Point {
    ld x, y;
    Point(ld x=0.0, ld y=0.0) : x(x), y(y) {}
    string to_s() {
        ostringstream os;
        os << x << ',' << y;
        return os.str();
    }
    Point operator+(const Point& o) const {
        return Point(x+o.x, y+o.y);
    }
};

struct Circle {
    ld x, y;
    ld r;
    Circle(ld x=0.0, ld y=0.0, ld r=1.0) : x(x), y(y), r(r) {}
    const static ld eps = 1e-8;
    bool isIntersect(const Circle& o) const {
        ld dx = o.x - x;
        ld dy = o.y - y;
        return (r + o.r) * (r + o.r) > (dx*dx + dy*dy);
    }
    pair<Point, Point> intersectionPoint(const Circle& o) const {
        ld dx = x - o.x;
        ld dy = y - o.y;
        ld a = (dx*dx + dy*dy + r*r - o.r*o.r) / 2.0;
        ld d = sqrt((dx*dx+dy*dy)*r*r-a*a);
        Point p, q;
        p.x = (a*dx + dy * d) / (dx*dx+dy*dy) + o.x;
        q.x = (a*dx - dy * d) / (dx*dx+dy*dy) + o.x;
        p.y = (a*dy - dx * d) / (dx*dx+dy*dy) + o.y;
        q.y = (a*dy + dx * d) / (dx*dx+dy*dy) + o.y;
        //cout << p.x << ',' << p.y << endl;
        //cout << q.x << ',' << q.y << endl;
        return make_pair(p, q);
    }
    bool contains(ld px, ld py) const {
        return r*r + eps >= (px - x) * (px - x) + (py - y) * (py - y);
    }
    string to_s() {
        ostringstream os;
        os << x << ',' << y << ':' << r << endl;
        return os.str();
    }
};

int countContainingPoints(const Circle& c, const vector<Circle>& ps) {
    int cnt = 0;
    for (int i = 0; i < ps.size(); i++) 
        if (c.contains(ps[i].x, ps[i].y)) 
            cnt++;
    return cnt;
}

int N;
int main() {
    //cout << Circle(0, 0, 5).intersectionPoint(Circle(8, 0, 5)).first.to_s() << endl;
    vector<Circle> cs;
    while (cin >> N, N) {
        cs.clear(); cs.resize(N);
        for (int i = 0; i < N; i++) cin >> cs[i].x >> cs[i].y;

        vector<Point> ips; // Intersection PointS
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                if (i == j) continue;
                if (!cs[i].isIntersect(cs[j])) continue;
                pair<Point, Point> ip = cs[i].intersectionPoint(cs[j]);
                ips.push_back(ip.first);
                ips.push_back(ip.second);
            }
        }

        if (ips.empty()) {
            cout << 1 << endl;
            continue;
        }

        int maxPoint = 0;
        for (int i = 0; i < ips.size(); i++) {
            maxPoint = max(maxPoint, countContainingPoints(Circle(ips[i].x, ips[i].y, 1.0), cs));
        }
        cout << maxPoint << endl;
    }
    return 0;
}