Refactor ChineseRemainderTheorem: reuse PrimeFactorization, add tests

Remove duplicated primeFactorization, pollardRho, isPrime, and gcf methods
in favor of the shared PrimeFactorization class. Simplify prime-power grouping
loop and clean up Javadoc. Add 17 JUnit 5 tests covering eliminateCoefficient,
crt, reduce, egcd, and full reduce+crt integration.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: williamfiset

src/main/java/com/williamfiset/algorithms/math/ChineseRemainderTheorem.java

@@ -1,20 +1,14 @@
 /**
- * Use the chinese remainder theorem to solve a set of congruence equations.
+ * Solve a set of congruence equations using the Chinese Remainder Theorem.
  *
- * <p>The first method (eliminateCoefficient) is used to reduce an equation of the form cx≡a(mod
- * m)cx≡a(mod m) to the form x≡a_new(mod m_new)x≡anew(mod m_new), which gets rids of the
- * coefficient. A value of null is returned if the coefficient cannot be eliminated.
+ * <p>{@link #eliminateCoefficient(long, long, long)} reduces cx ≡ a (mod m) to x ≡ a' (mod m'),
+ * removing the coefficient. Returns null if unsolvable.
  *
- * <p>The second method (reduce) is used to reduce a set of equations so that the moduli become
- * pairwise co-prime (which means that we can apply the Chinese Remainder Theorem). The input and
- * output are of the form x≡a_0(mod m_0),...,x≡a_n−1(mod m_n−1)x≡a_0(mod m_0),...,x≡a_n−1(mod
- * m_n−1). Note that the number of equations may change during this process. A value of null is
- * returned if the set of equations cannot be reduced to co-prime moduli.
+ * <p>{@link #reduce(long[], long[])} rewrites a system so that moduli become pairwise coprime,
+ * enabling CRT. The number of equations may change. Returns null if the system is inconsistent.
  *
- * <p>The third method (crt) is the actual Chinese Remainder Theorem. It assumes that all pairs of
- * moduli are co-prime to one another. This solves a set of equations of the form x≡a_0(mod
- * m_0),...,x≡v_n−1(mod m_n−1)x≡a_0(mod m_0),...,x≡v_n−1(mod m_n−1). It's output is of the form
- * x≡a_new(mod m_new)x≡a_new(mod m_new).
+ * <p>{@link #crt(long[], long[])} is the core CRT step. It assumes all moduli are pairwise coprime
+ * and returns x ≡ a' (mod M) where M is the product of all moduli.
  *
  * @author Micah Stairs
  */
@@ -24,12 +18,16 @@
 
 public class ChineseRemainderTheorem {
 
-  // eliminateCoefficient() takes cx≡a(mod m) and gives x≡a_new(mod m_new).
+  /**
+   * Reduces cx ≡ a (mod m) to x ≡ a' (mod m').
+   *
+   * @return {a', m'} or null if unsolvable.
+   */
   public static long[] eliminateCoefficient(long c, long a, long m) {
-
     long d = egcd(c, m)[0];
 
-    if (a % d != 0) return null;
+    if (a % d != 0)
+      return null;
 
     c /= d;
     a /= d;
@@ -42,36 +40,35 @@ public static long[] eliminateCoefficient(long c, long a, long m) {
     return new long[] {a, m};
   }
 
-  // reduce() takes a set of equations and reduces them to an equivalent
-  // set with pairwise co-prime moduli (or null if not solvable).
+  /**
+   * Reduces a system x ≡ a[i] (mod m[i]) to an equivalent system with pairwise coprime moduli.
+   *
+   * @return {a[], m[]} with coprime moduli, or null if the system is inconsistent.
+   */
   public static long[][] reduce(long[] a, long[] m) {
+    List<Long> aNew = new ArrayList<>();
+    List<Long> mNew = new ArrayList<>();
 
-    List<Long> aNew = new ArrayList<Long>();
-    List<Long> mNew = new ArrayList<Long>();
-
-    // Split up each equation into prime factors
+    // Split each modulus into prime-power factors
     for (int i = 0; i < a.length; i++) {
-      List<Long> factors = primeFactorization(m[i]);
+      List<Long> factors = PrimeFactorization.primeFactorization(m[i]);
       Collections.sort(factors);
-      ListIterator<Long> iterator = factors.listIterator();
-      while (iterator.hasNext()) {
-        long val = iterator.next();
-        long total = val;
-        while (iterator.hasNext()) {
-          long nextVal = iterator.next();
-          if (nextVal == val) {
-            total *= val;
-          } else {
-            iterator.previous();
-            break;
-          }
+
+      int j = 0;
+      while (j < factors.size()) {
+        long p = factors.get(j);
+        long pk = p;
+        j++;
+        while (j < factors.size() && factors.get(j) == p) {
+          pk *= p;
+          j++;
         }
-        aNew.add(a[i] % total);
-        mNew.add(total);
+        aNew.add(a[i] % pk);
+        mNew.add(pk);
       }
     }
 
-    // Throw away repeated information and look for conflicts
+    // Remove redundant equations and detect conflicts
     for (int i = 0; i < aNew.size(); i++) {
       for (int j = i + 1; j < aNew.size(); j++) {
         if (mNew.get(i) % mNew.get(j) == 0 || mNew.get(j) % mNew.get(i) == 0) {
@@ -81,111 +78,58 @@ public static long[][] reduce(long[] a, long[] m) {
               mNew.remove(j);
               j--;
               continue;
-            } else return null;
+            } else
+              return null;
           } else {
             if ((aNew.get(j) % mNew.get(i)) == aNew.get(i)) {
               aNew.remove(i);
               mNew.remove(i);
               i--;
               break;
-            } else return null;
+            } else
+              return null;
           }
         }
       }
     }
 
-    // Put result into an array
     long[][] res = new long[2][aNew.size()];
     for (int i = 0; i < aNew.size(); i++) {
       res[0][i] = aNew.get(i);
       res[1][i] = mNew.get(i);
     }
-
     return res;
   }
 
+  /**
+   * Solves x ≡ a[i] (mod m[i]) assuming all moduli are pairwise coprime.
+   *
+   * @return {x, M} where M is the product of all moduli.
+   */
   public static long[] crt(long[] a, long[] m) {
-
     long M = 1;
-    for (int i = 0; i < m.length; i++) M *= m[i];
-
-    long[] inv = new long[a.length];
-    for (int i = 0; i < inv.length; i++) inv[i] = egcd(M / m[i], m[i])[1];
+    for (long mi : m)
+      M *= mi;
 
     long x = 0;
     for (int i = 0; i < m.length; i++) {
-      x += (M / m[i]) * a[i] * inv[i]; // Overflow could occur here
+      long Mi = M / m[i];
+      long inv = egcd(Mi, m[i])[1];
+      x += Mi * a[i] * inv;
       x = ((x % M) + M) % M;
     }
 
     return new long[] {x, M};
   }
 
-  private static ArrayList<Long> primeFactorization(long n) {
-    ArrayList<Long> factors = new ArrayList<Long>();
-    if (n <= 0) throw new IllegalArgumentException();
-    else if (n == 1) return factors;
-    PriorityQueue<Long> divisorQueue = new PriorityQueue<Long>();
-    divisorQueue.add(n);
-    while (!divisorQueue.isEmpty()) {
-      long divisor = divisorQueue.remove();
-      if (isPrime(divisor)) {
-        factors.add(divisor);
-        continue;
-      }
-      long next_divisor = pollardRho(divisor);
-      if (next_divisor == divisor) {
-        divisorQueue.add(divisor);
-      } else {
-        divisorQueue.add(next_divisor);
-        divisorQueue.add(divisor / next_divisor);
-      }
-    }
-    return factors;
-  }
-
-  private static long pollardRho(long n) {
-    if (n % 2 == 0) return 2;
-    // Get a number in the range [2, 10^6]
-    long x = 2 + (long) (999999 * Math.random());
-    long c = 2 + (long) (999999 * Math.random());
-    long y = x;
-    long d = 1;
-    while (d == 1) {
-      x = (x * x + c) % n;
-      y = (y * y + c) % n;
-      y = (y * y + c) % n;
-      d = gcf(Math.abs(x - y), n);
-      if (d == n) break;
-    }
-    return d;
-  }
-
-  // Extended euclidean algorithm
-  private static long[] egcd(long a, long b) {
-    if (b == 0) return new long[] {a, 1, 0};
-    else {
-      long[] ret = egcd(b, a % b);
-      long tmp = ret[1] - ret[2] * (a / b);
-      ret[1] = ret[2];
-      ret[2] = tmp;
-      return ret;
-    }
-  }
-
-  private static long gcf(long a, long b) {
-    return b == 0 ? a : gcf(b, a % b);
-  }
-
-  private static boolean isPrime(long n) {
-    if (n < 2) return false;
-    if (n == 2 || n == 3) return true;
-    if (n % 2 == 0 || n % 3 == 0) return false;
-
-    int limit = (int) Math.sqrt(n);
-
-    for (int i = 5; i <= limit; i += 6) if (n % i == 0 || n % (i + 2) == 0) return false;
-
-    return true;
+  /** Extended Euclidean algorithm. Returns {gcd(a,b), x, y} where ax + by = gcd(a,b). */
+  static long[] egcd(long a, long b) {
+    if (b == 0)
+      return new long[] {a, 1, 0};
+    long[] ret = egcd(b, a % b);
+    long tmp = ret[1] - ret[2] * (a / b);
+    ret[1] = ret[2];
+    ret[2] = tmp;
+    return ret;
   }
 }

src/test/java/com/williamfiset/algorithms/math/BUILD

@@ -47,3 +47,14 @@ java_test(
     runtime_deps = JUNIT5_RUNTIME_DEPS,
     deps = TEST_DEPS,
 )
+
+# bazel test //src/test/java/com/williamfiset/algorithms/math:ChineseRemainderTheoremTest
+java_test(
+    name = "ChineseRemainderTheoremTest",
+    srcs = ["ChineseRemainderTheoremTest.java"],
+    main_class = "org.junit.platform.console.ConsoleLauncher",
+    use_testrunner = False,
+    args = ["--select-class=com.williamfiset.algorithms.math.ChineseRemainderTheoremTest"],
+    runtime_deps = JUNIT5_RUNTIME_DEPS,
+    deps = TEST_DEPS,
+)