"""
Jane Street Puzzle
November 2024: Beside the Point
By Kauri Beckmann
kauri.b@outlook.co.nz
03/11/2024
"""

# I should really upload my hand calculations...
# This is just the python code that executes the mathematical solution
# As the resulting integral is too complex even for the likes of Wolfram Alpha!

# Mathematical approach:
    # One point spawns closer to the edge.
        # We assume this point spawns in a triangle 1/8 the area of the square, and extrapolate on the basis of symmetry.
    # From this point, we find two 1/4-circles, with radius drawn from point to corner of square.
        # If the second point lands in the area of these two circles, but NOT the overlap, then an equidistant path to border can be found.
    # Write a formula defining this area
    # Integrate this formula over the area that the first point can spawn.
    # equals Probability an equidistant path can be drawn.

import math
from scipy.integrate import dblquad

def integrate_me(y,x):
    # This calculates the area of two 1/4-circles, minus the overlap.

    d = 1 # simply the size of the square

    # Radius
    Ra = math.sqrt(x**2 + y**2)
    # Area
    Aa = Ra**2 * math.pi / 4


    Rb = math.sqrt((1-x)**2 + y**2)
    Ab = Rb**2 * math.pi / 4

    alpha = math.acos((Ra**2 - Rb**2 + d**2)/(2 * Ra * d))
    beta = math.acos((Rb**2 - Ra**2 + d**2)/(2 * Rb * d))

    A_overlap = 0.5 * alpha * Ra**2 + 0.5 * beta * Rb**2 - 0.25 * Ra**2 * math.sin(2 * alpha) - 0.25 * Rb**2 * math.sin(2 * beta)

    A = Aa + Ab - 2*A_overlap
    return A

# We set limits the first point to spawn in a triangle, 1/8 the area of the square.
x_lower = 0
x_upper = 0.5
y_lower = 0
y_upper = lambda x: x  # y goes from 0 to x

# Compute the double integral
integral, error = dblquad(integrate_me, x_lower, x_upper, lambda x: y_lower, lambda x: y_upper(x))

# Confirming computational error well below required precision
print(f"Computational error = ", error)

# The area of the triangular region
area_triangle = 0.125

# We only integrated within the triangle; divide by this area.
mean_value = integral / area_triangle
print(f"Probability = ", mean_value)

Computational error =  2.485643838466147e-15
Probability =  0.491407578838308