277 lines · plain
1// Copyright John Maddock 2006, 2010.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_SP_FACTORIALS_HPP7#define BOOST_MATH_SP_FACTORIALS_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <boost/math/tools/config.hpp>14#include <boost/math/tools/type_traits.hpp>15#include <boost/math/tools/precision.hpp>16#include <boost/math/policies/error_handling.hpp>17#include <boost/math/special_functions/gamma.hpp>18#include <boost/math/special_functions/detail/unchecked_factorial.hpp>19#include <boost/math/special_functions/math_fwd.hpp>20 21#ifdef _MSC_VER22#pragma warning(push) // Temporary until lexical cast fixed.23#pragma warning(disable: 4127 4701)24#endif25#ifdef _MSC_VER26#pragma warning(pop)27#endif28 29namespace boost { namespace math30{31 32template <class T, class Policy>33BOOST_MATH_GPU_ENABLED inline T factorial(unsigned i, const Policy& pol)34{35 static_assert(!boost::math::is_integral<T>::value, "Type T must not be an integral type");36 // factorial<unsigned int>(n) is not implemented37 // because it would overflow integral type T for too small n38 // to be useful. Use instead a floating-point type,39 // and convert to an unsigned type if essential, for example:40 // unsigned int nfac = static_cast<unsigned int>(factorial<double>(n));41 // See factorial documentation for more detail.42 43 BOOST_MATH_STD_USING // Aid ADL for floor.44 45 if(i <= max_factorial<T>::value)46 return unchecked_factorial<T>(i);47 T result = boost::math::tgamma(static_cast<T>(i+1), pol);48 if(result > tools::max_value<T>())49 return result; // Overflowed value! (But tgamma will have signalled the error already).50 return floor(result + 0.5f);51}52 53template <class T>54BOOST_MATH_GPU_ENABLED inline T factorial(unsigned i)55{56 return factorial<T>(i, policies::policy<>());57}58/*59// Can't have these in a policy enabled world?60template<>61inline float factorial<float>(unsigned i)62{63 if(i <= max_factorial<float>::value)64 return unchecked_factorial<float>(i);65 return tools::overflow_error<float>(BOOST_CURRENT_FUNCTION);66}67 68template<>69inline double factorial<double>(unsigned i)70{71 if(i <= max_factorial<double>::value)72 return unchecked_factorial<double>(i);73 return tools::overflow_error<double>(BOOST_CURRENT_FUNCTION);74}75*/76template <class T, class Policy>77BOOST_MATH_GPU_ENABLED T double_factorial(unsigned i, const Policy& pol)78{79 static_assert(!boost::math::is_integral<T>::value, "Type T must not be an integral type");80 BOOST_MATH_STD_USING // ADL lookup of std names81 if(i & 1)82 {83 // odd i:84 if(i < max_factorial<T>::value)85 {86 unsigned n = (i - 1) / 2;87 return ceil(unchecked_factorial<T>(i) / (ldexp(T(1), (int)n) * unchecked_factorial<T>(n)) - 0.5f);88 }89 //90 // Fallthrough: i is too large to use table lookup, try the91 // gamma function instead.92 //93 T result = boost::math::tgamma(static_cast<T>(i) / 2 + 1, pol) / sqrt(constants::pi<T>());94 if(ldexp(tools::max_value<T>(), -static_cast<int>(i+1) / 2) > result)95 return ceil(result * ldexp(T(1), static_cast<int>(i+1) / 2) - 0.5f);96 }97 else98 {99 // even i:100 unsigned n = i / 2;101 T result = factorial<T>(n, pol);102 if(ldexp(tools::max_value<T>(), -(int)n) > result)103 return result * ldexp(T(1), (int)n);104 }105 //106 // If we fall through to here then the result is infinite:107 //108 return policies::raise_overflow_error<T>("boost::math::double_factorial<%1%>(unsigned)", 0, pol);109}110 111template <class T>112BOOST_MATH_GPU_ENABLED inline T double_factorial(unsigned i)113{114 return double_factorial<T>(i, policies::policy<>());115}116 117// TODO(mborland): We do not currently have support for tgamma_delta_ratio118#ifndef BOOST_MATH_HAS_GPU_SUPPORT119 120namespace detail{121 122template <class T, class Policy>123T rising_factorial_imp(T x, int n, const Policy& pol)124{125 static_assert(!boost::math::is_integral<T>::value, "Type T must not be an integral type");126 if(x < 0)127 {128 //129 // For x less than zero, we really have a falling130 // factorial, modulo a possible change of sign.131 //132 // Note that the falling factorial isn't defined133 // for negative n, so we'll get rid of that case134 // first:135 //136 bool inv = false;137 if(n < 0)138 {139 x += n;140 n = -n;141 inv = true;142 }143 T result = ((n&1) ? -1 : 1) * falling_factorial(-x, n, pol);144 if(inv)145 result = 1 / result;146 return result;147 }148 if(n == 0)149 return 1;150 if(x == 0)151 {152 if(n < 0)153 return static_cast<T>(-boost::math::tgamma_delta_ratio(x + 1, static_cast<T>(-n), pol));154 else155 return 0;156 }157 if((x < 1) && (x + n < 0))158 {159 const auto val = static_cast<T>(boost::math::tgamma_delta_ratio(1 - x, static_cast<T>(-n), pol));160 return (n & 1) ? T(-val) : val;161 }162 //163 // We don't optimise this for small n, because164 // tgamma_delta_ratio is already optimised for that165 // use case:166 //167 return 1 / static_cast<T>(boost::math::tgamma_delta_ratio(x, static_cast<T>(n), pol));168}169 170template <class T, class Policy>171inline T falling_factorial_imp(T x, unsigned n, const Policy& pol)172{173 static_assert(!boost::math::is_integral<T>::value, "Type T must not be an integral type");174 BOOST_MATH_STD_USING // ADL of std names175 if(x == 0)176 return 0;177 if(x < 0)178 {179 //180 // For x < 0 we really have a rising factorial181 // modulo a possible change of sign:182 //183 return (n&1 ? -1 : 1) * rising_factorial(-x, n, pol);184 }185 if(n == 0)186 return 1;187 if(x < 0.5f)188 {189 //190 // 1 + x below will throw away digits, so split up calculation:191 //192 if(n > max_factorial<T>::value - 2)193 {194 // If the two end of the range are far apart we have a ratio of two very large195 // numbers, split the calculation up into two blocks:196 T t1 = x * boost::math::falling_factorial(x - 1, max_factorial<T>::value - 2, pol);197 T t2 = boost::math::falling_factorial(x - max_factorial<T>::value + 1, n - max_factorial<T>::value + 1, pol);198 if(tools::max_value<T>() / fabs(t1) < fabs(t2))199 return boost::math::sign(t1) * boost::math::sign(t2) * policies::raise_overflow_error<T>("boost::math::falling_factorial<%1%>", 0, pol);200 return t1 * t2;201 }202 return x * boost::math::falling_factorial(x - 1, n - 1, pol);203 }204 if(x <= n - 1)205 {206 //207 // x+1-n will be negative and tgamma_delta_ratio won't208 // handle it, split the product up into three parts:209 //210 T xp1 = x + 1;211 unsigned n2 = itrunc((T)floor(xp1), pol);212 if(n2 == xp1)213 return 0;214 auto result = static_cast<T>(boost::math::tgamma_delta_ratio(xp1, -static_cast<T>(n2), pol));215 x -= n2;216 result *= x;217 ++n2;218 if(n2 < n)219 result *= falling_factorial(x - 1, n - n2, pol);220 return result;221 }222 //223 // Simple case: just the ratio of two224 // (positive argument) gamma functions.225 // Note that we don't optimise this for small n,226 // because tgamma_delta_ratio is already optimised227 // for that use case:228 //229 return static_cast<T>(boost::math::tgamma_delta_ratio(x + 1, -static_cast<T>(n), pol));230}231 232} // namespace detail233 234template <class RT>235inline typename tools::promote_args<RT>::type236 falling_factorial(RT x, unsigned n)237{238 typedef typename tools::promote_args<RT>::type result_type;239 return detail::falling_factorial_imp(240 static_cast<result_type>(x), n, policies::policy<>());241}242 243template <class RT, class Policy>244inline typename tools::promote_args<RT>::type245 falling_factorial(RT x, unsigned n, const Policy& pol)246{247 typedef typename tools::promote_args<RT>::type result_type;248 return detail::falling_factorial_imp(249 static_cast<result_type>(x), n, pol);250}251 252template <class RT>253inline typename tools::promote_args<RT>::type254 rising_factorial(RT x, int n)255{256 typedef typename tools::promote_args<RT>::type result_type;257 return detail::rising_factorial_imp(258 static_cast<result_type>(x), n, policies::policy<>());259}260 261template <class RT, class Policy>262inline typename tools::promote_args<RT>::type263 rising_factorial(RT x, int n, const Policy& pol)264{265 typedef typename tools::promote_args<RT>::type result_type;266 return detail::rising_factorial_imp(267 static_cast<result_type>(x), n, pol);268}269 270#endif // BOOST_MATH_HAS_GPU_SUPPORT271 272} // namespace math273} // namespace boost274 275#endif // BOOST_MATH_SP_FACTORIALS_HPP276 277