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

---
 src/Magnum/Math/Functions.h            | 28 ++++++++++++++++++++++++++
 src/Magnum/Math/Test/FunctionsTest.cpp | 13 ++++++++++++
 2 files changed, 41 insertions(+)

diff --git a/src/Magnum/Math/Functions.h b/src/Magnum/Math/Functions.h
index 92d8e9cb1..754691183 100644
--- a/src/Magnum/Math/Functions.h
+++ b/src/Magnum/Math/Functions.h
@@ -34,6 +34,7 @@
 #include <cstdlib> /* std::div() */
 #include <type_traits>
 #include <utility>
+#include <limits>
 #include <Corrade/Utility/StlMath.h>
 
 #include "Magnum/visibility.h"
@@ -410,6 +411,33 @@ template<std::size_t size, class T> inline Vector<size, T> ceil(const Vector<siz
     return out;
 }
 
+
+/**
+@brief Computes a Binomial Coefficient
+
+Given n elements computes the number of ways to choose k elements
+without replacement from n elements. The number of ways to do this is
+given by
+@f[
+    \frac{n}{k} = \frac{n! \cdot (n-k)!}{k!} = \frac{n \cdot (n-1) \cdots (n - k + 1)}{k \cdots 1}.
+@f]
+*/
+template<class T> inline typename std::enable_if<IsIntegral<T>::value, T>::type binomialCoefficient(T n, T k) {
+    CORRADE_ASSERT(k >= 0, "k must be positive to compute the Binomial Coefficient", T{0});
+    if (k > n) return 0;
+    if (k * 2 > n)
+        k = n-k;
+    if (k == 0) return 1;
+
+    T result = n;
+    for(T i = 2; i <= k; ++i ) {
+        CORRADE_ASSERT(result < std::numeric_limits<T>::max()/(n-i+1), "Overflow In Binomial Coefficient", T{0});
+        result *= (n-i+1);
+        result /= i;
+    }
+    return result;
+}
+
 /**
 @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..2718f183f 100644
--- a/src/Magnum/Math/Test/FunctionsTest.cpp
+++ b/src/Magnum/Math/Test/FunctionsTest.cpp
@@ -53,6 +53,8 @@ struct FunctionsTest: Corrade::TestSuite::Tester {
     void ceil();
     void fmod();
 
+    void binomialCoefficient();
+
     void sqrt();
     void sqrtInverted();
     void lerp();
@@ -115,6 +117,8 @@ FunctionsTest::FunctionsTest() {
               &FunctionsTest::ceil,
               &FunctionsTest::fmod,
 
+              &FunctionsTest::binomialCoefficient,
+
               &FunctionsTest::sqrt,
               &FunctionsTest::sqrtInverted,
               &FunctionsTest::lerp,
@@ -293,6 +297,15 @@ void FunctionsTest::ceil() {
     CORRADE_COMPARE(Math::ceil(2.7_degf), 3.0_degf);
 }
 
+void FunctionsTest::binomialCoefficient() {
+    CORRADE_COMPARE(Math::binomialCoefficient(1, 1), 1);
+    CORRADE_COMPARE(Math::binomialCoefficient(1, 0), 1);
+    CORRADE_COMPARE(Math::binomialCoefficient(19, 11), 75582);
+    CORRADE_COMPARE(Math::binomialCoefficient(1000, 999), 1000);
+    CORRADE_COMPARE(Math::binomialCoefficient(0, 0), 1);
+    CORRADE_COMPARE(Math::binomialCoefficient(32, 11), 129024480);
+}
+
 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));