works in julia now poggers
This commit is contained in:
parent
162bc5102f
commit
3a3c4a3e28
GPG Key ID:
BDA47A31A3C8EE6B
|
|
@ -0,0 +1,72 @@
|
|||
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)
|
||||
round(sum(l) / length(l), digits=8)
|
||||
113
727.py
113
727.py
|
|
@ -1,113 +0,0 @@
|
|||
import joblib
|
||||
import tqdm
|
||||
import multiprocessing
|
||||
from math import *
|
||||
|
||||
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 sqrt(dx * dx + dy * dy)
|
||||
|
||||
def circumcenter(A: Point, B: Point, C: Point) -> Point:
|
||||
D = 2.0 * (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 = 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):
|
||||
p = (a + b + c) / 2.0
|
||||
return sqrt(p * (p - a) * (p - b) * (p - c))
|
||||
|
||||
def inv_law_cosines(a, b, c):
|
||||
# given a, b, c return C
|
||||
return acos((a * a + b * b - c * c) / (2.0 * a * b))
|
||||
|
||||
def compute_d(A: Point, B: Point, C: Point, r1: float, r2: float, r3: float) -> float:
|
||||
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
|
||||
e1 = k1 + k2 + k3 + 2.0 * sqrt(k1 * k2 + k2 * k3 + k1 * k3)
|
||||
e2 = k1 + k2 + k3 - 2.0 * sqrt(k1 * k2 + k2 * k3 + k1 * k3)
|
||||
result = 1.0 / max(e1, e2)
|
||||
# re = sympy.symbols("re")
|
||||
|
||||
# re = result
|
||||
# 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)
|
||||
# print(t1 + t2 + t3, tfull, "=", tfull - t1 - t2 - t3)
|
||||
# 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) * cos(EAB)
|
||||
Ey = (r1 + result) * sin(EAB)
|
||||
E = Point(Ex, Ey)
|
||||
|
||||
# print(r1, r2, r3, D, E)
|
||||
return D.distance_to(E)
|
||||
|
||||
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 = r2 + r3
|
||||
b = r1 + r3
|
||||
c = r1 + r2
|
||||
|
||||
A = Point(0, 0)
|
||||
B = Point(c, 0)
|
||||
|
||||
CAB = inv_law_cosines(b, c, a)
|
||||
Cx = b * cos(CAB)
|
||||
Cy = b * sin(CAB)
|
||||
C = Point(Cx, Cy)
|
||||
|
||||
d = compute_d(A, B, C, r1, r2, r3)
|
||||
return d
|
||||
# return (f"d({r1}, {r2}, {r3}) = {d}")
|
||||
|
||||
total = len(list(create_loop()))
|
||||
L = joblib.Parallel(n_jobs=multiprocessing.cpu_count())(joblib.delayed(calc)(tup) for tup in tqdm.tqdm(create_loop(), total=total))
|
||||
# print(L)
|
||||
print(round(sum(L) / len(L), 8))
|
||||
|
|
@ -0,0 +1,12 @@
|
|||
import Base: *, +, -
|
||||
export Point, magnitude, +, -, distance
|
||||
|
||||
struct Point{T}
|
||||
x :: T
|
||||
y :: T
|
||||
end
|
||||
|
||||
magnitude(p) = p.x * p.x + p.y * p.y
|
||||
+(a, b) = Point(a.x + b.x, a.y + b.y)
|
||||
-(a, b) = Point(a.x - b.x, a.y - b.y)
|
||||
distance(a, b) = sqrt(magnitude(a - b))
|
||||
Loading…
Reference in New Issue