include("utils/point.jl")

# list of things we want to iterate over
domain = Iterators.filter(
  a -> a[3] > a[2] && a[2] > a[1] && gcd(collect(a)) == 1,
  Iterators.product(1:100, 1:100, 1:100)
)

# law of cosines, determines angle based on the side lengths
inv_law_cosines(a, b, c) = acos((a * a + b * b - c * c) / (2.0 * a * b))

function circumcenter(A, B, C)
  magA = magnitude(A)
  magB = magnitude(B)
  magC = magnitude(C)
  D = 2.0 * (A.x * (B.y - C.y) + B.x * (C.y - A.y) + C.x * (A.y - B.y))
  Ux = (magA * (B.y - C.y) + magB * (C.y - A.y) + magC * (A.y - B.y)) / D
  Uy = (magA * (C.x - B.x) + magB * (A.x - C.x) + magC * (B.x - A.x)) / D
  Point(Ux, Uy)
end

function lerp(src, dst, dist)
  total = distance(src, dst)
  ratio = dist / total
  dx = src.x + (dst.x - src.x) * ratio
  dy = src.y + (dst.y - src.y) * ratio
  Point(dx, dy)
end

function compute_d(A, B, C, r1, r2, r3)
  ABmid = lerp(A, B, r1)
  BCmid = lerp(B, C, r2)
  CAmid = lerp(C, A, r3)
  D = circumcenter(ABmid, BCmid, CAmid)

  k1, k2, k3 = 1.0 / r1, 1.0 / r2, 1.0 / r3
  s = k1 + k2 + k3 
  s2 = k1 * k2 + k2 * k3 + k1 * k3
  e1 = s + 2.0 * sqrt(s2)
  e2 = s - 2.0 * sqrt(s2)
  re = 1.0 / max(e1, e2)

  EAB = inv_law_cosines(r1 + re, r1 + r2, r2 + re)
  Ex = (r1 + re) * cos(EAB)
  Ey = (r1 + re) * sin(EAB)
  E = Point(Ex, Ey)
  
  distance(D, E)
end

function calc(args)
  r1 = args[1]
  r2 = args[2]
  r3 = args[3]
  
  a = r2 + r3
  b = r1 + r3
  c = r1 + r2

  A = Point(0.0, 0.0)
  B = Point(Float64(c), 0.0)

  CAB = inv_law_cosines(b, c, a)
  Cx = b * cos(CAB)
  Cy = b * sin(CAB)
  C = Point(Cx, Cy)
  
  compute_d(A, B, C, r1, r2, r3)
end

l = map(calc, domain)
print(round(sum(l) / length(l), digits=8))