Calculation of binomial coefficients ... still programmatically

Earlier, in three articles, the topic of calculating binomial coefficients was raised.

Calculation of binomial coefficients in C (C ++)
Calculation of binomial coefficients using Fourier transforms
Calculating the binomial coefficients ... manually

In the last article, the author demonstrated a purely mathematical optimization of the calculation of binomial coefficients. It turned out that for their calculation and the computer is not particularly needed. Enough paper, pens, calculator and beer . However, if one of the above things is not at hand, but an idle computer is in the area of ​​reach, you can do the same on the computer.

Algorithm

In the comments on the article, the author was blamed for the lack of a formalized algorithm. After some thought, the following was drawn.
')
When calculating the binomial coefficient after reducing the factorial formula to the form:

C (n, k) = P (n, k + 1) / P (1, nk)

We get a fraction that needs to be reduced. What is the easiest way to reduce fractions in integers? Spread out the numerator and denominator into prime factors. The matter is facilitated by the fact that in this case, both the numerator and the denominator are products of small numbers. Those. each factor can be decomposed into prime factors separately, accumulating the result in the form of a list of prime numbers with corresponding exponents.

Then, we go through the list obtained for the denominator, take each number and from the exponent for that number in the numerator we subtract the exponent for this number in the denominator. Based on the fact that the binomial coefficient is an integer number, the denominator should be reduced completely, that is, become equal to 1.

The last step: we go through the list obtained for the numerator in the previous step, and multiply all the numbers raised to the power corresponding to this number.

Implementation

It was decided to store the “prime number - exponent” pairs in the map container, where the prime number is the key and the exponent is the value.

For the factorization of numbers was taken the most primitive algorithm. To optimize future use, an external container is passed to the factorization function by reference.

#include <iostream> #include <map> using namespace std; typedef unsigned long long ullong; //      . //       //    map,      : //  -  ,  -     . template <typename integral_type> void factorize(integral_type num, map<integral_type, unsigned int>& res) { integral_type j = 2; while (num / integral_type(2) >= j) { if (num % j == 0) { res[j]++; num /= j; j = integral_type(2); } else { ++j; } } res[num]++; } typedef map<ullong, unsigned int> factmap; //    . ullong C(int n, int k) { ullong result = 1uLL; //        // C(n, k) = (n, k+1) / (1, nk) int lim_numer_min = k + 1; //       int lim_numer_max = n; //       int lim_denom_min = 2; //       int lim_denom_max = n - k; //       //          factmap numerator, denominator; //        for (int i = lim_numer_min; i <= lim_numer_max; ++i) { factorize(ullong(i), numerator); } //        for (int i = lim_denom_min; i <= lim_denom_max; ++i) { factorize(ullong(i), denominator); } //      for (auto i = denominator.begin(); i != denominator.end(); i++) { numerator[i->first] -= i->second; } //      //     for (auto i = numerator.begin(); i != numerator.end(); i++) { //        for (ullong p = 0; p < i->second; ++p) { result *= i->first; } } return result; } 

Testing

The head program was made like the program from the first article :

 int main() { int n; for (n = 34; n < 68; ++n) { cout << "C(" << n << ", " << n/2 << ") = " << C(n, n/2) << endl; } return 0; }