Problem:
 Score : 600 points    Problem Statement Given are N points (x_i, y_i) in a two-dimensional plane.  Find the minimum radius of a circle such that all the points are inside or on it.      Constraints  2 \leq N \leq 50  0 \leq x_i \leq 1000  0 \leq y_i \leq 1000  The given N points are all different.  The values in input are all integers.         Input Input is given from Standard Input in the following format:  N  x_1  y_1  :  x_N  y_N       Output Print the minimum radius of a circle such that all the N points are inside or on it.  Your output will be considered correct if the absolute or relative error from our answer is at most 10^{-6} .        Sample Input 1 2
0 0
1 0      Sample Output 1 0.500000000000000000  Both points are contained in the circle centered at (0.5,0) with a radius of 0.5 .       Sample Input 2 3
0 0
0 1
1 0      Sample Output 2 0.707106781186497524       Sample Input 3 10
10 9
5 9
2 0
0 0
2 7
3 3
2 5
10 0
3 7
1 9      Sample Output 3 6.726812023536805158  If the absolute or relative error from our answer is at most 10^{-6} , the output will be considered correct.

Solution (Python):
from math import sqrt

N = int(input())
XY = [tuple(map(int, input().split())) for _ in range(N)]

def calc(r): 
    lst = [] 
    for i in range(N - 1):
        x1 = XY[i][0]
        y1 = XY[i][1]
        for j in range(i + 1, N):
            x2 = XY[j][0]
            y2 = XY[j][1]
            diff = sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)

            if diff > 2 * r:
                return False
            
            h = sqrt(r ** 2 - (diff / 2) ** 2)
            
            if x1 == x2:
                ny = (y1 + y2) / 2
                lst.append([x1 - h, ny])
                lst.append([x1 + h, ny]) 
            elif y1 == y2:
                nx = (x1 + x2) / 2
                lst.append([nx, y1 - h])
                lst.append([nx, y1 + h])
            else:
                a = (y2 - y1) / (x2 - x1) 

                b = -1 / a

                size = sqrt(1 + b ** 2)

                nx = h / size
                ny = nx * b

                
                xc = (x1 + x2) / 2
                yc = (y1 + y2) / 2

                lst.append([xc + nx, yc + ny])
                lst.append([xc - nx, yc - ny])
    
    
    nr = r + eps
    for x, y in lst: 
        flag = True
        for X, Y in XY:
            tmp = (x - X) ** 2 + (y - Y) ** 2
            if tmp > nr ** 2:
                flag = False
                break
        if flag:
            return True
    return False
    
l = 0
r = 1000

eps = 10 ** (-9)

for i in range(100):
    mid = (l + r) / 2
    if calc(mid):
        r = mid
    else:
        l = mid

print (r)

Solution (Python):
import sys; input = sys.stdin.buffer.readline
from collections import defaultdict
import math
con = 10 ** 9 + 7; INF = float("inf")

def getlist():
	return list(map(int, input().split()))

def ccw(ax, ay, bx, by, cx, cy):
	x1 = bx - ax; y1 = by - ay
	x2 = cx - ax; y2 = cy - ay
	cross = x1 * y2 - x2 * y1 
	if cross > 0:
		return 1
	elif cross < 0:
		return - 1
	else:
		return 0

class triangle(object):
	def __init__(self, Ax, Ay, Bx, By, Cx, Cy):
		self.a = ((Cx - Bx) ** 2 + (Cy - By) ** 2) ** 0.5
		self.b = ((Ax - Cx) ** 2 + (Ay - Cy) ** 2) ** 0.5
		self.c = ((Bx - Ax) ** 2 + (By - Ay) ** 2) ** 0.5
		self.alpha = math.acos((self.b ** 2 + self.c ** 2 - self.a ** 2) / (2 * self.b * self.c))
		self.preA = (self.a ** 2) * (self.b ** 2 + self.c ** 2 - self.a ** 2)
		self.preB = (self.b ** 2) * (self.c ** 2 + self.a ** 2 - self.b ** 2)
		self.preC = (self.c ** 2) * (self.a ** 2 + self.b ** 2 - self.c ** 2)
		self.Rx = (self.preA * Ax + self.preB * Bx + self.preC * Cx) / (self.preA + self.preB + self.preC)
		self.Ry = (self.preA * Ay + self.preB * By + self.preC * Cy) / (self.preA + self.preB + self.preC)
		self.R = self.a / (2 * math.sin(self.alpha))

def dist(Ax, Ay, Bx, By):
	return (((Ax - Bx) ** 2 + (Ay - By) ** 2) ** 0.5)

def judge(N, node, R, Rx, Ry):
	jud = "Yes"
	for i in range(N):
		x = node[i][0]; y = node[i][1]
		if dist(Rx, Ry, x, y) > R + 0.00001:
			jud = "No"

	if jud == "Yes":
		return R
	else:
		return INF


def main():
	N = int(input())
	node = []
	for i in range(N):
		x, y = getlist()
		node.append([x, y])

	if N == 2:
		R = dist(node[0][0], node[0][1], node[1][0], node[1][1]) / 2
		print(R)
		return

	ans = INF
	
	for i in range(N):
		for j in range(i + 1, N):
			Rx = (node[i][0] + node[j][0]) / 2
			Ry = (node[i][1] + node[j][1]) / 2
			R = dist(node[i][0], node[i][1], node[j][0], node[j][1]) / 2
			ans = min(ans, judge(N, node, R, Rx, Ry))

	
	for i in range(N):
		for j in range(i + 1, N):
			for k in range(j + 1, N):
				x1 = node[i][0]; y1 = node[i][1]
				x2 = node[j][0]; y2 = node[j][1]
				x3 = node[k][0]; y3 = node[k][1]
				if ccw(x1, y1, x2, y2, x3, y3) == 0:
					continue
				T = triangle(x1, y1, x2, y2, x3, y3)
				Rx = T.Rx; Ry = T.Ry; R = T.R
				ans = min(ans, judge(N, node, R, Rx, Ry))

	print(ans)

if __name__ == '__main__':
	main()

Solution (Python):
import bisect
import copy
import heapq
import math
import sys
from collections import *
from itertools import accumulate, combinations, permutations, product

def input():
    return sys.stdin.readline()[:-1]
def ruiseki(lst):
    return [0]+list(accumulate(lst))
mod=pow(10,9)+7
al=[chr(ord('a') + i) for i in range(26)]
direction=[[1,0],[0,1],[-1,0],[0,-1]]

n=int(input())
xy=[list(map(int,input().split())) for i in range(n)]

mi,ma=0,2000
com=list(combinations(range(n),2))

def dist(p,q):
    a,b=p
    c,d=q
    return math.sqrt((a-c)**2+(b-d)**2)

while ma-mi>10**(-9):
    half=(ma+mi)/2
    
    lst=[]
    for co in com:
        a,b=co
        ax,ay=xy[a]
        bx,by=xy[b]
        kyori=dist(xy[a],xy[b])
        if kyori/2>half:
            mi=half
            break
        else:
            h=math.sqrt(half**2-(dist(xy[a],xy[b])/2)**2)
            tyuuten=((ax+bx)/2,(ay+by)/2)
            tmp=h/kyori
            lst.append((tyuuten[0]+(ay-by)*tmp,tyuuten[1]-(ax-bx)*tmp))
            lst.append((tyuuten[0]-(ay-by)*tmp,tyuuten[1]+(ax-bx)*tmp))
    else:
        for co in lst:
            bx,by=co
            flag=1
            for i in range(n):
                ax,ay=xy[i]
                if (ax-bx)**2+(ay-by)**2>half**2+9e-7:
                    flag=0
                    break
            if flag==1:
                ma=half
                break
        else:
            mi=half
print(ma)

Solution (Python):
import bisect
import copy
import heapq
import math
import sys
from collections import *
from itertools import accumulate, combinations, permutations, product

def input():
    return sys.stdin.readline()[:-1]
def ruiseki(lst):
    return [0]+list(accumulate(lst))
mod=pow(10,9)+7
al=[chr(ord('a') + i) for i in range(26)]
direction=[[1,0],[0,1],[-1,0],[0,-1]]

n=int(input())
xy=[list(map(int,input().split())) for i in range(n)]

mi,ma=0,2000
com=list(combinations(range(n),2))

def dist(p,q):
    a,b=p
    c,d=q
    return math.sqrt((a-c)**2+(b-d)**2)

while ma-mi>10**(-9):
    half=(ma+mi)/2
    
    lst=[]
    for co in com:
        a,b=co
        ax,ay=xy[a]
        bx,by=xy[b]
        kyori=dist(xy[a],xy[b])/2
        if kyori>half:
            mi=half
            break
        else:
            h=math.sqrt(half**2-kyori**2)
            tyuuten=((ax+bx)/2,(ay+by)/2)
            tmp=h/2/kyori
            lst.append((tyuuten[0]+(ay-by)*tmp,tyuuten[1]-(ax-bx)*tmp))
            lst.append((tyuuten[0]-(ay-by)*tmp,tyuuten[1]+(ax-bx)*tmp))
    else:
        for co in lst:
            bx,by=co
            flag=1
            for i in range(n):
                ax,ay=xy[i]
                if (ax-bx)**2+(ay-by)**2>half**2+9e-7:
                    flag=0
                    break
            if flag==1:
                ma=half
                break
        else:
            mi=half
print(ma)

Solution (Python):
import math
N=int(input())
x=[0 for i in range(N)]
y=[0 for i in range(N)]
p=[]
for i in range(N):
    x[i],y[i]=map(int,input().split())
    p.append((x[i],y[i]))
def dist(p,q):
    a,b=p
    c,d=q
    return math.sqrt((a-c)**2+(b-d)**2)
def check(R):
    Q=[]
    for i in range(N):
        for j in range(i+1,N):
            d=dist(p[i],p[j])/2
            if R<d:
                return False
            h=math.sqrt(R*R-d*d)
            
            nx,ny=(x[j]-x[i],y[j]-y[i])
            nr=math.sqrt(nx**2+ny**2)
            dx,dy=h*nx/nr,h*ny/nr
            mx,my=(x[i]+x[j])/2,(y[i]+y[j])/2
            
            
            Q.append((mx+dy,my-dx))
            Q.append((mx-dy,my+dx))
    
    for pa in Q:
        flag=1
        for i in range(N):
            
            if dist(pa,p[i])>R+9e-7:
                flag=0
                break
        if flag==1:
            return True
    return False

low=0
high=2000
for line in range(50):
    mid=(low+high)/2
    if check(mid):
        high=mid
    else:
        low=mid
print(mid)

Solution (Python):
import sys
INF = 1 << 60
MOD = 10**9 + 7 
sys.setrecursionlimit(2147483647)
input = lambda:sys.stdin.readline().rstrip()
def resolve():
    n = int(input())
    C = [complex(*map(int, input().split())) for _ in range(n)]

    res = INF
    for k in range(n):
        p1 = C[k]
        for j in range(k):
            p2 = C[j]
            for i in range(j):
                p3 = C[i]
                a = abs(p2 - p3) ** 2
                b = abs(p3 - p1) ** 2
                c = abs(p1 - p2) ** 2
                numerator = a * (b + c - a) * p1 + b * (c + a - b) * p2 + c * (a + b - c) * p3
                denominator = a * (b + c - a) + b * (c + a - b) + c * (a + b - c)
                if denominator == 0.0:
                    continue
                center = numerator / denominator
                radius = -INF
                for p in C:
                    radius = max(radius, abs(center - p))
                res = min(radius, res)

    for j in range(n):
        for i in range(n):
            center = (C[i] + C[j]) / 2
            radius = -INF
            for p in C:
                radius = max(radius, abs(center - p))
            res = min(radius, res)

    print(res)
resolve()

Solution (Python):
import sys
INF = 1 << 60
MOD = 10**9 + 7 
sys.setrecursionlimit(2147483647)
input = lambda:sys.stdin.readline().rstrip()
def resolve():
    n = int(input())
    C = [complex(*map(int, input().split())) for _ in range(n)]

    res = INF
    for k in range(n):
        p1 = C[k]
        for j in range(k):
            p2 = C[j]
            for i in range(j):
                p3 = C[i]
                a = abs(p2 - p3)
                b = abs(p3 - p1)
                c = abs(p1 - p2)
                numerator = a**2 * (b**2 + c**2 - a**2) * p1 + b**2 * (c**2 + a**2 - b**2) * p2 + c**2 * (a**2 + b**2 - c**2) * p3
                denominator = a**2 * (b**2 + c**2 - a**2) + b**2 * (c**2 + a**2 - b**2) + c**2 * (a**2 + b**2 - c**2)
                if denominator < 0.001:
                    continue
                center = numerator / denominator
                radius = -INF
                for p in C:
                    radius = max(radius, abs(center - p))
                res = min(radius, res)

    for j in range(n):
        for i in range(n):
            center = (C[i] + C[j]) / 2
            radius = -INF
            for p in C:
                radius = max(radius, abs(center - p))
            res = min(radius, res)

    print(res)
resolve()

Solution (Python):
import numpy as np

EPS = 1e-12

class MinEnclosingCircle:
  def __init__(self, pts):
    self.pts = pts
    self.center = None
    self.radius = None

  def fit(self):
    self.center = (self.pts[0] + self.pts[1]) / 2
    self.radius = np.linalg.norm(self.pts[0] - self.pts[1]) / 2 + EPS
    for i in range(2, len(self.pts)):
      distance = np.linalg.norm(self.pts[i] - self.center)
      if distance >= self.radius:
        self.center, self.radius = self.find_second_point(i)

  def find_second_point(self, i):
    center = (self.pts[0] + self.pts[i]) / 2
    radius = np.linalg.norm(self.pts[0] - self.pts[i]) / 2 + EPS
    for j in range(1, i):
      distance = np.linalg.norm(self.pts[j] - center)
      if distance >= radius:
        new_center, new_radius = self.find_third_point(i, j)
        if new_radius > 0:
          center = new_center
          radius = new_radius
    return center, radius

  def find_third_point(self, i, j):
    center = (self.pts[j] + self.pts[i]) / 2
    radius = np.linalg.norm(self.pts[j] - self.pts[i]) / 2 + EPS
    for k in range(j):
      distance = np.linalg.norm(self.pts[k] - center)
      if distance >= radius:
        new_center, new_radius = self.find_circle_3pts(i, j, k)
        if new_radius > 0:
          center = new_center
          radius = new_radius
    return center, radius

  def find_circle_3pts(self, i, j, k):
    v1 = self.pts[j] - self.pts[i]
    v2 = self.pts[k] - self.pts[i]
    mid_point1 = (self.pts[j] + self.pts[i]) / 2
    c1 = np.dot(mid_point1, v1)
    mid_point2 = (self.pts[k] + self.pts[i]) / 2
    c2 = np.dot(mid_point2, v2)
    det = np.cross(v1, v2)
    if abs(det) <= EPS:
      d1 = np.sum(v1 * v1)
      d2 = np.sum(v2 * v2)
      d3 = np.sum((v1 - v2) ** 2)
      radius = np.sqrt(max(d1, d2, d3)) / 2 + EPS
      if d1 >= d2 and d1 >= d3:
        center = (self.pts[i] + self.pts[j]) / 2
      elif d2 >= d1 and d2 >= d3:
        center = (self.pts[i] + self.pts[k]) / 2
      else:
        center = (self.pts[j] + self.pts[k]) / 2
    else:
      center = np.zeros_like(v1)
      center[0] = c1 * v2[1] - c2 * v1[1]
      center[1] = v1[0] * c2 - v2[0] * c1
      center = center / det
      radius = np.linalg.norm(center - self.pts[i]) + EPS
    return center, radius

n = int(input())
point = [np.array(input().split(), dtype=int) for _ in range(n)]
mec = MinEnclosingCircle(point)
mec.fit()

print(mec.radius)

Solution (Python):
from math import sqrt

def colinear(x1, y1, x2, y2, x3, y3):
	a = x1 * (y2 - y3) +  x2 * (y3 - y1) +  x3 * (y1 - y2)
	if a == 0:
		return True
	else:
		return False 



def distance(x1, y1, x2, y2):
	xd = (x1 - x2) ** 2
	yd = (y1 - y2) ** 2
	return sqrt(xd + yd)


def check(cx, cy, r, x, y, n):
	for i  in range(n):
		d = distance(cx, cy, x[i], y[i])
		if d - r > 00.1:
			return False
	return True
	


def circlebydiaendpoints(x1, y1, x2, y2):
	cx = (x1 + x2) / 2
	cy = (y1 + y2) / 2
	r = distance(cx, cy, x1, y1)
	return cx, cy, r


def makeCircumcircle(x1, y1, x2, y2, x3, y3):

	ox = (min(min(x1, x2), x3) + max(min(x1, x2), x3)) / 2
	oy = (min(min(y1, y2), y3) + max(min(y1, y2), y3)) / 2
	ax = x1 - ox
	ay = y1 - oy
	bx = x2 - ox
	by = y2 - oy
	cx = x3 - ox
	cy = y3 - oy
	d = (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) * 2
	x = ((ax*ax + ay*ay) * (by - cy) + (bx*bx + by*by) * (cy - ay) + (cx*cx + cy*cy) * (ay - by)) / d
	y = ((ax*ax + ay*ay) * (cx - bx) + (bx*bx + by*by) * (ax - cx) + (cx*cx + cy*cy) * (bx - ax)) / d
	cx = ox + x
	cy = oy + y
	r = max(max(distance(cx, cy, x1, y1), distance(cx, cy, x2, y2)), distance(cx, cy, x3, y3))
	return cx, cy, r


def main():
	n = int(input())
	x = []
	y = []
	for _ in range(n):
		a, b = input().split(' ')
		a = float(a)
		b = float(b)
		x.append(a)
		y.append(b)
	ans = -1
	for i in range(n - 1):
		for j in range(i + 1, n):
			cx, cy, r = circlebydiaendpoints(x[i], y[i], x[j], y[j])
			if check(cx, cy, r, x, y, n):
				if ans == -1:
					ans = r
				elif ans > r:
					ans = r



	for i in range(n - 2):
		for j in range(i + 1, n - 1):
			for k in range(j + 1, n):
				if colinear(x[i], y[i], x[j], y[j], x[k], y[k]):
					continue
				else:
					cx, cy, r = makeCircumcircle(x[i], y[i], x[j], y[j], x[k], y[k])
					if check(cx, cy, r, x, y, n):
						if ans == -1:
							ans = r
						elif ans > r:
							ans = r


	print(ans)
	
	
if __name__ == "__main__":
	main()

Solution (Python):
n = int(input())
C = []
for i in range(n):
    x, y = map(int, input().split())
    C.append(complex(x, y))


def calc_kyori(O, C):
    res = 0
    for i in range(n):
        temp = abs(O - C[i])
        res = max(temp, res)
    return res


def sannbunn_first(x, my, MY):
    cnt = 100
    while cnt > 0:
        y1 = (2 * my + MY) / 3
        y2 = (my + 2 * MY) / 3
        v1 = calc_kyori(complex(x, y1), C)
        v2 = calc_kyori(complex(x, y2), C)
        if v2 < v1:
            my = y1
        else:
            MY = y2
        cnt -= 1
    return my


def sannbunn_second(mx, MX, my, MY):
    cnt = 100
    while cnt > 0:
        x1 = (2 * mx + MX) / 3
        x2 = (mx + 2 * MX) / 3
        yy1 = sannbunn_first(x1, my, MY)
        v1 = calc_kyori(complex(x1, yy1), C)
        yy2 = sannbunn_first(x2, my, MY)
        v2 = calc_kyori(complex(x2, yy2), C)
        if v2 < v1:
            mx = x1
        else:
            MX = x2
        cnt -= 1
    return mx, yy1, v1


if n == 2:
    ans = abs(C[0] - C[1]) / 2
else:
    rx, ry, ans = sannbunn_second(0, 1000, 0, 1000)
    
print(ans)

Solution (Python):
n = int(input())
XY = []
for i in range(n):
    x, y = map(int, input().split())
    XY.append((x, y))



def cross(a, b):
    return a.real * b.imag - b.real * a.imag



def dot(a, b):
    return a.real * b.real + a.imag * b.imag



def circumcenter(a, b, c):
    bunshi = (a - b) * dot(c, c) + (b - c) * dot(a, a) + (c - a) * dot(b, b)
    bunbo = (a - b) * c.conjugate() + (b - c) * a.conjugate() + (c - a) * b.conjugate()
    if bunbo == 0:
        return False, 0
    return True, bunshi / bunbo


kouho = []
for i in range(n - 1):
    a, b = XY[i]
    z1 = complex(a, b)
    for j in range(i + 1, n):
        c, d = XY[j]
        z2 = complex(c, d)
        kouho.append((z1 + z2) / 2)
if n >= 3:
    for i in range(n - 2):
        a, b = XY[i]
        z1 = complex(a, b)
        for j in range(i + 1, n - 1):
            c, d = XY[j]
            z2 = complex(c, d)
            for k in range(j + 1, n):
                e, f = XY[k]
                z3 = complex(e, f)
                flag, gaishin = circumcenter(z1, z2, z3)
                if flag:
                    kouho.append(gaishin)
ans = 10 ** 10
for z in kouho:
    temp = 0
    for i in range(n):
        x, y = XY[i]
        z1 = complex(x, y)
        temp = max(temp, abs(z - z1))
    ans = min(ans, temp)
print(ans)

Solution (Python):
from itertools import combinations
from copy import copy

class circle:
    def __init__(self, pts):
        assert(0 <= len(pts) and len(pts) <= 3)
        self.__pts = copy(pts)
        if len(pts) == 0:
            self.__cx, self.__cy, self.__r = float('inf'), float('inf'), 0
        elif len(pts) == 1:
            self.__cx, self.__cy, self.__r = pts[0], 0
        elif len(pts) == 2:
            self.__cx, self.__cy, self.__r = self.__c2()
        else:
            self.__cx, self.__cy, self.__r = self.__c3()

    def __c2(self):
        assert(len(self.__pts) == 2)
        (x0, y0), (x1, y1) = self.__pts
        return ((x0+x1)/2, (y0+y1)/2,
                0.5*((x0-x1)*(x0-x1) + (y0-y1)*(y0-y1))**0.5)

    def __c3(self):
        assert(len(self.__pts) == 3)
        return self.__o3() if self.__obtuse() else self.__a3()

    def __o3(self):
        p, q = max(((xi-xj)**2 + (yi-yj)**2, ((xi, yi), (xj, yj)))
                   for (xi, yi), (xj, yj) in combinations(self.__pts, 2))[1]
        self.__pts = [p, q]
        return self.__c2()

    def __obtuse(self):
        """ Determine whether given tripe forms obtuse, cf.,
                https://math.stackexchange.com/questions/3200316/
        """
        (x0, y0), (x1, y1), (x2, y2) = self.__pts
        a = (x0-x1)**2 + (y0-y1)**2
        b = (x1-x2)**2 + (y1-y2)**2
        c = (x2-x0)**2 + (y2-y0)**2
        return a > b+c or b > c+a or c > a+b

    def __a3(self):
        """ Equation of circle when three points on the circle are given.
            https://www.geeksforgeeks.org/equation-of-circle-when-three-points-on-the-circle-are-given/
        """
        (x0, y0), (x1, y1), (x2, y2) = self.__pts

        x01, x02 = x0 - x1, x0 - x2
        y01, y02 = y0 - y1, y0 - y2

        y20, y10 = y2 - y0, y1 - y0
        x20, x10 = x2 - x0, x1 - x0

        sx02, sy02 = x0*x0 - x2*x2, y0*y0 - y2*y2
        sx10, sy10 = x1*x1 - x0*x0, y1*y1 - y0*y0

        f = (sx02*x01+sy02*x01 + sx10*x02+sy10*x02) / (2*(y20*x01 - y10*x02))
        g = (sx02*y01+sy02*y01 + sx10*y02+sy10*y02) / (2*(x20*y01 - x10*y02))
        c = -x0*x0 - y0*y0 - 2*g*x0 - 2*f*y0

        return (-g, -f, (f*f + g*g - c)**0.5)

    def draw(self):
        return self.__r


    def contains(self, x, y):
        return ((x-self.__cx)**2 + (y-self.__cy)**2)**0.5 <= self.__r + 1e-9

    def pts(self):
        return set(self.__pts)


def welzl(pts):
    """ Find the smallest enclosing circle for given a set of points in plane.
        from https://en.wikipedia.org/wiki/Smallest-circle_problem
    """
    return _welzl(set(pts), set())


def _welzl(p, r):
    """ Find the circle in recursively
        input:  Finite sets P and R of points in the plane |R| <= 3.
        output: Minimal disk enclosing P with R on the boundary.
    """
    if len(p) == 0 or len(r) == 3:
        return circle(list(r))
    q = sample(p, 1)[0]   
    d = _welzl(p - {q}, r)
    if d.contains(*q):
        return d
    return _welzl(p - {q}, r | {q})


def msw(pts):
    return _msw(set(pts), set(sample(pts, 3)))


def _msw(p, r):
    """ Find the circle in recursively
        input:  Finite sets P and R of points in the plane |R| = 3.
        output: Minimal disk enclosing P U R.
        from wikipedia.
    """
    if len(p) == 0:
        return circle(list(r))
    q = sample(p, 1)[0]   
    d = _msw(p - {q}, r)
    if d.contains(*q):
        return d
    x = (r | {q}) - welzl(r | {q}).pts()
    return _msw((p - {q}) | x, (r | {q}) - x)

if __name__ == '__main__':

    n=int(input())
    e=[]
    for i in range(n):
        x,y=map(int,input().split())
        e.append((x,y))

if n==2:
    d = circle(e)
    print(d.draw())
    exit()
ans=0

for i in range(n-2):
    for j in range(i+1,n-1):
        for k in range(j+1,n):
            pts=(e[i],e[j],e[k])
            d=circle(pts)
            ans=max(ans,d.draw())

print(ans)

Solution (Python):

import sys
import cmath
sys.setrecursionlimit(10**9)
INF=10**18
MOD=10**9+7
eps=1e-8
input=lambda: sys.stdin.readline().rstrip()
YesNo=lambda b: bool([print('Yes')] if b else print('No'))
YESNO=lambda b: bool([print('YES')] if b else print('NO'))
int1=lambda x:int(x)-1

def main():
    N=int(input())
    x,y=[0]*N,[0]*N
    for i in range(N):
        x[i],y[i]=map(int,input().split())
    
    def get_Intersection(c1,c2,R):
        mid=(c1+c2)/2
        if R<abs(c1-mid):
            return None,None
        d=(R**2-abs(c1-mid)**2)**0.5
        p=cmath.polar((c1-c2)*1j)
        v=cmath.rect(d,p[1])
        i1=mid+v
        i2=mid-v
        return i1,i2
    
    def nibutan(ok,ng):
        for _ in range(100):
            mid = (ok + ng) / 2
            if solve(mid):
                ok = mid
            else:
                ng = mid
        return ok
    
    def solve(mid):
        for i in range(N):
            for j in range(N):
                if i==j:
                    continue
                c1=x[i]+y[i]*1j
                c2=x[j]+y[j]*1j
                i1,i2=get_Intersection(c1,c2,mid)
                if i1 is None:
                    continue
                for k in range(N):
                    i3=x[k]+y[k]*1j
                    if abs(i3-i1)>mid+eps:
                        break
                else:
                    return True
                for k in range(N):
                    i3=x[k]+y[k]*1j
                    if abs(i3-i2)>mid+eps:
                        break
                else:
                    return True
        return False
    
    print(nibutan(2000,0))
    


    

if __name__ == '__main__':
    main()