boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Merge branch 'mydevelop' of https://github.com/cohomology/math into integrate_1251

Browse Source
This commit is contained in:
John Maddock
2025-05-23 16:27:12 +01:00
parent aecd335492 0d6ff719b1
commit 036cf85c12
3 changed files with 165 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

40
include/boost/math/special_functions/gamma.hpp
View File

@@ -934,7 +934,7 @@ BOOST_MATH_GPU_ENABLED T full_igamma_prefix(T a, T z, const Policy& pol)
{
BOOST_MATH_STD_USING
if (z > tools::max_value<T>())
if (z > tools::max_value<T>() || (a > 0 && z == 0))
return 0;
T alz = a * log(z);
@@ -992,7 +992,7 @@ template <class T, class Policy, class Lanczos>
BOOST_MATH_GPU_ENABLED T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
{
BOOST_MATH_STD_USING
if (z >= tools::max_value<T>())
if (z >= tools::max_value<T>() || (a > 0 && z == 0))
return 0;
T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);
T prefix{};
@@ -1321,7 +1321,11 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp_final(T a, T x, bool normalised, b
int eval_method;
if(is_int && (x > 0.6))
if (x == 0)
{
eval_method = 2;
}
else if(is_int && (x > 0.6))
{
// calculate Q via finite sum:
invert = !invert;
@@ -1501,14 +1505,14 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp_final(T a, T x, bool normalised, b
#ifdef BOOST_MATH_HAS_NVRTC
if (boost::math::is_same_v<T, float>)
{
init_value = (normalised ? 1 : ::tgammaf(a));
init_value = (normalised ? T(1) : ::tgammaf(a));
}
else
{
init_value = (normalised ? 1 : ::tgamma(a));
init_value = (normalised ? T(1) : ::tgamma(a));
}
#else
init_value = (normalised ? 1 : boost::math::tgamma(a, pol));
init_value = (normalised ? T(1) : boost::math::tgamma(a, pol));
#endif
if(normalised || (result >= 1) || (tools::max_value<T>() * result > init_value))
@@ -1640,32 +1644,26 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp_final(T a, T x, bool normalised, b
T gam;
if (boost::math::is_same_v<T, float>)
{
gam = normalised ? 1 : ::tgammaf(a);
gam = normalised ? T(1) : ::tgammaf(a);
}
else
{
gam = normalised ? 1 : ::tgamma(a);
gam = normalised ? T(1) : ::tgamma(a);
}
#else
T gam = normalised ? 1 : boost::math::tgamma(a, pol);
T gam = normalised ? T(1) : boost::math::tgamma(a, pol);
#endif
result = gam - result;
}
if(p_derivative)
{
/*
* We can never reach this:
if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
if((x == 0) || ((x < 1) && (tools::max_value<T>() * x < *p_derivative)))
{
// overflow, just return an arbitrarily large value:
*p_derivative = tools::max_value<T>() / 2;
}
*/
BOOST_MATH_ASSERT((x >= 1) || (tools::max_value<T>() * x >= *p_derivative));
//
// Need to convert prefix term to derivative:
//
*p_derivative /= x;
else
*p_derivative /= x;
}
return result;
@@ -1687,7 +1685,7 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in
T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
if(a >= max_factorial<T>::value && !normalised)
if(x > 0 && a >= max_factorial<T>::value && !normalised)
{
//
// When we're computing the non-normalized incomplete gamma
@@ -2143,8 +2141,8 @@ BOOST_MATH_GPU_ENABLED T gamma_p_derivative_imp(T a, T x, const Policy& pol)
//
if(x == 0)
{
return (a > 1) ? 0 :
(a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", nullptr, pol);
return (a > 1) ? T(0) :
(a == 1) ? T(1) : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", nullptr, pol);
}
//
// Normal case:

1
test/Jamfile.v2
View File

@@ -196,6 +196,7 @@ test-suite special_fun :
[ run git_issue_1139.cpp ]
[ run git_issue_1175.cpp ]
[ run git_issue_1194.cpp ]
[ run git_issue_1249.cpp /boost/test//boost_unit_test_framework ]
[ run git_issue_1255.cpp ]
[ run git_issue_1247.cpp ]
[ run special_functions_test.cpp /boost/test//boost_unit_test_framework ]

145
test/git_issue_1249.cpp Normal file
View File

@@ -0,0 +1,145 @@
// (C) Copyright Kilian Kilger 2025.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/test/results_collector.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
using namespace std;
using namespace boost::math;
using namespace boost::math::policies;
using namespace boost::multiprecision;
typedef policy<
policies::domain_error<errno_on_error>,
policies::pole_error<errno_on_error>,
policies::overflow_error<errno_on_error>,
policies::evaluation_error<errno_on_error>
> c_policy;
template<typename T>
struct test_lower
{
T operator()(T a, T x) const
{
return tgamma_lower(a, x, c_policy());
}
T expected(T a) const
{
return T(0.0);
}
};
template<typename T>
struct test_upper
{
T operator()(T a, T x) const
{
return tgamma(a, x, c_policy());
}
T expected(T a) const
{
return tgamma(a, c_policy());
}
};
template<typename T>
struct test_gamma_p
{
T operator()(T a, T x) const
{
return gamma_p(a, x, c_policy());
}
T expected(T) const
{
return T(0.0);
}
};
template<typename T>
struct test_gamma_q
{
T operator()(T a, T x) const
{
return gamma_q(a, x, c_policy());
}
T expected(T) const
{
return T(1.0);
}
};
template<typename T, template<typename> class Fun>
void test_impl(T a)
{
Fun<T> fn;
errno = 0;
T x = T(0.0);
T result = fn(a, x);
int saveErrno = errno;
errno = 0;
T expected = fn.expected(a);
BOOST_CHECK(errno == saveErrno);
BOOST_CHECK_EQUAL(result, expected);
}
template<template<typename> class Fun, typename T>
void test_type_dispatch(T val)
{
if (val <= (std::numeric_limits<float>::max)())
test_impl<float, Fun>(static_cast<float>(val));
if (val <= (std::numeric_limits<double>::max)())
test_impl<double, Fun>(static_cast<double>(val));
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
test_impl<long double, Fun>(static_cast<long double>(val));
#endif
test_impl<cpp_bin_float_50, Fun>(static_cast<cpp_bin_float_50>(val));
}
template<template<typename> class Fun>
void test_impl()
{
test_type_dispatch<Fun, double>(1.0);
test_type_dispatch<Fun, double>(0.1);
test_type_dispatch<Fun, double>(0.5);
test_type_dispatch<Fun, double>(0.6);
test_type_dispatch<Fun, double>(1.3);
test_type_dispatch<Fun, double>(1.5);
test_type_dispatch<Fun, double>(2);
test_type_dispatch<Fun, double>(100);
test_type_dispatch<Fun, double>((std::numeric_limits<float>::max)());
test_type_dispatch<Fun, double>((std::numeric_limits<double>::max)());
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
test_type_dispatch<Fun, long double>((std::numeric_limits<long double>::max)());
#endif
}
void test_derivative()
{
using namespace boost::math::detail;
double derivative = 0;
double result = gamma_incomplete_imp(1.0, 0.0, true, false, c_policy(), &derivative);
BOOST_CHECK(errno == 0);
BOOST_CHECK_EQUAL(derivative, tools::max_value<double>() / 2);
BOOST_CHECK_EQUAL(result, 0);
}
BOOST_AUTO_TEST_CASE( test_main )
{
test_impl<test_lower>();
test_impl<test_upper>();
test_impl<test_gamma_p>();
test_impl<test_gamma_q>();
test_derivative();
}