add constexpr methods for factorials and binomials

dune/common/math.hh

@@ -80,6 +80,55 @@ namespace Dune
     enum { factorial = 1 };
   };
 
+
+  //! calculate the factorial of n as a constexpr
+  // T has to be an integral type
+  template<class T,
+           std::enable_if_t<std::is_integral<std::decay_t<T>>::value, int> = 0>
+  constexpr inline static auto factorial(T&& n) noexcept
+    -> std::decay_t< T >
+  {
+    std::decay_t<T> fac = 1;
+    for(std::decay_t<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 and U have to be the same integral type
+  template<class T, class U,
+           std::enable_if_t<std::is_same<std::decay_t<T>, std::decay_t<U>>::value, int> = 0,
+           std::enable_if_t<std::is_integral<std::decay_t<T>>::value, int> = 0>
+  constexpr inline static auto binomial (T&& n, U&& k) noexcept
+    -> std::decay_t< T >
+  {
+    if( k < 0 || k > n )
+      return 0;
+
+    if (2*k > n)
+      return binomial(n, n-k);
+
+    std::decay_t<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)>{};
+  }
+
   //! 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();
+}