From 3a3c4a3e287a2943545c0ba82d3db14b98a34315 Mon Sep 17 00:00:00 2001 From: Michael Zhang Date: Mon, 7 Dec 2020 21:53:32 -0600 Subject: [PATCH] works in julia now poggers --- 001.py | 1 - 727.jl | 72 +++++++++++++++++++++++++++++++ 727.py | 113 ------------------------------------------------- utils/point.jl | 12 ++++++ 4 files changed, 84 insertions(+), 114 deletions(-) delete mode 100644 001.py create mode 100644 727.jl delete mode 100644 727.py create mode 100644 utils/point.jl diff --git a/001.py b/001.py deleted file mode 100644 index 2264aa5..0000000 --- a/001.py +++ /dev/null @@ -1 +0,0 @@ -print(sum(filter(lambda c: c % 3 == 0 or c % 5 == 0, range(1000)))) diff --git a/727.jl b/727.jl new file mode 100644 index 0000000..0ce6c89 --- /dev/null +++ b/727.jl @@ -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) diff --git a/727.py b/727.py deleted file mode 100644 index 88d8fdd..0000000 --- a/727.py +++ /dev/null @@ -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)) diff --git a/utils/point.jl b/utils/point.jl new file mode 100644 index 0000000..3be2d1f --- /dev/null +++ b/utils/point.jl @@ -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))