From f3685ed245b7b55efb0e1fed77c37481af87d936 Mon Sep 17 00:00:00 2001
From: mjacobse <44684927+mjacobse@users.noreply.github.com>
Date: Fri, 9 Nov 2018 13:21:48 +0100
Subject: [PATCH] Use boost factorial for integer arguments in cpp version.

---
 cpp/CMakeLists.txt  | 1 +
 cpp/src/lambert.cpp | 9 +++++----
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt
index f6818e3..ac6ff73 100644
--- a/cpp/CMakeLists.txt
+++ b/cpp/CMakeLists.txt
@@ -4,6 +4,7 @@ project(ICATT)
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -O2")
 set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
 find_package(Eigen3 REQUIRED)
+find_package(Boost REQUIRED)
 
 add_subdirectory(../lib/dopri ${CMAKE_BINARY_DIR}/dopri)
 
diff --git a/cpp/src/lambert.cpp b/cpp/src/lambert.cpp
index c12ab7a..eb62563 100644
--- a/cpp/src/lambert.cpp
+++ b/cpp/src/lambert.cpp
@@ -3,6 +3,7 @@
 #include <stdexcept>
 #include <Eigen/Dense>
 #include <cmath>
+#include <boost/math/special_functions/factorials.hpp>
 
 using Eigen::Vector3d;
 using std::runtime_error;
@@ -22,12 +23,12 @@ namespace lambert {
             res = (cosh(sqrt(-psi)) - 1) / (-psi);
         } else {
             res = 1.0 / 2.0;
-            auto delta = (-psi) / tgamma(2 + 2 + 1);
+            auto delta = (-psi) / boost::math::factorial<double>(2 + 2 + 1);
             auto k = 1;
             while (res + delta != res) {
                 res += delta;
                 k += 1;
-                delta = pow(-psi, k) / tgamma(2*k + 2 + 1);
+                delta = pow(-psi, k) / boost::math::factorial<double>(2*k + 2 + 1);
             }
         }
         return res;
@@ -42,12 +43,12 @@ namespace lambert {
             res = (sinh(sqrt(-psi)) - sqrt(-psi)) / (-psi * sqrt(-psi));
         } else {
             res = 1.0 / 6.0;
-            auto delta = (-psi) / tgamma(2 + 3 + 1);
+            auto delta = (-psi) / boost::math::factorial<double>(2 + 3 + 1);
             int k = 1;
             while (res + delta != res) {
                 res += delta;
                 k += 1;
-                delta = pow(-psi, k) / tgamma(2*k + 3 + 1);
+                delta = pow(-psi, k) / boost::math::factorial<double>(2*k + 3 + 1);
             }
         }
         return res;