From e16232c3a35be238664fc7179ab8942eb3888385 Mon Sep 17 00:00:00 2001
From: Seongbeom Park <sparkamita90@gmail.com>
Date: Tue, 11 Jan 2022 23:45:23 +0900
Subject: [PATCH] finish level3/doomsday-fuel

---
 README.md                        |  5 ++
 level3/doomsday-fuel/solution.py | 97 ++++++++++++++++++++++++++++++++
 2 files changed, 102 insertions(+)
 create mode 100644 level3/doomsday-fuel/solution.py

diff --git a/README.md b/README.md
index 73b9a82..646964c 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,7 @@
 # Google Foobar Challenge
 
+## Level3
+### doomsday-fuel
+* Reference
+  * L26.6 Absorption Probabilities, MIT OpenCourseWare, https://www.youtube.com/watch?v=vEsUsaK1HBk
+
diff --git a/level3/doomsday-fuel/solution.py b/level3/doomsday-fuel/solution.py
new file mode 100644
index 0000000..78b5094
--- /dev/null
+++ b/level3/doomsday-fuel/solution.py
@@ -0,0 +1,97 @@
+import numpy as np
+
+def solution(m):
+    stable_ores, unstable_ores = categorize_ores(m)
+
+    if 0 in stable_ores:
+        return [1]+[0]*(len(stable_ores)-1)+[1]
+
+    p = np.array(m)
+    q = np.zeros(p.shape, dtype=int)
+    for i in range(len(m)):
+        q[i][i] = sum(p[i])
+    r = q - p
+
+    t = gaussian_elimination(r, stable_ores, unstable_ores)
+
+    numerators = -t[0][stable_ores]
+    denominator = sum(numerators)
+
+    result = np.hstack((numerators, denominator))
+    g = abs_gcd(*list(result))
+
+    return list(result / g)
+
+def categorize_ores(m):
+    stable_ores = []
+    unstable_ores = []
+    for i in range(len(m)):
+        if sum(m[i]) == 0:
+            stable_ores += [i]
+        else:
+            unstable_ores += [i]
+    return stable_ores, unstable_ores
+
+def gaussian_elimination(r, s, u):
+    arr = r
+    for i in range(len(u)):
+        for j in range(i+1, len(u)):
+            row_up, row_down = arr[u[i]], arr[u[j]]
+            if row_down[u[i]] != 0:
+                arr[u[j]] = row_down*arr[u[i]][u[i]] - row_up*arr[u[j]][u[i]]
+                arr[u[j]] = arr[u[j]] / abs_gcd(*list(arr[u[j]]))
+    for i in reversed(range(len(u))):
+        for j in reversed(range(i)):
+            row_down, row_up = arr[u[i]], arr[u[j]]
+            if row_up[u[i]] != 0:
+                arr[u[j]] = row_up*row_down[u[i]] - row_down*row_up[u[i]]
+                arr[u[j]] = arr[u[j]] / abs_gcd(*list(arr[u[j]]))
+    return arr
+
+def abs_gcd(*args):
+    if len(args) == 1:
+        return abs(args[0])
+    if len(args) == 2:
+        m, n = abs(args[0]), abs(args[1])
+        while n != 0:
+            t = m % n
+            m, n = n, t
+        return abs(m)
+    return abs_gcd(args[0], abs_gcd(*args[1:]))
+
+tests = [
+        [
+            [[0, 2, 1, 0, 0], [0, 0, 0, 3, 4], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]],
+            [7, 6, 8, 21]
+        ],
+        [
+            [[0, 1, 0, 0, 0, 1], [4, 0, 0, 3, 2, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]],
+            [0, 3, 2, 9, 14]
+        ],
+        [
+            [[0, 1, 0, 0, 0, 1], [4, 0, 0, 3, 2, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 0]],
+            [0, 3, 11, 14]
+        ],
+        [
+            [[0, 1, 0, 0, 0, 1], [4, 0, 0, 3, 2, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]],
+            [0, 1, 6, 7]
+        ],
+        [
+            [[0]],
+            [1, 1]
+        ],
+        [
+            [[0, 1], [0, 0]],
+            [1, 1]
+        ],
+        [
+            [[0, 5, 0, 3], [0, 0, 4, 1], [0, 0, 0, 0], [0, 0, 0, 0]],
+            [1, 1, 2]
+        ],
+]
+
+#print(solution(tests[1][0]), tests[1][1])
+
+for test in tests:
+    result = solution(test[0])
+    print(result == test[1], result, test[1])