From 8013afc28cd2f372bada51bd4f1d5f57209834e0 Mon Sep 17 00:00:00 2001
From: janos <janos.meny@googlemail.com>
Date: Sat, 18 Jul 2020 16:00:48 +0200
Subject: [PATCH] Math: add binomial coefficient function.

---
 src/Magnum/CMakeLists.txt              |  2 +-
 src/Magnum/Math/Functions.cpp          | 18 +++++++++++++++++
 src/Magnum/Math/Functions.h            | 12 ++++++++++++
 src/Magnum/Math/Test/FunctionsTest.cpp | 27 ++++++++++++++++++++++++++
 4 files changed, 58 insertions(+), 1 deletion(-)

diff --git a/src/Magnum/CMakeLists.txt b/src/Magnum/CMakeLists.txt
index 49be7a307..f38bf29d2 100644
--- a/src/Magnum/CMakeLists.txt
+++ b/src/Magnum/CMakeLists.txt
@@ -111,11 +111,11 @@ set(MagnumMath_SRCS
     Math/Angle.cpp
     Math/Color.cpp
     Math/Half.cpp
-    Math/Functions.cpp
     Math/Packing.cpp
     Math/instantiation.cpp)
 
 set(MagnumMath_GracefulAssert_SRCS
+    Math/Functions.cpp
     Math/PackingBatch.cpp)
 
 # Objects shared between main and math test library
diff --git a/src/Magnum/Math/Functions.cpp b/src/Magnum/Math/Functions.cpp
index b547cc12a..a52819446 100644
--- a/src/Magnum/Math/Functions.cpp
+++ b/src/Magnum/Math/Functions.cpp
@@ -23,6 +23,8 @@
     DEALINGS IN THE SOFTWARE.
 */
 
+#include <Corrade/Utility/Assert.h>
+
 #include "Functions.h"
 
 namespace Magnum { namespace Math {
@@ -41,4 +43,20 @@ UnsignedInt log2(UnsignedInt number) {
     return log;
 }
 
+
+UnsignedLong binomialCoefficient(UnsignedInt n, UnsignedInt k) {
+    if (k > n) return 0;
+    if (k * 2 > n)
+        k = n-k;
+    if (k == 0) return 1;
+
+    UnsignedLong result = n;
+    for(UnsignedInt i = 2; i <= k; ++i ) {
+        CORRADE_ASSERT(result < ~UnsignedLong{} / (n-i+1), "Math::binomialCoefficient(): overflow for (" << Corrade::Utility::Debug::nospace << n << "choose" << k << Corrade::Utility::Debug::nospace << ")", 0ul);
+        result *= (n-i+1);
+        result /= i;
+    }
+    return result;
+}
+
 }}
diff --git a/src/Magnum/Math/Functions.h b/src/Magnum/Math/Functions.h
index 92d8e9cb1..2fe2cf1fc 100644
--- a/src/Magnum/Math/Functions.h
+++ b/src/Magnum/Math/Functions.h
@@ -410,6 +410,18 @@ template<std::size_t size, class T> inline Vector<size, T> ceil(const Vector<siz
     return out;
 }
 
+
+/**
+@brief [Binomial Coefficient](https://en.wikipedia.org/wiki/Binomial_coefficient)
+
+Returns the number of combinations of @f$ n @f$ things taken @f$ k @f$ at a time.
+The number of ways to do this is given by
+@f[
+    \begin{pmatrix} n \\ k \end{pmatrix} = \frac{n! \cdot (n-k)!}{k!} = \frac{n \cdot (n-1) \cdots (n - k + 1)}{k \cdots 1}.
+@f]
+*/
+UnsignedLong MAGNUM_EXPORT binomialCoefficient(UnsignedInt n, UnsignedInt k);
+
 /**
 @brief Floating point division remainder
 @param a     Numerator
diff --git a/src/Magnum/Math/Test/FunctionsTest.cpp b/src/Magnum/Math/Test/FunctionsTest.cpp
index 8f874e2e4..465400145 100644
--- a/src/Magnum/Math/Test/FunctionsTest.cpp
+++ b/src/Magnum/Math/Test/FunctionsTest.cpp
@@ -53,6 +53,9 @@ struct FunctionsTest: Corrade::TestSuite::Tester {
     void ceil();
     void fmod();
 
+    void binomialCoefficient();
+    void binomialCoefficientOverflow();
+
     void sqrt();
     void sqrtInverted();
     void lerp();
@@ -115,6 +118,9 @@ FunctionsTest::FunctionsTest() {
               &FunctionsTest::ceil,
               &FunctionsTest::fmod,
 
+              &FunctionsTest::binomialCoefficient,
+              &FunctionsTest::binomialCoefficientOverflow,
+
               &FunctionsTest::sqrt,
               &FunctionsTest::sqrtInverted,
               &FunctionsTest::lerp,
@@ -293,6 +299,27 @@ void FunctionsTest::ceil() {
     CORRADE_COMPARE(Math::ceil(2.7_degf), 3.0_degf);
 }
 
+void FunctionsTest::binomialCoefficient() {
+    CORRADE_COMPARE(Math::binomialCoefficient(1, 1), 1ul);
+    CORRADE_COMPARE(Math::binomialCoefficient(1, 0), 1ul);
+    CORRADE_COMPARE(Math::binomialCoefficient(19, 11), 75582ul);
+    CORRADE_COMPARE(Math::binomialCoefficient(1000, 999), 1000ul);
+    CORRADE_COMPARE(Math::binomialCoefficient(0, 0), 1ul);
+    CORRADE_COMPARE(Math::binomialCoefficient(32, 11), 129024480ul);
+    CORRADE_COMPARE(Math::binomialCoefficient(62, 31), 465428353255261088ul);
+}
+
+void FunctionsTest::binomialCoefficientOverflow() {
+    #ifdef CORRADE_NO_ASSERT
+    CORRADE_SKIP("CORRADE_NO_ASSERT defined, can't test assertions");
+    #endif
+
+    std::ostringstream out;
+    Error redirectError{&out};
+    Math::binomialCoefficient(63, 31);
+    CORRADE_COMPARE(out.str(), "Math::binomialCoefficient(): overflow for (63 choose 31)\n");
+}
+
 void FunctionsTest::fmod() {
     CORRADE_COMPARE(Math::fmod(5.1f, 3.0f), 2.1f);
     CORRADE_COMPARE(Math::fmod(Vector3(5.1f, -5.1f, 6.8f), Vector3(3.0f, 3.0f, 1.1f)), Vector3(2.1f, -2.1f, 0.2f));