Parallel computation of multiple-precision binomal coefficients C(a,b)

comb.cpp

#include <iostream>
#include <gmp.h>
#include <gmpxx.h>
#include <omp.h>

int main(int argc, char *argv[])
{
    if (argc != 1 + 2) {
        std::cerr << "Specify a and b for comb(a, b)\n";
        return 1;
    }

    mpz_class a, b;
    try {
        a = argv[1];
        b = argv[2];
    } catch (std::invalid_argument &e) {
        std::cerr << e.what() << ": Invalid argument\n";
        return 1;
    }

    if (a < 0 || b < 0) {
        std::cerr << "Invalid a and/or b\n";
        return 1;
    }
    if (b > a / 2)
        b = a - b;

    mpz_class num_prod_all = 1, den_prod_all = 1;

#pragma omp parallel shared(num_prod_all, den_prod_all) firstprivate(a, b)
    {
        const unsigned thread_num = omp_get_thread_num();
        const unsigned num_threads = omp_get_num_threads();

        mpz_class num = a - b + 1 + thread_num;
        mpz_class num_prod = 1, den_prod = 1;

        for (mpz_class i = 1 + thread_num; i <= b;
                i += num_threads, num += num_threads) {
            num_prod *= num;
            den_prod *= i;
        }

#pragma omp critical
        num_prod_all *= num_prod;
#pragma omp critical
        den_prod_all *= den_prod;
    }

    mpz_class c;
    mpz_divexact(c.get_mpz_t(),
            num_prod_all.get_mpz_t(), den_prod_all.get_mpz_t());

    std::cout << "comb(" << argv[1] << ", " << argv[2] << ") =\n";
    if (b != mpz_class(argv[2]))
        std::cout << "comb(" << a.get_str() << ", " << b.get_str() << ") =\n";
    std::cout << c.get_str() << "\n";

    return 0;
}