diff --git a/math/ncr_modulo_p.cpp b/math/ncr_modulo_p.cpp index 03aecf67e88..cbf5e9e1a00 100644 --- a/math/ncr_modulo_p.cpp +++ b/math/ncr_modulo_p.cpp @@ -10,7 +10,7 @@ */ #include /// for assert -#include /// for io operations +#include /// for std::cout #include /// for std::vector /** @@ -25,71 +25,95 @@ namespace math { * implementation. */ namespace ncr_modulo_p { + /** - * @brief Class which contains all methods required for calculating nCr mod p + * @namespace utils + * @brief this namespace contains the definitions of the functions called from + * the class math::ncr_modulo_p::NCRModuloP */ -class NCRModuloP { - private: - std::vector fac{}; /// stores precomputed factorial(i) % p value - uint64_t p = 0; /// the p from (nCr % p) - - public: - /** Constructor which precomputes the values of n! % mod from n=0 to size - * and stores them in vector 'fac' - * @params[in] the numbers 'size', 'mod' - */ - NCRModuloP(const uint64_t& size, const uint64_t& mod) { - p = mod; - fac = std::vector(size); - fac[0] = 1; - for (int i = 1; i <= size; i++) { - fac[i] = (fac[i - 1] * i) % p; - } +namespace utils { +/** + * @brief finds the values x and y such that a*x + b*y = gcd(a,b) + * + * @param[in] a the first input of the gcd + * @param[in] a the second input of the gcd + * @param[out] x the Bézout coefficient of a + * @param[out] y the Bézout coefficient of b + * @return the gcd of a and b + */ +int64_t gcdExtended(const int64_t& a, const int64_t& b, int64_t& x, + int64_t& y) { + if (a == 0) { + x = 0; + y = 1; + return b; } - /** Finds the value of x, y such that a*x + b*y = gcd(a,b) - * - * @params[in] the numbers 'a', 'b' and address of 'x' and 'y' from above - * equation - * @returns the gcd of a and b - */ - uint64_t gcdExtended(const uint64_t& a, const uint64_t& b, int64_t* x, - int64_t* y) { - if (a == 0) { - *x = 0, *y = 1; - return b; - } + int64_t x1 = 0, y1 = 0; + const int64_t gcd = gcdExtended(b % a, a, x1, y1); - int64_t x1 = 0, y1 = 0; - uint64_t gcd = gcdExtended(b % a, a, &x1, &y1); + x = y1 - (b / a) * x1; + y = x1; + return gcd; +} - *x = y1 - (b / a) * x1; - *y = x1; - return gcd; +/** Find modular inverse of a modulo m i.e. a number x such that (a*x)%m = 1 + * + * @param[in] a the number for which the modular inverse is queried + * @param[in] m the modulus + * @return the inverce of a modulo m, if it exists, -1 otherwise + */ +int64_t modInverse(const int64_t& a, const int64_t& m) { + int64_t x = 0, y = 0; + const int64_t g = gcdExtended(a, m, x, y); + if (g != 1) { // modular inverse doesn't exist + return -1; + } else { + return ((x + m) % m); } +} +} // namespace utils +/** + * @brief Class which contains all methods required for calculating nCr mod p + */ +class NCRModuloP { + private: + const int64_t p = 0; /// the p from (nCr % p) + const std::vector + fac; /// stores precomputed factorial(i) % p value - /** Find modular inverse of a with m i.e. a number x such that (a*x)%m = 1 - * - * @params[in] the numbers 'a' and 'm' from above equation - * @returns the modular inverse of a + /** + * @brief computes the array of values of factorials reduced modulo mod + * @param max_arg_val argument of the last factorial stored in the result + * @param mod value of the divisor used to reduce factorials + * @return vector storing factorials of the numbers 0, ..., max_arg_val + * reduced modulo mod */ - int64_t modInverse(const uint64_t& a, const uint64_t& m) { - int64_t x = 0, y = 0; - uint64_t g = gcdExtended(a, m, &x, &y); - if (g != 1) { // modular inverse doesn't exist - return -1; - } else { - int64_t res = ((x + m) % m); - return res; + static std::vector computeFactorialsMod(const int64_t& max_arg_val, + const int64_t& mod) { + auto res = std::vector(max_arg_val + 1); + res[0] = 1; + for (int64_t i = 1; i <= max_arg_val; i++) { + res[i] = (res[i - 1] * i) % mod; } + return res; } - /** Find nCr % p - * - * @params[in] the numbers 'n', 'r' and 'p' - * @returns the value nCr % p + public: + /** + * @brief constructs an NCRModuloP object allowing to compute (nCr)%p for + * inputs from 0 to size */ - int64_t ncr(const uint64_t& n, const uint64_t& r, const uint64_t& p) { + NCRModuloP(const int64_t& size, const int64_t& p) + : p(p), fac(computeFactorialsMod(size, p)) {} + + /** + * @brief computes nCr % p + * @param[in] n the number of objects to be chosen + * @param[in] r the number of objects to choose from + * @return the value nCr % p + */ + int64_t ncr(const int64_t& n, const int64_t& r) const { // Base cases if (r > n) { return 0; @@ -101,50 +125,71 @@ class NCRModuloP { return 1; } // fac is a global array with fac[r] = (r! % p) - int64_t denominator = modInverse(fac[r], p); - if (denominator < 0) { // modular inverse doesn't exist - return -1; - } - denominator = (denominator * modInverse(fac[n - r], p)) % p; - if (denominator < 0) { // modular inverse doesn't exist + const auto denominator = (fac[r] * fac[n - r]) % p; + const auto denominator_inv = utils::modInverse(denominator, p); + if (denominator_inv < 0) { // modular inverse doesn't exist return -1; } - return (fac[n] * denominator) % p; + return (fac[n] * denominator_inv) % p; } }; } // namespace ncr_modulo_p } // namespace math /** - * @brief Test implementations - * @param ncrObj object which contains the precomputed factorial values and - * ncr function - * @returns void + * @brief tests math::ncr_modulo_p::NCRModuloP */ -static void tests(math::ncr_modulo_p::NCRModuloP ncrObj) { - // (52323 C 26161) % (1e9 + 7) = 224944353 - assert(ncrObj.ncr(52323, 26161, 1000000007) == 224944353); - // 6 C 2 = 30, 30%5 = 0 - assert(ncrObj.ncr(6, 2, 5) == 0); - // 7C3 = 35, 35 % 29 = 8 - assert(ncrObj.ncr(7, 3, 29) == 6); +static void tests() { + struct TestCase { + const int64_t size; + const int64_t p; + const int64_t n; + const int64_t r; + const int64_t expected; + + TestCase(const int64_t size, const int64_t p, const int64_t n, + const int64_t r, const int64_t expected) + : size(size), p(p), n(n), r(r), expected(expected) {} + }; + const std::vector test_cases = { + TestCase(60000, 1000000007, 52323, 26161, 224944353), + TestCase(20, 5, 6, 2, 30 % 5), + TestCase(100, 29, 7, 3, 35 % 29), + TestCase(1000, 13, 10, 3, 120 % 13), + TestCase(20, 17, 1, 10, 0), + TestCase(45, 19, 23, 1, 23 % 19), + TestCase(45, 19, 23, 0, 1), + TestCase(45, 19, 23, 23, 1), + TestCase(20, 9, 10, 2, -1)}; + for (const auto& tc : test_cases) { + assert(math::ncr_modulo_p::NCRModuloP(tc.size, tc.p).ncr(tc.n, tc.r) == + tc.expected); + } + + std::cout << "\n\nAll tests have successfully passed!\n"; } /** - * @brief Main function - * @returns 0 on exit + * @brief example showing the usage of the math::ncr_modulo_p::NCRModuloP class */ -int main() { - // populate the fac array - const uint64_t size = 1e6 + 1; - const uint64_t p = 1e9 + 7; - math::ncr_modulo_p::NCRModuloP ncrObj = - math::ncr_modulo_p::NCRModuloP(size, p); - // test 6Ci for i=0 to 7 +void example() { + const int64_t size = 1e6 + 1; + const int64_t p = 1e9 + 7; + + // the ncrObj contains the precomputed values of factorials modulo p for + // values from 0 to size + const auto ncrObj = math::ncr_modulo_p::NCRModuloP(size, p); + + // having the ncrObj we can efficiently query the values of (n C r)%p + // note that time of the computation does not depend on size for (int i = 0; i <= 7; i++) { - std::cout << 6 << "C" << i << " = " << ncrObj.ncr(6, i, p) << "\n"; + std::cout << 6 << "C" << i << " mod " << p << " = " << ncrObj.ncr(6, i) + << "\n"; } - tests(ncrObj); // execute the tests - std::cout << "Assertions passed\n"; +} + +int main() { + tests(); + example(); return 0; }