109 lines
3.1 KiB
Python
109 lines
3.1 KiB
Python
|
import sympy
|
||
|
|
import joblib
|
||
|
|
import tqdm
|
||
|
|
|
||
|
|
def gcd(a, b):
|
||
|
|
if b == 0: return a
|
||
|
|
return gcd(b, a % b)
|
||
|
|
|
||
|
|
class Point:
|
||
|
|
def __init__(self, x, y):
|
||
|
|
self.x = x
|
||
|
|
self.y = y
|
||
|
|
|
||
|
|
def __str__(self):
|
||
|
|
return f"({self.x}, {self.y})"
|
||
|
|
|
||
|
|
def magnitude(self):
|
||
|
|
return self.x * self.x + self.y * self.y
|
||
|
|
|
||
|
|
def distance_to(self, other):
|
||
|
|
dx = other.x - self.x
|
||
|
|
dy = other.y - self.y
|
||
|
|
return sympy.sqrt(dx * dx + dy * dy)
|
||
|
|
|
||
|
|
def circumcenter(A: Point, B: Point, C: Point) -> Point:
|
||
|
|
D = sympy.S(2) * (A.x * (B.y - C.y) + B.x * (C.y - A.y) + C.x * (A.y - B.y))
|
||
|
|
Ux = (A.magnitude() * (B.y - C.y) + B.magnitude() * (C.y - A.y) + C.magnitude() * (A.y - B.y)) / D
|
||
|
|
Uy = (A.magnitude() * (C.x - B.x) + B.magnitude() * (A.x - C.x) + C.magnitude() * (B.x - A.x)) / D
|
||
|
|
return Point(Ux, Uy)
|
||
|
|
|
||
|
|
def lerp(src: Point, dst: Point, dist: float) -> Point:
|
||
|
|
total = src.distance_to(dst)
|
||
|
|
ratio = sympy.Rational(dist, total)
|
||
|
|
dx = src.x + (dst.x - src.x) * ratio
|
||
|
|
dy = src.y + (dst.y - src.y) * ratio
|
||
|
|
return Point(dx, dy)
|
||
|
|
|
||
|
|
def heron(a, b, c):
|
||
|
|
a, b, c = sympy.S(a), sympy.S(b), sympy.S(c)
|
||
|
|
p = (a + b + c) / sympy.S(2)
|
||
|
|
return sympy.sqrt(p * (p - a) * (p - b) * (p - c))
|
||
|
|
|
||
|
|
def inv_law_cosines(a, b, c):
|
||
|
|
# given a, b, c return C
|
||
|
|
return sympy.acos((a * a + b * b - c * c) / (sympy.S(2) * a * b))
|
||
|
|
|
||
|
|
|
||
|
|
def compute_d(A: Point, B: Point, C: Point, r1: float, r2: float, r3: float) -> float:
|
||
|
|
r1, r2, r3 = sympy.S(r1), sympy.S(r2), sympy.S(r3)
|
||
|
|
|
||
|
|
ABmid = lerp(A, B, r1)
|
||
|
|
BCmid = lerp(B, C, r2)
|
||
|
|
CAmid = lerp(C, A, r3)
|
||
|
|
D = circumcenter(ABmid, BCmid, CAmid)
|
||
|
|
|
||
|
|
re = sympy.symbols("re")
|
||
|
|
t1 = heron(r1 + r3, r1 + re, r3 + re)
|
||
|
|
t2 = heron(r1 + r2, r1 + re, r2 + re)
|
||
|
|
t3 = heron(r2 + r3, r2 + re, r3 + re)
|
||
|
|
tfull = heron(r1 + r2, r1 + r3, r2 + r3)
|
||
|
|
eq = sympy.Eq(t1 + t2 + t3, tfull)
|
||
|
|
result = sympy.solve(eq, re)
|
||
|
|
result = result[0].evalf()
|
||
|
|
# print(f"radius of E is {result}")
|
||
|
|
|
||
|
|
EAB = inv_law_cosines(r1 + result, r1 + r2, r2 + result)
|
||
|
|
Ex = (r1 + result) * sympy.cos(EAB)
|
||
|
|
Ey = (r1 + result) * sympy.sin(EAB)
|
||
|
|
E = Point(Ex, Ey)
|
||
|
|
|
||
|
|
# print(D, E)
|
||
|
|
return D.distance_to(E).evalf()
|
||
|
|
|
||
|
|
total = 0
|
||
|
|
count = 0
|
||
|
|
def create_loop():
|
||
|
|
for r1 in range(1, 101 - 2):
|
||
|
|
for r2 in range(r1 + 1, 101 - 1):
|
||
|
|
for r3 in range(r2 + 1, 101):
|
||
|
|
if gcd(gcd(r1, r2), r3) != 1: continue
|
||
|
|
yield (r1, r2, r3)
|
||
|
|
|
||
|
|
def calc(x):
|
||
|
|
r1, r2, r3 = x
|
||
|
|
# build a triangle out of the side lengths
|
||
|
|
# A is at (0, 0), a = r2 + r3
|
||
|
|
# B is at (c, 0), b = r1 + r3
|
||
|
|
# c = r1 + r2
|
||
|
|
# C is at ()
|
||
|
|
a = sympy.S(r2 + r3)
|
||
|
|
b = sympy.S(r1 + r3)
|
||
|
|
c = sympy.S(r1 + r2)
|
||
|
|
|
||
|
|
A = Point(sympy.S(0), sympy.S(0))
|
||
|
|
B = Point(c, sympy.S(0))
|
||
|
|
|
||
|
|
CAB = sympy.acos(sympy.Rational(b * b + c * c - a * a, 2.0 * b * c))
|
||
|
|
Cx = b * sympy.cos(CAB)
|
||
|
|
Cy = b * sympy.sin(CAB)
|
||
|
|
C = Point(Cx, Cy)
|
||
|
|
|
||
|
|
d = compute_d(A, B, C, r1, r2, r3)
|
||
|
|
return d
|
||
|
|
print(f"d({r1}, {r2}, {r3}) = {d}")
|
||
|
|
|
||
|
|
total = len(list(create_loop()))
|
||
|
|
L = joblib.Parallel(n_jobs=8)(joblib.delayed(calc)(tup) for tup in tqdm.tqdm(create_loop(), total=total))
|
||
|
|
print(sum(L) / len(L))
|