Euler #587

from numpy import pi
from numpy import sqrt
full = (4 - pi)/4


def intersection(circles):
    a = (circles**2 + 1)/circles**2
    b = -2/circles - 2
    c = 1
    return (-b - sqrt(b**2 - 4*a*c))/(2*a)


def integrate(f, a, b, dx=0.00001):
    area = 0
    while a <= b:
        area += f(a)*dx
        a += dx
    return area


def height(circles):
    return intersection(circles)/circles

for circleNum in range(1, 10**5):
    xCoord = intersection(circleNum)
    rightHalf = integrate(lambda x: -sqrt(1-(x-1)**2)+1, xCoord, 1)
    leftHalf = 0.5*xCoord*height(circleNum)
    fractionArea = (rightHalf + leftHalf)/full
    print(fractionArea)
    if fractionArea < 0.001:
        print(circleNum)
        quit()