Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: Add ncr mod p code #1325

Merged
merged 18 commits into from Nov 22, 2020
Merged
Show file tree
Hide file tree
Changes from 7 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions DIRECTORY.md
Expand Up @@ -143,6 +143,7 @@
* [Least Common Multiple](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/least_common_multiple.cpp)
* [Miller Rabin](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/miller_rabin.cpp)
* [Modular Inverse Fermat Little Theorem](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/modular_inverse_fermat_little_theorem.cpp)
* [Ncr Modulo P](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/ncr_modulo_p.cpp)
* [Number Of Positive Divisors](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/number_of_positive_divisors.cpp)
* [Power For Huge Numbers](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/power_for_huge_numbers.cpp)
* [Prime Factorization](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/prime_factorization.cpp)
Expand Down
138 changes: 138 additions & 0 deletions math/ncr_modulo_p.cpp
@@ -0,0 +1,138 @@
/**
* @file
* @brief This program aims at calculating nCr modulo p
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add a detailed description using @details.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Provide a Wikipedia link in Markdown format (if available/possible).
If there's no Wikipedia link, add another web source for algorithm explanation. 🙂

* @details nCr is defined as n! / (r! * (n-r)!) where n! represents factorial
* of n. In many cases, the value of nCr is too large to fit in a 64 bit
* integer. Hence, in competitive programming, there are many problems or
* subproblems to compute nCr modulo p where p is a given number.
* @author [Kaustubh Damania](https://github.com/KaustubhDamania)
*/

#include <cassert> /// for assert
#include <iostream> /// for io operations
#include <vector> /// for std::vector

/**
* @namespace math
* @brief Mathematical algorithms
*/
namespace math {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
namespace math {
namespace math {
/**
* @namespace ncr_modulo_p
* @brief Functions for nCr modulo p implementation.
*/
namespace ncr_modulo_p {

/**
* @brief Class which contains all methods required for calculating nCr mod p
*/
class NCRModuloP {
private:
std::vector<uint64_t> fac;
uint64_t p;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
std::vector<uint64_t> fac;
uint64_t p;
std::vector<uint64_t> fac{};
uint64_t p = 0;


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(uint64_t size, uint64_t mod) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
NCRModuloP(uint64_t size, uint64_t mod) {
NCRModuloP(const uint64_t& size, const uint64_t& mod) {

p = mod;
fac = std::vector<uint64_t>(size);
fac[0] = 1;
for (int i = 1; i <= size; i++) {
fac[i] = (fac[i - 1] * i) % p;
}
}

/** 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(uint64_t a, uint64_t b, int64_t *x, int64_t *y) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
uint64_t gcdExtended(uint64_t a, uint64_t b, int64_t *x, int64_t *y) {
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;
uint64_t gcd = gcdExtended(b % a, a, &x1, &y1);

*x = y1 - (b / a) * x1;
*y = x1;
return gcd;
}

/** 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
*/
int64_t modInverse(uint64_t a, uint64_t m) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
int64_t modInverse(uint64_t a, uint64_t m) {
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;
}
}

/** Find nCr % p
*
* @params[in] the numbers 'n', 'r' and 'p'
* @returns the value nCr % p
*/
int64_t ncr(uint64_t n, uint64_t r, uint64_t p) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
int64_t ncr(uint64_t n, uint64_t r, uint64_t p) {
int64_t ncr(const uint64_t& n, const uint64_t& r, const uint64_t& p) {

// Base cases
if (r > n) {
return 0;
}
if (r == 1) {
return n % p;
}
if (r == 0 || r == n) {
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
return -1;
}
return (fac[n] * denominator) % p;
}
};
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
};
};
} // namespace math

} // namespace math
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
} // namespace math
} // namespace ncr_modulo_p
} // namespace math


/**
* @brief Test implementations
* @param ncrObj object which contains the precomputed factorial values and
* ncr function
* @returns void
*/
void tests(math::NCRModuloP ncrObj) {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
void tests(math::NCRModuloP ncrObj) {
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);
}

/**
* @brief Main function
* @returns 0 on exit
*/
int main() {
KaustubhDamania marked this conversation as resolved.
Show resolved Hide resolved
// populate the fac array
const uint64_t size = 1e6 + 1;
const uint64_t p = 1e9 + 7;
math::NCRModuloP ncrObj = math::NCRModuloP(size, p);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
math::NCRModuloP ncrObj = math::NCRModuloP(size, p);
math::ncr_modulo_p::NCRModuloP ncrObj = math::ncr_modulo_p::NCRModuloP(size, p);

// test 6Ci for i=0 to 7
for (int i = 0; i <= 7; i++) {
std::cout << 6 << "C" << i << " = " << ncrObj.ncr(6, i, p) << "\n";
}
tests(ncrObj); // execute the tests
std::cout << "Assertions passed\n";
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
std::cout << "Assertions passed\n";
std::cout << "Assertions passed\n";
return 0;

}