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: #!/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
        pos = []
        for loop in range(N):
            S = list(map(float, lst.pop().split()))
            pos.append(S)
        ans = 0
        for i in range(N):
            for j in range(N):
                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 i == j:
                    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 range(N):
                    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: import math


def find_point(x1, y1, x2, y2):
    # 二頂点(x1,y1),(x2,y2)を円周上にもつ時の半径1の円の中心
    if (x1-x2)**2+(y1-y2)**2 > 4:
        return False
    mx, my = (x1+x2)/2, (y1+y2)/2
    L = math.sqrt((y2-y1)**2+(x2-x1)**2)
    k = math.sqrt(1-((y2-y1)**2+(x2-x1)**2)/4)
    # (mx,my)+-(k/L)(-y2+y1,x2-x1)
    X1, Y1 = mx+(k/L)*(-y2+y1), my+(k/L)*(x2-x1)
    X2, Y2 = mx-(k/L)*(-y2+y1), my-(k/L)*(x2-x1)
    return ((X1, Y1), (X2, Y2))


def solve():
    N = int(input())
    if N == 0:
        return
    que = [tuple(map(float, input().split())) for _ in range(N)]
    ans = 1
    for i in range(N):
        for j in range(i+1, N):
            A = find_point(que[i][0], que[i][1], que[j][0], que[j][1])
            if A != False:
                a, b = A
                suba = 0
                subb = 0
                for k in range(N):
                    if math.hypot(que[k][0]-a[0], que[k][1]-a[1]) < 1+10**(-7):
                        suba += 1
                    if math.hypot(que[k][0]-b[0], que[k][1]-b[1]) < 1+10**(-7):
                        subb += 1
                ans = max(ans, suba, subb)
    print(ans)
    return solve()


def main():
    solve()


if __name__ == "__main__":
    main()

Python: import cmath
EPS = 1e-5

while True:
	n = int(input())
	if n == 0:
		break
	d = []
	for _ in range(n):
		x, y = list(map(float, input().split()))
		d.append(complex(x, y))

	res = 1
	for i in range(n):
		o1 = d[i]
		cnt0 = 1
		for j in range(n):
			if j == i:
				continue
			o2 = d[j]
			if abs(o1 - o2) + EPS <= 1:
				cnt0 += 1

		res = max(res, cnt0)
		for j in range(i+1, n):
			o2 = d[j]
			if abs(o1 - o2) - EPS >= 2:
				continue
			else:
				co = (abs(o1 - o2)**2) / (2 * abs(o1 - o2))
				si = cmath.sqrt(1 - co**2)
				p1 = o1 + (o2 - o1) * (co + 1j * si) / abs(o1 - o2)
				p2 = o1 + (o2 - o1) * (co - 1j * si) / abs(o1 - o2)

			cnt1, cnt2 = 2, 2
			for l in range(n):
				if l == i or l == j:
					continue
				else:
					o3 = d[l]
					if abs(p1 - o3) + EPS <= 1:
						cnt1 += 1
					if abs(p2 - o3) + EPS <= 1:
						cnt2 += 1
			res = max(res, cnt1, cnt2)
	print(res)

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)
        if n == 0:
            return None,None
        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()
    elif n == 1:
        print(1)
    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
            p_1,p_2 = cent(ps[i],ps[j])
            if p_1 != None:
                for m in range(n):
                    if abs(p_1-ps[m]) <= 1.001:
                        ans1 += 1
                    if abs(p_2-ps[m]) <= 1.001:
                        ans2 += 1
                ans = max(ans,ans1,ans2)
    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)
        if n == 0:
            return None,None
        t = r * 1j * ((1-n**2)**0.5/n)
        return(q + t, q - t)

def abs(x):
    return(((x.real) ** 2 + (x.imag) ** 2) ** 0.5)

while(True):
    ans = 1
    n = int(input())
    if n == 0:
        quit()
    elif n == 1:
        print(1)
    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
            p_1,p_2 = cent(ps[i],ps[j])
            if p_1 != None:
                for m in range(n):
                    if abs(p_1-ps[m]) <= 1.001:
                        ans1 += 1
                    if abs(p_2-ps[m]) <= 1.001:
                        ans2 += 1
                ans = max(ans,ans1,ans2)
    print(ans)

Python: #複素数でやってみる
from itertools import combinations
import math

while True:
    N=int(input())
    if N==0:
        break
    mp=[]
    for _ in range(N):
        a,b=map(float,input().split())
        mp.append(complex(a,b))
    center_lst=[]
    for p,q in combinations(mp,2):
        mid=(p+q)/2
        if abs(p-q)>2:
            continue
        d_vec=complex(-(p-q).imag,(p-q).real)/abs(p-q)
        t=1-abs(p-mid)**2
        center_vec1=mid+d_vec*math.sqrt(t)
        center_vec2=mid-d_vec*math.sqrt(t)
        center_lst.append(center_vec1)
        center_lst.append(center_vec2)
    
    ans=1
    for center in center_lst:
        tmp=0
        for x in mp:
            x_center=x-center
            if abs(x_center)<=1+1e-7:
                tmp+=1
        if tmp>ans:
            ans=tmp
    print(ans)