Add simultaneous congruence solving using the Chinese remainder theorem (#17)

Author: vil02

.spec/math/solve_simultaneous_congruences_spec.lua

@@ -0,0 +1,35 @@
+describe("Solve simultaneous congruences", function()
+	local ssc = require("math.solve_simultaneous_congruences")
+
+	it("should handle cases with solutions", function()
+		assert.equal(0, ssc({ { 0, 3 } }))
+		assert.equal(65, ssc({ { 1, 8 }, { 2, 9 } }))
+		assert.equal(1, ssc({ { 1, 54 }, { 1, 73 }, { 1, 997 }, { 1, 102353 } }))
+		assert.equal(39, ssc({ { 0, 3 }, { 3, 4 }, { 4, 5 } }))
+		assert.equal(23, ssc({ { 2, 3 }, { 3, 5 }, { 2, 7 } }))
+		assert.equal(34, ssc({ { 1, 3 }, { 4, 5 }, { 6, 7 } }))
+		assert.equal(388, ssc({ { 3, 7 }, { 3, 5 }, { 4, 12 } }))
+		assert.equal(87, ssc({ { 2, 5 }, { 3, 7 }, { 10, 11 } }))
+		assert.equal(29, ssc({ { 2, 3 }, { 1, 4 }, { 7, 11 } }))
+		assert.equal(125, ssc({ { 6, 7 }, { 8, 9 }, { 4, 11 }, { 8, 13 } }))
+		assert.equal(89469, ssc({ { 6, 11 }, { 13, 16 }, { 9, 21 }, { 19, 25 } }))
+	end)
+
+	it("should handle cases without solution", function()
+		assert.equal(nil, ssc({ { 5, 17 }, { 4, 17 } }))
+		assert.equal(nil, ssc({ { 3, 4 }, { 0, 6 } }))
+		assert.equal(nil, ssc({ { 3, 13 }, { 7, 1 } }))
+	end)
+
+	it("should throw error when a modulus is zero", function()
+		assert.has_error(function()
+			ssc({ { 0, 0 } })
+		end)
+	end)
+
+	it("should throw error when a modulus is negative", function()
+		assert.has_error(function()
+			ssc({ { 0, -1 } })
+		end)
+	end)
+end)

src/math/solve_simultaneous_congruences.lua

@@ -0,0 +1,42 @@
+local modular_inverse = require("math.modular_inverse")
+
+local function all_moduli_prod(congruences)
+	local res = 1
+	for _, congruence in ipairs(congruences) do
+		res = res * congruence[2]
+	end
+	return res
+end
+
+-- Solves system of simultaneous congruences, i.e. the system of the form:
+-- ```
+-- x = a_0 mod m_0
+-- x = a_1 mod m_1
+-- ...
+-- ```
+-- where `a_0`, `a_1`, ... and `m_0`, `m_1`, ... are given.
+--
+-- The system is represented by a list of pairs.
+-- The pair `{a, m}` represents a congruence `x = a mod m`.
+-- For the input `{{a_0, m_0}, {a_1, m_1}, ...}`,
+-- it finds a number `x`, such that
+-- `x = a_i mod m_i` and `0 < x < m_0 * m_1 * ...`.
+-- The implementation is based on the Chinese remainder theorem.
+return function(congruences)
+	local all_prod = all_moduli_prod(congruences)
+
+	local res = 0
+	for _, congruence in ipairs(congruences) do
+		local residue = congruence[1]
+		local modulus = congruence[2]
+		local cur_prod = math.floor(all_prod / modulus)
+		local cur_inv = modular_inverse(cur_prod, modulus)
+		if cur_inv == nil then
+			-- moduli of the congruences are not co-prime
+			return nil
+		end
+		res = (res + residue * cur_inv * cur_prod) % all_prod
+	end
+
+	return res
+end