Permutation function for a math library in C++

Question

Continuing the code for combinatorics library, adding the definition of ⁿPᵣ. This is an extract from a larger codebase.

1. Is there any way to limit the code of gcd() to this file? Should I convert the namespace Combinatorics to a class?
2. Any other advice you would want to give me?

Code

// MathLibrary.cpp : Defines the functions for the static library.
//
#include "pch.h"
#include "framework.h"
#include <vector>
#include <iostream>
#include <limits>
#include <cstddef>

namespace Combinatorics
{
    const std::vector<unsigned long long> representableFactors = {1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800,
        479001600,6227020800,87178291200,1307674368000,20922789888000,355687428096000,6402373705728000,121645100408832000,2432902008176640000};

    unsigned long long factorial(unsigned int n) 
    {
        if (n >= representableFactors.size())
        {
            throw std::invalid_argument("Combinatorics:factorial - The argument is too large");
        }
        return representableFactors[n];
    }

    inline long long gcd(long long a, long long b)
    {
        if (a % b == 0)
        {
            return b;
        }
        return gcd(b, (a % b));
    }

    long long combinations(int n, int r)
    {
        if (n < 0 || r < 0)
        {
            throw std::invalid_argument("Combinatorics::combinations - N and R cannot be less than 0");
        }
        if (n < r)
        {
            throw std::invalid_argument("Combinatorics::combinations - The value of r cannot be greater than n");
        }
        if (n - r < r)
        {
            r = n - r;
        }
        long long int denominatorProduct = 1;
        long long int numeratorProduct = 1;
        for (long long int denomCount = r, numCount = n ; denomCount >= 1; denomCount--, numCount--)
        {
            //TODO : Convert to limits
            if((LLONG_MAX / n) < numeratorProduct || (LLONG_MAX / n) < denominatorProduct)
            {
                throw std::invalid_argument("Combinatorics::combinations - Overflow detected aborting");
            }
            denominatorProduct *= denomCount;
            numeratorProduct *=  numCount;
            long long gcdCommonFactor = gcd(denominatorProduct, numeratorProduct);
            denominatorProduct /= gcdCommonFactor;
            numeratorProduct /= gcdCommonFactor;
        }
        return (numeratorProduct / denominatorProduct);
    }
    //Code to be reviewed.
    long long permutations(int n, int r)
    {
        if (n < 0 || r < 0)
        {
            throw std::invalid_argument("Combinatorics::permutations - N and R cannot be less than 0");
        }
        if (n < r)
        {
            throw std::invalid_argument("Combinatorics::permutations - The value of r cannot be greater than n");
        }
        long long permutations = 1;
        for (long long int numeratorCount = n; numeratorCount >= (n - r + 1); numeratorCount--)
        {
            //TODO: Convert to limits
            if ((LLONG_MAX / n) < permutations)
            {
                throw std::invalid_argument("Combinatorics::permutations - Overflow detected aborting");
            }
            permutations *= numeratorCount;
        }
        return permutations;
    }
}

Test code

#include "pch.h"
#include <iostream>
#include "../MathLibrary/MathLibrary.cpp"

TEST(Combinatorial_Factorial, small_ints)
{
    EXPECT_EQ(Combinatorics::factorial(0), 1);
    EXPECT_EQ(Combinatorics::factorial(1), 1);
    EXPECT_EQ(Combinatorics::factorial(5), 120);
    EXPECT_EQ(Combinatorics::factorial(20), 2432902008176640000);
}

TEST(Combinatorial_Factorial, too_big)
{
    EXPECT_THROW(Combinatorics::factorial(500), std::invalid_argument);
}

TEST(Combinatorial_Combinations, small_ints)
{
    EXPECT_EQ(Combinatorics::combinations(5,5), 1);
    EXPECT_EQ(Combinatorics::combinations(5, 0), 1);
    EXPECT_EQ(Combinatorics::combinations(5, 1), 5);
    EXPECT_EQ(Combinatorics::combinations(20,10),184756);
    EXPECT_EQ(Combinatorics::combinations(40, 35),658008);
}

TEST(Combinatorial_Combinations, negative_n)
{
    EXPECT_THROW(Combinatorics::combinations(-5, 5), std::invalid_argument);
}

TEST(Combinatorial_Combinations, r_greater_than_n)
{
    EXPECT_THROW(Combinatorics::combinations(4, 5), std::invalid_argument);
}

TEST(Combinatorial_Combinations, overflow)
{
    EXPECT_THROW(Combinatorics::combinations(100, 50), std::invalid_argument);
}

TEST(Combinatorial_Permutations, small_ints)
{
    EXPECT_EQ(Combinatorics::permutations(5, 5), 120);
    EXPECT_EQ(Combinatorics::permutations(5, 0), 1);
    EXPECT_EQ(Combinatorics::permutations(5, 2), 20);
    EXPECT_EQ(Combinatorics::permutations(10, 2), 90);
    EXPECT_EQ(Combinatorics::permutations(40, 3), 59280);
    EXPECT_EQ(Combinatorics::permutations(40, 7), 93963542400);
    EXPECT_EQ(Combinatorics::permutations(50, 4), 5527200);
}

TEST(Combinatorial_Permutations, r_negative)
{
    EXPECT_THROW(Combinatorics::permutations(5, -5), std::invalid_argument);
}

TEST(Combinatorial_Permutations, n_negative)
{
    EXPECT_THROW(Combinatorics::permutations(-5, 5), std::invalid_argument);
}

TEST(Combinatorial_Permutations,r_greater)
{
    EXPECT_THROW(Combinatorics::permutations(5, 6), std::invalid_argument);
}

TEST(Combinatorial_Permutations,overflow)
{
    EXPECT_THROW(Combinatorics::permutations(50,46), std::invalid_argument);
}

Answer 1

Unnecessary includes:

#include "pch.h"  // in both files
#include "framework.h"
#include <iostream>  // in both files
#include <limits>

Missing includes:

// in implementation
#include <climits>
#include <stdexcept>

// in tests
#include <gtest/gtest.h>

---

There is no need to call gcd() in the implementation of combinations(), or to keep separate track of denominator. Consider the expansion of ⁿCᵣ:

\$ ^nC_r = \frac{(r+1) ✕ (r+2) ✕ (r+3) ✕ ... ✕ n}{1✕2✕3✕...✕(n-r)} \$

At any time, consider the partial result formed by taking i terms from each of the numerator and denominator. The denominator part is obviously just i!; we can show that the numerator part must be an exact multiple of i! because it must contain at least one multiple of each of 0,...,i.

So if we change our loop to count upwards in the denominator, we can show that we always end up with a denominator of 1 (see the assert()):

for (long long int denomCount = 1, numCount = n ; denomCount <= r; ++denomCount, --numCount)
{
    //TODO : Convert to limits
    if((LLONG_MAX / n) < numeratorProduct || (LLONG_MAX / n) < denominatorProduct)
    {
        throw std::invalid_argument("Combinatorics::combinations - Overflow detected aborting");
    }
    denominatorProduct *= denomCount;
    numeratorProduct *=  numCount;
    long long gcdCommonFactor = gcd(denominatorProduct, numeratorProduct);
    denominatorProduct /= gcdCommonFactor;
    numeratorProduct /= gcdCommonFactor;
    assert(denominatorProduct == 1);
}

That observation allows us to remove denominatorProduct altogether:

long long int product = 1;
for (long long int denomCount = 1, numCount = n;  denomCount <= r;  ++denomCount, --numCount)
{
    if (LLONG_MAX / n < product)
    {
        throw std::range_error("Combinatorics::permutations - overflow detected");
    }
    product *= numCount;
    product /= denomCount;
}
return product;

---

Counting down in permutations() means that this test might catch some false-positives:

    for (long long int numeratorCount = n; numeratorCount >= (n - r + 1); numeratorCount--)
    {
        //TODO: Convert to limits
        if ((LLONG_MAX / n) < permutations)

If we count upwards, then the final multiplication is by n, and the test is exact.

The same problem is present in combinations().

---

Since the combinatoric functions are undefined for negative values, we should be throwing std::domain_error to represent that, or simply accept unsigned types (enable your compiler's "signed conversion" warnings to avoid problems with implicit conversions).

And we should return an unsigned type too, giving us a little extra range to represent the larger values.

---

To answer your specific question, the way to make gcd() internal to the translation unit is to declare it in the anonymous namespace (in the implementation file, not in a header!) or to use the old C-style static keyword.

But we should be using std::gcd() instead of reinventing that wheel. (If we had to reimplement it, I'd prefer to see it written iteratively rather than recursively, since many C++ compilers don't eliminate tail-calls).

---

It's probably simpler to use a std::array for the table of factorials, and that should be a local static constant within the factorial() function - it doesn't need to be visible anywhere else. I'm not convinced that a table is a good idea - it inhibits changing this to a template function (which might support larger result types).