Day 13

chinese_remainder/src/lib.rs

@@ -0,0 +1,64 @@
+/// Represents the GCD computed as part of the Extended Euclidean Algorithm.
+/// This invariant should always be valid : ax + by = gcd
+#[derive(Debug,PartialEq)]
+pub struct EGCD {
+    pub a: i64,
+    pub b: i64,
+    pub x: i64,
+    pub y: i64,
+    pub gcd: i64,
+}
+
+/// Computes the Extended Euclidean Algorithm's solution to finding a & b's GCD
+/// Additionally to the standard Euclidean Algorithm, it provides x and y so
+/// that ax + by = GCD
+pub fn egcd(a: i64, b: i64) -> EGCD {
+    // If b == 0, then a == gcd(a,b)
+    // Then gcd = 1 * a + 0 * b follows
+    if b == 0 {
+        EGCD {
+            a,
+            b,
+            x: 1,
+            // Anything other than 0 gives a valid solution but
+            // With different end x and y
+            y: 0,
+            gcd: a,
+        }
+    } else {
+        let rec = egcd(b, a%b);
+        EGCD {
+            a,
+            b,
+            x: rec.y,
+            y: rec.x - rec.y * (a / b),
+            gcd: rec.gcd,
+        }
+    }
+}
+
+/// Solves the chinese remainder problem, stated as follows:
+/// ```noexec
+/// for i in residues.len(),
+///     (chinese_remainder(residues, moduli) - residues[i]) / moduli[i] == 0
+/// ``` 
+pub fn chinese_remainder(residues: &[i64], moduli: &[i64]) -> Option<i64> {
+    // TODO: Filter out 0s in moduli (meaning big_n is 0 -> div by 0)
+    let big_n: i64 = moduli.iter().product();
+    Some(moduli.iter()
+        .map(|&ni| {
+            egcd(ni, big_n / ni)
+        })
+        .zip(residues)
+        .map(|(egcd,ai)| {
+            if egcd.gcd != 1 {
+                // Fail in case moduli are not co-prime
+                None
+            } else {
+                Some(ai * egcd.y * egcd.b)
+            }
+        })
+        .sum::<Option<i64>>()?
+        .rem_euclid(big_n)) // Sum on Option<_> returns None if the Iter
+                            // contains a None
+}