Problem:
 Score : 600 points    Problem Statement Takahashi wants to grill N pieces of meat on a grilling net, which can be seen as a two-dimensional plane. The coordinates of the i -th piece of meat are \left(x_i, y_i\right) , and its hardness is c_i .  Takahashi can use one heat source to grill the meat. If he puts the heat source at coordinates \left(X, Y\right) , where X and Y are real numbers, the i -th piece of meat will be ready to eat in c_i \times \sqrt{\left(X - x_i\right)^2 + \left(Y-y_i\right)^2} seconds.  Takahashi wants to eat K pieces of meat. Find the time required to have K or more pieces of meat ready if he put the heat source to minimize this time.      Constraints  All values in input are integers.  1 \leq N \leq 60  1 \leq K \leq N  -1000 \leq x_i , y_i \leq 1000  \left(x_i, y_i\right) \neq \left(x_j, y_j\right) \left(i \neq j \right)  1 \leq c_i \leq 100         Input Input is given from Standard Input in the following format:  N  K  x_1  y_1  c_1  \vdots  x_N  y_N  c_N       Output Print the answer.  It will be considered correct if its absolute or relative error from our answer is at most 10^{-6} .        Sample Input 1 4 3
-1 0 3
0 0 3
1 0 2
1 1 40      Sample Output 1 2.4  If we put the heat source at \left(-0.2, 0\right) , the 1 -st, 2 -nd, and 3 -rd pieces of meat will be ready to eat within 2.4 seconds. This is the optimal place to put the heat source.       Sample Input 2 10 5
-879 981 26
890 -406 81
512 859 97
362 -955 25
128 553 17
-885 763 2
449 310 57
-656 -204 11
-270 76 40
184 170 16      Sample Output 2 7411.2252

Solution (Python):
import sys

input = sys.stdin.buffer.readline
import math
from itertools import combinations
EPS = 10**-7


class P2():
    def __init__(self, x: float, y: float):
        self.x = x
        self.y = y
        self.norm = (self.x**2 + self.y**2)**0.5

    def __add__(self, other):
        return P2(self.x + other.x, self.y + other.y)

    def __sub__(self, other):
        return P2(self.x - other.x, self.y - other.y)

    def __mul__(self, c: float):
        return P2(self.x * c, self.y * c)

    def __truediv__(self, c: float):
        return P2(self.x / c, self.y / c)

    def __eq__(self, other):
        return (abs(self.x - other.x) < EPS) and (abs(self.y - other.y) < EPS)

    def dot(self, other) -> float:
        return self.x * other.x + self.y * other.y

    def det(self, other) -> float:
        return self.x * other.y - self.y * other.x

    def rot90(self):
        return P2(-self.y, self.x)


class Circle():
    def __init__(self, p: P2, r: float):
        self.p = p  
        self.r = r  

    def contains(self, p: P2, allow_on_edge: bool) -> bool:
        if allow_on_edge:
            return (self.p - p).norm < self.r + EPS
        else:
            return (self.p - p).norm < self.r - EPS

    def cct(self, other) -> int:
        
        if self.p == other.p:
            
            if abs(self.r - other.r) < EPS:
                return float('inf')  
            else:
                return 0  

        d = other.p - self.p
        if d.norm > self.r + other.r + EPS:
            return 4  
        elif d.norm > self.r + other.r - EPS:
            return 3  
        elif d.norm > abs(self.r - other.r) + EPS:
            return 2  
        elif d.norm > abs(self.r - other.r) - EPS:
            return 1  
        else:
            return 0  

    def intersection_points_with_circle(self, other):
        cct = self.cct(other)
        d = other.p - self.p

        if cct == 1:
            
            if self.r > other.r:
                return [self.p + d * (self.r / d.norm)]
            else:
                return [self.p - d * (self.r / d.norm)]
        elif cct == 3:
            
            return [self.p + d * (self.r / d.norm)]
        elif cct == 2:
            
            x = (self.r**2 - other.r**2 + d.norm**2) / (2 * d.norm)
            h = math.sqrt(max(self.r**2 - x**2, 0.0))
            v = d.rot90() * (h / d.norm)
            dx = self.p + d * (x / d.norm)
            return [dx + v, dx - v]

        return []


n, k = map(int, input().split())
p = []
c = []
for _ in range(n):
    x, y, cc = map(float, input().split())
    p.append(P2(x, y))
    c.append(cc)


def check(r):
    circles = [Circle(p[i], r / c[i]) for i in range(n)]

    
    cand = [pi for pi in p]
    for c_i, c_j in combinations(circles, r=2):
        cand_ij = c_i.intersection_points_with_circle(c_j)
        cand.extend(cand_ij)

    for cand_p in cand:
        cnt = 0
        for c_i in circles:
            if (c_i.p - cand_p).norm < c_i.r + EPS:
                cnt += 1
                if cnt >= k:
                    return True

    return False


lb = 0.0  
ub = 3000 * 100  
while ub - lb > EPS:
    mid = (ub + lb) / 2
    if check(mid):
        ub = mid
    else:
        lb = mid

print(ub)

Solution (Python):
import sys

input = sys.stdin.buffer.readline
import math
from itertools import combinations
EPS = 10**-7


class P2():
    def __init__(self, x: float, y: float):
        self.x = x
        self.y = y
        self.norm = (self.x**2 + self.y**2)**0.5

    def __add__(self, other):
        return P2(self.x + other.x, self.y + other.y)

    def __sub__(self, other):
        return P2(self.x - other.x, self.y - other.y)

    def __mul__(self, c: float):
        return P2(self.x * c, self.y * c)

    def __truediv__(self, c: float):
        return P2(self.x / c, self.y / c)

    def __eq__(self, other):
        return (abs(self.x - other.x) < EPS) and (abs(self.y - other.y) < EPS)

    def dot(self, other) -> float:
        return self.x * other.x + self.y * other.y

    def det(self, other) -> float:
        return self.x * other.y - self.y * other.x

    def rot90(self):
        return P2(-self.y, self.x)


class Circle():
    def __init__(self, p: P2, r: float):
        self.p = p  
        self.r = r  

    def contains(self, p: P2, allow_on_edge: bool) -> bool:
        if allow_on_edge:
            return (self.p - p).norm < self.r + EPS
        else:
            return (self.p - p).norm < self.r - EPS

    def cct(self, other) -> int:
        
        if self.p == other.p:
            
            return float('inf') if abs(self.r - other.r) < EPS else 0

        d = other.p - self.p
        if d.norm > self.r + other.r + EPS:
            
            return 4
        elif d.norm > self.r + other.r - EPS:
            
            return 3
        elif d.norm > abs(self.r - other.r) + EPS:
            
            return 2
        elif d.norm > abs(self.r - other.r) - EPS:
            
            return 1
        else:
            
            return 0

    def intersection_points_with_circle(self, other):
        cct = self.cct(other)
        d = other.p - self.p
        intersections = []
        if cct == 1:
            
            if self.r > other.r:
                intersections.append(self.p + d * (self.r / d.norm))
            else:
                intersections.append(self.p - d * (self.r / d.norm))
        elif cct == 3:
            
            intersections.append(self.p + d * (self.r / d.norm))
        elif cct == 2:
            
            x = (self.r**2 - other.r**2 + d.norm**2) / (2 * d.norm)
            h = math.sqrt(max(self.r**2 - x**2, 0.0))
            v = d.rot90() * (h / d.norm)
            dx = self.p + d * (x / d.norm)
            intersections.append(dx + v)
            intersections.append(dx - v)

        return intersections


n, k = map(int, input().split())
p = []
c = []
for _ in range(n):
    x, y, cc = map(float, input().split())
    p.append(P2(x, y))
    c.append(cc)

dist = [[(p[i] - p[j]).norm for j in range(n)] for i in range(n)]


def check(r):
    circles = [Circle(p[i], r / c[i]) for i in range(n)]

    
    cand = [pi for pi in p]
    for c_i, c_j in combinations(circles, r=2):
        cand_ij = c_i.intersection_points_with_circle(c_j)
        cand.extend(cand_ij)

    for cand_p in cand:
        cnt = 0
        for c_i in circles:
            if (c_i.p - cand_p).norm < c_i.r + EPS:
                cnt += 1
                if cnt >= k:
                    return True

    return False


lb = 0.0  
ub = 3000 * 100  
while ub - lb > EPS:
    mid = (ub + lb) / 2
    if check(mid):
        ub = mid
    else:
        lb = mid

print(ub)

Solution (Python):
import sys
from itertools import combinations
from math import atan2, cos, sin


def apollonius(p1, c1, p2, c2):
    m1 = p1 + (p2 - p1) * c2 / (c1 + c2)
    m2 = p1 + (p2 - p1) * c2 / (c2 - c1)
    m = (m1 + m2) / 2
    r = abs(m1 - m)
    return m, r


def apollonius_intersections(p1, c1, p2, c2, p3, c3):
    m12, r12 = apollonius(p1, c1, p2, c2)
    m13, r13 = apollonius(p1, c1, p3, c3)
    v = m13 - m12
    d = abs(v)
    if d > r12 + r13 or d < abs(r12 - r13):
        return None, None
    theta = atan2(v.imag, v.real)
    xx = (r12 ** 2 - r13 ** 2 + d ** 2) / (2 * d)
    s = (r12 + r13 + d) / 2
    yy = 2 * (s * (s - r12) * (s - r13) * (s - d)) ** 0.5 / d
    st = sin(theta)
    ct = cos(theta)
    e1 = (xx * ct - yy * st) + (xx * st + yy * ct) * 1j + m12
    e2 = (xx * ct + yy * st) + (xx * st - yy * ct) * 1j + m12

    t1 = abs(e1 - p1) * c1
    t2 = abs(e2 - p1) * c1
    if t1 < t2:
        return e1, t1
    return e2, t2


def get_circumscribed_center(p1, p2, p3):
    A = abs(p2 - p3) ** 2
    B = abs(p3 - p1) ** 2
    C = abs(p1 - p2) ** 2
    T = A * (B + C - A)
    U = B * (C + A - B)
    W = C * (A + B - C)
    try:
        return (T * p1 + U * p2 + W * p3) / (T + U + W)
    except:
        return 0j


def solve(n, k, xyc):
    if k == 1:
        return 0

    xxx = xyc[0::3]
    yyy = xyc[1::3]
    ccc = xyc[2::3]
    beefs = [(i, x + y * 1j, c) for i, (x, y, c) in enumerate(zip(xxx, yyy, ccc))]

    ans = 1e18
    for (i1, p1, c1), (i2, p2, c2) in combinations(beefs, 2):
        
        
        px = p1 + (p2 - p1) * (c2 / (c1 + c2))
        t = abs(px - p1) * c1 + 1e-8

        ok = 0
        ng = []
        for i, p, c in beefs:
            if abs(px - p) * c <= t:
                ok += 1
            else:
                ng.append(i)
        if ok >= k:
            ans = min(ans, t)
            continue
        
        if ans <= t:
            continue

        
        
        
        
        
        
        for i3 in ng:
            if i3 < i2:
                continue
            _, p3, c3 = beefs[i3]
            if c1 == c2:
                if c2 == c3:
                    e = get_circumscribed_center(p1, p2, p3)
                    t = abs(e - p1) * c1
                else:
                    e, t = apollonius_intersections(p3, c3, p1, c1, p2, c2)
            elif c1 == c3:
                e, t = apollonius_intersections(p2, c2, p1, c1, p3, c3)
            else:
                e, t = apollonius_intersections(p1, c1, p2, c2, p3, c3)

            if e is None:
                continue

            t += 1e-8
            ok2 = 0
            for i, p, c in beefs:
                if abs(e - p) * c <= t:
                    ok2 += 1
            if ok2 >= k:
                ans = min(ans, t)

    return ans


n, k = map(int, input().split())
xyc = list(map(int, sys.stdin.read().split()))
print(solve(n, k, xyc))

Solution (Python):
import math
def intersectionCC0(x1, y1, r1, r2):
    r = x1**2 + y1**2
    t = (r + r1**2 - r2**2) / (2*r)
    dd = r1**2/r - t**2
    if - (10**-8) < dd < 0: dd = 0
    if dd < 0: return []
    x, y = t * x1, t * y1
    if dd == 0: return [(x, y)]
    sq = math.sqrt(dd)
    dx, dy = y1 * sq, -x1 * sq
    return [(x+dx, y+dy), (x-dx, y-dy)]

def intersectionCC(x1, y1, x2, y2, r1, r2):
    return [(x1+x, y1+y) for x, y in intersectionCC0(x2-x1, y2-y1, r1, r2)]

def chk(t):
    L = [a for a in L0]
    for i, (x1, y1, c1) in enumerate(X):
        for x2, y2, c2 in X[:i]:
            L += intersectionCC(x1, y1, x2, y2, t/c1, t/c2)
    for x0, y0 in L:
        cnt = 0
        for x, y, c in X:
            if math.sqrt((x-x0) ** 2 + (y-y0) ** 2) * c <= t + 10**-8:
                cnt += 1
        if cnt >= K: return 1
    return 0

N, K = map(int, input().split())
X = []
L0 = []
for _ in range(N):
    x, y, c = map(int, input().split())
    X.append((x, y, c))
    L0.append((x, y))

l, r = 0, 150000
while (r-l) * (1<<20) > max(1, l):
    m = (l+r) / 2
    if chk(m):
        r = m
    else:
        l = m
    
print((l+r)/2)

Solution (Python):
import math
def intersectionCC0(x1, y1, r1, r2):
    r = x1**2 + y1**2
    t = (r + r1**2 - r2**2) / (2*r)
    dd = r1**2/r - t**2
    if - (10**-8) < dd < 0: dd = 0
    if dd < 0: return []
    x, y = t * x1, t * y1
    if dd == 0: return [(x, y)]
    sq = math.sqrt(dd)
    dx, dy = y1 * sq, -x1 * sq
    return [(x+dx, y+dy), (x-dx, y-dy)]

def intersectionCC(x1, y1, x2, y2, r1, r2):
    return [(x1+x, y1+y) for x, y in intersectionCC0(x2-x1, y2-y1, r1, r2)]

def chk(t):
    L = [a for a in L0]
    for i, (x1, y1, c1) in enumerate(X):
        for x2, y2, c2 in X[:i]:
            L += intersectionCC(x1, y1, x2, y2, t/c1, t/c2)
    for x0, y0 in L:
        cnt = 0
        for x, y, c in X:
            if math.sqrt((x-x0) ** 2 + (y-y0) ** 2) * c <= t + 10**-8:
                cnt += 1
        if cnt >= K: return 1
    return 0

N, K = map(int, input().split())
X = []
L0 = []
for _ in range(N):
    x, y, c = map(int, input().split())
    X.append((x, y, c))
    L0.append((x, y))

l, r = 0, 150000
while (r-l) * (1<<20) > max(1, l):
    m = (l+r) / 2
    if chk(m):
        r = m
    else:
        l = m
    
print((l+r)/2)

Solution (Python):
import math
def intersectionCC0(x1, y1, r1, r2):
    r = x1**2 + y1**2
    t = (r + r1**2 - r2**2) / (2*r)
    dd = r1**2/r - t**2
    if - (10**-16) < dd < 0: dd = 0
    if dd < 0: return []
    x, y = t * x1, t * y1
    if dd == 0: return [(x, y)]
    sq = math.sqrt(dd)
    dx, dy = y1 * sq, -x1 * sq
    return [(x+dx, y+dy), (x-dx, y-dy)]

def intersectionCC(x1, y1, x2, y2, r1, r2):
    return [(x1+x, y1+y) for x, y in intersectionCC0(x2-x1, y2-y1, r1, r2)]

def chk(t):
    L = [a for a in L0]
    for i, (x1, y1, c1) in enumerate(X):
        for x2, y2, c2 in X[:i]:
            L += intersectionCC(x1, y1, x2, y2, t/c1, t/c2)
    for x0, y0 in L:
        cnt = 0
        for x, y, c in X:
            if math.sqrt((x-x0) ** 2 + (y-y0) ** 2) * c <= t + 10**-8:
                cnt += 1
        if cnt >= K: return 1
    return 0

N, K = map(int, input().split())
X = []
L0 = []
for _ in range(N):
    x, y, c = map(int, input().split())
    X.append((x, y, c))
    L0.append((x, y))

l, r = 0, 1 << 30
while (r-l) * (1<<25) > max(1, l):
    m = (l+r) / 2
    if chk(m):
        r = m
    else:
        l = m
    
print(l)

Solution (Python):


def Line_Line(p0x, p0y, p1x, p1y, p2x, p2y, p3x, p3y):
    retx = ((p0y-p2y)*(p0x-p1x)*(p2x-p3x) - p0x*(p2x-p3x)*(p0y-p1y) + p2x*(p2y-p3y)*(p0x-p1x)) \
        / ((p2y-p3y)*(p0x-p1x)-(p2x-p3x)*(p0y-p1y))
    if p0x != p1x:
        rety = (p0y-p1y)*(retx-p0x)/(p0x-p1x) + p0y
    else:
        rety = (p2y-p3y)*(retx-p2x)/(p2x-p3x) + p2y
    return retx, rety



def Circle_Line(p0x, p0y, p1x, p1y, cx, cy, cr):
    tan_p = (-p1y+p0y, p1x-p0x)
    nearestx, nearesty = Line_Line(p0x, p0y, p1x, p1y, cx, cy, cx+tan_p[0], cy+tan_p[1])
    distance = ((nearestx-cx)**2+(nearesty-cy)**2)**0.5
    p_e = ((p1x-p0x)/(tan_p[0]**2+tan_p[1]**2)**0.5, (p1y-p0y)/(tan_p[0]**2+tan_p[1]**2)**0.5)
    near_to_ret = (cr**2-distance**2)**0.5
    return [
        (nearestx + near_to_ret*p_e[0], nearesty + near_to_ret*p_e[1]),
        (nearestx - near_to_ret*p_e[0], nearesty - near_to_ret*p_e[1]),        
        ]



def Circle_Circle(c0x, c0y, c0r, c1x, c1y, c1r):
    d = ((c0x-c1x)**2+(c0y-c1y)**2)**0.5
    if not abs(c0r-c1r) <= d <= c0r+c1r:
        return []
    mid_p = (
        ((c0r**2-c1r**2+d**2)*c1x+(c1r**2-c0r**2+d**2)*c0x)/(2*d**2),
        ((c0r**2-c1r**2+d**2)*c1y+(c1r**2-c0r**2+d**2)*c0y)/(2*d**2)
    )
    vx = c0y - c1y
    vy = c1x - c0x
    return Circle_Line(mid_p[0], mid_p[1], mid_p[0]+vx, mid_p[1]+vy, c0x, c0y, c0r)


import sys
input = sys.stdin.readline
from math import sqrt
eps = 10**(-9)

N, K = map(int, input().split())
XYC = [list(map(int, input().split())) for _ in range(N)]

l = 0
r = 10**8
while r-l > 0.00000001:
    D = (l+r)/2
    Points = []
    for i, (x0, y0, c0) in enumerate(XYC):
        Points.append((x0, y0))
        for j, (x1, y1, c1) in enumerate(XYC):
            if i >= j: continue
            for px, py in Circle_Circle(x0, y0, D/c0, x1, y1, D/c1):
                Points.append((px, py))

    ok = False
    for px, py in Points:
        count = 0
        for x, y, c in XYC:
            if c*((x-px)**2+(y-py)**2)**0.5 <= D+eps:
                count += 1
        if count >= K:
            ok = True
            break
        
    if ok:
        r = D
    else:
        l = D

print(r)

Solution (Python):
N,K = map(int,input().split())
XYC = [list(map(int,input().split())) for _ in [0]*N]
eps = 1e-7
def f(x,y,t):
    ret = 0
    for X,Y,C in XYC:
        if ((x-X)**2+(y-Y)**2)<=(t/C)**2+eps: ret += 1
    return ret

def get_cross_section(x1,y1,r1,x2,y2,r2):
    d = ((x1-x2)**2+(y1-y2)**2)**0.5
    X = (d*d+r1*r1-r2*r2)/(2*d)
    Y = max(0,(r1*r1-X*X))**0.5
    dx = (x2-x1)/d
    dy = (y2-y1)/d
    ret = [(x1+X*dx+Y*dy,y1+X*dy-Y*dx)]
    if Y>0 : ret.append((x1+X*dx-Y*dy,y1+X*dy+Y*dx))
    return ret
def check(t):
    kouho = []
    for i in range(N):
        x1,y1,c1 = XYC[i]
        kouho.append((x1,y1))
        for j in range(i+1,N):
            x2,y2,c2 = XYC[j]
            tmp = get_cross_section(x1,y1,t/c1,x2,y2,t/c2)
            for tmp2 in tmp:
                kouho.append(tmp2)
    for tmp in kouho:
        if f(tmp[0],tmp[1],t) >= K:return True
    return False

l=0
r=1000000

while abs(r-l)>1e-7:
    i = (l+r)/2
    if check(i):
        r = i
    else:
        l = i
print(r)

Solution (Python):
N,K = map(int,input().split())
XYC = [list(map(int,input().split())) for _ in [0]*N]
eps = 1e-7
def f(x,y,t):
    ret = 0
    for X,Y,C in XYC:
        if ((x-X)**2+(y-Y)**2)<=(t/C)**2+eps: ret += 1
    return ret

def get_cross_section(x1,y1,r1,x2,y2,r2):
    d = ((x1-x2)**2+(y1-y2)**2)**0.5
    X = (d*d+r1*r1-r2*r2)/(2*d)
    Y = max(0,(r1*r1-X*X))**0.5
    dx = (x2-x1)/d
    dy = (y2-y1)/d
    ret = [(x1+X*dx+Y*dy,y1+X*dy-Y*dx)]
    if Y>0 : ret.append((x1+X*dx-Y*dy,y1+X*dy+Y*dx))
    return ret
def check(t):
    kouho = []
    for i in range(N):
        x1,y1,c1 = XYC[i]
        kouho.append((x1,y1))
        for j in range(i+1,N):
            x2,y2,c2 = XYC[j]
            tmp = get_cross_section(x1,y1,t/c1,x2,y2,t/c2)
            for tmp2 in tmp:
                kouho.append(tmp2)
    for tmp in kouho:
        if f(tmp[0],tmp[1],t) >= K:return True
    return False

l=0
r=1000000
for _ in range(100):
    i = (l+r)/2
    if check(i):
        r = i
    else:
        l = i
print(r)

Solution (Python):

from collections import defaultdict,deque
from heapq import heappush, heappop
from itertools import permutations
import sys
import math
import bisect
def LI(): return [int(x) for x in sys.stdin.readline().split()]
def I(): return int(sys.stdin.readline())
def LS():return [list(x) for x in sys.stdin.readline().split()]
def S():
    res = list(sys.stdin.readline())
    if res[-1] == "\n":
        return res[:-1]
    return res
def IR(n):
    return [I() for i in range(n)]
def LIR(n):
    return [LI() for i in range(n)]
def SR(n):
    return [S() for i in range(n)]
def LSR(n):
    return [LS() for i in range(n)]

sys.setrecursionlimit(1000000)
mod = 1000000007

def solve():
    
    def D2(x1,y1,x2,y2):
        return (x1-x2)**2+(y1-y2)**2

    
    def check(R1,R2):
        x1,y1,r1 = R1
        x2,y2,r2 = R2
        d2 = D2(x1,y1,x2,y2)
        if (r1-r2)**2 < d2 <= (r1+r2)**2:
            return 1
        else:
            return 0

    
    def point(R1,R2):
        if not check(R1,R2):
            return []
        else:
            x1,y1,r1 = R1
            x2,y2,r2 = R2
            x2 -= x1
            y2 -= y1
            r = x2**2+y2**2
            a = (r+r1**2-r2**2)/2
            d = (r*r1**2-a**2)**0.5
            X,Y = a*x2/r, a*y2/r
            k = d/r
            X1,Y1 = X+y2*k+x1, Y-x2*k+y1
            X2,Y2 = X-y2*k+x1, Y+x2*k+y1
            return [(X1,Y1),(X2,Y2)]

    
    def inR(p,R):
        x2,y2 = p
        x1,y1,r1 = R
        d2 = D2(x1,y1,x2,y2)
        if d2 <= r1**2+1e-3:
            return 1
        else:
            return 0

    n,k = LI()
    p = LIR(n)
    l = 0
    r = 1000000
    while r-l > 1e-6:
        t = (l+r) / 2
        P = []
        for i in range(n):
            xi,yi,ci = p[i]
            ri = t/ci
            P.append((xi,yi))
            for j in range(i):
                xj,yj,cj = p[j]
                rj = t/cj
                ps = point((xi,yi,ri),(xj,yj,rj))
                for pi in ps:
                    P.append(pi)
        for pi in P:
            s = 0
            for x,y,c in p:
                r1 = t/c
                if inR(pi,(x,y,r1)):
                    s += 1
                if s >= k:
                    break
            if s >= k:
                r = t
                break
        else:
            l = t
    print(l)
    return


if __name__ == "__main__":
    solve()

Solution (Python):
import sys
from itertools import combinations
from math import atan2, cos, sin


def apollonius(p1, c1, p2, c2):
    m1 = p1 + (p2 - p1) * c2 / (c1 + c2)
    m2 = p1 + (p2 - p1) * c2 / (c2 - c1)
    m = (m1 + m2) / 2
    r = abs(m1 - m)
    return m, r


def apollonius_intersections(p1, c1, p2, c2, p3, c3):
    m12, r12 = apollonius(p1, c1, p2, c2)
    m13, r13 = apollonius(p1, c1, p3, c3)
    v = m13 - m12
    d = abs(v)
    if d > r12 + r13 or d < abs(r12 - r13):
        return None, None
    theta = atan2(v.imag, v.real)
    xx = (r12 ** 2 - r13 ** 2 + d ** 2) / (2 * d)
    s = (r12 + r13 + d) / 2
    yy = 2 * (s * (s - r12) * (s - r13) * (s - d)) ** 0.5 / d
    st = sin(theta)
    ct = cos(theta)
    e1 = (xx * ct - yy * st) + (xx * st + yy * ct) * 1j + m12
    e2 = (xx * ct + yy * st) + (xx * st - yy * ct) * 1j + m12

    t1 = abs(e1 - p1) * c1
    t2 = abs(e2 - p1) * c1
    if t1 < t2:
        return e1, t1
    return e2, t2


def get_circumscribed_center(p1, p2, p3):
    A = abs(p2 - p3) ** 2
    B = abs(p3 - p1) ** 2
    C = abs(p1 - p2) ** 2
    T = A * (B + C - A)
    U = B * (C + A - B)
    W = C * (A + B - C)
    
    try:
        return (T * p1 + U * p2 + W * p3) / (T + U + W)
    except:
        return 0j


def solve(n, k, xyc):
    if k == 1:
        return 0

    xxx = xyc[0::3]
    yyy = xyc[1::3]
    ccc = xyc[2::3]
    beefs = [(i, x + y * 1j, c) for i, (x, y, c) in enumerate(zip(xxx, yyy, ccc))]

    ans = 1e18
    for (i1, p1, c1), (i2, p2, c2) in combinations(beefs, 2):
        
        
        px = p1 + (p2 - p1) * (c2 / (c1 + c2))
        t = abs(px - p1) * c1 + 1e-8

        ok = 0
        ng = []
        for i, p, c in beefs:
            if abs(px - p) * c <= t:
                ok += 1
            else:
                ng.append(i)
        if ok >= k:
            ans = min(ans, t)
            continue

        

        
        
        
        
        
        
        for i3 in ng:
            if i3 < i2:
                continue
            _, p3, c3 = beefs[i3]
            if c1 == c2:
                if c2 == c3:
                    e = get_circumscribed_center(p1, p2, p3)
                    t = abs(e - p1) * c1
                else:
                    e, t = apollonius_intersections(p3, c3, p1, c1, p2, c2)
            elif c1 == c3:
                e, t = apollonius_intersections(p2, c2, p1, c1, p3, c3)
            else:
                e, t = apollonius_intersections(p1, c1, p2, c2, p3, c3)

            if e is None:
                continue

            t += 1e-8
            ok2 = 0
            for i, p, c in beefs:
                if abs(e - p) * c <= t:
                    ok2 += 1
            if ok2 >= k:
                ans = min(ans, t)

    return ans


n, k = map(int, input().split())
xyc = list(map(int, sys.stdin.read().split()))
print(solve(n, k, xyc))

Solution (Python):
import sys

stdin = sys.stdin

ns = lambda: stdin.readline().rstrip()
ni = lambda: int(stdin.readline().rstrip())
nm = lambda: map(int, stdin.readline().split())
nl = lambda: list(map(int, stdin.readline().split()))

def intersection(p, q, eps = 10**-7):
    x1, y1, r1 = p
    x2, y2, r2 = q
    X = x2 - x1; Y = y2 - y1
    d = X**2 + Y**2
    c1 = r1**2; c2 = r2**2

    rc = (d + c1 - c2) / 2
    rs = d * c1 - rc ** 2
    
    if rs < eps: return list()
    
    rs **= .5
    p1x = x1 + (rc * X + rs * Y) / d
    p1y = y1 + (rc * Y - rs * X) / d
    p2x = x1 + (rc * X + rs * Y) / d
    p2y = y1 + (rc * Y + rs * X) / d
    return [(p1x, p1y), (p2x, p2y)]

eps = 10**-7
n,k = nm()
if k == 1: print(0); exit(0)
niku = [tuple(nm()) for _ in range(n)]
l = [c*(x**2+y**2)**.5 for x,y,c in niku]
l.sort()
hi = l[k-1]
lo = 0
for _ in range(50):
  mid = (hi+lo)/2
  l = [(x,y,mid/c) for x,y,c in niku]
  p = [(x,y) for x,y,c in niku]
  for i in range(n-1):
    for j in range(i+1,n):
      p.extend(intersection(l[i], l[j]))
  for x,y in p:
    f = sum((v-x)**2 + (w-y)**2 <= r**2 + eps for v,w,r in l)
    if f >= k:
      hi = mid
      break
  else:
    lo = mid
print(hi)

Solution (Python):
import math

n, k = map(int, input().split())
x, y, c = zip(*[map(int, input().split()) for _ in range(n)])
que = [(0, 0, 1000)]
res = 1e9
i = 0
while i < len(que):
    px, py, r = que[i]
    i += 1
    lb = [c[j] * math.hypot(max(abs(px - x[j]) - r, 0), max(abs(py - y[j]) - r, 0)) for j in range(n)]
    if res * (1 - 5e-7) > sorted(lb)[k - 1]:
        dis = [c[j] * math.hypot(px - x[j], py - y[j]) for j in range(n)]
        res = min(res, sorted(dis)[k - 1])
        for dx in (-0.5, 0.5):
            for dy in (-0.5, 0.5):
                que.append((px + dx * r, py + dy * r, r / 2))
print(res)

Solution (Python):
import cmath,sys
input = sys.stdin.readline


"""
cdist(z1,z2) : z1とz2の距離
cInDiv(z1,z2) : z1とz2をp:qに内分する点
cOutDiv(z1,z2) : z1とz2をp:qに外分する点
"""


EPS = 0.000000001

def cdist(z1,z2=0):
    return abs(z1-z2)

def cInDiv(z1,z2,p,q):
    return (q*z1+p*z2)/(p+q)

def cOutDiv(z1,z2,p,q):
    return (-q*z1+p*z2)/(p-q)

"""
cInCircle(z,c,r) : 点zは中心c半径rの円の中にあるか(0:円周上,-1:円内,1:円外)
cIntersectionPnt(z1,z2,r1,r2) : 中心z1半径r1,中心z2半径r2の円の交点の数、座標1,座標2
"""

def cInCircle(z,c,r):
    d = abs(z-c)
    if abs(r-d) < EPS:
        return 0
    if r-d > 0:
        return -1
    return 1

def cIntersectionPnt(z1,z2,r1,r2):
    d = abs(z1-z2)

    if d < EPS:
        if abs(r1-r2) < EPS:
            return float("inf"),-1,-1
        else:
            return 0,-1,-1
    
    
    if abs(d - (r1+r2)) < EPS:
        diff = (z2-z1)/d
        return 1,z1+diff*r1,-1
    
    if abs(d - (r1-r2)) < EPS:
        diff = (z2-z1)/d
        return 1,z1+diff*r1,-1
    
    if abs(d - (r2-r1)) < EPS:
        diff = (z1-z2)/d
        return 1,z2+diff*r2,-1

    if d > r1+r2:
        return 0,-1,-1
    if d < abs(r1-r2):
        return 0,-1,-1


    rc = (r1*r1+d*d-r2*r2)/(2*d)
    rs = cmath.sqrt(r1*r1-rc*rc)
    diff = (z2-z1)/d
    return 2,z1+diff*rc+diff*1j*rs,z1+diff*rc+diff*(-1j)*rs


n,k = map(int,input().split())
xy = []
c = []
for i in range(n):
    xx,yy,cc = map(int,input().split())
    xy.append(xx+yy*1j)
    c.append(cc)


def judge(t):
    points = []
    for i in range(n):
        for j in range(i+1,n):
            flag,pt1,pt2 = cIntersectionPnt(xy[i],xy[j],t/c[i],t/c[j])
            if flag == 1:
                points.append(pt1)
            elif flag == 2:
                points.append(pt1)
                points.append(pt2)
    for pc in xy:
        points.append(pc)
    for z in points:
        cnt = 0
        for i in range(n):
            flag = cInCircle(z,xy[i],t/c[i])
            if flag <= 0:
                cnt += 1
        if cnt >= k:
            return True
    return False


ok = 10**10
ng = 0
while abs(ok-ng) > EPS:
    mid = (ok+ng)/2
    if judge(mid):
        ok = mid
    else:
        ng = mid

print(ok)

Solution (Python):
import sys
from itertools import combinations
from math import sqrt
import numpy as np

read = sys.stdin.read
readline = sys.stdin.readline


def main():
    N, K, *xyc = map(int, read().split())
    x = xyc[::3]
    y = xyc[1::3]
    xy = list(zip(x, y))
    c = xyc[2::3]
    left = 0
    right = 10 ** 6
    delta = 0.1 ** 7
    while right - left > delta:
        mid = (left + right) / 2
        intersections = []
        r = [(mid / i) ** 2 for i in c]
        for i, j in combinations(range(N), 2):
            x1, y1 = xy[i]
            x2, y2 = xy[j]
            X = x2 - x1
            Y = y2 - y1
            XY = X ** 2 + Y ** 2
            r1 = r[i]
            r2 = r[j]

            a = (XY + r1 - r2) / 2
            b = XY * r1 - a ** 2
            if b < delta:
                continue
            b = sqrt(b)
            xi1 = (a * X + Y * b) / XY + x1
            yi1 = (a * Y - X * b) / XY + y1
            xi2 = (a * X - Y * b) / XY + x1
            yi2 = (a * Y + X * b) / XY + y1
            intersections.append((xi1, yi1))
            intersections.append((xi2, yi2))

        intersections.extend(xy)
        intersections = np.array(intersections, np.float)
        cnt = np.zeros(len(intersections), np.int64)

        for (i, j), k in zip(xy, r):
            XY = (intersections[:, 0] - i) ** 2 + (intersections[:, 1] - j) ** 2
            cnt[XY < k + delta] += 1

        if np.any(cnt >= K):
            right = mid
        else:
            left = mid

    print(right)


if __name__ == '__main__':
    main()

Solution (Python):
import sys

sys.setrecursionlimit(10 ** 6)
int1 = lambda x: int(x) - 1
p2D = lambda x: print(*x, sep="\n")
def II(): return int(sys.stdin.readline())
def MI(): return map(int, sys.stdin.readline().split())
def LI(): return list(map(int, sys.stdin.readline().split()))
def LLI(rows_number): return [LI() for _ in range(rows_number)]
def SI(): return sys.stdin.readline()[:-1]

class Vector:
    def __init__(self, x, y):
        self.x, self.y = x, y
        self.sqrnorm = pow(self.x, 2) + pow(self.y, 2)
        self.norm = self.sqrnorm ** 0.5

    def __add__(self, other):
        return Vector(self.x + other.x, self.y + other.y)

    def __sub__(self, other):
        return Vector(self.x - other.x, self.y - other.y)

    def __mul__(self, other):
        return Vector(self.x * other, self.y * other)

    __rmul__ = __mul__

class Circle:
    def __init__(self, ox, oy, r):
        self.o, self.r = Vector(ox, oy), r

def rot90(v): return Vector(-v.y, v.x)

def comv(c1, c2):
    if c1.r < c2.r: c1, c2 = c2, c1
    r1, o1 = c1.r, c1.o
    r2, o2 = c2.r, c2.o
    d = (o1 - o2).norm
    s, mx = r1 + r2 + d, max(r1, r2, d)
    if s + eps < 2 * mx: return []
    elif s > 2 * mx:
        a = (r1 * r1 - r2 * r2 + d * d) / d / 2
        h = (r1 * r1 - a * a) ** 0.5
        o2o1 = o2 - o1
        v1 = o1 + o2o1 * (a / d) + rot90(o2o1) * (h / d)
        v2 = o1 + o2o1 * (a / d) - rot90(o2o1) * (h / d)
        return [v1, v2]
    else: return [(o2 - o1) * (r1 / d)]

def isinside(c, v): return (c.o - v).norm <= c.r + eps

eps = 10 ** -7
def main():
    def ok(t):
        
        cc = []
        vv = []
        for x, y, c in xyc:
            cc.append(Circle(x, y, t / c))
            vv.append(Vector(x, y))
        
        for j in range(n):
            for i in range(j):
                vv += comv(cc[i], cc[j])
        for v in vv:
            cnt=0
            for c in cc:
                if isinside(c,v):cnt+=1
                if cnt==k:return True
        return False

    n, k = MI()
    xyc = LLI(n)
    l, r = 0, 200000
    for _ in range(50):
        t = (l + r) / 2
        if ok(t): r = t
        else: l = t
    print(t)

main()