[!471] add dynamic and static methods for factorials and binomials

Merge branch 'feature/dynamic-static-factorial-binomial' into 'master'

ref:core/dune-common Added some dynamic and static methods for calculating
factorials and binomials. There are several different implementation in
several modules of varying accessibility and age which can be replaced by
this.

See merge request [core/dune-common!471]

  [core/dune-common!471]: gitlab.dune-project.org/core/dune-common/merge_requests/471

dune/common/math.hh

@@ -9,6 +9,7 @@
 
 #include <cmath>
 #include <complex>
+#include <limits>
 #include <type_traits>
 
 #include <dune/common/typeutilities.hh>
@@ -80,6 +81,60 @@ namespace Dune
     enum { factorial = 1 };
   };
 
+
+  //! calculate the factorial of n as a constexpr
+  // T has to be an integral type
+  template<class T>
+  constexpr inline static T factorial(const T& n) noexcept
+  {
+    static_assert(std::numeric_limits<T>::is_integer, "`factorial(n)` has to be called with an integer type.");
+    T fac = 1;
+    for(T k = 0; k < n; ++k)
+      fac *= k+1;
+    return fac;
+  }
+
+  //! calculate the factorial of n as a constexpr
+  template<class T, T n>
+  constexpr inline static auto factorial (std::integral_constant<T, n>) noexcept
+  {
+    return std::integral_constant<T, factorial(n)>{};
+  }
+
+
+  //! calculate the binomial coefficient n over k as a constexpr
+  // T has to be an integral type
+  template<class T>
+  constexpr inline static T binomial (const T& n, const T& k) noexcept
+  {
+    static_assert(std::numeric_limits<T>::is_integer, "`binomial(n, k)` has to be called with an integer type.");
+
+    if( k < 0 || k > n )
+      return 0;
+
+    if (2*k > n)
+      return binomial(n, n-k);
+
+    T bin = 1;
+    for(auto i = n-k; i < n; ++i)
+      bin *= i+1;
+    return bin / factorial(k);
+  }
+
+  //! calculate the binomial coefficient n over k as a constexpr
+  template<class T, T n, T k>
+  constexpr inline static auto binomial (std::integral_constant<T, n>, std::integral_constant<T, k>) noexcept
+  {
+    return std::integral_constant<T, binomial(n, k)>{};
+  }
+
+  template<class T, T n>
+  constexpr inline static auto binomial (std::integral_constant<T, n>, std::integral_constant<T, n>) noexcept
+  {
+    return std::integral_constant<T, (n >= 0 ? 1 : 0)>{};
+  }
+
+
   //! compute conjugate complex of x
   // conjugate complex does nothing for non-complex types
   template<class K>

dune/common/test/CMakeLists.txt

@@ -364,6 +364,11 @@ dune_add_test(SOURCES versiontest.cc
               LINK_LIBRARIES dunecommon
               LABELS quick)
 
+dune_add_test(SOURCES mathtest.cc
+              LINK_LIBRARIES dunecommon
+              LABELS quick)
+
+
 install(
   FILES
   arithmetictestsuite.hh

dune/common/test/mathtest.cc

+#include "config.h"
+
+#include <iostream>
+
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/indices.hh>
+#include <dune/common/math.hh>
+#include <dune/common/test/testsuite.hh>
+
+#include <dune/common/math.hh>
+
+
+using namespace Dune::Hybrid;
+using namespace Dune::Indices;
+using Dune::TestSuite;
+
+template<class T, T n>
+constexpr inline static auto next(std::integral_constant<T, n>)
+  -> std::integral_constant<T, n+1>
+{
+  return {};
+}
+
+template<class T, T k>
+auto testStaticFactorial (std::integral_constant<T, k> _k = {}) -> TestSuite
+{
+  TestSuite t;
+
+  std::cout << "test static factorial\n{";
+
+  forEach(integralRange(_k), [&](auto _i) {
+    auto value = Dune::factorial(_i);
+
+    t.check(decltype(value)::value == Dune::Factorial<decltype(_i)::value>::factorial);
+
+    std::cout<< ' ' << value() << ',';
+  });
+
+  std::cout << "};\n\n";
+
+  return t;
+}
+
+template<class T, T k>
+auto testStaticBinomial (std::integral_constant<T, k> _k = {}) -> TestSuite
+{
+  TestSuite t;
+
+  std::cout << "test static binomial\n";
+
+  forEach(integralRange(_k), [&](auto _i) {
+    std::cout << "{";
+    forEach(integralRange(next(_i)), [&](auto _j) {
+      const auto value = Dune::binomial(_i,_j);
+
+      auto control = Dune::Factorial<decltype(_i)::value>::factorial
+                     / Dune::Factorial<decltype(_j)::value>::factorial
+                     / Dune::Factorial<decltype(_i)::value - decltype(_j)::value>::factorial;
+      t.check(decltype(value)::value == control);
+
+      std::cout<< ' ' << value() << ',';
+    });
+
+    std::cout << "};\n";
+  });
+
+  std::cout << "\n";
+
+  return t;
+}
+
+int main(int argc, char** argv)
+{
+  TestSuite t;
+
+  t.subTest(testStaticFactorial(_5));
+  t.subTest(testStaticBinomial(_5));
+
+  return t.exit();
+}