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

from math import acos
from cmath import exp
from itertools import combinations
N, K = map(int, input().split())
XYC = []
for _ in range(N):
    x, y, c = map(int, input().split())
    XYC.append((complex(x, y), c))
if K==1:
    print(0)
    exit()
ans = 1e300
eps = 1e-7
for (p1, c1), (p2, c2) in combinations(XYC, 2):
    p_center = (p1 * c1 + p2 * c2) / (c1 + c2)
    t = abs(p_center-p1) * c1
    k = 0
    for xy, c in XYC:
        if abs(p_center-xy) * c < t + eps:
            k += 1
    if k >= K and t < ans:
        ans = t

def get_center(p1, p2, p3, c1, c2, c3):
    if c1 == c2 == c3:
        t1, t2, t3 = p2-p3, p3-p1, p1-p2
        denom = (t1*p1.conjugate()+t2*p2.conjugate()+t3*p3.conjugate())
        if denom == 0:
            return []
        return [(t1*abs(p1)**2+t2*abs(p2)**2+t3*abs(p3)**2) / denom]
    elif c1 == c3:
        c2, c3 = c3, c2
        p2, p3 = p3, p2
    elif c2 == c3:
        c1, c3 = c3, c1
        p1, p3 = p3, p1
    o1 = (c1*c1*p1-c3*c3*p3)/(c1*c1-c3*c3)
    r1 = (c1*c3*abs(p1-p3))/abs(c1*c1-c3*c3)
    o2 = (c2*c2*p2-c3*c3*p3)/(c2*c2-c3*c3)
    r2 = (c2*c3*abs(p2-p3))/abs(c2*c2-c3*c3)
    l = abs(o1-o2)
    if r1+r2 < l or l+r1<r2 or l+r2<r1:
        return []
    cos_a = (r1*r1+l*l-r2*r2)/(2*r1*l)
    a = acos(cos_a)
    res1 = o1 + (o2-o1) * exp(a*1j) * r1 / l
    res2 = o1 + (o2-o1) * exp(-a*1j) * r1 / l
    return [res1, res2]

for (p1, c1), (p2, c2), (p3, c3) in combinations(XYC, 3):
    p_centers = get_center(p1, p2, p3, c1, c2, c3)
    for p_center in p_centers:
        t = abs(p1 - p_center) * c1
        
        
        
        k = 0
        for p, c in XYC:
            if abs(p_center - p) * c < t + eps:
                k += 1
        if k >= K and t < ans:
            ans = t

print(ans)

Solution (Python):
import numpy as np
from scipy.optimize import fmin
from sys import stdin as s
n,k=map(int,s.buffer.readline().split())
Z=np.array(s.buffer.read().split(), np.float64)
X,Y,C=Z[::3],Z[1::3],Z[2::3]
Z=[Z[i*3:i*3+3]for i in range(n)]
def f(x):
  x,y=x
  t=np.sort(C*np.sqrt((X-x)*(X-x)+(Y-y)*(Y-y)))
  return t[k-1]
A=[fmin(f,(x,y),disp=False,ftol=1e-7,initial_simplex=[(x+0.5,y),(x-0.5,y+0.5),(x-0.5,y-0.5)])for x,y,c in Z]
B=list(map(lambda x:f(x),A))
print(min(B))

Solution (Python):
import numpy as np
from scipy.optimize import fmin
import sys
from sys import stdin as s
n,k=map(int,s.buffer.readline().split())
Z=np.array(s.buffer.read().split(), np.float64)
X,Y,C=Z[::3],Z[1::3],Z[2::3]
Z=[Z[i*3:i*3+3]for i in range(n)]
def f(x):
  x,y=x
  t=np.sort(C*np.sqrt((X-x)*(X-x)+(Y-y)*(Y-y)))
  return t[k-1]
A=[fmin(f,(x,y),disp=False,ftol=1e-7,initial_simplex=[(x+0.5,y),(x-0.5,y+0.5),(x-0.5,y-0.5)])for x,y,c in Z]
B=list(map(lambda x:f(x),A))
print(min(B))

Solution (Python):
import numpy as np
from scipy.optimize import fmin

n,k=map(int,input().split())
XYC=[list(map(int,input().split())) for i in range(n)]
X,Y,C=[],[],[]
for x,y,c in XYC:
    X.append(x)
    Y.append(y)
    C.append(c)
X=np.array(X,dtype=np.float)
Y=np.array(Y,dtype=np.float)
C=np.array(C,dtype=np.float)
def f(xy):
    x,y=xy
    times=(X-x)*(X-x)+(Y-y)*(Y-y)
    times=C*np.sqrt(times)
    times=np.sort(times)
    return times[k-1]
deff=0.5
ansx=[fmin(f,(x,y),disp=False,ftol=10**(-7),initial_simplex=[(x+deff,y),(x-deff,y+deff),(x-deff,y-deff)]) for x,y,c in XYC]
ansy=list(map(lambda x:f(x),ansx))
print(min(ansy))

Solution (Python):
from math import sqrt, fabs


class Circle:
    def __init__(self, x, y, r):
        self.x = x
        self.y = y
        self.r = r
        self.eps = 1e-8

    def contain(self, x, y):
        d2 = (self.x - x) ** 2 + (self.y - y) ** 2
        return d2 <= self.r ** 2 + self.eps

    def intersect(self, other):
        eps = self.eps

        if self.r > other.r:
            x1, y1, R = self.x, self.y, self.r
            x2, y2, r = other.x, other.y, other.r
        else:
            x1, y1, R = other.x, other.y, other.r
            x2, y2, r = self.x, self.y, self.r
        d = sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)

        
        if fabs(d) < eps and fabs(R - r) < eps:
            return '相同', None

        
        if fabs(d - (R - r)) < eps:
            x3 = (R * x2 - r * x1) / (R - r)
            y3 = (R * y2 - r * y1) / (R - r)
            return '內切', [(x3, y3)]

        
        if fabs(d - (R + r)) < eps:
            x3 = (R * x1 + r * x2) / (R + r)
            y3 = (R * y1 + r * y2) / (R + r)
            return '外切', [(x3, y3)]

        
        if d < R - r - eps:
            return '內離', None

        
        if d > R + r + eps:
            return '外離', None

        
        a = (R ** 2 - r ** 2 + d ** 2) / (2 * d)
        h = sqrt(R ** 2 - a ** 2)
        cx = x1 + a * (x2 - x1) / d
        cy = y1 + a * (y2 - y1) / d
        x3 = cx + h * (y2 - y1) / d
        y3 = cy - h * (x2 - x1) / d
        x4 = cx - h * (y2 - y1) / d
        y4 = cy + h * (x2 - x1) / d
        return '相交', [(x3, y3), (x4, y4)]


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


def check(t):
    circles = [Circle(x, y, t / c) for x, y, c in zip(X, Y, C)]

    cands = []
    for i in range(N):
        cands.append((X[i], Y[i]))
        for j in range(i + 1, N):
            kind, points = circles[i].intersect(circles[j])
            if kind in ['相交', '內切', '外切']:
                cands.extend(points)

    for (x, y) in cands:
        cnt = sum(1 for circle in circles if circle.contain(x, y))
        if cnt >= K:
            return True

    return False



lb, ub = 0, 1e9
for _ in range(100):
    mid = (lb + ub) / 2
    if check(mid):
        ub = mid
    else:
        lb = mid
print(ub)

Solution (Python):
from math import sqrt, fabs


class Circle:
    def __init__(self, x, y, r):
        self.x = x
        self.y = y
        self.r = r
        self.eps = 1e-8

    def contain(self, x, y):
        d2 = (self.x - x) ** 2 + (self.y - y) ** 2
        return d2 <= self.r ** 2 + self.eps

    def intersect(self, other):
        eps = self.eps

        if self.r > other.r:
            x1, y1, R = self.x, self.y, self.r
            x2, y2, r = other.x, other.y, other.r
        else:
            x1, y1, R = other.x, other.y, other.r
            x2, y2, r = self.x, self.y, self.r
        d = sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)

        
        if fabs(d) < eps and fabs(R - r) < eps:
            return '相同', None

        
        if fabs(d - (R - r)) < eps:
            x3 = (R * x2 - r * x1) / (R - r)
            y3 = (R * y2 - r * y1) / (R - r)
            return '內切', [(x3, y3)]

        
        if fabs(d - (R + r)) < eps:
            x3 = (R * x1 + r * x2) / (R + r)
            y3 = (R * y1 + r * y2) / (R + r)
            return '外切', [(x3, y3)]

        
        if d < R - r - eps:
            return '內離', None

        
        if d > R + r + eps:
            return '外離', None

        
        a = (R ** 2 - r ** 2 + d ** 2) / (2 * d)
        h = sqrt(R ** 2 - a ** 2)
        cx = x1 + a * (x2 - x1) / d
        cy = y1 + a * (y2 - y1) / d
        x3 = cx + h * (y2 - y1) / d
        y3 = cy - h * (x2 - x1) / d
        x4 = cx - h * (y2 - y1) / d
        y4 = cy + h * (x2 - x1) / d
        return '相交', [(x3, y3), (x4, y4)]


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


def check(t):
    circles = [Circle(x, y, t / c) for x, y, c in zip(X, Y, C)]

    cands = []
    for i in range(N):
        cands.append((X[i], Y[i]))
        for j in range(i + 1, N):
            kind, points = circles[i].intersect(circles[j])
            if kind in ['相交', '內切', '外切']:
                cands.extend(points)

    for (x, y) in cands:
        cnt = sum(1 for circle in circles if circle.contain(x, y))
        if cnt >= K:
            return True

    return False



lb, ub = 0, 1e9
for _ in range(100):
    mid = (lb + ub) / 2
    if check(mid):
        ub = mid
    else:
        lb = mid
print(ub)

Solution (Python):
from math import sqrt, fabs


class Circle:
    def __init__(self, x, y, r):
        self.x = x
        self.y = y
        self.r = r
        self.eps = 1e-8

    def contain(self, x, y):
        d2 = (self.x - x) ** 2 + (self.y - y) ** 2
        return d2 <= self.r ** 2 + self.eps

    def intersect(self, other):
        eps = self.eps

        if self.r > other.r:
            x1, y1, R = self.x, self.y, self.r
            x2, y2, r = other.x, other.y, other.r
        else:
            x1, y1, R = other.x, other.y, other.r
            x2, y2, r = self.x, self.y, self.r
        d = sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)

        
        if fabs(d) < eps and fabs(R - r) < eps:
            return '相同', None

        
        if fabs(d - (R - r)) < eps:
            x3 = (R * x2 - r * x1) / (R - r)
            y3 = (R * y2 - r * y1) / (R - r)
            return '內切', [(x3, y3)]

        
        if fabs(d - (R + r)) < eps:
            x3 = (R * x1 + r * x2) / (R + r)
            y3 = (R * y1 + r * y2) / (R + r)
            return '外切', [(x3, y3)]

        
        if d < R - r + eps:
            return '內離', None

        
        if d > R + r + eps:
            return '外離', None

        
        a = (R ** 2 - r ** 2 + d ** 2) / (2 * d)
        h = sqrt(R ** 2 - a ** 2)
        cx = x1 + a * (x2 - x1) / d
        cy = y1 + a * (y2 - y1) / d
        x3 = cx + h * (y2 - y1) / d
        y3 = cy - h * (x2 - x1) / d
        x4 = cx - h * (y2 - y1) / d
        y4 = cy + h * (x2 - x1) / d
        return '相交', [(x3, y3), (x4, y4)]


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


def check(t):
    circles = [Circle(x, y, t / c) for x, y, c in zip(X, Y, C)]

    cands = []
    for i in range(N):
        cands.append((X[i], Y[i]))
        for j in range(i + 1, N):
            kind, points = circles[i].intersect(circles[j])
            if kind in ['相交', '內切', '外切']:
                cands.extend(points)

    for (x, y) in cands:
        cnt = sum(1 for circle in circles if circle.contain(x, y))
        if cnt >= K:
            return True

    return False



lb, ub = 0, 1e9
for _ in range(100):
    mid = (lb + ub) / 2
    if check(mid):
        ub = mid
    else:
        lb = mid
print(ub)

Solution (Python):
import math
import itertools


def intersections(T, x1, y1, c1, x2, y2, c2):
    
    left = ((c1 * c2)**2) * ((x1 - x2)**2 + (y1 - y2)**2)
    right = (T**2) * ((c1 + c2)**2)
    
    if left > right:
        return list()
    
    if left == right:
        return [((c1*x1 + c2*x2) / (c1 + c2), (c1*y1 + c2*y2) / (c1 + c2))]
    
    r1, r2 = T / c1, T / c2
    p12 = (x2 - x1, y2 - y1)
    coef = (r1**2 + p12[0]**2 + p12[1]**2 - r2**2) / (2 * (p12[0]**2 + p12[1]**2))
    p1h = (coef * p12[0], coef * p12[1])
    coef = (r1**2 - p1h[0]**2 - p1h[1]**2) / (p12[0]**2 + p12[1]**2)
    coef = math.sqrt(max(0, coef))  
    hq1 = (- coef * p12[1], coef * p12[0])
    hq2 = (coef * p12[1], - coef * p12[0])
    return [
        (x1 + p1h[0] + hq1[0], y1 + p1h[1] + hq1[1]),
        (x1 + p1h[0] + hq2[0], y1 + p1h[1] + hq2[1])
    ]


def check(T, K, points):
    N = len(points)
    target_points = [(p[0], p[1]) for p in points]
    for pair in itertools.combinations(range(N), 2):
        p1 = points[pair[0]]
        p2 = points[pair[1]]
        target_points += intersections(T, p1[0], p1[1], p1[2], p2[0], p2[1], p2[2])
    for tp in target_points:
        cnt = 0
        for p in points:
            d2 = (tp[0] - p[0])**2 + (tp[1] - p[1])**2
            if d2 <= (T/p[2])**2 + 10**(-6):  
                cnt += 1
            if cnt >= K:
                return True
    return False


def main():
    eps = 10 ** (-10)
    N, K = list(map(int, input().split()))
    points = [list(map(int, input().split())) for _ in range(N)]
    ng, ok = 0, 3 * (10**5)  
    while ok - ng > eps:
        mid = (ok + ng) / 2
        if check(mid, K, points):
            ok = mid
        else:
            ng = mid
    print(ok)


if __name__ == '__main__':
    main()

Solution (Python):
import math

n, k = map(int, input().split())

xyc = []
for _ in range(n):
    x, y, c = map(int, input().split())
    xyc.append((x, y, c))

if k == 1:
    print(0)
    exit()


def get_intersecs(x0, y0, r0, x1, y1, r1):
    
    

    d=math.sqrt((x1-x0)**2 + (y1-y0)**2)

    
    if d > r0 + r1 :
        return []
    
    if d < abs(r0-r1):
        return []
    
    if d == 0 and r0 == r1:
        return []
    else:
        a=(r0**2-r1**2+d**2)/(2*d)
        if r0 < a:
            return []
        h=math.sqrt(r0**2-a**2)
        x2=x0+a*(x1-x0)/d
        y2=y0+a*(y1-y0)/d
        x3=x2+h*(y1-y0)/d
        y3=y2-h*(x1-x0)/d

        x4=x2-h*(y1-y0)/d
        y4=y2+h*(x1-x0)/d

        return [(x3, y3), (x4, y4)]


def calc(x1, y1, x2, y2, c2):
    return c2 * math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)


def count(x, y, t):
    cnt = 0
    for i in range(n):
        xi, yi, ci = xyc[i]
        if calc(x, y, xi, yi, ci) < t + 0.0000001:
            cnt += 1
            if cnt >= k:
                return cnt
    return cnt


def check(t):
    
    for i in range(n):
        x, y, _ = xyc[i]
        if count(x, y, t) >= k:
            return True
    
    for i in range(n - 1):
        for j in range(i + 1, n):
            x1, y1, c1 = xyc[i]
            x2, y2, c2 = xyc[j]
            intersecs = get_intersecs(x1, y1, t/c1, x2, y2, t/c2)
            for (x, y) in intersecs:
                
                if count(x, y, t) >= k:
                    return True
    return False


low = 0
high = 10 ** 7
for _ in range(100):
    mid = (high + low) / 2
    if check(mid):
        high = mid
    else:
        low = mid
print((high + low) / 2)

Solution (Python):
import sys
from cmath import rect
import numpy as np
import numba
from numba import njit
i8 = numba.int64

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

eps = 1e-9

@njit((numba.complex128, numba.float64) * 2, cache=True)
def cross_circles(c1, r1, c2, r2):
    d = abs(c1 - c2)
    if d <= abs(r1 - r2) or r1 + r2 <= d:
        return None
    cos_theta = (r1 * r1 + d * d - r2 * r2) / (2 * r1 * d)
    theta = np.arccos(cos_theta)
    z1 = (c2 - c1) * rect(r1 / d, theta)
    z2 = (c2 - c1) * rect(r1 / d, -theta)
    return (c1 + z1, c1 + z2)

@njit
def test(Z, C, t):
    R = t / C
    cand = list(Z)
    for j in range(len(Z)):
        for i in range(j):
            points = cross_circles(Z[i], R[i], Z[j], R[j])
            if points is not None:
                for p in points:
                    cand.append(p)
    cand = np.array(cand)
    cnt = np.zeros(len(cand), np.int64)
    for i in range(len(Z)):
        cnt += np.abs(cand - Z[i]) <= R[i] + eps
    return np.any(cnt >= K)

N, K = map(int, readline().split())
nums = np.array(read().split(), np.float64)

Z = nums[::3] + 1j * nums[1::3]
C = nums[2::3]

l = 0
r = 10**9
for _ in range(100):
    m = (l + r) / 2
    if test(Z, C, m):
        r = m
    else:
        l = m

print(r)

Solution (Python):
import math
import sys
from itertools import combinations

stdin = sys.stdin


def ni(): return int(ns())


def na(): return list(map(int, stdin.readline().split()))


def naa(N): return [na() for _ in range(N)]


def ns(): return stdin.readline().rstrip()  


def circumcircle(ax, ay, bx, by, cx, cy):
    if (bx - ax) * (cy - ay) - (by - ay) * (cx - ax) == 0:
        return None
    D = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    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
    R = math.sqrt((X - ax) * (X - ax) + (Y - ay) * (Y - ay))
    return X, Y, R


def diametercircle(ax, ay, bx, by):
    X = (ax + bx) / 2
    Y = (ay + by) / 2
    R = math.sqrt((X - ax) * (X - ax) + (Y - ay) * (Y - ay))
    return X, Y, R


def distance(ax, ay, bx, by):
    return math.sqrt((ax - bx) ** 2 + (ay - by) ** 2)


def intersectioncircle(ax, ay, ra, bx, by, rb):
    d = distance(ax, ay, bx, by)
    if d > ra + rb:
        return None
    if d < abs(ra - rb):
        return None
    if d == 0:
        return None
    t = (ra**2 - rb**2 + d**2) / (2 * d)
    h = math.sqrt(ra**2 - t**2 + sys.float_info.epsilon)
    x2 = ax + t * (bx - ax) / d
    y2 = ay + t * (by - ay) / d
    x3 = x2 + h * (by - ay) / d
    y3 = y2 - h * (bx - ax) / d

    x4 = x2 - h * (by - ay) / d
    y4 = y2 + h * (bx - ax) / d

    return x3, y3, x4, y4


def eat_count(cx, cy, t):
    eat = 0
    for xyc in XYC_array:
        dx, dy, c = xyc
        if distance(cx, cy, dx, dy) <= (t / c + 10 ** -8):
            eat += 1
    return eat


N, K = na()
XYC_array = naa(N)


mint = 0
maxt = 10 ** 6
eps = 10 ** -6

while(maxt - mint > eps):
    t = (mint + maxt) / 2
    ok_flag = 0
    for xyc in XYC_array:
        cx, cy, _ = xyc
        if eat_count(cx, cy, t) >= K:
            ok_flag = 1
            break
    if ok_flag == 1:
        maxt = t
        continue
    for p in combinations(range(N), 2):
        ax, ay, ca = XYC_array[p[0]]
        bx, by, cb = XYC_array[p[1]]
        if intersectioncircle(ax, ay, t / ca, bx, by, t / cb):
            X1, Y1, X2, Y2 = intersectioncircle(ax, ay, t / ca, bx, by, t / cb)
            if eat_count(X1, Y1, t) >= K or eat_count(X2, Y2, t) >= K:
                ok_flag = 1
                break
    if ok_flag == 1:
        maxt = t
    else:
        mint = t

print(t)

Solution (Python):
from math import hypot

def solve():
    N, K = map(int, input().split())
    xs, ys, cs = [], [], []
    for _ in range(N):
        x, y, c = map(int, input().split())
        xs.append(x)
        ys.append(y)
        cs.append(c)

    EPS = 10**-10
    def isLE(a, b):
        return not (a > b+EPS)

    def getIntersectionLCs(a, b, c, xC, yC, r):
        e = a*xC + b*yC + c
        ab2 = a**2 + b**2
        v = r**2 * ab2 - e**2
        if v < 0:
            return []
        elif v == 0:
            xP = xC + -a*e / ab2
            yP = yC + -b*e / ab2
            return [(xP, yP)]
        else:
            v = v**0.5
            ae, be = a*e, b*e
            av, bv = a*v, b*v
            xP1 = xC + (-ae + bv) / ab2
            yP1 = yC + (-be - av) / ab2
            xP2 = xC + (-ae - bv) / ab2
            yP2 = yC + (-be + av) / ab2
            return [(xP1, yP1), (xP2, yP2)]

    def getIntersectionCCs(x1, y1, r1, x2, y2, r2):
        xC, yC = x2-x1, y2-y1
        c = -(xC**2 + yC**2 + r1**2 - r2**2) / 2
        intCCs = getIntersectionLCs(xC, yC, c, 0, 0, r1)
        intCCs = [(x1 + x, y1 + y) for x, y in intCCs]
        return intCCs

    def countNum(xP, yP, xs, ys, rs, K):
        num = 0
        for i, (x, y, r) in enumerate(zip(xs, ys, rs), start=1):
            dist = hypot(x-xP, y-yP)
            if isLE(dist, r):
                num += 1
            if num >= K:
                return True
            elif num+N-i < K:
                return False
        return num >= K

    def isOK(t):
        rs = [t/c for c in cs]

        for xP, yP in zip(xs, ys):
            if countNum(xP, yP, xs, ys, rs, K):
                return True

        for i in range(N):
            x1, y1, r1 = xs[i], ys[i], rs[i]
            for j in range(i+1, N):
                x2, y2, r2 = xs[j], ys[j], rs[j]
                intCCs = getIntersectionCCs(x1, y1, r1, x2, y2, r2)
                for xP, yP in intCCs:
                    if countNum(xP, yP, xs, ys, rs, K):
                        return True

        return False

    ng, ok = 0, 2*10**5
    while abs(ok-ng) > 10**-7:
        mid = (ng+ok) / 2
        if isOK(mid):
            ok = mid
        else:
            ng = mid

    print(ok)


solve()

Solution (Python):

import sys
import cmath
import math
import copy
INF = float("inf")

EPS = 1e-9


def ctc(p1, r1, p2, r2):
    
    
    
    
    
    
    if r2 < r1:
        p1, p2 = p2, p1
        r1, r2 = r2, r1

    d = abs(p1-p2)

    if p1 == p2 and d < EPS:    
        return INF
    elif (d - (r2+r1)) > EPS:   
        return 4
    elif (d - (r2+r1)) > -EPS:  
        return 3
    elif (d-(r2-r1)) > EPS:     
        return 2
    elif (d-(r2-r1)) > -EPS:    
        return 1
    else:                       
        return 0


def inner(p1, c1, p2, c2):
    coeff1 = c1/(c1+c2)
    coeff2 = c2/(c1+c2)
    return p1*coeff1+p2*coeff2


def intersection(p1, r1, p2, r2):
    status = ctc(p1, r1, p2, r2)
    d = abs(p1-p2)
    if not (1 <= status <= 3):
        return []
    elif status == 3:           
        return [inner(p1, r1, p2, r2)]
    elif status == 1:           
        return [p1 + r1*(p1-p2)/d]
    else:                       
        if r1 > r2:
            p1, p2 = p2, p1
            r1, r2 = r2, r1
        d = abs(p1-p2)
        cos = (r1**2+d**2-r2**2)/(2*r1*d)
        angle = math.acos(cos)
        phi = cmath.phase(p2-p1)
        return [
            p1 + cmath.rect(r1, phi+angle),
            p1 + cmath.rect(r1, phi-angle)
        ]


def circle_in(p, o, r1, allow_on_edge=True):
    if allow_on_edge:
        return abs(o-p)-r1 < EPS
    else:
        return abs(o-p)-r1 < -EPS


def solve(N: int, K: int, x: "List[int]", y: "List[int]", c: "List[int]"):

    
    comp = [complex(xx, yy) for xx, yy in zip(x, y)]

    left = 0
    right = 2000*max(c)
    while right - left > EPS:
        mid = (left + right)/2
        yakiniku = 0
        b = copy.copy(comp)
        for i in range(N):
            for j in range(i+1, N):
                b.extend(intersection(comp[i], mid/c[i], comp[j], mid/c[j]))

        for p in b:
            count = [circle_in(p, comp[k], mid/c[k])
                     for k in range(N)]
            yakiniku = max(yakiniku, sum(count))
        if yakiniku >= K:
            right = mid
        else:
            left = mid

    print(mid)
    return


def main():

    def iterate_tokens():
        for line in sys.stdin:
            for word in line.split():
                yield word
    tokens = iterate_tokens()
    N = int(next(tokens))  
    K = int(next(tokens))  
    x = [int()] * (N)  
    y = [int()] * (N)  
    c = [int()] * (N)  
    for i in range(N):
        x[i] = int(next(tokens))
        y[i] = int(next(tokens))
        c[i] = int(next(tokens))
    solve(N, K, x, y, c)


if __name__ == '__main__':
    main()

Solution (Python):
def main():
    import sys
    sys.setrecursionlimit(10**9)
    input = sys.stdin.readline
    from math import hypot
    from collections import deque

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

    dq = deque([(0, 0, 1000.0)])

    ans = 10**18
    while dq:
        X, Y, r = dq.popleft()
        c_hypot = [c * hypot(max(abs(X - x) - r, 0), max(abs(Y - y) - r, 0)) for x, y, c in XYC]
        c_hypot.sort()
        if ans * 0.9999995 > c_hypot[K-1]:
            c_hypot = [c * hypot(X - x, Y - y) for x, y, c in XYC]
            c_hypot.sort()
            if c_hypot[K-1] < ans:
                ans = c_hypot[K-1]

            r /= 2
            dq += [
                (X + r, Y + r, r),
                (X - r, Y + r, r),
                (X + r, Y - r, r),
                (X - r, Y - r, r),
            ]
    print(ans)

main()

Solution (Python):

import sys
INF = float("inf")

tol = 1e-9


def ctc(x1, y1, r1, x2, y2, r2):
    
    
    
    
    
    tol = 1e-9
    
    if r2 < r1:
        x1, x2 = x2, x1
        y1, y2 = y2, y1
        r1, r2 = r2, r1

    d = ((x1-x2)**2 + (y1-y2)**2)**0.5

    if x1 == x2 and y1 == y2 and d < tol:
        
        return INF
    elif (d - (r2+r1)) > tol:
        
        return 4
    elif (d - (r2+r1)) > -tol:
        
        return 3
    elif (d-(r2-r1)) > tol:
        
        return 2
    elif (d-(r2-r1)) > -tol:
        
        return 1
    else:
        
        return 0


def inner(x1, y1, c1, x2, y2, c2):
    coeff1 = c1/(c1+c2)
    coeff2 = c2/(c1+c2)
    b = (x1*coeff1+x2*coeff2, y1*coeff1+y2*coeff2)
    return b


def intersection(x1, y1, r1, x2, y2, r2):
    status = ctc(x1, y1, r1, x2, y2, r2)

    if not (1 <= status <= 3):
        return []
    elif status == 3:
        
        return [inner(x1, y1, r1, x2, y2, r2)]
    elif status == 1:
        
        return [(x1 + r1*(x1-x2)/d, y1 + r1*(y1-y2)/d)]
    else:
        
        L = ((x1-x2)**2 + (y1-y2)**2)**0.5
        cos = ((L**2+r1**2)-r2**2)/(2*r1*L)
        e1 = ((x2-x1)/L, (y2-y1)/L)
        e2 = (-(y2-y1)/L, (x2-x1)/L)
        P1H = ((r1*cos/L)*(x2-x1), (r1*cos/L)*(y2-y1))
        ah = (r1**2-(r1*cos)**2)**0.5
        HA = (e2[0]*ah, e2[1]*ah)
        A = (x1 + P1H[0]+HA[0], y1+P1H[1]+HA[1])
        B = (x1 + P1H[0]-HA[0], y1+P1H[1]-HA[1])
        return [A, B]


def circle_in(x, y, x1, y1, r1):
    tol = 1e-9
    L = (x1-x)**2 + (y1-y)**2
    R = r1**2
    LR = L-R

    if LR > tol:
        return False
    elif LR > -tol:
        return True
    else:
        return True


def solve(N: int, K: int, x: "List[int]", y: "List[int]", c: "List[int]"):

    

    left = 0
    right = 2000*max(c)
    for _ in range(300):
        if right - left < tol:
            break
        mid = (left + right)/2
        yakiniku = 0
        for i in range(N):
            for j in range(i+1, N):
                b = intersection(x[i], y[i], mid/c[i], x[j], y[j], mid/c[j])
                for xx, yy in b:
                    count = [circle_in(xx, yy, x[k], y[k], mid/c[k])
                             for k in range(N)]
                    if yakiniku < sum(count):
                        yakiniku = max(yakiniku, sum(count))
        for i in range(N):
            count = [circle_in(x[i], y[i], x[k], y[k], mid/c[k])
                     for k in range(N)]
            if yakiniku < sum(count):
                yakiniku = max(yakiniku, sum(count))
        
        
        
        
        if yakiniku >= K:
            right = mid
        else:
            left = mid

    print(mid)
    return


def main():

    def iterate_tokens():
        for line in sys.stdin:
            for word in line.split():
                yield word
    tokens = iterate_tokens()
    N = int(next(tokens))  
    K = int(next(tokens))  
    x = [int()] * (N)  
    y = [int()] * (N)  
    c = [int()] * (N)  
    for i in range(N):
        x[i] = int(next(tokens))
        y[i] = int(next(tokens))
        c[i] = int(next(tokens))
    solve(N, K, x, y, c)


if __name__ == '__main__':
    main()

Solution (Python):
import sys
INF = 1 << 60
MOD = 10**9 + 7 
sys.setrecursionlimit(2147483647)
input = lambda:sys.stdin.readline().rstrip()
from cmath import phase, rect
from math import acos
def resolve():
    n, k = map(int, input().split())
    ZC = [None] * n
    for i in range(n):
        x, y, c = map(int, input().split())
        ZC[i] = (complex(x, y), c)

    def check(t):
        
        circle = [(z, t / c) for z, c in ZC]
        candidate = [c for c, r in circle]
        for c1, r1 in circle:
            for c2, r2 in circle:
                if c1 == c2:
                    continue
                d = abs(c2 - c1)
                alpha = phase(c2 - c1)
                s = (r1**2 + d**2 - r2**2) / (2 * r1 * d)
                if abs(s) >= 1:
                    continue
                theta = acos(s)
                candidate.append(c1 + rect(r1, alpha + theta))
                candidate.append(c1 + rect(r1, alpha - theta))

        return any(sum(abs(z - c) <= r + 0.00001 for c, r in circle) >= k for z in candidate)

    l = 0
    r = 1 << 30
    for _ in range(100):
        m = (l + r) / 2
        if check(m):
            r = m
        else:
            l = m

    print(m)
resolve()

Solution (Python):
import sys

input = sys.stdin.buffer.readline

n, k = map(int, input().split())
XY = []
C = []
kouho = []
for i in range(n):
    x, y, c = map(int, input().split())
    XY.append((x, y))
    C.append(c)
    kouho.append(complex(x, y))
ct = 100
eps = 10 ** (-7)


class circle:
    def __init__(self, center, r):
        self.center = center
        self.r = r


def circle_cross(c1, c2):
    r1 = c1.r
    o1 = c1.center
    r2 = c2.r
    o2 = c2.center
    d = abs(o1 - o2)
    if d > r1 + r2:
        return 0, 0, 0
    elif r1 + r2 < d + eps and r1 + r2 > d - eps:
        
        return 1, (r2 * o1 + r1 * o2) / (r1 + r2), 0
    elif d < r1 + r2 and d > abs(r1 - r2):
        
        o21 = o2 - o1
        oh_nagasa = (d ** 2 + r1 ** 2 - r2 ** 2) / (2 * d)
        if oh_nagasa < eps:
            return 1, (r2 * o1 + r1 * o2) / (r1 + r2), 0
        oh = oh_nagasa * (o21 / d)
        housen = (o21 / d) * complex(0, 1)
        suisen = max(0, (r1 ** 2 - oh_nagasa ** 2) ** (1 / 2))
        ans1 = o1 + oh + housen * suisen
        ans2 = o1 + oh - housen * suisen
        return 2, ans1, ans2
    elif abs(r1 - r2) > d - eps and abs(r1 - r2) < d + eps:
        
        ans = (-r2 * o1 + r1 * o2) / (r1 - r2)
        return 1, ans, 0
    else:
        
        return 0, 0, 0


if n == 1:
    print(0)
    exit()
if k == 1:
    print(0)
    exit()

ok = 10 ** 10
ng = 0

while ct > 0:
    mid = (ok + ng) / 2
    kouho_temp = kouho[::]
    for i in range(n - 1):
        a, b = XY[i]
        p = C[i]
        A = circle(complex(a, b), mid / p)
        for j in range(i + 1, n):
            c, d = XY[j]
            q = C[j]
            B = circle(complex(c, d), mid / q)
            t1, t2, t3 = circle_cross(A, B)
            if t1 == 0:
                continue
            if t1 == 1:
                kouho_temp.append(t2)
            if t1 == 2:
                kouho_temp.append(t2)
                kouho_temp.append(t3)
    res = 0
    for z in kouho_temp:
        cnt = 0
        for i in range(n):
            x, y = XY[i]
            p = C[i]
            z1 = complex(x, y)
            if abs(z - z1) < mid / p + eps:
                cnt += 1
        res = max(cnt, res)
    if res >= k:
        ok = mid
    else:
        ng = mid
    ct -= 1


print(mid)