Problem:


Problem H: Don't Burst the Balloon

An open-top box having a square bottom is placed on the floor. You see a number of needles vertically planted on its bottom.

You want to place a largest possible spheric balloon touching the box bottom, interfering with none of the side walls nor the needles.

Java Specific: Submitted Java programs may not use "java.awt.geom.Area". You may use it for your debugging purposes.

Figure H.1 shows an example of a box with needles and the corresponding largest spheric balloon. It corresponds to the first dataset of the sample input below.

Figure H.1. The upper shows an example layout and the lower shows the largest spheric balloon that can be placed.

Input

The input is a sequence of datasets. Each dataset is formatted as follows.

n w
x1 y1 h1
:
xn yn hn

The first line of a dataset contains two positive integers, n and w, separated by a space. n represents the number of needles, and w represents the height of the side walls.

The bottom of the box is a 100 &times; 100 square. The corners of the bottom face are placed at positions (0, 0, 0), (0, 100, 0), (100, 100, 0), and (100, 0, 0).

Each of the n lines following the first line contains three integers, xi, yi, and hi. (xi, yi, 0) and hi represent the base position and the height of the i-th needle. No two needles stand at the same position.

You can assume that 1 ≤ n ≤ 10, 10 ≤ w ≤ 200, 0 < xi < 100, 0 < yi < 100 and 1 ≤ hi ≤ 200. You can ignore the thicknesses of the needles and the walls.

The end of the input is indicated by a line of two zeros. The number of datasets does not exceed 1000.

Output

For each dataset, output a single line containing the maximum radius of a balloon that can touch the bottom of the box without interfering with the side walls or the needles. The output should not contain an error greater than 0.0001.

Sample Input

5 16
70 66 40
38 52 20
40 35 10
70 30 10
20 60 10
1 100
54 75 200
1 10
90 10 1
1 11
54 75 200
3 10
53 60 1
61 38 1
45 48 1
4 10
20 20 10
20 80 10
80 20 10
80 80 10
0 0

Output for the Sample Input

26.00000
39.00000
130.00000
49.49777
85.00000
95.00000



Sample input:
5 16
70 66 40
38 52 20
40 35 10
70 30 10
20 60 10
1 100
54 75 200
1 10
90 10 1
1 11
54 75 200
3 10
53 60 1
61 38 1
45 48 1
4 10
20 20 10
20 80 10
80 20 10
80 80 10
0 0


Sample output:
26.00000
39.00000
130.00000
49.49777
85.00000
95.00000

Submission status: Accepted

Solution code:
from math import sqrt

def getcircle(nee, R): #針にぶつかる範囲の円
    (x, y, h) = nee
    if R <= h:
        return (x, y, R)
    r = sqrt(R**2 - (R-h)**2)
    return (x, y, r)
    
def crosscircle(c1, c2): #円の交点
    (x1, y1, r1) = c1
    (x2, y2, r2) = c2
    if (x1-x2)**2 + (y1-y2)**2 > (r1+r2)**2:
        return ()
    elif x1 == x2 and y1 == y2:
        return ()
    elif (x1-x2)**2 + (y1-y2)**2 == (r1+r2)**2:
        x = (x1*r2+x2*r1)/(r1+r2)
        y = (y1*r2+y2*r1)/(r1+r2)
        return ((x, y),)
    elif y1 == y2:
        x = (r1**2-r2**2-x1**2+x2**2)/(2*(x2-x1))
        if y1**2-(y1**2-r1**2+(x-x1)**2) >= 0:
            yans1 = y1+sqrt(y1**2-(y1**2-r1**2+(x-x1)**2))
            yans2 = y1-sqrt(y1**2-(y1**2-r1**2+(x-x1)**2))
            return ((x, yans1), (x, yans2))
        else:
            return ()
    elif x1 == x2:
        y = (r1**2-r2**2-y1**2+y2**2)/(2*(y2-y1))
        if x1**2-(x1**2-r1**2+(y-y1)**2) >= 0:
            xans1 = x1+sqrt(x1**2-(x1**2-r1**2+(y-y1)**2))
            xans2 = x1-sqrt(x1**2-(x1**2-r1**2+(y-y1)**2))
            return ((xans1, y), (xans2, y))
        else:
            return ()
    else:
        A = ((x1-x2)/(y1-y2))**2 + 1
        B = (x1-x2)*(r1**2-r2**2-x1**2+x2**2+(y1-y2)**2)/(y1-y2)**2 - 2*x1
        C = ((r1**2-r2**2-x1**2+x2**2+(y1-y2)**2)/(2*(y1-y2)))**2 + x1**2 - r1**2
        if B**2 - 4*A*C >= 0:
            xans1 = (-B+sqrt(B**2-4*A*C))/(2*A)
            xans2 = (-B-sqrt(B**2-4*A*C))/(2*A)
            yans1 = -((r1**2-r2**2-x1**2+x2**2-y1**2+y2**2))/(2*(y1-y2))-xans1*(x1-x2)/(y1-y2)
            yans2 = -((r1**2-r2**2-x1**2+x2**2-y1**2+y2**2))/(2*(y1-y2))-xans2*(x1-x2)/(y1-y2)
            return ((xans1, yans1), (xans2, yans2))  
        else:
            return ()
        

while True:
    n, w = map(int, input().split())
    if n == w == 0:
        break
    nees = [tuple(map(float, input().split())) for i in range(n)]
    eps = 0.000001
    U = 130
    L = 0
    while U - L > 0.0001:
        ok = False
        R = (U+L)/2
        cirs = []
        points = []
        if R > w:
            l = sqrt(R**2 - (R-w)**2) #球が壁にぶつからない範囲の正方形(l<=x<=100-l, l<=y<=100-l)
        else:
            l = R
        if l >= 100-l:
            U = R
            continue
        for nee in nees:
            c = getcircle(nee, R)
            cirs.append(c)
        for c in cirs: #正方形と円の交点をリストに入れる
            (x, y, r) = c
            if abs(l-x) <= r:
                if l <= y+sqrt(r**2-(x-l)**2) <= 100-l:
                    points.append((l,y+sqrt(r**2-(x-l)**2)))
                if l <= y-sqrt(r**2-(x-l)**2) <= 100-l:
                    points.append((l,y-sqrt(r**2-(x-l)**2)))
            if abs(l-y) <= r:
                if l <= x+sqrt(r**2-(y-l)**2) <= 100-l:
                    points.append((x+sqrt(r**2-(y-l)**2),y))
                if l <= x-sqrt(r**2-(y-l)**2) <= 100-l:
                    points.append((x-sqrt(r**2-(y-l)**2),y))
            if abs(100-l-x) <= r:
                if l <= y+sqrt(r**2-(x+l-100)**2) <= 100-l:
                    points.append((100-l,y+sqrt(r**2-(x+l-100)**2)))
                if l <= y-sqrt(r**2-(x+l-100)**2) <= 100-l:
                    points.append((100-l,y-sqrt(r**2-(x+l-100)**2)))
            if abs(100-l-y) <= r:
                if l <= x+sqrt(r**2-(100-l-y)**2) <= 100-l:
                    points.append((x+sqrt(r**2-(100-l-y)**2),y))
                if l <= x-sqrt(r**2-(100-l-y)**2) <= 100-l:
                    points.append((x-sqrt(r**2-(100-l-y)**2),y))
        for i in range(n): #円と円の交点をリストに入れる
            if n != 1 and i != n-1:
                for j in range(i+1, n):
                    p = crosscircle(cirs[i], cirs[j])
                    if p == ():
                        continue
                    for a in p:
                        if l <= a[0] <= 100-l and l <= a[1] <= 100-l:
                            points.append(a)
        points.append((l,l)) #正方形の四隅をリストに入れる
        points.append((100-l,l))
        points.append((l,100-l))
        points.append((100-l,100-l))
        for p in points: #リスト内の点が円の内部にないか判定する
            if not l <= p[0] <=100-l:
                continue
            if not l <= p[1] <=100-l:
                continue
            pok = True
            for c in cirs:
                if (p[0]-c[0])**2 + (p[1]-c[1])**2 <= (c[2]-eps)**2:
                    pok = False
                    break
            if pok:
                ok = True
                break
        if ok:
            L = R
        else:
            U = R
    print(R)

Submission status: Accepted

Solution code:
from math import sqrt


def get_cross_point_of_circles(x1, y1, r1, x2, y2, r2):
    rr0 = (x2 - x1) ** 2 + (y2 - y1) ** 2
    xd = x2 - x1
    yd = y2 - y1
    rr1 = r1 ** 2
    rr2 = r2 ** 2
    cv = (rr0 + rr1 - rr2)
    sv2 = 4 * rr0 * rr1 - cv ** 2
    if sv2 < 0:
        return None

    sv = sqrt(sv2)
    return (
        (x1 + (cv * xd - sv * yd) / (2. * rr0), y1 + (cv * yd + sv * xd) / (2. * rr0)),
        (x1 + (cv * xd + sv * yd) / (2. * rr0), y1 + (cv * yd - sv * xd) / (2. * rr0)),
    )


def get_solution_of_quadratic(a, b, c):
    if b**2 >= 4*a*c:
        return {(-b + sqrt(b**2 - 4*a*c))/(2*a), (-b - sqrt(b**2 - 4*a*c))/(2*a)}
    else:
        return {}


def get_nearest_wall(x, y):
    if y <= x:
        if y <= 100 - x:
            return x, 0
        else:
            return 100, y
    else:
        if y <= 100 - x:
            return 0, y
        else:
            return x, 100


while True:
    n, w = map(int, input().split())
    if n == w == 0:
        break

    xyh = [tuple(map(float, input().split())) for _ in range(n)]

    # 目いっぱいでもOKか
    m = (2500+w**2) / (2*w) if w < 50 else 50
    P = (50, 50, m)
    flag0 = True
    for k in range(n):
        if (50 - xyh[k][0]) ** 2 + (50 - xyh[k][1]) ** 2 + max(.0, m - xyh[k][2]) ** 2 <= m ** 2:
            flag0 = False
    if flag0:
        print(m)
        continue

    l = 0
    r = 130
    while abs(r - l) >= 0.0001:
        m = (l + r) / 2

        flag = False

        for i in range(n):
            if n != 1 and i != n - 1:
                # needle二つで抑えられる点
                for j in range(i + 1, n):
                    h1 = min(xyh[i][2], m)
                    h2 = min(xyh[j][2], m)

                    P = get_cross_point_of_circles(xyh[i][0], xyh[i][1], sqrt(2 * h1 * m - h1 ** 2),
                                                   xyh[j][0], xyh[j][1], sqrt(2 * h2 * m - h2 ** 2))
                    if P is None:
                        continue

                    for p in P:
                        if 0 <= p[0] <= 100 and 0 <= p[1] <= 100:

                            wx, wy = get_nearest_wall(p[0], p[1])

                            if (wx - p[0]) ** 2 + (wy - p[1]) ** 2 + max(0, m - w) ** 2 >= m ** 2:
                                flag2 = True
                                for k in range(n):
                                    if k != i and k != j:
                                        if (p[0] - xyh[k][0]) ** 2 + (p[1] - xyh[k][1]) ** 2 \
                                                + max(.0, m - xyh[k][2]) ** 2 <= m ** 2:
                                            flag2 = False
                                if flag2:
                                    flag = True

            # needle一つで抑える
            x, y, h = xyh[i]
            h = min(h, m)

            S1 = get_solution_of_quadratic(2, -2 * (x + y), x ** 2 + y ** 2 - 2 * m * h + h ** 2)
            S2 = get_solution_of_quadratic(2, -2 * (x + 100 - y), x ** 2 + (100 - y) ** 2 - 2 * m * h + h ** 2)
            for s in S1:
                if 0 <= s <= 100:
                    wx, wy = get_nearest_wall(s, s)

                    if (wx - s) ** 2 + (wy - s) ** 2 + max(0, m - w) ** 2 >= m ** 2:
                        flag2 = True
                        for k in range(n):
                            if k != i:
                                if (s - xyh[k][0]) ** 2 + (s - xyh[k][1]) ** 2 \
                                        + max(.0, m - xyh[k][2]) ** 2 <= m ** 2:
                                    flag2 = False
                        if flag2:
                            flag = True
            for s in S2:
                if 0 <= s <= 100:
                    wx, wy = get_nearest_wall(s, 100 - s)

                    if (wx - s) ** 2 + (wy - (100 - s)) ** 2 + max(0, m - w) ** 2 >= m ** 2:
                        flag2 = True
                        for k in range(n):
                            if k != i:
                                if (s - xyh[k][0]) ** 2 + ((100 - s) - xyh[k][1]) ** 2 \
                                        + max(.0, m - xyh[k][2]) ** 2 <= m ** 2:
                                    flag2 = False
                        if flag2:
                            flag = True
        if flag:
            l = m
        else:
            r = m

    print(m)