Skip to content

Commit

Permalink
factorial and lfactorial don't assert anymore
Browse files Browse the repository at this point in the history
  • Loading branch information
@jtlap
jtlap authored Feb 19, 2025
1 parent 2d6b1fe commit 8323b82
Show file tree
Hide file tree
Showing 8 changed files with 86 additions and 50 deletions.
50 changes: 35 additions & 15 deletions include/eve/module/special/regular/factorial.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,11 @@
#include <eve/arch.hpp>
#include <eve/traits/overload.hpp>
#include <eve/module/core/decorator/core.hpp>

#include <eve/module/special/regular/tgamma.hpp>
namespace eve
{
template<typename Options>
struct factorial_t : elementwise_callable<factorial_t, Options>
struct factorial_t : elementwise_callable<factorial_t, Options, pedantic_option>
{
template<eve::integral_value T>
EVE_FORCEINLINE constexpr
Expand Down Expand Up @@ -45,25 +45,30 @@ namespace eve
//! namespace eve
//! {
//! // Regular overload
//! template <value T> constexpr as_wide_as_t<double,T> factorial(T x) noexcept; // 1
//! template <value T> constexpr as_wide_as_t<double,T> factorial(T x) noexcept; // 1
//!
//! // Semantic options
//! template <value T> constexpr as_wide_as_t<double,T> factorial[pedantic](T x) noexcept; // 2
//!
//! // Lanes masking
//! constexpr auto factorial[conditional_expr auto c](value auto n) noexcept; // 2
//! constexpr auto factorial[logical_value auto m](value auto n) noexcept; // 2
//! constexpr auto factorial[conditional_expr auto c](value auto n) noexcept; // 3
//! constexpr auto factorial[logical_value auto m](value auto n) noexcept; // 3
//! }
//! @endcode
//!
//! **Parameters**
//!
//! * `n`: must be integral or flint.
//! * `n`: input value.
//! * `c`: [Conditional expression](@ref eve::conditional_expr) masking the operation.
//! * `m`: [Logical value](@ref eve::logical_value) masking the operation.
//!
//! **Return value**
//!
//! 1. The value of \f$ n!\f$ is returned. If the entry is of integral type a double based floating_value
//! is returned.
//! 2. [The operation is performed conditionnaly](@ref conditional)
//! is returned. The call return a Nan for any entry which is not a flint or not positive.
//! 3. With the pedantic option \f$\Gamma(x+1)\f$ is returned. (more expansive)
//! 4. [The operation is performed conditionnaly](@ref conditional)
//!
//!
//! @warning
//! This function will overflow as soon as the input is greater than 171 for integral or double
Expand All @@ -90,8 +95,14 @@ namespace eve
{
if constexpr(signed_integral_value<T>)
{
EVE_ASSERT(eve::all(is_gez(n)), "factorial : some entry elements are not positive");
return factorial(convert(n, uint_from<T>()));
if constexpr(O::contains(pedantic))
return tgamma(inc(convert(n, as<double>())));
else
{
auto ltzn = is_ltz(n);
auto nn = if_else(ltzn, zero, convert(eve::abs(n), uint_from<T>()));
return if_else(ltzn, allbits, factorial(nn));
}
}
else
{
Expand Down Expand Up @@ -281,11 +292,20 @@ namespace eve
T factorial_(EVE_REQUIRES(cpu_), O const&, T n) noexcept
{
using elt_t = element_type_t<T>;
EVE_ASSERT(eve::all(is_flint(n)), "factorial : some entry elements are not flint");
EVE_ASSERT(eve::all(is_gez(n)), "factorial : some entry elements are not positive");
auto r = factorial(convert(n, uint_from<T>()));
if constexpr( std::same_as<elt_t, double> ) return r;
else return convert(r, as<float>());
if constexpr(O::contains(pedantic))
{
auto r = tgamma(inc(convert(n, as<double>())));
if constexpr( std::same_as<elt_t, double> ) return r;
else return convert(r, as<float>());
}
else
{
auto bad = is_not_flint(n) || is_ltz(n) || n > 171;
auto nn = if_else(bad, zero, convert(eve::min(eve::abs(n), T(172)), uint_from<T>())); //never ub
auto r = if_else(bad, allbits, factorial(nn));
if constexpr( std::same_as<elt_t, double> ) return r;
else return convert(r, as<float>());
}
}
}
}
44 changes: 24 additions & 20 deletions include/eve/module/special/regular/lfactorial.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,11 +11,12 @@
#include <eve/module/core/decorator/core.hpp>
#include <eve/module/special/regular/factorial.hpp>
#include <eve/module/special/regular/log_abs_gamma.hpp>
#include <eve/module/special/regular/log_gamma.hpp>

namespace eve
{
template<typename Options>
struct lfactorial_t : elementwise_callable<lfactorial_t, Options>
struct lfactorial_t : elementwise_callable<lfactorial_t, Options, pedantic_option>
{
template<eve::integral_value T>
EVE_FORCEINLINE constexpr
Expand Down Expand Up @@ -50,9 +51,12 @@ namespace eve
//! // Regular overload
//! template <value T> constexpr as_wide_as_t<double,T> lfactorial(T x) noexcept; // 1
//!
//! // Semantic options
//! template <value T> constexpr as_wide_as_t<double,T> lfactorial[pedantic](T x) noexcept; // 2
//!
//! // Lanes masking
//! constexpr auto factorial[conditional_expr auto c](value auto n) noexcept; // 2
//! constexpr auto factorial[logical_value auto m](value auto n) noexcept; // 2
//! constexpr auto factorial[conditional_expr auto c](value auto n) noexcept; // 3
//! constexpr auto factorial[logical_value auto m](value auto n) noexcept; // 3
//! }
//! @endcode
//!
Expand All @@ -67,8 +71,9 @@ namespace eve
//! element type is always double to try to avoid overflow as possible.
//! * If the entry is a floating point value which must be a flint,
//! the result is of the same type as the entry.
//! * If `n` elements are nor integer nor flint the result is undefined.
//! 2. [The operation is performed conditionnaly](@ref conditional)
//! * If `n` elements are nor integer nor flint the result is NaN.
//! 2. With the pedantic option \f$\log(\Gamma(x+1))\f$ is returned.
//! 3. [The operation is performed conditionnaly](@ref conditional)
//!
//! @groupheader{External references}
//! * [Wolfram MathWorld: Erf](https://mathworld.wolfram.com/Factorial.html)
Expand All @@ -86,24 +91,23 @@ namespace eve
{
template<typename T, callable_options O>
constexpr EVE_FORCEINLINE decltype(eve::factorial(T()))
lfactorial_(EVE_REQUIRES(cpu_), O const&, T n) noexcept
lfactorial_(EVE_REQUIRES(cpu_), O const& , T n) noexcept
{
EVE_ASSERT(eve::all(is_flint(n)), "lfactorial : some entry elements are not flint");
EVE_ASSERT(eve::all(is_gez(n)), "lfactorial : some entry elements are not positive");
constexpr auto max = std::same_as<element_type_t<T>, double> ? 171 : 35;
auto r = eve::log(factorial(n));

if( eve::all(n < max) ) return r;
if constexpr(O::contains(pedantic))
{
using elt_t = element_type_t<T>;
auto r = eve::log_gamma(eve::convert(n, as<double>()));
if constexpr( std::same_as<elt_t, double> ) return r;
else return convert(r, as<float>());
}
else
{
auto np = [](auto x)
{
if constexpr(integral_value<T>) return convert(inc(x), as<double>());
else return inc(x);
}(n);

return if_else(n < max, r, log_abs_gamma(np));

auto bad = is_not_flint(n) || is_ltz(n);
auto np = [](auto x) {
if constexpr(integral_value<T>) return convert(inc(x), as<double>());
else return inc(x);
}(n);
return if_else(bad, allbits, log_abs_gamma(np));
}
}
}
Expand Down
2 changes: 1 addition & 1 deletion include/eve/module/special/regular/log_gamma.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ namespace eve
constexpr T
log_gamma_(EVE_REQUIRES(cpu_), O const&, T a0) noexcept
{
auto aa0 = if_else(a0 == minf(as(a0)) || is_lez(signgam(a0)), allbits, a0);
auto aa0 = if_else(a0 == minf(as(a0)) || is_lez(signgam(a0)) ||(is_ltz(a0) && is_flint(a0)), allbits, a0);
return log_abs_gamma(aa0);
}
}
Expand Down
5 changes: 5 additions & 0 deletions test/doc/special/factorial.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,14 +5,19 @@
int main()
{
eve::wide<float> wf([](auto i, auto c)->float{ return 2*(i+c/2);});
eve::wide wf1{-3.0, 0.5, 3.0, 4.0};
eve::wide wi = {93, 25, 32, 180, 1, 2, 3, 4};
eve::wide<std::uint32_t> wu([](auto i, auto )->std::uint32_t{ return i;});

std::cout << "<- wf = " << wf << "\n";
std::cout << "<- wf1 = " << wf1<< "\n";
std::cout << "<- wi = " << wi << "\n";
std::cout << "<- wu = " << wu << "\n";

std::cout << "-> factorial(wf) = " << eve::factorial(wf) << "\n";
std::cout << "-> factorial(wi) = " << eve::factorial(wi) << "\n";
std::cout << "-> factorial(wu) = " << eve::factorial(wu) << "\n";
std::cout << "-> factorial(wf1) = " << eve::factorial(wf1) << "\n";
std::cout << "-> factorial[pedantic](wf1) = " << eve::factorial[eve::pedantic](wf1) << "\n";

}
4 changes: 4 additions & 0 deletions test/doc/special/lfactorial.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,14 +5,18 @@
int main()
{
eve::wide<float> wf([](auto i, auto c)->float{ return 2*(i+c/2);});
eve::wide wf1{-3.0, 0.5, 3.0, 4.0};
eve::wide wi = {93, 25, 32, 180, 1, 2, 3, 4};
eve::wide<std::uint32_t> wu([](auto i, auto )->std::uint32_t{ return i;});

std::cout << "<- wf = " << wf << "\n";
std::cout << "<- wf1 = " << wf1 << "\n";
std::cout << "<- wi = " << wi << "\n";
std::cout << "<- wu = " << wu << "\n";

std::cout << "-> lfactorial(wf) = " << eve::lfactorial(wf) << "\n";
std::cout << "-> lfactorial(wi) = " << eve::lfactorial(wi) << "\n";
std::cout << "-> lfactorial(wu) = " << eve::lfactorial(wu) << "\n";
std::cout << "-> lfactorial(wf1) = " << eve::lfactorial(wf1) << "\n";
std::cout << "-> lfactorial[pedantic](wf1) = " << eve::lfactorial[eve::pedantic](wf1) << "\n";
}
8 changes: 5 additions & 3 deletions test/doc/special/log_gamma.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,12 @@
int main()
{
eve::wide<float> wf([](auto i, auto c)->float{ return 2*(i+c/2);});

std::cout << "<- wf = " << wf << "\n";
eve::wide bads{0.0, -0.0, -1.0, -2.5};
std::cout << "<- wf = " << wf << "\n";
std::cout << "<- bads = " << bads << "\n";

std::cout << "-> log_gamma(wf) = " << eve::log_gamma(wf) << "\n";
std::cout << "-> log_gamma[ignore_last(2)](wf)= " << eve::log_gamma[eve::ignore_last(2)](wf) << "\n";
std::cout << "-> log_gamma[wf != -2.0f](wf) = " << eve::log_gamma[wf != -2.0f](wf) << "\n";
std::cout << "-> log_gamma[wf != 12.0f](wf) = " << eve::log_gamma[wf != 12.0f](wf) << "\n";
std::cout << "-> log_gamma(bads) = " << eve::log_gamma(bads)<< "\n";
}
20 changes: 10 additions & 10 deletions test/unit/module/special/lfactorial.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,18 +45,18 @@ TTS_CASE_TPL("Check corner-cases behavior of eve::lfactorial on wide", eve::test
if constexpr( eve::integral_value<T> )
{
using v_t = eve::element_type_t<T>;
TTS_EQUAL(eve::lfactorial(T(10)), eve::log(d_t(3628800)));
TTS_EQUAL(eve::lfactorial(T(5)), eve::log(d_t(120)));
TTS_EQUAL(eve::lfactorial(T(0)), eve::log(d_t(1)));
TTS_EQUAL(eve::lfactorial(T(1)), eve::log(d_t(1)));
if constexpr( sizeof(v_t) > 1 ) TTS_EQUAL(eve::lfactorial(T(200)), eve::log_abs_gamma(d_t(201)));
TTS_ULP_EQUAL(eve::lfactorial(T(10)), eve::log(d_t(3628800)), 1.0);
TTS_ULP_EQUAL(eve::lfactorial(T(5)), eve::log(d_t(120)), 1.0);
TTS_ULP_EQUAL(eve::lfactorial(T(0)), eve::log(d_t(1)), 1.0);
TTS_ULP_EQUAL(eve::lfactorial(T(1)), eve::log(d_t(1)), 1.0);
if constexpr( sizeof(v_t) > 1 ) TTS_ULP_EQUAL(eve::lfactorial(T(200)), eve::log_abs_gamma(d_t(201)), 1.0);
}
else
{
TTS_EQUAL(eve::lfactorial(T(10)), eve::log(T(3628800)));
TTS_EQUAL(eve::lfactorial(T(5)), eve::log(T(120)));
TTS_EQUAL(eve::lfactorial(T(0)), eve::log(T(1)));
TTS_EQUAL(eve::lfactorial(T(1)), eve::log(T(1)));
TTS_EQUAL(eve::lfactorial(T(200)), eve::log_abs_gamma(T(201)));
TTS_ULP_EQUAL(eve::lfactorial(T(10)), eve::log(T(3628800)), 1.0);
TTS_ULP_EQUAL(eve::lfactorial(T(5)), eve::log(T(120)), 1.0);
TTS_ULP_EQUAL(eve::lfactorial(T(0)), eve::log(T(1)), 1.0);
TTS_ULP_EQUAL(eve::lfactorial(T(1)), eve::log(T(1)), 1.0);
TTS_ULP_EQUAL(eve::lfactorial(T(200)), eve::log_abs_gamma(T(201)), 1.0);
}
};
3 changes: 2 additions & 1 deletion test/unit/module/special/log_gamma.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,8 @@ TTS_CASE_WITH("Check behavior of log_gamma on wide",
eve::nan(eve::as(e)) : x; }, a0), 5e-2);

TTS_ULP_EQUAL(log_gamma(T(0.5)), T(std::lgamma(v_t(0.5))), 1.);
TTS_ULP_EQUAL(log_gamma(T(-35)), T(std::lgamma(v_t(-35))), 0.5);
TTS_ULP_EQUAL(log_gamma(T(-35)), eve::nan(eve::as<T>()), 0.5);
TTS_ULP_EQUAL(log_gamma(T(-2.5)), eve::nan(eve::as<T>()), 0.5);
TTS_ULP_EQUAL(log_gamma(T(-75.3)), T(std::lgamma(v_t(-75.3))), 1);

if constexpr( eve::platform::supports_invalids )
Expand Down

0 comments on commit 8323b82

Please sign in to comment.