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 copy import copy
from itertools import combinations
from random import random, sample, seed


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, c='g', ax=None):
        if ax is None:
            ax = plt.gca()
        c = plt.Circle((self.__cx, self.__cy), self.__r, fill=False, color=c)
        ax.add_artist(c)

    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)

    @property
    def center(self):
        return (self.__cx, self.__cy)

    @property
    def radius(self):
        return self.__r


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

        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(set(p) - {q}, r)
    if d.contains(*q):
        return d
    return welzl(set(p) - {q}, r | {q})


def msw(pts):
    """ Matousek, Sharir, Welzl's Algorithm for the minimum-circle problem """
    if len(pts) <= 3:
        return circle(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)


def gen(n=20, s=1):
    seed(s)
    return [(random(), random()) for i in range(n)]


if __name__ == '__main__':
    N = int(input())
    pts = []
    for _ in range(N):
        x, y = map(int, input().split())
        pts.append((x,y))

    print(msw(pts).radius)

Solution (Python):
def gaishin(p,q,r):
  a = ((q[0]-r[0])**2+(q[1]-r[1])**2)**(0.5) 
  b = ((r[0]-p[0])**2+(r[1]-p[1])**2)**(0.5) 
  c = ((p[0]-q[0])**2+(p[1]-q[1])**2)**(0.5) 
  cosp = (b**2+c**2-a**2)/(2*b*c)
  cosq = (c**2+a**2-b**2)/(2*c*a)
  cosr = (a**2+b**2-c**2)/(2*a*b)
  if cosp >= 1: cosp = 1
  if cosp <= -1: cosp = -1
  if cosq >= 1: cosq = 1
  if cosq <= -1: cosq = -1
  if cosr >= 1: cosr = 1
  if cosr <= -1: cosr = -1
  sinp = (1-cosp**2)**(0.5)
  sinq = (1-cosq**2)**(0.5)
  sinr = (1-cosr**2)**(0.5)
  sin2p = 2*sinp*cosp
  sin2q = 2*sinq*cosq
  sin2r = 2*sinr*cosr
  if sin2p+sin2q+sin2r == 0:
    infinite = float("inf")
    return (infinite, infinite)
  x = (sin2p*p[0]+sin2q*q[0]+sin2r*r[0])/(sin2p+sin2q+sin2r)
  y = (sin2p*p[1]+sin2q*q[1]+sin2r*r[1])/(sin2p+sin2q+sin2r)
  
  return (x,y)

N = int(input())
A = []
for i in range(N):
  x,y = map(int,input().split())
  A.append((x,y))
INF = float("inf")
ans = INF
for i in range(N):
  for j in range(i+1,N):
    z = ((A[i][0]+A[j][0])/2,(A[i][1]+A[j][1])/2) 
    dx = A[i][0]-A[j][0]
    dy = A[i][1]-A[j][1]
    hankei = ((dx**2+dy**2)**(0.5))/2
    Flag = True
    for t in range(N):
      if t == i or t == j: 
        continue
      temp = ((z[0]-A[t][0])**2+(z[1]-A[t][1])**2)**(0.5)
      hankei = max(hankei,temp)
    ans = min(ans,hankei)

if N == 2:
  print(ans)
  exit()
for i in range(N):
  for j in range(i+1,N):
    for k in range(j+1,N):
      p = A[i]; q = A[j]; r = A[k]
      z = gaishin(p,q,r)
      
      hankei = 0
      for t in range(N):
        temp = ((z[0]-A[t][0])**2+(z[1]-A[t][1])**2)**(0.5)
        hankei = max(hankei,temp)
      ans = min(ans,hankei)
print(ans)

Solution (Python):
def gaishin(p,q,r):
  a = ((q[0]-r[0])**2+(q[1]-r[1])**2)**(0.5) 
  b = ((r[0]-p[0])**2+(r[1]-p[1])**2)**(0.5) 
  c = ((p[0]-q[0])**2+(p[1]-q[1])**2)**(0.5) 
  cosp = (b**2+c**2-a**2)/(2*b*c)
  cosq = (c**2+a**2-b**2)/(2*c*a)
  cosr = (a**2+b**2-c**2)/(2*a*b)
  if cosp >= 1: cosp = 1
  if cosp <= -1: cosp = -1
  if cosq >= 1: cosq = 1
  if cosq <= -1: cosq = -1
  if cosr >= 1: cosr = 1
  if cosr <= -1: cosr = -1
  sinp = (1-cosp**2)**(0.5)
  sinq = (1-cosq**2)**(0.5)
  sinr = (1-cosr**2)**(0.5)
  sin2p = 2*sinp*cosp
  sin2q = 2*sinq*cosq
  sin2r = 2*sinr*cosr
  
  
  
  if sin2p+sin2q+sin2r < pow(10,-7):
    infinite = float("inf")
    return (infinite, infinite)
  x = (sin2p*p[0]+sin2q*q[0]+sin2r*r[0])/(sin2p+sin2q+sin2r)
  y = (sin2p*p[1]+sin2q*q[1]+sin2r*r[1])/(sin2p+sin2q+sin2r)
  R = a/(2*sinp) 
  return (x,y)

N = int(input())
A = []
for i in range(N):
  x,y = map(int,input().split())
  A.append((x,y))
INF = float("inf")
ans = INF
for i in range(N):
  for j in range(i+1,N):
    z = ((A[i][0]+A[j][0])/2,(A[i][1]+A[j][1])/2) 
    dx = A[i][0]-A[j][0]
    dy = A[i][1]-A[j][1]
    hankei = ((dx**2+dy**2)**(0.5))/2
    Flag = True
    for t in range(N):
      if t == i or t == j: 
        continue
      temp = ((z[0]-A[t][0])**2+(z[1]-A[t][1])**2)**(0.5)
      hankei = max(hankei,temp)
    ans = min(ans,hankei)

if N == 2:
  print(ans)
  exit()
for i in range(N):
  for j in range(i+1,N):
    for k in range(j+1,N):
      p = A[i]; q = A[j]; r = A[k]
      z = gaishin(p,q,r)
      
      hankei = 0
      for t in range(N):
        temp = ((z[0]-A[t][0])**2+(z[1]-A[t][1])**2)**(0.5)
        hankei = max(hankei,temp)
      ans = min(ans,hankei)
print(ans)

Solution (Python):
from fractions import Fraction
from math import sqrt

n = int(input())
p = []

for i in range(n):
    p.append(list(map(int, input().split())))

ans = 0

for i in range(n-1):
    for j in range(i+1,n):
        flag = True
        x1,y1 = p[i][0],p[i][1]
        x2,y2 = p[j][0],p[j][1]
        d = (x2 - x1)**2 + (y2 - y1)**2
        for k in range(n):
            if (i == k or j == k):
                continue
            x3,y3 = p[k][0],p[k][1]
            e = (x3-x1)**2 + (y3-y1)**2 + (x3-x2)**2 + (y3-y2)**2
            if (d < e):
                flag = False
                break
        if flag:
            ans = sqrt(d)/2.0

if (ans != 0):
    print(ans)
    exit()

ans = 100000000

for i in range(n-2):
    for j in range(i+1,n-1):
        for k in range(j+1,n):
            a,b = p[i][0],p[i][1]
            c,d = p[j][0],p[j][1]
            e,f = p[k][0],p[k][1]
            if ((a == c == e) or (b == d == f)):
                continue
            ny = (e-a) * (a*a + b*b - c*c - d*d) - (c-a) * (a*a + b*b - e*e - f*f)
            dy = 2 * (e-a) * (b-d) - 2 * (c-a) * (b-f)
            if (dy == 0):
                continue
            py = Fraction(ny,dy)
            if (c != a):
                px = 2*(b-d) * py - a*a - b*b + c*c + d*d
                px /= 2*(c-a)
            else:
                px = 2*(b-f) * py - a*a - b*b + e*e + f*f
                px /= 2*(e-a)
            r = (c-px)**2 + (d-py)**2
            flag = True
            for l in range(n):
                if (l == i or l == j or l == k):
                    continue
                x,y = p[l][0],p[l][1]
                s = (x-px)**2 + (y-py)**2
                if (r < s):
                    flag = False
                    break
            if (flag):
                q = r.numerator / r.denominator
                q = sqrt(q)
                ans = min(ans,q)

print(ans)

Solution (Python):
import itertools
import math


def main():
    N = int(input())
    coords = [list(map(int, input().split(' '))) for _ in range(N)]
    ans_squared = 10 ** 12
    
    for i, j in itertools.combinations(range(N), 2):
        x1, y1 = coords[i]
        x2, y2 = coords[j]
        d_squared = (x1 - x2)**2 + (y1 - y2)**2
        cover = True
        for n in range(N):
            x, y = coords[n]
            if (2 * x - (x1 + x2))**2 + (2 * y - (y1 + y2))**2 > d_squared:
                cover = False
                break
        if cover:
            ans_squared = min(ans_squared, d_squared / 4)
    if N == 2:
        print(math.sqrt(ans_squared))
        return
    
    for i, j, k in itertools.combinations(range(N), 3):
        x1, y1 = coords[i]
        x2, y2 = coords[j]
        x3, y3 = coords[k]
        dx12, dx23, dy12, dy23 = x1 - x2, x2 - x3, y1 - y2, y2 - y3
        if dx12 * dy23 == dy12 * dx23:
            
            continue
        v1, v2, v3 = x1 ** 2 + y1 ** 2, x2 ** 2 + y2 ** 2, x3 ** 2 + y3 ** 2
        scale = 2 * (dx12 * dy23 - dy12 * dx23)
        
        scaled_cx = dy23 * v1 - (dy12 + dy23) * v2 + dy12 * v3
        scaled_cy = - dx23 * v1 + (dx12 + dx23) * v2 - dx12 * v3
        r_squared = (x1 - scaled_cx / scale)**2 + (y1 - scaled_cy / scale)**2
        cover = True
        for n in range(N):
            x, y = coords[n]
            if (scale * x - scaled_cx)**2 + (scale * y - scaled_cy)**2 \
                    > (scale * x1 - scaled_cx)**2 + (scale * y1 - scaled_cy)**2:
                cover = False
        if cover:
            ans_squared = min(ans_squared, r_squared)
    print(math.sqrt(ans_squared))


if __name__ == '__main__':
    main()

Solution (Python):
import sys,bisect,string,math,time,functools,random,fractions
from heapq import heappush,heappop,heapify
from collections import deque,defaultdict,Counter
from itertools import permutations,combinations,groupby
rep=range;R=range
def Golf():n,*t=map(int,open(0).read().split())
def I():return int(input())
def S_():return input()
def IS():return input().split()
def LS():return [i for i in input().split()]
def MI():return map(int,input().split())
def LI():return [int(i) for i in input().split()]
def LI_():return [int(i)-1 for i in input().split()]
def NI(n):return [int(input()) for i in range(n)]
def NI_(n):return [int(input())-1 for i in range(n)]
def StoLI():return [ord(i)-97 for i in input()]
def ItoS(n):return chr(n+97)
def LtoS(ls):return ''.join([chr(i+97) for i in ls])
def RA():return map(int,open(0).read().split())
def GI(V,E,ls=None,Directed=False,index=1):
    org_inp=[];g=[[] for i in range(V)]
    FromStdin=True if ls==None else False
    for i in range(E):
        if FromStdin:
            inp=LI()
            org_inp.append(inp)
        else:
            inp=ls[i]
        if len(inp)==2:
            a,b=inp;c=1
        else:
            a,b,c=inp
        if index==1:a-=1;b-=1
        aa=(a,c);bb=(b,c);g[a].append(bb)
        if not Directed:g[b].append(aa)
    return g,org_inp
def GGI(h,w,search=None,replacement_of_found='.',mp_def={'
    
    mp=[boundary]*(w+2);found={}
    for i in R(h):
        s=input()
        for char in search:
            if char in s:
                found[char]=((i+1)*(w+2)+s.index(char)+1)
                mp_def[char]=mp_def[replacement_of_found]
        mp+=[boundary]+[mp_def[j] for j in s]+[boundary]
    mp+=[boundary]*(w+2)
    return h+2,w+2,mp,found
def TI(n):return GI(n,n-1)
def accum(ls):
    rt=[0]
    for i in ls:rt+=[rt[-1]+i]
    return rt
def bit_combination(n,base=2):
    rt=[]
    for tb in R(base**n):s=[tb//(base**bt)%base for bt in R(n)];rt+=[s]
    return rt
def gcd(x,y):
    if y==0:return x
    if x%y==0:return y
    while x%y!=0:x,y=y,x%y
    return y
def YN(x):print(['NO','YES'][x])
def Yn(x):print(['No','Yes'][x])
def show(*inp,end='\n'):
    if show_flg:print(*inp,end=end)

mo=10**9+7
inf=float('inf')
FourNb=[(-1,0),(1,0),(0,1),(0,-1)];EightNb=[(-1,0),(1,0),(0,1),(0,-1),(1,1),(-1,-1),(1,-1),(-1,1)];compas=dict(zip('WENS',FourNb));cursol=dict(zip('LRUD',FourNb))
l_alp=string.ascii_lowercase

read=sys.stdin.buffer.read;readline=sys.stdin.buffer.readline;input=lambda:sys.stdin.readline().rstrip()

show_flg=False
show_flg=True

ans=0


def check(x,y,r):
    for s,t in c:
        if (x-s)**2+(y-t)**2>r+1e-6:
            return False
    return True
    
 
def center_of_3points(x1,y1,x2,y2,x3,y3):
    cx = ((y1-y3)*(y1**2 -y2**2 +x1**2 -x2**2) -(y1-y2)*(y1**2 -y3**2 +x1**2 -x3**2)) / (2*(y1-y3)*(x1-x2)-2*(y1-y2)*(x1-x3))
    cy = ((x1-x3)*(x1**2 -x2**2 +y1**2 -y2**2) -(x1-x2)*(x1**2 -x3**2 +y1**2 -y3**2)) / (2*(x1-x3)*(y1-y2)-2*(x1-x2)*(y1-y3))
    r=((x1-cx)**2+(y1-cy)**2)
    return cx,cy,r
 
def rotate2D(x,y,q):
    s=x*math.cos(q)-y*math.sin(q)
    t=x*math.sin(q)+y*math.cos(q)
    return (s,t)


def cir2(x1,y1,x2,y2):
    p,q=(x1+x2)/2,(y1+y2)/2
    r=((x1-p)**2+(y1-q)**2)
    return p,q,r
 
ans=inf
 
n=I()
c=[]
CC=[]
rx=random.randint(1,10**3)
ry=random.randint(1,10**3)
for i in range(n):
    x,y=LI()
    x+=rx
    y+=ry
    CC.append((x,y))
 
 
def solve():
    ans=inf
    for i in range(n):
        x,y=CC[i]
        c.append((x,y))
    
    for i in range(n):
        for j in range(i+1,n):
                po=c[i]+c[j]
                x1,y1,x2,y2=po
                p,q,r=cir2(x1,y1,x2,y2)
                
                if check(p,q,r):
                    ans=min(ans,r)
    
    for i in range(n):
        for j in range(i+1,n):
            for k in range(j+1,n):
                x1,y1,x2,y2,x3,y3=CC[i]+CC[j]+CC[k]
                a,b= (x2-x1),(y2-y1)
                C,D= (x3-x1),(y3-y1)
                if a*D==b*C:
                    continue
                x1,y1,x2,y2,x3,y3=c[i]+c[j]+c[k]
                p,q,r=center_of_3points(x1,y1,x2,y2,x3,y3)
                
                if check(p,q,r):
                    ans=min(ans,r)

    return ans

print(solve()**0.5)

Solution (Python):


import numpy as np

N=int(input())
x=np.zeros([N,2])
for i in range(N):
  x[i,:] = list(map(float,input().split()))
  
def norm(x,y):return np.sqrt(((x - y)**2).sum())

x0 = x.mean(0)
old = x0
rate = 0.06
decay = 0.999
epsilon=1e-15

while 1:
  d = np.sqrt(((x - old)**2).sum(1))
  max_, argmax = d.max(), d.argmax()
  direction = x[argmax, :] - old
  new = old + direction*rate
  rate *= decay
  if norm(new,old) <= epsilon:break
  old = new

print(norm(direction, 0))

Solution (Python):

import sys
import itertools
import numpy as np
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

N = int(readline())
XY = []
for i in range(N):
    x,y = map(int,readline().split())
    XY.append(x+y*1j)
    
center = []    
for z0,z1 in itertools.combinations(XY,2):
    center.append((z0+z1)/2)
    
for z0,z1,z2 in itertools.combinations(XY,3):
    a = z1-z0
    b = z2-z0
    nume = a*b*(a.conjugate()-b.conjugate())
    deno = a.conjugate()*b - a*b.conjugate()
    if deno == 0:
        continue
    c = nume / deno
    c += z0
    center.append(c)

XY = np.array(XY)
center = np.array(center)

ans = abs(XY[:,None] - center[None,:]).max(axis=0).min()
print(ans)