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  x 1  y 1  h 1 : x n  y n  h n  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 × 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, x i , y i , and h i . ( x i , y i , 0) and h i 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 < x i < 100, 0 < y i < 100 and 1 ≤ h i ≤ 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

Solution (Python):
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) 
        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)

Solution (Python):
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)]

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

            
            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)