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)
	}
}

Python: from math import sqrt

def center(a,b):
    m = (a+b)/2
    d = abs(a-m)
    vec= complex((b-a).imag,(a-b).real)
    length = sqrt(1-d**2)
    c1,c2 = m+vec*(length/(2*d)),m-vec*(length/(2*d))
    return c1,c2



def check(c):
    cnt = 0
    for loc in location:
        if abs(c-loc) <= 1.0+EPS:
            cnt += 1
    return cnt 

while True:
    n = int(input())
    if n==0:
        break
    location = []
    ans = 1
    EPS = 1e-7
    for _ in range(n):
        a,b = map(float ,input().split())
        location.append(complex(a,b))

    for i in range(n):
        for j in range(i+1,n):
            if abs(location[i]-location[j]) >=1.0*2:
                continue
            c1,c2 = center(location[i],location[j])
            ans = max(ans,check(c1),check(c2))
            
    print(ans)

Python: def cent(p1,p2):
    if abs(p1-p2) >= 2:
        return None,None
    else:
        q = (p1+p2)/2
        r = (p1-p2)/2
        n = abs(r)
        t = r * 1j * ((1-n**2)**0.5/n)
        return(q + t, q - t)

while(True):
    ans = 1
    n = int(input())
    if n == 0:
        quit()
    ps = []
    for _ in range(n):
        x,y = map(float,input().split())
        ps.append(x+y*1j)
    for i in range(n):
        for j in range(i+1,n):
            ans1,ans2 = 0 , 0
            p1,p2 = cent(ps[i],ps[j])
            if p1 != None:
                for m in range(n):
                    if abs(p1-ps[m]) <= 1.001:
                        ans1 += 1
                    if abs(p2-ps[m]) <= 1.001:
                        ans2 += 1
                ans = max(ans,ans1,ans2)
    print(ans)

Python: #!/usr/bin/env python

from collections import deque
import itertools as it
import sys
import math

sys.setrecursionlimit(1000000)

def func():
    lst = []
    while True:
        S = input()
        lst.append(S)
        if S == '0':
            break
    lst = list(reversed(lst))
    while True:
        N = int(lst.pop())
        if N == 0:
            break
        R = range(N)
        pos = []
        for loop in R:
            S = list(map(float, lst.pop().split()))
            pos.append(S)
        ans = 0
        for i in R:
            for j in R:
                v = [pos[j][0] - pos[i][0], pos[j][1] - pos[i][1]]
                norm = math.sqrt(v[0] * v[0] + v[1] * v[1])
                if norm > 2 or norm == 0:
                    continue
                vn = [v[0] / norm, v[1] / norm]
                norm /= 2
                norm2 = math.sqrt(1 - norm * norm)
                v90 = [-vn[1], vn[0]]
                cent = [pos[i][0] + v[0] / 2 + v90[0] * norm2, 
                pos[i][1] + v[1] / 2 + v90[1] * norm2]
                S = 2
                for k in R:
                    if i == k or j == k:
                        continue
                    x = cent[0] - pos[k][0]
                    y = cent[1] - pos[k][1]
                    if x * x + y * y <= 1:
                        S += 1
                ans = max(S, ans)
        print(max(1, ans))

if __name__ == '__main__':
    func()

Python: #!/usr/bin/env python

from collections import deque
import itertools as it
import sys
import math

sys.setrecursionlimit(1000000)

def func():
    while True:
        N = int(input())
        if N == 0:
            break
        R = range(N)
        pos = []
        for loop in R:
            S = list(map(float, input().split()))
            pos.append(S)
        ans = 0
        for i in R:
            for j in R:
                v = [pos[j][0] - pos[i][0], pos[j][1] - pos[i][1]]
                norm = math.sqrt(v[0] * v[0] + v[1] * v[1])
                if norm > 2 or norm == 0:
                    continue
                vn = [v[0] / norm, v[1] / norm]
                norm /= 2
                norm2 = math.sqrt(1 - norm * norm)
                v90 = [-vn[1], vn[0]]
                cent = [pos[i][0] + v[0] / 2 + v90[0] * norm2, 
                pos[i][1] + v[1] / 2 + v90[1] * norm2]
                S = 2
                for k in R:
                    if i == k or j == k:
                        continue
                    x = cent[0] - pos[k][0]
                    y = cent[1] - pos[k][1]
                    if x * x + y * y <= 1:
                        S += 1
                ans = max(S, ans)
        print(max(1, ans))

if __name__ == '__main__':
    func()

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

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.0 + eps) {
					cnt++;
				}
			}
			res = max(res, cnt);
		}
		cout << res << endl;
	}
	return 0;
}