diff --git a/.drone.star b/.drone.star index 691246e337..17b2006d7a 100644 --- a/.drone.star +++ b/.drone.star @@ -47,22 +47,22 @@ def main(ctx): result.append(linux_cxx("Ubuntu g++-8 " + cxx + " " + suite, "g++-8", packages="g++-8", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-8', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubuntu g++-9 " + cxx + " " + suite, "g++-9", packages="g++-9", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-9', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in clang_6_stds: - result.append(linux_cxx("Ubuntu clang++-6 " + cxx + " " + suite, "clang++-6.0", packages="clang-6.0", llvm_os="xenial", llvm_ver="6.0", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-6.0', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - result.append(linux_cxx("Ubuntu clang++-7 " + cxx + " " + suite, "clang++-7", packages="clang-7", llvm_os="xenial", llvm_ver="7", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-7', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - result.append(linux_cxx("Ubuntu clang++-8 " + cxx + " " + suite, "clang++-8", packages="clang-8", llvm_os="xenial", llvm_ver="8", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-8', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - result.append(linux_cxx("Ubuntu clang++-9 " + cxx + " " + suite, "clang++-9", packages="clang-9", llvm_os="xenial", llvm_ver="9", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-9', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu clang++-6 " + cxx + " " + suite, "clang++-6.0", packages="clang-6.0", llvm_os="bionic", llvm_ver="6.0", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-6.0', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu clang++-7 " + cxx + " " + suite, "clang++-7", packages="clang-7", llvm_os="bionic", llvm_ver="7", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-7', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu clang++-8 " + cxx + " " + suite, "clang++-8", packages="clang-8", llvm_os="bionic", llvm_ver="8", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-8', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu clang++-9 " + cxx + " " + suite, "clang++-9", packages="clang-9", llvm_os="focal", llvm_ver="9", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-9', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in gnu_9_stds: - result.append(linux_cxx("Ubuntu g++-10 " + cxx + " " + suite, "g++-10", packages="g++-10", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-10', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - result.append(linux_cxx("Ubuntu g++-11 " + cxx + " " + suite, "g++-11", packages="g++-11", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-11', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu g++-10 " + cxx + " " + suite, "g++-10", packages="g++-10", buildtype="boost", image="cppalliance/droneubuntu2204:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-10', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu g++-11 " + cxx + " " + suite, "g++-11", packages="g++-11", buildtype="boost", image="cppalliance/droneubuntu2204:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-11', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubuntu g++-12 " + cxx + " " + suite, "g++-12", packages="g++-12", buildtype="boost", image="cppalliance/droneubuntu2204:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-12', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubuntu clang++-14 " + cxx + " " + suite, "clang++-14", packages="clang-14", llvm_os="jammy", llvm_ver="14", buildtype="boost", image="cppalliance/droneubuntu2204:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-14', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubuntu clang++-15 " + cxx + " " + suite, "clang++-15", packages="clang-15", llvm_os="jammy", llvm_ver="15", buildtype="boost", image="cppalliance/droneubuntu2204:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-15', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in clang_10_stds: result.append(linux_cxx("Ubuntu clang++-10 " + cxx + " " + suite, "clang++-10", packages="clang-10", llvm_os="xenial", llvm_ver="10", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-10', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in gnu_non_native: - result.append(linux_cxx("Ubuntu g++ s390s " + cxx + " " + suite, "g++", packages="g++", buildtype="boost", image="cppalliance/droneubuntu2204:multiarch", arch="s390x", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu g++ s390s " + cxx + " " + suite, "g++", packages="g++", buildtype="boost", image="cppalliance/droneubuntu2404:multiarch", arch="s390x", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in gnu_non_native: - result.append(linux_cxx("Ubuntu g++ ARM64" + cxx + " " + suite, "g++", packages="g++", buildtype="boost", image="cppalliance/droneubuntu2204:multiarch", arch="arm64", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu g++ ARM64" + cxx + " " + suite, "g++", packages="g++", buildtype="boost", image="cppalliance/droneubuntu2404:multiarch", arch="arm64", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in gnu_non_native: result.append(osx_cxx("M1 Clang " + cxx + " " + suite, "clang++", buildscript="drone", buildtype="boost", xcode_version="14.1", environment={'TOOLSET': 'clang', 'CXXSTD': cxx, 'TEST_SUITE': suite, 'DEFINE': 'BOOST_MATH_NO_REAL_CONCEPT_TESTS,BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS,BOOST_MATH_MULTI_ARCH_CI_RUN', }, globalenv=globalenv)) for suite in gcc13_things_to_test: diff --git a/.drone/drone.sh b/.drone/drone.sh index 86dbf129f5..7fc9d61b11 100644 --- a/.drone/drone.sh +++ b/.drone/drone.sh @@ -51,8 +51,8 @@ echo '==================================> BEFORE_SCRIPT' echo '==================================> SCRIPT' echo "using $TOOLSET : : $COMPILER : -std=$CXXSTD $OPTIONS ;" > ~/user-config.jam -(cd libs/config/test && ../../../b2 config_info_travis_install toolset=$TOOLSET && ./config_info_travis) -(cd libs/math/test && ../../../b2 toolset=$TOOLSET $TEST_SUITE define=$DEFINE) +(cd libs/config/test && ../../../b2 print_config_info print_math_info toolset=$TOOLSET) +(cd libs/math/test && ../../../b2 -j3 toolset=$TOOLSET $TEST_SUITE) echo '==================================> AFTER_SUCCESS' diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index fee18d22c7..8a6bbcf460 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -28,7 +28,7 @@ jobs: standard: [ c++20 ] suite: [ github_ci_block_1, github_ci_block_2 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Set TOOLSET @@ -86,7 +86,7 @@ jobs: compiler: [ g++-13, clang++-19 ] standard: [ c++14, c++17, c++20, c++23 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Set TOOLSET @@ -145,7 +145,7 @@ jobs: standard: [ 14, 17, 20 ] suite: [ github_ci_block_1, github_ci_block_2 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Checkout main boost @@ -188,7 +188,7 @@ jobs: standard: [ 14 ] suite: [ github_ci_block_1, github_ci_block_2 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Checkout main boost @@ -231,7 +231,7 @@ jobs: standard: [ 14, 17, latest ] suite: [ github_ci_block_1, github_ci_block_2 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Checkout main boost @@ -271,7 +271,7 @@ jobs: env: TOOLSET: gcc steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Install Cygwin @@ -305,7 +305,7 @@ jobs: matrix: compiler: [ g++-13 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Add repository @@ -345,7 +345,7 @@ jobs: matrix: compiler: [ clang++-19 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Add repository @@ -387,7 +387,7 @@ jobs: standard: [ c++14, c++17, c++20, c++23 ] suite: [ github_ci_block_1, github_ci_block_2 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Set TOOLSET @@ -448,7 +448,7 @@ jobs: runs-on: ${{matrix.os}} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Install packages if: matrix.install @@ -486,6 +486,227 @@ jobs: cd ../boost-root/__build__ cmake --build . --target tests + posix-cmake-subdir: + strategy: + fail-fast: false + matrix: + include: + - os: ubuntu-22.04 + - os: ubuntu-24.04 + - os: macos-15 + + runs-on: ${{matrix.os}} + + steps: + - uses: actions/checkout@v6 + + - name: Use library with add_subdirectory (header-only) + run: | + cd test/cmake_subdir_test + mkdir __build__ && cd __build__ + cmake .. + cmake --build . + ctest --output-on-failure --no-tests=error + + posix-cmake-subdir-legacy: + strategy: + fail-fast: false + matrix: + include: + - os: ubuntu-22.04 + - os: ubuntu-24.04 + - os: macos-15 + + runs-on: ${{matrix.os}} + + steps: + - uses: actions/checkout@v6 + + - name: Setup Boost + run: | + echo GITHUB_REPOSITORY: $GITHUB_REPOSITORY + LIBRARY=${GITHUB_REPOSITORY#*/} + echo LIBRARY: $LIBRARY + echo "LIBRARY=$LIBRARY" >> $GITHUB_ENV + echo GITHUB_BASE_REF: $GITHUB_BASE_REF + echo GITHUB_REF: $GITHUB_REF + REF=${GITHUB_BASE_REF:-$GITHUB_REF} + REF=${REF#refs/heads/} + echo REF: $REF + BOOST_BRANCH=develop && [ "$REF" == "master" ] && BOOST_BRANCH=master || true + echo BOOST_BRANCH: $BOOST_BRANCH + cd .. + git clone -b $BOOST_BRANCH --depth 1 https://github.com/boostorg/boost.git boost-root + cd boost-root + mkdir -p libs/$LIBRARY + cp -r $GITHUB_WORKSPACE/* libs/$LIBRARY + git submodule update --init tools/boostdep + python tools/boostdep/depinst/depinst.py --git_args "--jobs 3" $LIBRARY + + - name: Use library with add_subdirectory (legacy compiled) + run: | + cd ../boost-root/libs/$LIBRARY/test/cmake_subdir_test_legacy + mkdir __build__ && cd __build__ + cmake .. + cmake --build . + ctest --output-on-failure --no-tests=error + + posix-cmake-install: + strategy: + fail-fast: false + matrix: + include: + - os: macos-15 + + runs-on: ${{matrix.os}} + + steps: + - uses: actions/checkout@v6 + + - name: Setup Boost + run: | + echo GITHUB_REPOSITORY: $GITHUB_REPOSITORY + LIBRARY=${GITHUB_REPOSITORY#*/} + echo LIBRARY: $LIBRARY + echo "LIBRARY=$LIBRARY" >> $GITHUB_ENV + echo GITHUB_BASE_REF: $GITHUB_BASE_REF + echo GITHUB_REF: $GITHUB_REF + REF=${GITHUB_BASE_REF:-$GITHUB_REF} + REF=${REF#refs/heads/} + echo REF: $REF + BOOST_BRANCH=develop && [ "$REF" == "master" ] && BOOST_BRANCH=master || true + echo BOOST_BRANCH: $BOOST_BRANCH + cd .. + git clone -b $BOOST_BRANCH --depth 1 https://github.com/boostorg/boost.git boost-root + cd boost-root + mkdir -p libs/$LIBRARY + cp -r $GITHUB_WORKSPACE/* libs/$LIBRARY + git submodule update --init tools/boostdep + python tools/boostdep/depinst/depinst.py --git_args "--jobs 3" $LIBRARY + + - name: Configure + run: | + cd ../boost-root + mkdir __build__ && cd __build__ + cmake -DBOOST_INCLUDE_LIBRARIES=$LIBRARY -DCMAKE_INSTALL_PREFIX=~/.local .. + + - name: Install + run: | + cd ../boost-root/__build__ + cmake --build . --target install + + - name: Use the installed library (header-only) + run: | + cd ../boost-root/libs/$LIBRARY/test/cmake_install_test + mkdir __build__ && cd __build__ + cmake -DCMAKE_INSTALL_PREFIX=~/.local .. + cmake --build . + ctest --output-on-failure --no-tests=error + + posix-cmake-install-legacy-static: + strategy: + fail-fast: false + matrix: + include: + - os: macos-15 + + runs-on: ${{matrix.os}} + + steps: + - uses: actions/checkout@v6 + + - name: Setup Boost + run: | + echo GITHUB_REPOSITORY: $GITHUB_REPOSITORY + LIBRARY=${GITHUB_REPOSITORY#*/} + echo LIBRARY: $LIBRARY + echo "LIBRARY=$LIBRARY" >> $GITHUB_ENV + echo GITHUB_BASE_REF: $GITHUB_BASE_REF + echo GITHUB_REF: $GITHUB_REF + REF=${GITHUB_BASE_REF:-$GITHUB_REF} + REF=${REF#refs/heads/} + echo REF: $REF + BOOST_BRANCH=develop && [ "$REF" == "master" ] && BOOST_BRANCH=master || true + echo BOOST_BRANCH: $BOOST_BRANCH + cd .. + git clone -b $BOOST_BRANCH --depth 1 https://github.com/boostorg/boost.git boost-root + cd boost-root + mkdir -p libs/$LIBRARY + cp -r $GITHUB_WORKSPACE/* libs/$LIBRARY + git submodule update --init tools/boostdep + python tools/boostdep/depinst/depinst.py --git_args "--jobs 3" $LIBRARY + + - name: Configure + run: | + cd ../boost-root + mkdir __build__ && cd __build__ + cmake -DBOOST_INCLUDE_LIBRARIES=$LIBRARY -DBOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS=ON -DCMAKE_INSTALL_PREFIX=~/.local .. + + - name: Install + run: | + cd ../boost-root/__build__ + cmake --build . --target install + + - name: Use the installed library (legacy compiled) + run: | + cd ../boost-root/libs/$LIBRARY/test/cmake_install_test_legacy + mkdir __build__ && cd __build__ + cmake -DCMAKE_INSTALL_PREFIX=~/.local .. + cmake --build . + ctest --output-on-failure --no-tests=error + + posix-cmake-install-legacy-shared: + strategy: + fail-fast: false + matrix: + include: + - os: macos-15 + + runs-on: ${{matrix.os}} + + steps: + - uses: actions/checkout@v6 + + - name: Setup Boost + run: | + echo GITHUB_REPOSITORY: $GITHUB_REPOSITORY + LIBRARY=${GITHUB_REPOSITORY#*/} + echo LIBRARY: $LIBRARY + echo "LIBRARY=$LIBRARY" >> $GITHUB_ENV + echo GITHUB_BASE_REF: $GITHUB_BASE_REF + echo GITHUB_REF: $GITHUB_REF + REF=${GITHUB_BASE_REF:-$GITHUB_REF} + REF=${REF#refs/heads/} + echo REF: $REF + BOOST_BRANCH=develop && [ "$REF" == "master" ] && BOOST_BRANCH=master || true + echo BOOST_BRANCH: $BOOST_BRANCH + cd .. + git clone -b $BOOST_BRANCH --depth 1 https://github.com/boostorg/boost.git boost-root + cd boost-root + mkdir -p libs/$LIBRARY + cp -r $GITHUB_WORKSPACE/* libs/$LIBRARY + git submodule update --init tools/boostdep + python tools/boostdep/depinst/depinst.py --git_args "--jobs 3" $LIBRARY + + - name: Configure + run: | + cd ../boost-root + mkdir __build__ && cd __build__ + cmake -DBOOST_INCLUDE_LIBRARIES=$LIBRARY -DBOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS=ON -DCMAKE_INSTALL_PREFIX=~/.local -DBUILD_SHARED_LIBS=ON .. + + - name: Install + run: | + cd ../boost-root/__build__ + cmake --build . --target install + + - name: Use the installed library (legacy compiled) + run: | + cd ../boost-root/libs/$LIBRARY/test/cmake_install_test_legacy + mkdir __build__ && cd __build__ + cmake -DCMAKE_INSTALL_PREFIX=~/.local .. + cmake --build . + ctest --output-on-failure --no-tests=error + sycl-cmake-test: strategy: fail-fast: false @@ -512,7 +733,7 @@ jobs: printenv >> $GITHUB_ENV - name: checkout project code - uses: actions/checkout@v4 + uses: actions/checkout@v6 - name: Install Packages run: | @@ -570,7 +791,7 @@ jobs: echo "Installed cuda version is: ${{steps.cuda-toolkit.outputs.cuda}}"+ echo "Cuda install location: ${{steps.cuda-toolkit.outputs.CUDA_PATH}}" nvcc -V - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Install Packages run: | @@ -599,7 +820,7 @@ jobs: run: | cd ../boost-root mkdir __build__ && cd __build__ - cmake -DBOOST_INCLUDE_LIBRARIES=$LIBRARY -DBUILD_TESTING=ON -DCMAKE_CUDA_COMPILER=/usr/local/cuda/bin/nvcc -DBOOST_MATH_ENABLE_CUDA=1 -DCMAKE_CUDA_ARCHITECTURES=86 -DCUDA_TOOLKIT_ROOT_DIR=/usr/local/cuda-12.8 .. + cmake -DBOOST_INCLUDE_LIBRARIES=$LIBRARY -DBUILD_TESTING=ON -DCMAKE_CUDA_COMPILER=/usr/local/cuda/bin/nvcc -DBOOST_MATH_ENABLE_CUDA=1 -DCMAKE_CUDA_ARCHITECTURES="75;86" -DCMAKE_CUDA_STANDARD=17 .. - name: Build tests run: | cd ../boost-root/__build__ diff --git a/.github/workflows/codecov.yml b/.github/workflows/codecov.yml index 6c76d1f238..f867c38311 100644 --- a/.github/workflows/codecov.yml +++ b/.github/workflows/codecov.yml @@ -92,7 +92,7 @@ jobs: fi git config --global pack.threads 0 - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: # For coverage builds fetch the whole history, else only 1 commit using a 'fake ternary' fetch-depth: ${{ matrix.coverage && '0' || '1' }} @@ -106,7 +106,7 @@ jobs: restore-keys: ${{matrix.os}}-${{matrix.container}}-${{matrix.compiler}}- - name: Fetch Boost.CI - uses: actions/checkout@v4 + uses: actions/checkout@v6 with: repository: boostorg/boost-ci ref: master diff --git a/CMakeLists.txt b/CMakeLists.txt index b789eb22aa..616c21e27f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,7 +5,7 @@ cmake_minimum_required(VERSION 3.8...3.16) -project(boost_math VERSION 1.91.0 LANGUAGES CXX) +project(boost_math VERSION 1.92.0 LANGUAGES CXX) add_library(boost_math INTERFACE) @@ -46,6 +46,136 @@ endif() target_compile_features(boost_math INTERFACE cxx_std_14) +# Legacy C99/TR1 compiled libraries +option(BOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS "Build the C99 and TR1 compiled libraries" OFF) + +if(BOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS) + +include(CheckCXXSourceCompiles) +set(_config_inc "") +if(NOT BOOST_MATH_STANDALONE) + get_target_property(_config_type Boost::config TYPE) + if(_config_type STREQUAL "INTERFACE_LIBRARY") + get_target_property(_config_inc Boost::config INTERFACE_INCLUDE_DIRECTORIES) + endif() +endif() +set(CMAKE_REQUIRED_INCLUDES "${CMAKE_CURRENT_SOURCE_DIR}/include" ${_config_inc}) +check_cxx_source_compiles("#include <${CMAKE_CURRENT_SOURCE_DIR}/config/has_long_double_support.cpp> \n int main() { return 0;}" BOOST_MATH_HAS_LONG_DOUBLE) +unset(CMAKE_REQUIRED_INCLUDES) + +set(C99_SOURCES + acosh + asinh + atanh + cbrt + copysign + erfc + erf + expm1 + fmax + fmin + fpclassify + hypot + lgamma + llround + log1p + lround + nextafter + nexttoward + round + tgamma + trunc +) + +set(TR1_SOURCES + assoc_laguerre + assoc_legendre + beta + comp_ellint_1 + comp_ellint_2 + comp_ellint_3 + cyl_bessel_i + cyl_bessel_j + cyl_bessel_k + cyl_neumann + ellint_1 + ellint_2 + ellint_3 + expint + hermite + laguerre + legendre + riemann_zeta + sph_bessel + sph_legendre + sph_neumann +) + +list(TRANSFORM C99_SOURCES PREPEND "src/tr1/") +list(TRANSFORM TR1_SOURCES PREPEND "src/tr1/") + +list(TRANSFORM C99_SOURCES APPEND "f.cpp" OUTPUT_VARIABLE C99_SOURCESf) +list(TRANSFORM TR1_SOURCES APPEND "f.cpp" OUTPUT_VARIABLE TR1_SOURCESf) + +set(types "" f) + +if(BOOST_MATH_HAS_LONG_DOUBLE) + list(TRANSFORM C99_SOURCES APPEND "l.cpp" OUTPUT_VARIABLE C99_SOURCESl) + list(TRANSFORM TR1_SOURCES APPEND "l.cpp" OUTPUT_VARIABLE TR1_SOURCESl) + list(APPEND types l) +endif() + +list(TRANSFORM C99_SOURCES APPEND ".cpp") +list(TRANSFORM TR1_SOURCES APPEND ".cpp") + +foreach(type IN LISTS types) + add_library(boost_math_tr1${type} ${TR1_SOURCES${type}}) + add_library(Boost::math_tr1${type} ALIAS boost_math_tr1${type}) + if(NOT BOOST_MATH_STANDALONE) + target_link_libraries(boost_math_tr1${type} PUBLIC Boost::config) + endif() + target_include_directories(boost_math_tr1${type} PRIVATE src/tr1) + target_include_directories(boost_math_tr1${type} PRIVATE include) + + add_library(boost_math_c99${type} ${C99_SOURCES${type}}) + add_library(Boost::math_c99${type} ALIAS boost_math_c99${type}) + if(NOT BOOST_MATH_STANDALONE) + target_link_libraries(boost_math_c99${type} PUBLIC Boost::config) + endif() + target_include_directories(boost_math_c99${type} PRIVATE src/tr1) + target_include_directories(boost_math_c99${type} PRIVATE include) + + if(BUILD_SHARED_LIBS) + target_compile_definitions(boost_math_tr1${type} PUBLIC BOOST_MATH_TR1_DYN_LINK=1) + target_compile_definitions(boost_math_c99${type} PUBLIC BOOST_MATH_TR1_DYN_LINK=1) + if(WIN32) + target_compile_definitions(boost_math_tr1${type} PRIVATE "BOOST_SYMBOL_EXPORT=__declspec(dllexport)" BOOST_ALL_NO_LIB) + target_compile_definitions(boost_math_c99${type} PRIVATE "BOOST_SYMBOL_EXPORT=__declspec(dllexport)" BOOST_ALL_NO_LIB) + else() + target_compile_definitions(boost_math_tr1${type} PRIVATE "BOOST_SYMBOL_EXPORT=__attribute__((visibility(\"default\")))" BOOST_ALL_NO_LIB) + target_compile_definitions(boost_math_c99${type} PRIVATE "BOOST_SYMBOL_EXPORT=__attribute__((visibility(\"default\")))" BOOST_ALL_NO_LIB) + endif() + endif() + + target_compile_features(boost_math_tr1${type} PUBLIC cxx_std_14) + target_compile_features(boost_math_c99${type} PUBLIC cxx_std_14) + + if(APPLE) + target_compile_definitions(boost_math_tr1${type} PRIVATE _DARWIN_C_SOURCE) + target_compile_definitions(boost_math_c99${type} PRIVATE _DARWIN_C_SOURCE) + endif() +endforeach() + +if(BOOST_SUPERPROJECT_VERSION AND NOT CMAKE_VERSION VERSION_LESS 3.13) + set(boost_math_legacy_targets "") + foreach(type IN LISTS types) + list(APPEND boost_math_legacy_targets boost_math_c99${type} boost_math_tr1${type}) + endforeach() + boost_install(TARGETS ${boost_math_legacy_targets} VERSION "${BOOST_SUPERPROJECT_VERSION}") +endif() + +endif() + if(BUILD_TESTING AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/test/CMakeLists.txt") add_subdirectory(test) diff --git a/README.md b/README.md index 39b0bddd63..b44755ec55 100644 --- a/README.md +++ b/README.md @@ -138,8 +138,12 @@ The Boost Math Library is located in `libs/math/`. ### Running tests -First, make sure you are in `libs/math/test`. -You can either run all the tests listed in `Jamfile.v2` or run a single test: +The Boost build program `b2` is used to run the tests. To build `b2`, from the +top-level Boost directory, run either `./bootstrap.sh` (Linux/MacOS) or +`bootstrap.bat` (Windows). + +Then move to `libs/math/test`. You can either run all the tests listed in +`Jamfile.v2` or run a single test: test$ ../../../b2 <- run all tests test$ ../../../b2 static_assert_test <- single test diff --git a/doc/distributions/students_t.qbk b/doc/distributions/students_t.qbk index 3396048c5f..2ff2812cbb 100644 --- a/doc/distributions/students_t.qbk +++ b/doc/distributions/students_t.qbk @@ -29,6 +29,11 @@ RealType beta, RealType sd, RealType hint = 100); + + // degrees of freedom inversion from a quantile and probability: + BOOST_MATH_GPU_ENABLED static RealType find_degrees_of_freedom( + RealType t, + RealType p); }; }} // namespaces @@ -106,6 +111,13 @@ For more information on this function see the [@http://www.itl.nist.gov/div898/handbook/prc/section2/prc222.htm NIST Engineering Statistics Handbook]. + BOOST_MATH_GPU_ENABLED static RealType find_degrees_of_freedom( + RealType t, + RealType p); + +Returns the degrees of freedom [nu] such that CDF(/x/; [nu]) = /p/. +Requires 0 < /p/ < 1 and /x/ [ne] 0, otherwise calls __domain_error. + [h4 Non-member Accessors] All the [link math_toolkit.dist_ref.nmp usual non-member accessor functions] that are generic to all @@ -164,6 +176,7 @@ without the subtraction implied above.]] [[skewness][if (v > 3) 0 else NaN ]] [[kurtosis][if (v > 4) 3 * (v - 2) \/ (v - 4) else NaN]] [[kurtosis excess][if (v > 4) 6 \/ (df - 4) else NaN]] +[[find_degrees_of_freedom(t, p)][Uses a 2nd-order Edgeworth expansion as initial guess for a numerical root finder]] ] If the moment index /k/ is less than /v/, then the moment is undefined. diff --git a/include/boost/math/constants/constants.hpp b/include/boost/math/constants/constants.hpp index bdc0a74e27..3daf2d266e 100644 --- a/include/boost/math/constants/constants.hpp +++ b/include/boost/math/constants/constants.hpp @@ -129,7 +129,7 @@ namespace boost{ namespace math { initializer() { - F(); + (F)(); } void force_instantiate()const{} }; @@ -151,7 +151,7 @@ namespace boost{ namespace math { initializer() { - F(); + (F)(); } void force_instantiate()const{} }; diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index 29feac48fe..4fe3c7d751 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -164,7 +164,104 @@ namespace boost return check_dist(function, x_min, x_max, result, pol) && check_prob(function, p, result, pol); } // bool check_dist_and_prob + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_mean(const RealType x_min, const RealType x_max) + { + return (x_min + x_max) / 2; + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_variance(const RealType x_min, const RealType x_max) + { + return (x_max - x_min) * (x_max - x_min) / 8; + } + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_pdf(const RealType x, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + using boost::math::constants::pi; + return 1 / (pi() * sqrt((x - x_min) * (x_max - x))); + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_cdf(const RealType x, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases: + if (x == x_min) + { + return 0; + } + else if (x == x_max) + { + return 1; + } + using boost::math::constants::pi; + return 2 * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_ccdf(const RealType x, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases + if (x == x_min) + { + return 1; + } + else if (x == x_max) + { + return 0; + } + using boost::math::constants::pi; + // Naive version x = 1 - x; + // result = static_cast(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); + // is less accurate, so use acos instead of asin for complement. + return 2 * acos(sqrt((x - x_min) / (x_max - x_min))) / pi(); + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_quantile( const RealType p, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases: + if (p == 0) + { + return 0; + } + if (p == 1) + { + return 1; + } + using boost::math::constants::half_pi; + RealType sin2hpip = sin(half_pi() * p); + RealType sin2hpip2 = sin2hpip * sin2hpip; + return -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; + } + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_cquantile(const RealType q, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases: + if (q == 1) + { + return 0; + } + if (q == 0) + { + return 1; + } + // Naive RealType p = 1 - q; result = sin(half_pi() * p); loses accuracy, so use a cos alternative instead. + //result = cos(half_pi() * q); // for arcsine(0,1) + //result = result * result; + // For generalized arcsine: + using boost::math::constants::half_pi; + RealType cos2hpip = cos(half_pi() * q); + RealType cos2hpip2 = cos2hpip * cos2hpip; + return -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; + } } // namespace arcsine_detail template > @@ -230,7 +327,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::mean(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -239,7 +336,9 @@ namespace boost { return result; } - return (x_min + x_max) / 2; + + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_mean(static_cast(x_min), static_cast(x_max))); } // mean template @@ -248,7 +347,7 @@ namespace boost RealType result; RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::variance(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -257,7 +356,9 @@ namespace boost { return result; } - return (x_max - x_min) * (x_max - x_min) / 8; + + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_variance(static_cast(x_min), static_cast(x_max))); } // variance template @@ -277,7 +378,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); RealType result; - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::median(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -286,7 +387,9 @@ namespace boost { return result; } - return (x_min + x_max) / 2; + // The median is the same as the mean + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_mean(static_cast(x_min), static_cast(x_max))); } template @@ -296,7 +399,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::skewness(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -315,7 +418,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::kurtosis_excess(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -324,8 +427,7 @@ namespace boost { return result; } - result = -3; - return result / 2; + return static_cast(-1.5f); } // kurtosis_excess template @@ -335,7 +437,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::kurtosis(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -344,75 +446,61 @@ namespace boost { return result; } - - return 3 + kurtosis_excess(dist); + return static_cast(1.5f); } // kurtosis template BOOST_MATH_GPU_ENABLED inline RealType pdf(const arcsine_distribution& dist, const RealType& xx) { // Probability Density/Mass Function arcsine. BOOST_FPU_EXCEPTION_GUARD - BOOST_MATH_STD_USING // For ADL of std functions. - - constexpr auto function = "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)"; - RealType lo = dist.x_min(); - RealType hi = dist.x_max(); + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); RealType x = xx; // Argument checks: - RealType result = 0; - if (false == arcsine_detail::check_dist_and_x( - function, - lo, hi, x, + RealType result; + if (!arcsine_detail::check_dist_and_x( + "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + x, &result, Policy())) { return result; } - using boost::math::constants::pi; - result = static_cast(1) / (pi() * sqrt((x - lo) * (hi - x))); - return result; + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_pdf(static_cast(x), + static_cast(x_min), + static_cast(x_max))); } // pdf template BOOST_MATH_GPU_ENABLED inline RealType cdf(const arcsine_distribution& dist, const RealType& x) { // Cumulative Distribution Function arcsine. - BOOST_MATH_STD_USING // For ADL of std functions. - - constexpr auto function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; - RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); // Argument checks: RealType result = 0; - if (false == arcsine_detail::check_dist_and_x( - function, - x_min, x_max, x, + if (!arcsine_detail::check_dist_and_x( + "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + x, &result, Policy())) { return result; } - // Special cases: - if (x == x_min) - { - return 0; - } - else if (x == x_max) - { - return 1; - } - using boost::math::constants::pi; - result = static_cast(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); - return result; + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_cdf(static_cast(x), + static_cast(x_min), + static_cast(x_max))); } // arcsine cdf template BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type, RealType>& c) { // Complemented Cumulative Distribution Function arcsine. - BOOST_MATH_STD_USING // For ADL of std functions. - constexpr auto function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; - RealType x = c.param; arcsine_distribution const& dist = c.dist; RealType x_min = dist.x_min(); @@ -420,27 +508,19 @@ namespace boost // Argument checks: RealType result = 0; - if (false == arcsine_detail::check_dist_and_x( - function, - x_min, x_max, x, + if (!arcsine_detail::check_dist_and_x( + "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + x, &result, Policy())) { return result; } - if (x == x_min) - { - return 0; - } - else if (x == x_max) - { - return 1; - } - using boost::math::constants::pi; - // Naive version x = 1 - x; - // result = static_cast(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); - // is less accurate, so use acos instead of asin for complement. - result = static_cast(2) * acos(sqrt((x - x_min) / (x_max - x_min))) / pi(); - return result; + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_ccdf(static_cast(x), + static_cast(x_min), + static_cast(x_max))); } // arcsine ccdf template @@ -454,37 +534,22 @@ namespace boost // and return a value such that the probability that a random variable x // will be less than or equal to that value // is whatever probability you supplied as an argument. - BOOST_MATH_STD_USING // For ADL of std functions. - - using boost::math::constants::half_pi; - - constexpr auto function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; - RealType result = 0; // of argument checks: RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist_and_prob( - function, - x_min, x_max, p, + if (!arcsine_detail::check_dist_and_prob( + "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + p, &result, Policy())) { return result; } - // Special cases: - if (p == 0) - { - return 0; - } - if (p == 1) - { - return 1; - } - - RealType sin2hpip = sin(half_pi() * p); - RealType sin2hpip2 = sin2hpip * sin2hpip; - result = -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; - - return result; + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_quantile(static_cast(p), + static_cast(x_min), + static_cast(x_max))); } // quantile template @@ -493,19 +558,14 @@ namespace boost // Complement Quantile or Percent Point arcsine function. // Return the number of expected x for a given // complement of the probability q. - BOOST_MATH_STD_USING // For ADL of std functions. - - using boost::math::constants::half_pi; - constexpr auto function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; - // Error checks: RealType q = c.param; const arcsine_distribution& dist = c.dist; RealType result = 0; RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist_and_prob( - function, + if (!arcsine_detail::check_dist_and_prob( + "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)", x_min, x_max, q, @@ -513,26 +573,11 @@ namespace boost { return result; } - // Special cases: - if (q == 1) - { - return 0; - } - if (q == 0) - { - return 1; - } - // Naive RealType p = 1 - q; result = sin(half_pi() * p); loses accuracy, so use a cos alternative instead. - //result = cos(half_pi() * q); // for arcsine(0,1) - //result = result * result; - // For generalized arcsine: - RealType cos2hpip = cos(half_pi() * q); - RealType cos2hpip2 = cos2hpip * cos2hpip; - result = -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; - - return result; + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_cquantile(static_cast(q), + static_cast(x_min), + static_cast(x_max))); } // Quantile Complement - } // namespace math } // namespace boost diff --git a/include/boost/math/distributions/binomial.hpp b/include/boost/math/distributions/binomial.hpp index b17893e422..262823170f 100644 --- a/include/boost/math/distributions/binomial.hpp +++ b/include/boost/math/distributions/binomial.hpp @@ -157,7 +157,7 @@ namespace boost template BOOST_MATH_CUDA_ENABLED inline bool check_dist_and_prob(const char* function, const RealType& N, RealType p, RealType prob, RealType* result, const Policy& pol) { - if((check_dist(function, N, p, result, pol) && detail::check_probability(function, prob, result, pol)) == false) + if(!(check_dist(function, N, p, result, pol) && detail::check_probability(function, prob, result, pol))) return false; return true; } @@ -226,7 +226,7 @@ namespace boost // but zero is the best we can do: return 0; } - if(p == 1 || success_fraction == 1) + if(q == 0 || success_fraction == 1) { // Probability of n or fewer successes is always one, // so n is the most sensible answer here: return trials; @@ -321,11 +321,11 @@ namespace boost BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_lower_bound_on_p"; // Error checks: RealType result = 0; - if(false == binomial_detail::check_dist_and_k( + if(!(binomial_detail::check_dist_and_k( function, trials, RealType(0), successes, &result, Policy()) && binomial_detail::check_dist_and_prob( - function, trials, RealType(0), probability, &result, Policy())) + function, trials, RealType(0), probability, &result, Policy()))) { return result; } if(successes == 0) @@ -346,11 +346,11 @@ namespace boost BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_upper_bound_on_p"; // Error checks: RealType result = 0; - if(false == binomial_detail::check_dist_and_k( + if(!(binomial_detail::check_dist_and_k( function, trials, RealType(0), successes, &result, Policy()) && binomial_detail::check_dist_and_prob( - function, trials, RealType(0), probability, &result, Policy())) + function, trials, RealType(0), probability, &result, Policy()))) { return result; } if(trials == successes) @@ -373,11 +373,11 @@ namespace boost BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_minimum_number_of_trials"; // Error checks: RealType result = 0; - if(false == binomial_detail::check_dist_and_k( + if(!(binomial_detail::check_dist_and_k( function, k, p, k, &result, Policy()) && binomial_detail::check_dist_and_prob( - function, k, p, alpha, &result, Policy())) + function, k, p, alpha, &result, Policy()))) { return result; } result = ibetac_invb(k + 1, p, alpha, Policy()); // returns n - k @@ -392,11 +392,11 @@ namespace boost BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_maximum_number_of_trials"; // Error checks: RealType result = 0; - if(false == binomial_detail::check_dist_and_k( + if(!(binomial_detail::check_dist_and_k( function, k, p, k, &result, Policy()) && binomial_detail::check_dist_and_prob( - function, k, p, alpha, &result, Policy())) + function, k, p, alpha, &result, Policy()))) { return result; } result = ibeta_invb(k + 1, p, alpha, Policy()); // returns n - k diff --git a/include/boost/math/distributions/cauchy.hpp b/include/boost/math/distributions/cauchy.hpp index 00a9ac9ac9..5a19f499a3 100644 --- a/include/boost/math/distributions/cauchy.hpp +++ b/include/boost/math/distributions/cauchy.hpp @@ -57,11 +57,11 @@ BOOST_MATH_GPU_ENABLED RealType cdf_imp(const cauchy_distribution((complement) ? 1 : 0); } #endif - if(false == detail::check_x(function, x, &result, Policy())) + if(!detail::check_x(function, x, &result, Policy())) { // Catches x == NaN return result; } @@ -111,15 +111,15 @@ BOOST_MATH_GPU_ENABLED RealType quantile_imp( RealType result = 0; RealType location = dist.location(); RealType scale = dist.scale(); - if(false == detail::check_location(function, location, &result, Policy())) + if(!detail::check_location(function, location, &result, Policy())) { return result; } - if(false == detail::check_scale(function, scale, &result, Policy())) + if(!detail::check_scale(function, scale, &result, Policy())) { return result; } - if(false == detail::check_probability(function, p, &result, Policy())) + if(!detail::check_probability(function, p, &result, Policy())) { return result; } @@ -224,11 +224,11 @@ BOOST_MATH_GPU_ENABLED inline RealType pdf(const cauchy_distribution::find constexpr auto function = "boost::math::chi_squared_distribution<%1%>::find_degrees_of_freedom(%1%,%1%,%1%,%1%,%1%)"; // Check for domain errors: RealType error_result; - if(false == - detail::check_probability(function, alpha, &error_result, Policy()) - && detail::check_probability(function, beta, &error_result, Policy())) + if(!(detail::check_probability(function, alpha, &error_result, Policy()) + && detail::check_probability(function, beta, &error_result, Policy()))) { // Either probability is outside 0 to 1. return error_result; } diff --git a/include/boost/math/distributions/extreme_value.hpp b/include/boost/math/distributions/extreme_value.hpp index 73454d29d4..7049ac6a15 100644 --- a/include/boost/math/distributions/extreme_value.hpp +++ b/include/boost/math/distributions/extreme_value.hpp @@ -173,8 +173,6 @@ BOOST_MATH_GPU_ENABLED inline RealType cdf(const extreme_value_distribution&, %1%)", x, &result, Policy())) return result; @@ -199,8 +197,6 @@ BOOST_MATH_GPU_ENABLED inline RealType logcdf(const extreme_value_distribution&, %1%)", x, &result, Policy())) return result; diff --git a/include/boost/math/distributions/find_location.hpp b/include/boost/math/distributions/find_location.hpp index 33c2e64b62..c23b41eb10 100644 --- a/include/boost/math/distributions/find_location.hpp +++ b/include/boost/math/distributions/find_location.hpp @@ -10,6 +10,7 @@ #include // for all distribution signatures. #include +#include #include #include #include @@ -51,10 +52,10 @@ namespace boost return policies::raise_domain_error( function, "z parameter was %1%, but must be finite!", z, pol); } - if(!(boost::math::isfinite)(scale)) + typename Dist::value_type result; + if(!(boost::math::detail::check_scale)(function, scale, &result, pol)) { - return policies::raise_domain_error( - function, "scale parameter was %1%, but must be finite!", scale, pol); + return result; } //cout << "z " << z << ", p " << p << ", quantile(Dist(), p) " @@ -94,11 +95,11 @@ namespace boost return policies::raise_domain_error( function, "z parameter was %1%, but must be finite!", z, policies::policy<>()); } + typename Dist::value_type result; typename Dist::value_type scale = c.param2; - if(!(boost::math::isfinite)(scale)) + if(!(boost::math::detail::check_scale)(function, scale, &result, policies::policy<>())) { - return policies::raise_domain_error( - function, "scale parameter was %1%, but must be finite!", scale, policies::policy<>()); + return result; } // cout << "z " << c.dist << ", quantile (Dist(), " << c.param1 << ") * scale " << c.param2 << endl; return z - quantile(Dist(), p) * scale; @@ -123,11 +124,11 @@ namespace boost return policies::raise_domain_error( function, "z parameter was %1%, but must be finite!", z, c.param3); } + typename Dist::value_type result; typename Dist::value_type scale = c.param2; - if(!(boost::math::isfinite)(scale)) + if(!(boost::math::detail::check_scale)(function, scale, &result, c.param3)) { - return policies::raise_domain_error( - function, "scale parameter was %1%, but must be finite!", scale, c.param3); + return result; } // cout << "z " << c.dist << ", quantile (Dist(), " << c.param1 << ") * scale " << c.param2 << endl; return z - quantile(Dist(), p) * scale; diff --git a/include/boost/math/distributions/inverse_gaussian.hpp b/include/boost/math/distributions/inverse_gaussian.hpp index 20d3b6bdd5..795ad23d4a 100644 --- a/include/boost/math/distributions/inverse_gaussian.hpp +++ b/include/boost/math/distributions/inverse_gaussian.hpp @@ -52,6 +52,7 @@ #include #include #include // for erf/erfc. +#include // for erf/erfc. #include #include #include @@ -298,7 +299,7 @@ struct inverse_gaussian_quantile_complement_functor namespace detail { template - BOOST_MATH_GPU_ENABLED inline RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1) + BOOST_MATH_GPU_ENABLED inline RealType guess_ig(RealType p, RealType q, RealType mu, RealType lambda) { // guess at random variate value x for inverse gaussian quantile. BOOST_MATH_STD_USING using boost::math::policies::policy; @@ -323,27 +324,16 @@ namespace detail x = mu * exp(quantile(n01, p) / sqrt(phi) - 1/(2 * phi)); } else - { // phi < 2 so much less symmetrical with long tail, - // so use gamma distribution as an approximation. - using boost::math::gamma_distribution; - - // Define the distribution, using gamma_nooverflow: - using gamma_nooverflow = gamma_distribution; - - gamma_nooverflow g(static_cast(0.5), static_cast(1.)); - - // R qgamma(0.2, 0.5, 1) = 0.0320923 - RealType qg = quantile(complement(g, p)); - x = lambda / (qg * 2); - // - if (x > mu/2) // x > mu /2? - { // x too large for the gamma approximation to work well. - //x = qgamma(p, 0.5, 1.0); // qgamma(0.270614, 0.5, 1) = 0.05983807 - RealType q = quantile(g, p); - // x = mu * exp(q * static_cast(0.1)); // Said to improve at high p - // x = mu * x; // Improves at high p? - x = mu * exp(q / sqrt(phi) - 1/(2 * phi)); - } + { + // Construct a gamma distribution with the same mean and standard deviation, and hope it + // get's us in more or less the right ballpark: + inverse_gaussian_distribution ig(mu, lambda); + RealType m = mean(ig); + RealType v = standard_deviation(ig); + RealType k = m * m / v; + RealType theta = v / m; + gamma_distribution n01(k, theta); + x = p < q ? quantile(n01, p) : quantile(complement(n01, q)); } return x; } // guess_ig @@ -379,7 +369,7 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const inverse_gaussian_distribut return result; // infinity; } - RealType guess = detail::guess_ig(p, dist.mean(), dist.scale()); + RealType guess = detail::guess_ig(p, RealType(1 - p), dist.mean(), dist.scale()); using boost::math::tools::max_value; RealType min = static_cast(0); // Minimum possible value is bottom of range of distribution. @@ -390,8 +380,19 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const inverse_gaussian_distribut int get_digits = policies::digits();// get digits from policy, boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); // and max iterations. using boost::math::tools::newton_raphson_iterate; + inverse_gaussian_quantile_functor func(dist, p); + // + // Our guess is often not very good, so lets bracket the root just in case: + // + RealType f0 = boost::math::get<0>(func(guess)); + if (f0 < 0) + tools::detail::bracket_root_towards_max(func, guess, f0, min, max, max_iter); + else + tools::detail::bracket_root_towards_min(func, guess, f0, min, max, max_iter); + if((guess < min) || (guess > max)) + guess = (min + max) / 2; result = - newton_raphson_iterate(inverse_gaussian_quantile_functor(dist, p), guess, min, max, get_digits, max_iter); + newton_raphson_iterate(func, guess, min, max, get_digits, max_iter); if (max_iter >= policies::get_max_root_iterations()) { return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE @@ -455,7 +456,7 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type(); boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); using boost::math::tools::newton_raphson_iterate; - result = newton_raphson_iterate(inverse_gaussian_quantile_complement_functor(c.dist, q), guess, min, max, get_digits, max_iter); + inverse_gaussian_quantile_complement_functor func(c.dist, q); + RealType f0 = boost::math::get<0>(func(guess)); + if (f0 > 0) + tools::detail::bracket_root_towards_max(func, guess, f0, min, max, max_iter); + else + tools::detail::bracket_root_towards_min(func, guess, f0, min, max, max_iter); + if ((guess < min) || (guess > max)) + guess = (min + max) / 2; + result = + newton_raphson_iterate(func, guess, min, max, get_digits, max_iter); + result = newton_raphson_iterate(func, guess, min, max, get_digits, max_iter); if (max_iter >= policies::get_max_root_iterations()) { return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 783c321980..9bae5c41ff 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -44,7 +45,8 @@ namespace boost }; template - BOOST_MATH_GPU_ENABLED RealType find_non_centrality_f(const RealType x, const RealType v1, const RealType v2, const RealType p, const RealType q, const RealType p_q_precision, const Policy& pol) { + BOOST_MATH_GPU_ENABLED RealType find_non_centrality_f(const RealType x, const RealType v1, const RealType v2, const RealType p, const RealType q, const RealType p_q_precision, const Policy& pol) + { constexpr auto function = "non_central_f<%1%>::find_non_centrality"; if ( p == 0 || q == 0) { @@ -80,6 +82,103 @@ namespace boost } return result; } + + template + struct f_degrees_of_freedom_finder + { + f_degrees_of_freedom_finder( + RealType x_, RealType v_, RealType nc_, bool find_v1_, RealType p_, bool c) + : x(x_), v(v_), nc(nc_), find_v1(find_v1_), p(p_), comp(c) {} + + RealType operator()(const RealType& input_v) + { + RealType v1 = find_v1 ? input_v : v; + RealType v2 = find_v1 ? v : input_v; + non_central_f_distribution d(v1, v2, nc); + return comp ? + p - cdf(complement(d, x)) + : cdf(d, x) - p; + } + private: + RealType x; + RealType v; + RealType nc; + bool find_v1; + RealType p; + bool comp; + }; + + template + RealType large_v2_approximation(RealType x, RealType v1, RealType p, RealType q, RealType nc) + { // For v2 -> inf approximate f_degreese_of_freedome_finder with chi squared distribution + // with degrees of freedom v1 at the cdf at x * v1 + bool comp = p < q ? false : true; + RealType pval = p < q ? p : q; + non_central_chi_squared_distribution d(v1, nc); + return comp ? pval - cdf(complement(d, x*v1)) : cdf(d, x*v1) - pval; + } + + template + inline RealType find_degrees_of_freedom_f( + const RealType x, const RealType v, const RealType nc, const bool find_v1, const RealType p, const RealType q, const Policy& pol) + { + BOOST_MATH_STD_USING + using std::fabs; + const char* function = "non_central_f<%1%>::find_degrees_of_freedom_f"; + if((p == 0) || (q == 0)) + { + // + // Can't find a thing if one of p and q is zero: + // + return policies::raise_domain_error(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } + if (x < tools::epsilon()) + { + return policies::raise_evaluation_error(function, "Can't find degrees of freedom when the abscissa value is very close to zero as all degrees of freedom generate the same CDF at x=0: try again further out in the tails!!", + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } + f_degrees_of_freedom_finder f(x, v, nc, find_v1, p < q ? p : q, p < q ? false : true); + + // There are times when f has two roots - thus, two degrees of freedom can + // be found. We find this case by checking if the sign of f for large and + // small values of v have the same sign. If the sign is the same, then there + // are an even number of roots. If the signs differ, there is only one root + // and we can safely find the root. + RealType vLarge = sqrt(boost::math::tools::max_value()); + RealType vSmall = 1 / vLarge; + + // As v2 -> infinity, noncentral f converges to chi-squared distribution. + // Rather than evaluating f(vLarge), we can use the more stable chi-squared approximation + RealType large_difference; + if (find_v1) + { + large_difference = large_v2_approximation(x, v, p, q, nc); + } + else + large_difference = f(vLarge); + + if ((large_difference < 0) == (f(vSmall) < 0)){ + return policies::raise_evaluation_error(function, "Can't find degrees of freedom because two degrees of freedom can be found using the given parameters", + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } + + tools::eps_tolerance tol(policies::digits()); + std::uintmax_t max_iter = policies::get_max_root_iterations(); + // + // Pick an initial guess: + // + RealType guess = 1; + std::pair ir = tools::bracket_and_solve_root( + f, guess, RealType(2), f(guess) < 0 ? true : false, tol, max_iter, pol); + RealType result = ir.first + (ir.second - ir.first) / 2; + if(max_iter >= policies::get_max_root_iterations()) + { + return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE + " or there is no answer to problem. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + } + return result; + } } // namespace detail template > @@ -163,6 +262,96 @@ namespace boost result, function); } + BOOST_MATH_GPU_ENABLED static RealType find_v1(const RealType x, const RealType v2, const RealType nc, const RealType p) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_v1"; + typedef typename policies::evaluation::type eval_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + eval_type result = detail::find_degrees_of_freedom_f( + static_cast(x), + static_cast(v2), + static_cast(nc), + true, + static_cast(p), + static_cast(1-p), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } + template + BOOST_MATH_GPU_ENABLED static RealType find_v1(const complemented4_type& c) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; + typedef typename policies::evaluation::type eval_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + eval_type result = detail::find_degrees_of_freedom_f( + static_cast(c.dist), + static_cast(c.param1), + static_cast(c.param2), + true, + static_cast(1-c.param3), + static_cast(c.param3), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } + BOOST_MATH_GPU_ENABLED static RealType find_v2(const RealType x, const RealType v2, const RealType nc, const RealType p) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_v1"; + typedef typename policies::evaluation::type eval_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + eval_type result = detail::find_degrees_of_freedom_f( + static_cast(x), + static_cast(v2), + static_cast(nc), + false, + static_cast(p), + static_cast(1-p), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } + template + BOOST_MATH_GPU_ENABLED static RealType find_v2(const complemented4_type& c) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; + typedef typename policies::evaluation::type eval_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + eval_type result = detail::find_degrees_of_freedom_f( + static_cast(c.dist), + static_cast(c.param1), + static_cast(c.param2), + false, + static_cast(1-c.param3), + static_cast(c.param3), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } private: // Data member, initialized by constructor. RealType v1; // alpha. diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 39f20d6e41..859118768e 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -19,9 +19,11 @@ #include #include // for ibeta(a, b, x). #include +#include #include #include -#include +#include +#include #include #ifdef _MSC_VER @@ -56,7 +58,11 @@ class students_t_distribution RealType alpha, RealType beta, RealType sd, - RealType hint = 100); + RealType hint = 0); + + BOOST_MATH_GPU_ENABLED static RealType find_degrees_of_freedom( + RealType t, + RealType p); private: // Data member: @@ -123,6 +129,11 @@ BOOST_MATH_GPU_ENABLED inline RealType pdf(const students_t_distribution sqrt(tools::max_value())) + { + // Exact limit: the pdf underflows to zero long before x * x can overflow. + return 0; + } RealType basem1 = x * x / df; if(basem1 < 0.125) { @@ -140,6 +151,8 @@ BOOST_MATH_GPU_ENABLED inline RealType pdf(const students_t_distribution BOOST_MATH_GPU_ENABLED inline RealType cdf(const students_t_distribution& dist, const RealType& x) { + BOOST_MATH_STD_USING // for ADL of std functions + RealType error_result; // degrees_of_freedom > 0 or infinity check: RealType df = dist.degrees_of_freedom(); @@ -195,6 +208,11 @@ BOOST_MATH_GPU_ENABLED inline RealType cdf(const students_t_distribution sqrt(tools::max_value())) + { + // Exact tail values: avoids a spurious intermediate overflow in x * x below. + return (x < 0) ? static_cast(0) : static_cast(1); + } RealType x2 = x * x; RealType probability; if(df > 2 * x2) @@ -281,6 +299,18 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type +BOOST_MATH_GPU_ENABLED inline bool analytical_df_if_cdf_matches(const Distribution& dist, RealType x, RealType p) +{ + RealType cdf_val = cdf(dist, x); + return epsilon_difference(cdf_val, p) < 4; +} + +// Minimum degrees-of-freedom used as the warm-start fallback when the +// Edgeworth approximation yields no valid positive root or is inaccurate +constexpr double df_hint_fallback = 0.01; + // // Functors for finding degrees of freedom: // @@ -297,8 +327,31 @@ struct sample_size_func return 1; } students_t_distribution t(df); - RealType qa = quantile(complement(t, alpha)); - RealType qb = quantile(complement(t, beta)); + RealType qa, qb; +#ifndef BOOST_MATH_NO_EXCEPTIONS + try { +#endif + qa = quantile(complement(t, alpha)); +#ifndef BOOST_MATH_NO_EXCEPTIONS + } + catch (const std::overflow_error&) + { + // Any arbitrary large value will do, must be positive since we calculate qa * qa below: + return tools::max_value(); + } +#endif +#ifndef BOOST_MATH_NO_EXCEPTIONS + try { +#endif + qb = quantile(complement(t, beta)); +#ifndef BOOST_MATH_NO_EXCEPTIONS + } + catch (const std::overflow_error&) + { + // Any arbitrary large value will do, will be negative if beta < 0: + return beta < 0 ? -tools::max_value() : tools::max_value(); + } +#endif qa += qb; qa *= qa; qa *= ratio; @@ -308,6 +361,117 @@ struct sample_size_func RealType alpha, beta, ratio; }; +template +struct cdf_to_df_func +{ + BOOST_MATH_GPU_ENABLED cdf_to_df_func(RealType x_val, RealType p_val, bool c) + : x(x_val), p(p_val), comp(c) {} + + BOOST_MATH_GPU_ENABLED RealType operator()(const RealType& df) + { + students_t_distribution t(df); + return comp ? + RealType(p - cdf(complement(t, x))) + : RealType(cdf(t, x) - p); + } + RealType x, p; + bool comp; +}; + +// +// Shared root-finding helper used by find_degrees_of_freedom. +// +template +BOOST_MATH_GPU_ENABLED RealType solve_for_degrees_of_freedom( + Func f, + RealType hint, + bool rising, + const char* function, + Policy const&) +{ + tools::eps_tolerance tol(policies::digits()); + boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); + boost::math::pair r = tools::bracket_and_solve_root( + f, hint, RealType(2), rising, tol, max_iter, Policy()); + RealType result = r.first + (r.second - r.first) / 2; + if (max_iter >= policies::get_max_root_iterations()) + { + return policies::raise_evaluation_error(function, // LCOV_EXCL_LINE + "Unable to locate solution in a reasonable time: either there is no answer to " // LCOV_EXCL_LINE + "the degrees of freedom or the answer is infinite. Current best guess is %1%", // LCOV_EXCL_LINE + result, Policy()); // LCOV_EXCL_LINE + } + return result; +} + +// +// Edgeworth warm-start for find_degrees_of_freedom(t, p). +// On success writes a df estimate into 'result' and returns true. +// Returns false when no positive root is found; 'result' is left unchanged. +// +template +BOOST_MATH_GPU_ENABLED bool approximate_df_with_edgeworth_expansion(RealType x, RealType p, RealType& result) +{ + BOOST_MATH_STD_USING + // F(x; nu) ~ cdf_normal(x) - pdf_normal(x)*(x + x^3)/(4*nu) + pdf_normal(x)*(3x + 5x^3 + 7x^5 - 3x^7)/(96*nu^2) + // Substituting u = 1/nu and setting F(x; nu) = p gives b*u^2 - a*u + c = 0, where: + // a = pdf_normal(x) * (x + x^3) / 4 + // b = pdf_normal(x) * (3x + 5x^3 + 7x^5 - 3x^7) / 96 + // c = cdf_normal(x) - p + normal_distribution std_normal(0, 1); + RealType pdf_val = pdf(std_normal, x); + RealType cdf_val = cdf(std_normal, x); + + RealType x2 = x * x; + RealType x3 = x2 * x; + RealType x5 = x3 * x2; + RealType x7 = x5 * x2; + + RealType a = pdf_val * (x + x3) / 4; + RealType b = pdf_val * (3 * x + 5 * x3 + 7 * x5 - 3 * x7) / 96; + RealType c = cdf_val - p; + + RealType discriminant = a * a - 4 * b * c; + if (discriminant >= 0 && b != 0) + { + RealType sqrt_disc = sqrt(discriminant); + // The two roots of b*u^2 - a*u + c = 0. Pick the smallest positive u (= largest df = 1/u). + RealType u = (a - sqrt_disc) / (2 * b); + if (u <= 0) + u = (a + sqrt_disc) / (2 * b); + if (u > 0) + { + result = 1 / u; + return true; + } + } + return false; +} + +template +BOOST_MATH_GPU_ENABLED RealType calculate_hill_df_guess( + RealType difference_from_mean, + RealType alpha, + RealType beta, + RealType sd) + + // Hill-corrected normal approximation used as an initial guess + // for find_degrees_of_freedom(difference_from_mean, alpha, beta, sd). + // t_p,nu approx z_p * (1 + (z_p^2 + 1) / (4*nu)) + // nu + 1 = (sd/diff)^2 * (t_alpha + t_beta)^2 + +{ + BOOST_MATH_STD_USING + normal_distribution n(0, 1); + RealType za = quantile(complement(n, alpha)); + RealType zb = quantile(complement(n, beta)); + + RealType v_norm = pow((za + zb) * sd / difference_from_mean, 2) - 1; + RealType hint = v_norm + (za * za - za * zb + zb * zb + 1) / 2; + + return (hint <= 0) ? RealType(1) : hint; +} + } // namespace detail template @@ -329,19 +493,64 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::find_ return error_result; if(hint <= 0) - hint = 1; + hint = detail::calculate_hill_df_guess(difference_from_mean, alpha, beta, sd); detail::sample_size_func f(alpha, beta, sd, difference_from_mean); - tools::eps_tolerance tol(policies::digits()); - boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); - boost::math::pair r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy()); - RealType result = r.first + (r.second - r.first) / 2; - if(max_iter >= policies::get_max_root_iterations()) + return detail::solve_for_degrees_of_freedom(f, hint, false, function, Policy()); +} + +template +BOOST_MATH_GPU_ENABLED RealType students_t_distribution::find_degrees_of_freedom( + RealType t, + RealType p) +{ + BOOST_MATH_STD_USING // for ADL of std functions + constexpr auto function = "boost::math::students_t_distribution<%1%>::find_degrees_of_freedom"; + // + // Check for domain errors: + // + RealType error_result; + if (false == detail::check_probability(function, p, &error_result, Policy())) + return error_result; + + if (t == 0) { - return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time: either there is no answer to how many degrees of freedom are required" // LCOV_EXCL_LINE - " or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + // CDF(0; df) = 0.5 for all df; only solvable when p == 0.5. + if (p == static_cast(0.5)) + return policies::raise_overflow_error(function, nullptr, Policy()); + return policies::raise_domain_error( + function, + "No degrees of freedom can satisfy CDF(0; df) == %1% (must be 0.5).", + p, Policy()); } - return result; + + // Analytical cases: df = infinity (normal), df = 1 (cauchy) + boost::math::normal_distribution norm(0, 1); + if (detail::analytical_df_if_cdf_matches(norm, t, p)) + return policies::raise_overflow_error(function, nullptr, Policy()); + boost::math::cauchy_distribution cauchy(0, 1); + if (detail::analytical_df_if_cdf_matches(cauchy, t, p)) + return RealType(1); + + // Edgeworth warm start: compute a df estimate; fall back to df_hint_fallback if it fails + // or is too inaccurate + RealType hint = static_cast(detail::df_hint_fallback); + if (detail::approximate_df_with_edgeworth_expansion(t, p, hint)) + { + // Check that approximation is at least somewhat close; + // for small degrees of freedom it does not fail but is very inaccurate. + RealType relative_error_threshold = static_cast(0.1); + students_t_distribution t_approx(hint); + RealType p_approx = cdf(t_approx, t); + if (relative_difference(p_approx, p) > relative_error_threshold) + hint = static_cast(detail::df_hint_fallback); + } + // Root-find on f(df) = CDF(t; df) - p. + // CDF(t; df) is strictly increasing in df for t > 0, decreasing for t < 0. + // Use the smaller of p and q=1-p to avoid cancellation near 1. + RealType q = 1 - p; + detail::cdf_to_df_func f(t, p < q ? p : q, p < q ? false : true); + return detail::solve_for_degrees_of_freedom(f, hint, t > 0, function, Policy()); } template diff --git a/include/boost/math/interpolators/pchip.hpp b/include/boost/math/interpolators/pchip.hpp index 6697166df0..a4ad5a3579 100644 --- a/include/boost/math/interpolators/pchip.hpp +++ b/include/boost/math/interpolators/pchip.hpp @@ -121,6 +121,11 @@ class pchip { } } + std::pair domain() const + { + return impl_->domain(); + } + private: std::shared_ptr> impl_; }; diff --git a/include/boost/math/octonion.hpp b/include/boost/math/octonion.hpp index 5d147df584..4239140564 100644 --- a/include/boost/math/octonion.hpp +++ b/include/boost/math/octonion.hpp @@ -942,6 +942,10 @@ namespace boost b /= rhs; \ c /= rhs; \ d /= rhs; \ + e /= rhs; \ + f /= rhs; \ + g /= rhs; \ + h /= rhs; \ \ return(*this); \ } diff --git a/include/boost/math/policies/error_handling.hpp b/include/boost/math/policies/error_handling.hpp index d7ba694564..adcb56f3e5 100644 --- a/include/boost/math/policies/error_handling.hpp +++ b/include/boost/math/policies/error_handling.hpp @@ -26,9 +26,9 @@ #include #include #include +#ifndef BOOST_MATH_NO_EXCEPTIONS #include #include -#ifndef BOOST_MATH_NO_EXCEPTIONS #include #endif #include @@ -95,6 +95,8 @@ T user_indeterminate_result_error(const char* function, const char* message, con namespace detail { +#ifndef BOOST_MATH_NO_EXCEPTIONS +// Only used to build exception messages; relies on template inline std::string prec_format(const T& val) { @@ -124,7 +126,9 @@ inline std::string prec_format(const std::float128_t& val) return std::string(buffer, r.ptr); } -#endif +#endif // BOOST_MATH_USE_CHARCONV_FOR_CONVERSION + +#endif // BOOST_MATH_NO_EXCEPTIONS inline void replace_all_in_string(std::string& result, const char* what, const char* with) { diff --git a/include/boost/math/policies/policy.hpp b/include/boost/math/policies/policy.hpp index ccc33bb9f8..6d80b2b292 100644 --- a/include/boost/math/policies/policy.hpp +++ b/include/boost/math/policies/policy.hpp @@ -766,6 +766,9 @@ struct evaluation using type = typename boost::math::conditional::type; }; +template +using evaluation_t = typename evaluation::type; + template struct precision { diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index f6327fa220..5c3a762ca0 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -19,7 +19,6 @@ #include #include #include -#include #include #include #include @@ -475,7 +474,6 @@ BOOST_MATH_GPU_ENABLED T ibeta_power_terms(T a, const char* = "boost::math::ibeta<%1%>(%1%, %1%, %1%)") { BOOST_MATH_STD_USING - if(!normalised) { return prefix * pow(x, a) * pow(y, b); @@ -1142,7 +1140,20 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no T x0 = a / (a + b); T y0 = b / (a + b); - T nu = x0 * log(x / x0) + y0 * log(y / y0); + + // Expand nu about x0 + T nu = 0; + for (int i=2; i<5; i++) + { + nu += pow(x-x0, i) / i * (pow(x0, -(i-1)) - pow(x0-1, -(i-1))) * pow(-1, i+1); + } + // Calculate the next term in the series + T remainder = pow(x-x0, 5) / 5 * (pow(x0, -4) - pow(x0-1, 4)); + + // If the remainder is large, then fall back to using the log formula + if (remainder >= tools::forth_root_epsilon()){ + nu = x0 * log(x / x0) + y0 * log(y / y0); + } // // Above compution is unstable, force nu to zero if // something went wrong: @@ -1181,7 +1192,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no T mul = 1; if (!normalised) mul = boost::math::beta(a, b, pol); - + // T log_erf_remainder = -0.5 * log(2 * constants::pi() * (a+b)) + a * log(x / x0) + b * log((1-x) / (1-x0)) + log(abs(1/nu - sqrt(x0 * (1-x0)) / (x-x0))); return mul * ((invert ? (1 + boost::math::erf(-nu * sqrt((a + b) / 2), pol)) / 2 : boost::math::erfc(-nu * sqrt((a + b) / 2), pol) / 2)); } @@ -1571,37 +1582,47 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b else { // a and b both large: - bool use_asym = false; T ma = BOOST_MATH_GPU_SAFE_MAX(a, b); - T xa = ma == a ? x : y; T saddle = ma / (a + b); - T powers = 0; - if ((ma > 1e-5f / tools::epsilon()) && (ma / BOOST_MATH_GPU_SAFE_MIN(a, b) < (xa < saddle ? 2 : 15))) + // If ma large and x is close to saddle, use `erf` approximation in `ibeta_large_ab`. + // In this case we don't need to invert + if ((ma > 1e-5f / tools::epsilon()) && (fabs(saddle - x) < 1e-12) && (1 - saddle > 0.125)) { - if (a == b) - use_asym = true; - else - { - powers = exp(log(x / (a / (a + b))) * a + log(y / (b / (a + b))) * b); - if (powers < tools::epsilon()) - use_asym = true; - } + fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol); + invert = false; } - if(use_asym) + else { - fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol); - if (fract * tools::epsilon() < powers) + // This is the same logic that is in `ibeta_fraction2`. We have repeated + // the implementation here because when x is close to saddle, the continued + // fraction will not converge. If this occurs, we fall back to the `erf` + // approximation. Simply using `ibeta_fraction2` will cause an evaluation + // error to be thrown and we won't be able to use the `erf` approximation. + typedef typename lanczos::lanczos::type lanczos_type; + T local_result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol); + if (p_derivative) { - // Erf approximation failed, correction term is too large, fall back: - fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); + *p_derivative = local_result; + BOOST_MATH_ASSERT(*p_derivative >= 0); + } + if (local_result != 0) + { + ibeta_fraction2_t f(a, b, x, y); + boost::math::uintmax_t max_terms = boost::math::policies::get_max_series_iterations(); + T local_fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon(), max_terms); + if (max_terms >= boost::math::policies::get_max_series_iterations()) + { + // Continued fraction failed, fall back to asymptotic expansion: + fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol); + invert = false; + } + else{ + fract = local_result / local_fract; + } } else - invert = false; + fract = 0; } - else - fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); - - BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } if(p_derivative) diff --git a/include/boost/math/special_functions/detail/fp_traits.hpp b/include/boost/math/special_functions/detail/fp_traits.hpp index 7cdb12b188..7362d0fada 100644 --- a/include/boost/math/special_functions/detail/fp_traits.hpp +++ b/include/boost/math/special_functions/detail/fp_traits.hpp @@ -270,13 +270,7 @@ template<> struct fp_traits_non_native // long double (64 bits) ------------------------------------------------------- -#if defined(BOOST_NO_INT64_T) || defined(BOOST_NO_INCLASS_MEMBER_INITIALIZATION)\ - || defined(BOOST_BORLANDC) || defined(__CODEGEAR__) || (defined(__APPLE__) && defined(__aarch64__)) || defined(_MSC_VER)\ - || (defined(__GNUC__) && defined(__aarch64__) && defined(_WIN32))\ - || (defined(__GNUC__) && (defined(__arm__) || defined(__thumb__)))\ - || defined(__SYCL_DEVICE_ONLY__) || (LDBL_MANT_DIG == 53) - -static_assert(LDBL_MANT_DIG == 53, "Oops, assumption that long double is a 64-bit quantity is incorrect!!"); +#if LDBL_MANT_DIG == 53 || defined(__SYCL_DEVICE_ONLY__) template<> struct fp_traits_non_native { @@ -305,14 +299,10 @@ template<> struct fp_traits_non_native // long double (>64 bits), x86 and x64 ----------------------------------------- -#elif defined(__i386) || defined(__i386__) || defined(_M_IX86) \ - || defined(__amd64) || defined(__amd64__) || defined(_M_AMD64) \ - || defined(__x86_64) || defined(__x86_64__) || defined(_M_X64) +#elif LDBL_MANT_DIG == 64 // Intel extended double precision format (80 bits) -static_assert(LDBL_MANT_DIG == 64, "Oops, assumption that long double is an 80-bit quantity is incorrect!!"); - template<> struct fp_traits_non_native { @@ -440,8 +430,6 @@ struct fp_traits_non_native // IEEE extended double precision format with 15 exponent bits (128 bits) -static_assert(LDBL_MANT_DIG == 113, "Oops, assumption that long double is a 128-bit quantity is incorrect!!"); - template<> struct fp_traits_non_native { diff --git a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp index 3adc86f74f..895c586482 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp @@ -23,6 +23,7 @@ unsigned set_crossover_locations(const Seq& aj, const Seq& bj, const Real& z, unsigned int* crossover_locations) { BOOST_MATH_STD_USING + using nothrow_policy = typename boost::math::policies::normalise, boost::math::policies::rounding_error>::type; unsigned N_terms = 0; if(aj.size() == 1 && bj.size() == 1) @@ -55,13 +56,13 @@ Real t = (-sqrt(sq) - b + z) / 2; if (t >= 0) { - crossover_locations[N_terms] = itrunc(t); + crossover_locations[N_terms] = itrunc(t, nothrow_policy()); ++N_terms; } t = (sqrt(sq) - b + z) / 2; if (t >= 0) { - crossover_locations[N_terms] = itrunc(t); + crossover_locations[N_terms] = itrunc(t, nothrow_policy()); ++N_terms; } } @@ -71,13 +72,13 @@ Real t = (-sqrt(sq) - b - z) / 2; if (t >= 0) { - crossover_locations[N_terms] = itrunc(t); + crossover_locations[N_terms] = itrunc(t, nothrow_policy()); ++N_terms; } t = (sqrt(sq) - b - z) / 2; if (t >= 0) { - crossover_locations[N_terms] = itrunc(t); + crossover_locations[N_terms] = itrunc(t, nothrow_policy()); ++N_terms; } } @@ -110,7 +111,7 @@ unsigned n = 0; for (auto bi = bj.begin(); bi != bj.end(); ++bi, ++n) { - crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi) + 1; + crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi, nothrow_policy()) + 1; } std::sort(crossover_locations, crossover_locations + bj.size(), std::less()); N_terms = (unsigned)bj.size(); @@ -122,6 +123,7 @@ std::pair hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, const Terminal& termination, long long& log_scale) { BOOST_MATH_STD_USING + using nothrow_policy = typename boost::math::policies::normalise>::type; Real result = 1; Real abs_result = 1; Real term = 1; @@ -129,8 +131,19 @@ Real tol = boost::math::policies::get_epsilon(); std::uintmax_t k = 0; Real upper_limit(sqrt(boost::math::tools::max_value())), diff; + if ((tools::max_value() / fabs(z) < upper_limit)) + { + upper_limit = tools::max_value() / fabs(z); + } + for (auto pa = aj.begin(); pa != aj.end(); ++pa) + { + if (tools::max_value() / fabs(*pa) < upper_limit) + { + upper_limit = tools::max_value() / fabs(*pa); + } + } Real lower_limit(1 / upper_limit); - long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value()) - 2; + long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value(), nothrow_policy()) - 2; Real scaling_factor = exp(Real(log_scaling_factor)); Real term_m1; long long local_scaling = 0; @@ -234,7 +247,10 @@ log_scale += log_scaling_factor; local_scaling += log_scaling_factor; } - if (fabs(abs_result) < lower_limit) + /* + * rescaling to avoid underflow is pointless here, given that the first term is 1. + * + else if ((fabs(abs_result) < lower_limit) && (fabs(abs_result) * scaling_factor < upper_limit)) { abs_result *= scaling_factor; result *= scaling_factor; @@ -242,7 +258,7 @@ log_scale -= log_scaling_factor; local_scaling -= log_scaling_factor; } - + */ if ((abs(result * tol) > abs(term)) && (abs(term0) > abs(term))) break; if (abs_result * tol > abs(result)) @@ -321,7 +337,7 @@ else { int ls = 1; - Real p = log_pochhammer(*ai, s, pol, &ls); + Real p = log_pochhammer(static_cast(*ai), s, pol, &ls); s1 *= ls; term += p; loop_error_scale = (std::max)(p, loop_error_scale); @@ -334,7 +350,7 @@ for (auto bi = bj.begin(); bi != bj.end(); ++bi) { int ls = 1; - Real p = log_pochhammer(*bi, s, pol, &ls); + Real p = log_pochhammer(static_cast(*bi), s, pol, &ls); s2 *= ls; term -= p; loop_error_scale = (std::max)(p, loop_error_scale); @@ -371,13 +387,13 @@ if (term <= tools::log_min_value()) { // rescale if we can: - long long scale = lltrunc(floor(term - tools::log_min_value()) - 2); + long long scale = lltrunc(floor(term - tools::log_min_value()) - 2, nothrow_policy()); term -= scale; loop_scale += scale; } if (term > 10) { - int scale = itrunc(floor(term)); + long long scale = lltrunc(floor(term), nothrow_policy()); term -= scale; loop_scale += scale; } @@ -402,7 +418,7 @@ term /= scaling_factor; loop_scale += log_scaling_factor; } - if (fabs(loop_result) < lower_limit) + else if ((fabs(loop_result) < lower_limit) && (fabs(loop_result) * scaling_factor < upper_limit)) { loop_result *= scaling_factor; loop_abs_result *= scaling_factor; diff --git a/include/boost/math/special_functions/ellint_2.hpp b/include/boost/math/special_functions/ellint_2.hpp index 0cc1fa0944..501d36b66b 100644 --- a/include/boost/math/special_functions/ellint_2.hpp +++ b/include/boost/math/special_functions/ellint_2.hpp @@ -122,7 +122,7 @@ BOOST_MATH_GPU_ENABLED T ellint_e_imp(T phi, T k, const Policy& pol) rphi = constants::half_pi() - rphi; } T k2 = k * k; - if(boost::math::pow<3>(rphi) * k2 / 6 < tools::epsilon() * fabs(rphi)) + if(boost::math::pow<3>(rphi) * k2 / 6 <= tools::epsilon() * fabs(rphi)) { // See http://functions.wolfram.com/EllipticIntegrals/EllipticE2/06/01/03/0001/ result = s * rphi; diff --git a/include/boost/math/special_functions/hypergeometric_pFq.hpp b/include/boost/math/special_functions/hypergeometric_pFq.hpp index 070b852dd2..96f68e510e 100644 --- a/include/boost/math/special_functions/hypergeometric_pFq.hpp +++ b/include/boost/math/special_functions/hypergeometric_pFq.hpp @@ -58,12 +58,22 @@ namespace boost { BOOST_MATH_STD_USING long long scale = 0; + static const char* function = "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)"; std::pair r = boost::math::detail::hypergeometric_pFq_checked_series_impl(aj, bj, value_type(z), pol, boost::math::detail::iteration_terminator(boost::math::policies::get_max_series_iterations()), scale); - r.first *= exp(Real(scale)); - r.second *= exp(Real(scale)); + // + // Overflow check: + // + if (static_cast(scale) > tools::log_max_value()) + return (r.first < 0 ? -1 : 1) * policies::raise_overflow_error(function, nullptr, pol); + Real mul = exp(Real(scale)); + if(fabs(r.first) > 1) + if(tools::max_value() / fabs(r.first) < mul) + return (r.first < 0 ? -1 : 1) * policies::raise_overflow_error(function, nullptr, pol); + r.first *= mul; + r.second *= mul; if (p_abs_error) *p_abs_error = static_cast(r.second) * boost::math::tools::epsilon(); - return policies::checked_narrowing_cast(r.first, "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)"); + return policies::checked_narrowing_cast(r.first, function); } template diff --git a/include/boost/math/special_functions/lanczos.hpp b/include/boost/math/special_functions/lanczos.hpp index 0ec24bddbf..409e14c5e1 100644 --- a/include/boost/math/special_functions/lanczos.hpp +++ b/include/boost/math/special_functions/lanczos.hpp @@ -2751,7 +2751,7 @@ struct lanczos } // namespace math } // namespace boost -#if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) +#if !defined(_CRAYC) && !defined(BOOST_MATH_ENABLE_CUDA) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) #if ((defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__) || defined(_M_AMD64) || defined(_M_X64)) && !defined(_MANAGED) && !defined(BOOST_MATH_HAS_GPU_SUPPORT) #include #endif diff --git a/include/boost/math/special_functions/next.hpp b/include/boost/math/special_functions/next.hpp index 74c34f06ad..f73cb2f6a5 100644 --- a/include/boost/math/special_functions/next.hpp +++ b/include/boost/math/special_functions/next.hpp @@ -25,7 +25,7 @@ #include -#if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) +#if !defined(_CRAYC) && !defined(BOOST_MATH_ENABLE_CUDA) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) #if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__) #include "xmmintrin.h" #define BOOST_MATH_CHECK_SSE2 diff --git a/include/boost/math/special_functions/sign.hpp b/include/boost/math/special_functions/sign.hpp index 4f76522654..5d2dfc2e23 100644 --- a/include/boost/math/special_functions/sign.hpp +++ b/include/boost/math/special_functions/sign.hpp @@ -14,7 +14,7 @@ #pragma once #endif -#ifndef __CUDACC_RTC__ +#ifndef BOOST_MATH_HAS_NVRTC #include #include @@ -234,7 +234,7 @@ BOOST_MATH_GPU_ENABLED T sign(T z) } // namespace math } // namespace boost -#endif // __CUDACC_RTC__ +#endif // BOOST_MATH_HAS_NVRTC #endif // BOOST_MATH_TOOLS_SIGN_HPP diff --git a/include/boost/math/tools/config.hpp b/include/boost/math/tools/config.hpp index 4829ebe36e..40346c9b0f 100644 --- a/include/boost/math/tools/config.hpp +++ b/include/boost/math/tools/config.hpp @@ -11,7 +11,7 @@ #pragma once #endif -#ifndef __CUDACC_RTC__ +#if !(defined(__CUDACC_RTC__) && defined(BOOST_MATH_ENABLE_NVRTC)) #include @@ -168,7 +168,7 @@ # define BOOST_MATH_NOINLINE __declspec(noinline) # elif defined(__GNUC__) && __GNUC__ > 3 // Clang also defines __GNUC__ (as 4) -# if defined(__CUDACC__) +# if defined(__CUDACC__) && defined(BOOST_MATH_ENABLE_CUDA) // nvcc doesn't always parse __noinline__, // see: https://svn.boost.org/trac/boost/ticket/9392 # define BOOST_MATH_NOINLINE __attribute__ ((noinline)) @@ -678,7 +678,7 @@ namespace boost{ namespace math{ // CUDA support: // -#ifdef __CUDACC__ +#if defined(__CUDACC__) && defined(BOOST_MATH_ENABLE_CUDA) // We have to get our include order correct otherwise you get compilation failures #include @@ -774,7 +774,7 @@ BOOST_MATH_GPU_ENABLED constexpr T gpu_safe_max(const T& a, const T& b) { return # define BOOST_MATH_STATIC_LOCAL_VARIABLE # else # define BOOST_MATH_INLINE_CONSTEXPR constexpr -# define BOOST_MATH_STATIC constexpr +# define BOOST_MATH_STATIC static # define BOOST_MATH_STATIC_LOCAL_VARIABLE static # endif #endif diff --git a/include/boost/math/tools/cubic_roots.hpp b/include/boost/math/tools/cubic_roots.hpp index 451c5de813..1097ad5ce4 100644 --- a/include/boost/math/tools/cubic_roots.hpp +++ b/include/boost/math/tools/cubic_roots.hpp @@ -6,11 +6,33 @@ #define BOOST_MATH_TOOLS_CUBIC_ROOTS_HPP #include #include +#include #include #include namespace boost::math::tools { +namespace detail { + +// Orders finite roots ascending and moves NaN entries to the back. +// Sorting NaNs is UB when using vanilla operator< (comparisons to NaN are always false), +// so we need to provide our own sort function that pushes NaN to the end +template +bool roots_less(const Real& lhs, const Real& rhs) +{ + if ((boost::math::isnan)(lhs)) + { + return false; + } + if ((boost::math::isnan)(rhs)) + { + return true; + } + return lhs < rhs; +} + +} // namespace detail + // Solves ax^3 + bx^2 + cx + d = 0. // Only returns the real roots, as types get weird for real coefficients and // complex roots. Follows Numerical Recipes, Chapter 5, section 6. NB: A better @@ -55,7 +77,7 @@ std::array cubic_roots(Real a, Real b, Real c, Real d) { roots[0] = x0; roots[1] = x1; roots[2] = 0; - std::sort(roots.begin(), roots.end()); + std::sort(roots.begin(), roots.end(), detail::roots_less); return roots; } Real p = b / a; @@ -115,7 +137,7 @@ std::array cubic_roots(Real a, Real b, Real c, Real d) { } } } - std::sort(roots.begin(), roots.end()); + std::sort(roots.begin(), roots.end(), detail::roots_less); return roots; } diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index b0b0fc246c..af9b76d998 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -456,8 +456,13 @@ namespace detail { if (count) { max = guess; + // + // We can't have this extra recursive tidy up step on CUDA: + // +#if !defined(BOOST_MATH_ENABLE_CUDA) && !defined(BOOST_MATH_ENABLE_SYCL) if (multiplier > 16) return (guess0 - guess) + bracket_root_towards_min(f, guess, f_current, min, max, count); +#endif } return guess0 - (max + min) / 2; } @@ -520,8 +525,13 @@ namespace detail { if (count) { min = guess; + // + // We can't have this extra recursive tidy up step on CUDA: + // +#if !defined(BOOST_MATH_ENABLE_CUDA) && !defined(BOOST_MATH_ENABLE_SYCL) if (multiplier > 16) return (guess0 - guess) + bracket_root_towards_max(f, guess, f_current, min, max, count); +#endif } return guess0 - (max + min) / 2; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4715cf379b..9d673adc7d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -10,13 +10,13 @@ if(HAVE_BOOST_TEST) message(STATUS "Building boost.math with CUDA") - find_package(CUDA REQUIRED) enable_language(CUDA) + find_package(CUDAToolkit REQUIRED) set(CMAKE_CUDA_EXTENSIONS OFF) enable_testing() - boost_test_jamfile(FILE cuda_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception Boost::unit_test_framework ${CUDA_LIBRARIES} INCLUDE_DIRECTORIES ${CUDA_INCLUDE_DIRS} ) + boost_test_jamfile(FILE cuda_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception Boost::unit_test_framework CUDA::cudart COMPILE_DEFINITIONS BOOST_MATH_ENABLE_CUDA=1 ) elseif (BOOST_MATH_ENABLE_NVRTC) diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index c5bde01635..18d5d9ae01 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -202,6 +202,8 @@ test-suite special_fun : [ run git_issue_1255.cpp ] [ run git_issue_1247.cpp ] [ run git_issue_1308.cpp /boost/test//boost_unit_test_framework ] + [ run git_issue_1412.cpp ] + [ run git_issue_1412_pt_2.cpp ] [ run special_functions_test.cpp /boost/test//boost_unit_test_framework ] [ run test_airy.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] [ run test_bessel_j.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] @@ -895,6 +897,7 @@ test-suite distribution_tests : : test_poisson_real_concept ] [ run test_rayleigh.cpp /boost/test//boost_unit_test_framework ] [ run test_students_t.cpp /boost/test//boost_unit_test_framework ] + [ run test_students_t_overflow.cpp /boost/test//boost_unit_test_framework ] [ run test_skew_normal.cpp /boost/test//boost_unit_test_framework ] [ run test_triangular.cpp pch /boost/test//boost_unit_test_framework ] [ run test_uniform.cpp pch /boost/test//boost_unit_test_framework ] diff --git a/test/cmake_install_test/CMakeLists.txt b/test/cmake_install_test/CMakeLists.txt new file mode 100644 index 0000000000..010d5e9fe3 --- /dev/null +++ b/test/cmake_install_test/CMakeLists.txt @@ -0,0 +1,17 @@ +# Copyright 2025 Matt Borland +# Distributed under the Boost Software License, Version 1.0. +# https://www.boost.org/LICENSE_1_0.txt + +cmake_minimum_required(VERSION 3.8...3.16) + +project(cmake_install_test LANGUAGES CXX) + +find_package(boost_math REQUIRED) + +add_executable(quick ../quick.cpp) +target_link_libraries(quick Boost::math) + +enable_testing() +add_test(quick quick) + +add_custom_target(check COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -C $) diff --git a/test/cmake_install_test_legacy/CMakeLists.txt b/test/cmake_install_test_legacy/CMakeLists.txt new file mode 100644 index 0000000000..b5db787104 --- /dev/null +++ b/test/cmake_install_test_legacy/CMakeLists.txt @@ -0,0 +1,19 @@ +# Copyright 2025 Matt Borland +# Distributed under the Boost Software License, Version 1.0. +# https://www.boost.org/LICENSE_1_0.txt + +cmake_minimum_required(VERSION 3.8...3.16) + +project(cmake_install_test_legacy LANGUAGES CXX) + +find_package(boost_math REQUIRED) +find_package(boost_math_c99 REQUIRED) +find_package(boost_math_tr1 REQUIRED) + +add_executable(quick_tr1 ../quick_tr1.cpp) +target_link_libraries(quick_tr1 Boost::math Boost::math_c99 Boost::math_tr1) + +enable_testing() +add_test(quick_tr1 quick_tr1) + +add_custom_target(check COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -C $) diff --git a/test/cmake_subdir_test/CMakeLists.txt b/test/cmake_subdir_test/CMakeLists.txt new file mode 100644 index 0000000000..09d57446d7 --- /dev/null +++ b/test/cmake_subdir_test/CMakeLists.txt @@ -0,0 +1,18 @@ +# Copyright 2025 Matt Borland +# Distributed under the Boost Software License, Version 1.0. +# https://www.boost.org/LICENSE_1_0.txt + +cmake_minimum_required(VERSION 3.8...3.20) + +project(cmake_subdir_test LANGUAGES CXX) + +# Math is header-only and standalone by default (no deps needed) +add_subdirectory(../.. boostorg/math) + +add_executable(quick ../quick.cpp) +target_link_libraries(quick Boost::math) + +enable_testing() +add_test(quick quick) + +add_custom_target(check COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -C $) diff --git a/test/cmake_subdir_test_legacy/CMakeLists.txt b/test/cmake_subdir_test_legacy/CMakeLists.txt new file mode 100644 index 0000000000..c606705fa0 --- /dev/null +++ b/test/cmake_subdir_test_legacy/CMakeLists.txt @@ -0,0 +1,23 @@ +# Copyright 2025 Matt Borland +# Distributed under the Boost Software License, Version 1.0. +# https://www.boost.org/LICENSE_1_0.txt + +cmake_minimum_required(VERSION 3.8...3.20) + +project(cmake_subdir_test_legacy LANGUAGES CXX) + +# config must be added before math (compiled targets depend on Boost::config) +add_subdirectory(../../../config boostorg/config) + +# Enable legacy compiled targets before adding math +set(BOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS ON CACHE BOOL "" FORCE) + +add_subdirectory(../.. boostorg/math) + +add_executable(quick_tr1 ../quick_tr1.cpp) +target_link_libraries(quick_tr1 Boost::math Boost::math_c99 Boost::math_tr1) + +enable_testing() +add_test(quick_tr1 quick_tr1) + +add_custom_target(check COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -C $) diff --git a/test/cuda_jamfile b/test/cuda_jamfile index e4d7505da8..f2f2c9aa82 100644 --- a/test/cuda_jamfile +++ b/test/cuda_jamfile @@ -9,6 +9,10 @@ project : requirements [ requires cxx14_decltype_auto cxx14_generic_lambdas cxx14_return_type_deduction cxx14_variable_templates cxx14_constexpr ] ; +# Github Issues +run github_issue_1383.cu ; +run github_issue_1383_pt_2.cu ; + # Quad run test_exp_sinh_quad_float.cu ; run test_exp_sinh_quad_double.cu ; diff --git a/test/git_issue_1412.cpp b/test/git_issue_1412.cpp new file mode 100644 index 0000000000..74735e695a --- /dev/null +++ b/test/git_issue_1412.cpp @@ -0,0 +1,23 @@ +// Copyright Matt Borland 2026 +// 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) +// +// Regression test for https://github.com/boostorg/math/issues/1412 + +#define BOOST_MATH_NO_EXCEPTIONS + +// Claim the include guard so any later #include expands to nothing and std::setprecision is unavailable. +// The original reproducer with _SSTREAM broke which iso outside our purview +#define _LIBCPP_IOMANIP // libc++ +#define _GLIBCXX_IOMANIP 1 // libstdc++ + +#include +#include + +int main() +{ + // Exercise a special function so the no-exceptions header path is + // instantiated, not just parsed. + return static_cast(boost::math::laguerre(1u, 0u, 0.5)); +} diff --git a/test/git_issue_1412_pt_2.cpp b/test/git_issue_1412_pt_2.cpp new file mode 100644 index 0000000000..a591fdf346 --- /dev/null +++ b/test/git_issue_1412_pt_2.cpp @@ -0,0 +1,26 @@ +// Copyright Matt Borland 2026 +// 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) +// +// Regression test for https://github.com/boostorg/math/issues/1412 +// +// Part 2: Test all of #include as that is what libc++ is going to use + +#define BOOST_MATH_NO_EXCEPTIONS + +// Claim the include guard so any later #include expands to nothing and std::setprecision is unavailable. +// The original reproducer with _SSTREAM broke which is outside our purview. +#define _LIBCPP_IOMANIP // libc++ +#define _GLIBCXX_IOMANIP 1 // libstdc++ + +#include "compile_test/instantiate.hpp" + +int main(int argc, char* []) +{ + if (argc > 10000) + { + instantiate(double(0)); + instantiate_mixed(double(0)); + } +} diff --git a/test/github_issue_1383.cu b/test/github_issue_1383.cu new file mode 100644 index 0000000000..119f5d1d5e --- /dev/null +++ b/test/github_issue_1383.cu @@ -0,0 +1,120 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024 - 2026. +// 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) + +#ifdef BOOST_MATH_ENABLE_CUDA +# undef BOOST_MATH_ENABLE_CUDA +#endif // BOOST_MATH_ENABLE_CUDA + +// Purposefully pull in headers that caused errors in the linked issue +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +#ifdef BOOST_MATH_ENABLE_CUDA +# error "We should not be enabling this ourselves" +#endif // BOOST_MATH_ENABLE_CUDA + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in, float_type *out, int numElements) +{ + using std::cos; + const int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = cos(in[i]); + } +} + +/** + * Host main routine + */ +int main() +{ + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + // Initialize the input vectors + // Check some of our numeric_limits for viability + std::mt19937_64 rng {42}; + std::uniform_real_distribution dist(0, boost::math::constants::pi()); + static_assert(boost::math::numeric_limits::is_specialized, "Should be since it's a double"); + static_assert(boost::math::numeric_limits::is_signed, "Should be since it's a double"); + + for (int i = 0; i < numElements; ++i) + { + input_vector[i] = dist(rng); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 256; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + + cuda_test<<>>(input_vector.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + + std::cout << "CUDA kernal done in: " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(std::cos(input_vector[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 10) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + + return 0; +} diff --git a/test/github_issue_1383_pt_2.cu b/test/github_issue_1383_pt_2.cu new file mode 100644 index 0000000000..8e746abfce --- /dev/null +++ b/test/github_issue_1383_pt_2.cu @@ -0,0 +1,116 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024 - 2026. +// 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) + +// Purposefully pull in headers that caused errors in the linked issue +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in, float_type *out, int numElements) +{ + using std::cos; + const int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = cos(in[i]); + if (out[i] > boost::math::numeric_limits::max() || !boost::math::numeric_limits::is_signed) + { + __trap(); + } + } +} + +/** + * Host main routine + */ +int main() +{ + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + // Initialize the input vectors + // Check some of our numeric_limits for viability + std::mt19937_64 rng {42}; + std::uniform_real_distribution dist(0, boost::math::constants::pi()); + static_assert(boost::math::numeric_limits::is_specialized, "Should be since it's a double"); + static_assert(boost::math::numeric_limits::is_signed, "Should be since it's a double"); + + for (int i = 0; i < numElements; ++i) + { + input_vector[i] = dist(rng); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 256; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + + cuda_test<<>>(input_vector.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + + std::cout << "CUDA kernal done in: " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(std::cos(input_vector[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 10) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + + return 0; +} diff --git a/test/octonion_test.cpp b/test/octonion_test.cpp index 96f43103c5..00ebf535ff 100644 --- a/test/octonion_test.cpp +++ b/test/octonion_test.cpp @@ -699,6 +699,34 @@ BOOST_TEST_CASE_TEMPLATE_FUNCTION(multiplication_test, T) } } +BOOST_TEST_CASE_TEMPLATE_FUNCTION(division_test, T) +{ + #if BOOST_WORKAROUND(__GNUC__, < 3) + #else /* BOOST_WORKAROUND(__GNUC__, < 3) */ + using ::std::numeric_limits; + + using ::boost::math::abs; + #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ + + + BOOST_TEST_MESSAGE("Testing division for " + << string_type_name::_() << "."); + + BOOST_REQUIRE_PREDICATE(::std::less_equal(), + (abs(::boost::math::octonion(1,0,0,0,0,0,0,0)/ + ::boost::math::octonion(1,0,0,0,0,0,0,0)- + static_cast(1))) + (numeric_limits::epsilon())); + + for (int idx = 1; idx < 8; ++idx) + { + ::boost::math::octonion toto = index_i_element(idx); + + BOOST_REQUIRE_PREDICATE(::std::less_equal(), + (abs(toto/toto-static_cast(1))) + (numeric_limits::epsilon())); + } +} BOOST_TEST_CASE_TEMPLATE_FUNCTION(exp_test, T) { @@ -752,6 +780,7 @@ boost::unit_test::test_suite * init_unit_test_suite(int, char *[]) #define BOOST_OCTONION_TEST \ BOOST_OCTONION_COMMON_GENERATOR(multiplication) \ + BOOST_OCTONION_COMMON_GENERATOR(division) \ BOOST_OCTONION_COMMON_GENERATOR_NEAR_EPS(exp) diff --git a/test/pchip_test.cpp b/test/pchip_test.cpp index 6e8ea235bf..5e8fb3c2ac 100644 --- a/test/pchip_test.cpp +++ b/test/pchip_test.cpp @@ -38,6 +38,10 @@ void test_constant() auto pchip_spline = pchip(std::move(x_copy), std::move(y_copy)); //std::cout << "Constant value pchip spline = " << pchip_spline << "\n"; + auto [lo, hi] = pchip_spline.domain(); + CHECK_ULP_CLOSE(lo, 0, 1); + CHECK_ULP_CLOSE(hi, 81, 1); + for (Real t = x[0]; t <= x.back(); t += Real(0.25)) { CHECK_ULP_CLOSE(Real(7), pchip_spline(t), 2); CHECK_ULP_CLOSE(Real(0), pchip_spline.prime(t), 2); diff --git a/test/quick.cpp b/test/quick.cpp new file mode 100644 index 0000000000..30ac1b896b --- /dev/null +++ b/test/quick.cpp @@ -0,0 +1,17 @@ +// Copyright 2026 Matt Borland +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include +#include + +int main() +{ + auto result = boost::math::tgamma(2.0); + // tgamma(2) == 1! == 1.0 + if (result < 0.99 || result > 1.01) + { + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} diff --git a/test/quick_tr1.cpp b/test/quick_tr1.cpp new file mode 100644 index 0000000000..33e21a1495 --- /dev/null +++ b/test/quick_tr1.cpp @@ -0,0 +1,26 @@ +// Copyright 2026 Matt Borland +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include +#include + +int main() +{ + // boost_cbrt from boost_math_c99 + auto cbrt_result = boost::math::tr1::boost_cbrt(8.0); + if (cbrt_result < 1.99 || cbrt_result > 2.01) + { + return EXIT_FAILURE; + } + + // boost_riemann_zeta from boost_math_tr1 + // zeta(2) == pi^2/6 ~= 1.6449 + auto zeta_result = boost::math::tr1::boost_riemann_zeta(2.0); + if (zeta_result < 1.64 || zeta_result > 1.65) + { + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/test/test_arcsine.cpp b/test/test_arcsine.cpp index eed8dc9481..dc2fba919e 100644 --- a/test/test_arcsine.cpp +++ b/test/test_arcsine.cpp @@ -105,6 +105,7 @@ void test_ignore_policy(RealType) //std::cout << "pdf(ignore_error_arcsine(-1, +1), std::numeric_limits::infinity()) = " << pdf(ignore_error_arcsine(-1, +1), std::numeric_limits::infinity()) << std::endl; // Outputs: pdf(ignore_error_arcsine(-1, +1), std::numeric_limits::infinity()) = 1.#QNAN } + // PDF BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(0, 1), std::numeric_limits::infinity()))); // x == infinity BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(-1, 1), std::numeric_limits::infinity()))); // x == infinity BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(0, 1), static_cast (-2)))); // x < xmin @@ -120,9 +121,31 @@ void test_ignore_policy(RealType) BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(0, 1), static_cast (+2)))); // x > x_max BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(-1, 1), static_cast (+2)))); // x > x_max + // CDF + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(0, 1), std::numeric_limits::infinity()))); // x == infinity + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(-1, 1), std::numeric_limits::infinity()))); // x == infinity + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(0, 1), static_cast (-2)))); // x < xmin + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(-1, 1), static_cast (-2)))); // x < xmin + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(0, 1), static_cast (+2)))); // x > x_max + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(-1, 1), static_cast (+2)))); // x > x_max + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_arcsine(0, 1), std::numeric_limits::infinity())))); // x == infinity + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_arcsine(0, 1), static_cast (-2))))); // x < xmin + + // Quantile + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_arcsine(0, 1), std::numeric_limits::infinity()))); // p == infinity + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_arcsine(0, 1), static_cast(-1)))); // p < 0 + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_arcsine(0, 1), static_cast(2)))); // p > 1 + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_arcsine(0, 1), std::numeric_limits::infinity())))); // q == infinity + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_arcsine(0, 1), static_cast(-1))))); // q < 0 + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_arcsine(0, 1), static_cast(2))))); // q > 1 + // Mean BOOST_CHECK((boost::math::isnan)(mean(ignore_error_arcsine(-nan, 0)))); BOOST_CHECK((boost::math::isnan)(mean(ignore_error_arcsine(+nan, 0)))); + + // Median + BOOST_CHECK((boost::math::isnan)(median(ignore_error_arcsine(-nan, 0)))); + BOOST_CHECK((boost::math::isnan)(median(ignore_error_arcsine(+nan, 0)))); if (std::numeric_limits::has_infinity) { @@ -276,7 +299,8 @@ void test_spots(RealType) BOOST_CHECK_EQUAL(variance(arcsine_01), 0.125); // 1/8 = 0.125 BOOST_CHECK_CLOSE_FRACTION(standard_deviation(arcsine_01), one_div_root_two() / 2, tolerance); // 1/ sqrt(s) = 0.35355339059327379 BOOST_CHECK_EQUAL(skewness(arcsine_01), 0); // - BOOST_CHECK_EQUAL(kurtosis_excess(arcsine_01), -1.5); // 3/2 + BOOST_CHECK_EQUAL(kurtosis_excess(arcsine_01), -1.5); // -3/2 + BOOST_CHECK_EQUAL(kurtosis(arcsine_01), 1.5); // 3/2 BOOST_CHECK_EQUAL(support(arcsine_01).first, 0); // BOOST_CHECK_EQUAL(range(arcsine_01).first, 0); // BOOST_CHECK_THROW(mode(arcsine_01), std::domain_error); // Two modes at x_min and x_max, so throw instead. @@ -374,8 +398,8 @@ void test_spots(RealType) BOOST_CHECK_EQUAL(variance(as_m11), 0.5); // 1 - (-1) = 2 ^ 2 = 4 /8 = 0.5 BOOST_CHECK_EQUAL(skewness(as_m11), 0); // - BOOST_CHECK_EQUAL(kurtosis_excess(as_m11), -1.5); // 3/2 - + BOOST_CHECK_EQUAL(kurtosis_excess(as_m11), -1.5); // -3/2 + BOOST_CHECK_EQUAL(kurtosis(as_m11), 1.5); // 3/2 BOOST_CHECK_CLOSE_FRACTION(pdf(as_m11, 0.05), static_cast(0.31870852113797122803869876869296281629727218095644L), tolerance); BOOST_CHECK_CLOSE_FRACTION(pdf(as_m11, 0.5), static_cast(0.36755259694786136634088433220864629426492432024443L), tolerance); @@ -428,7 +452,8 @@ void test_spots(RealType) BOOST_CHECK_EQUAL(median(as_m2m1), -1.5); // 1 / (1 + 1) = 1/2 exactly. BOOST_CHECK_EQUAL(variance(as_m2m1), 0.125); BOOST_CHECK_EQUAL(skewness(as_m2m1), 0); // - BOOST_CHECK_EQUAL(kurtosis_excess(as_m2m1), -1.5); // 3/2 + BOOST_CHECK_EQUAL(kurtosis(as_m2m1), 1.5); // 3/2 + BOOST_CHECK_EQUAL(kurtosis_excess(as_m2m1), -1.5); // -3/2 BOOST_CHECK_CLOSE_FRACTION(pdf(as_m2m1, -1.95), static_cast(1.4605059227421865250256574657088244053723856445614L), 4 * tolerance); BOOST_CHECK_CLOSE_FRACTION(pdf(as_m2m1, -1.5), static_cast(0.63661977236758134307553505349005744813783858296183L), tolerance); @@ -437,11 +462,19 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION(cdf(as_m2m1, -1.05), static_cast(0.85643370687129372924905811522494428117838480010259L), tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(as_m2m1, -1.5), static_cast(0.5L), tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(as_m2m1, -1.95), static_cast(0.14356629312870627075094188477505571882161519989741L), 8 * tolerance); // Not much less accurate. + BOOST_CHECK_EQUAL(cdf(as_m2m1, as_m2m1.x_min()), 0); + BOOST_CHECK_EQUAL(cdf(as_m2m1, as_m2m1.x_max()), 1); + BOOST_CHECK_EQUAL(cdf(complement(as_m2m1, as_m2m1.x_min())), 1); + BOOST_CHECK_EQUAL(cdf(complement(as_m2m1, as_m2m1.x_max())), 0); // Quantile BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.85643370687129372924905811522494428117838480010259L)), -static_cast(1.05L), 2 * tolerance); // BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.5L)), -static_cast(1.5L), 2 * tolerance); // BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.14356629312870627075094188477505571882161519989741L)), -static_cast(1.95L), 4 * tolerance); // + BOOST_CHECK_EQUAL(quantile(as_m2m1, 1), 1); + BOOST_CHECK_EQUAL(quantile(as_m2m1, 0), 0); + BOOST_CHECK_EQUAL(quantile(complement(as_m2m1, 1)), 0); + BOOST_CHECK_EQUAL(quantile(complement(as_m2m1, 0)), 1); BOOST_CHECK_CLOSE_FRACTION(quantile(complement(as_m2m1, static_cast(0.14356629312870627075094188477505571882161519989741L))), -static_cast(1.05L), 2 * tolerance); // BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.5L)), -static_cast(1.5L), 2 * tolerance); // @@ -607,7 +640,8 @@ void test_spots(RealType) BOOST_CHECK_EQUAL(variance(as), 0.125); //0.125 BOOST_CHECK_CLOSE_FRACTION(standard_deviation(as), one_div_root_two() / 2, std::numeric_limits::epsilon()); // 0.353553 BOOST_CHECK_EQUAL(skewness(as), 0); // - BOOST_CHECK_EQUAL(kurtosis_excess(as), -1.5); // 3/2 + BOOST_CHECK_EQUAL(kurtosis(as), 1.5); // 3/2 + BOOST_CHECK_EQUAL(kurtosis_excess(as), -1.5); // -3/2 BOOST_CHECK_EQUAL(support(as).first, 0); // BOOST_CHECK_EQUAL(range(as).first, 0); // BOOST_CHECK_THROW(mode(as), std::domain_error); // Two modes at x_min and x_max, so throw instead. diff --git a/test/test_bernoulli.cpp b/test/test_bernoulli.cpp index 5f78b0afc5..256a897783 100644 --- a/test/test_bernoulli.cpp +++ b/test/test_bernoulli.cpp @@ -114,6 +114,11 @@ void test_spots(RealType) bernoulli_distribution(static_cast(0.4L)), static_cast(0)), static_cast(0.6L), tolerance); // Expect 1- p. + BOOST_CHECK_CLOSE_FRACTION( + pdf( // OK k (or n) + bernoulli_distribution(static_cast(0.4L)), static_cast(1)), + static_cast(0.4L), tolerance); // Expect p. + BOOST_CHECK_EQUAL( mean(bernoulli_distribution(static_cast(0.5L))), static_cast(0.5L)); @@ -280,6 +285,46 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(quantile(complement(w, +inf)), std::domain_error); // p = + inf } // has_infinity #endif + + using boost::math::policies::policy; + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef bernoulli_distribution ignore_error_bernoulli; + + if (std::numeric_limits::has_quiet_NaN) + { + // PDF + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_bernoulli(0.5), std::numeric_limits::infinity()))); // k == infinity + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_bernoulli(0.5), 2))); // k > 1 + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_bernoulli(0.5), static_cast(0.5)))); // k != 0, 1 + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_bernoulli(0.5), static_cast(-0.5)))); // k != 0, 1 + + // CDF + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_bernoulli(0.5), std::numeric_limits::infinity()))); // k == infinity + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_bernoulli(0.5), 2))); // k > 1 + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_bernoulli(0.5), static_cast(0.5)))); // k != 0, 1 + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_bernoulli(0.5), static_cast(-0.5)))); // k != 0, 1 + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_bernoulli(0.5), std::numeric_limits::infinity())))); // k == infinity + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_bernoulli(0.5), 2)))); // k > 1 + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_bernoulli(0.5), static_cast(0.5))))); // k != 0, 1 + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_bernoulli(0.5), static_cast(-0.5))))); // k != 0, 1 + + // // Quantile + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_bernoulli(0.5), std::numeric_limits::infinity()))); // p == infinity + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_bernoulli(0.5), static_cast(-1)))); // p < 0 + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_bernoulli(0.5), static_cast(2)))); // p > 1 + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_bernoulli(0.5), std::numeric_limits::infinity())))); // q == infinity + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_bernoulli(0.5), static_cast(-1))))); // q < 0 + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_bernoulli(0.5), static_cast(2))))); // q > 1 + + } // has_quiet_NaN } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main ) diff --git a/test/test_beta_dist.cpp b/test/test_beta_dist.cpp index d989197b6b..3983f526ff 100644 --- a/test/test_beta_dist.cpp +++ b/test/test_beta_dist.cpp @@ -65,6 +65,18 @@ using std::numeric_limits; # include #endif +template +void test_find_alpha_beta(RealType mean, RealType variance){ + BOOST_MATH_CHECK_THROW(beta_distribution::find_alpha(mean, variance), std::domain_error); + BOOST_MATH_CHECK_THROW(beta_distribution::find_beta(mean, variance), std::domain_error); +} + +template +void test_find_alpha_beta(RealType param, RealType x, RealType p){ + BOOST_MATH_CHECK_THROW(beta_distribution::find_alpha(param, x, p), std::domain_error); + BOOST_MATH_CHECK_THROW(beta_distribution::find_beta(param, x, p), std::domain_error); +} + template void test_spot( RealType a, // alpha a @@ -175,7 +187,9 @@ void test_spots(RealType) using ::boost::math::pdf; // Tests that should throw: - BOOST_MATH_CHECK_THROW(mode(beta_distribution(static_cast(1), static_cast(1))), std::domain_error); + BOOST_MATH_CHECK_THROW(mode(beta_distribution(static_cast(1), static_cast(1))), std::domain_error); // alpha = beta = 1 + BOOST_MATH_CHECK_THROW(mode(beta_distribution(static_cast(2), static_cast(1))), std::domain_error); // beta = 1 + BOOST_MATH_CHECK_THROW(mode(beta_distribution(static_cast(1), static_cast(2))), std::domain_error); // alpha = 1 // mode is undefined, and throws domain_error! // BOOST_MATH_CHECK_THROW(median(beta_distribution(static_cast(1), static_cast(1))), std::domain_error); @@ -497,7 +511,19 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(quantile(dist, -1), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(complement(dist, -1)), std::domain_error); - // No longer allow any parameter to be NaN or inf, so all these tests should throw. + BOOST_MATH_CHECK_THROW(pdf(beta_distribution(static_cast(0.5), static_cast(1.5)), static_cast(0)), std::overflow_error); // alpha < 1, x = 0 + BOOST_MATH_CHECK_THROW(pdf(beta_distribution(static_cast(1.5), static_cast(0.5)), static_cast(1)), std::overflow_error); // beta < 1, x = 1 + + // There aren't any upper bounds on the mean/variance. Should we add these? + test_find_alpha_beta(static_cast(-1), static_cast(1)); // mean < 0 + test_find_alpha_beta(static_cast(-1), static_cast(-1)); // var < 0 + + test_find_alpha_beta(static_cast(1), static_cast(0.5), static_cast(-0.5)); // p < 0 + test_find_alpha_beta(static_cast(1), static_cast(0.5), static_cast(1.5)); // p > 0 + test_find_alpha_beta(static_cast(1), static_cast(1.5), static_cast(0.5)); // x > 0 + test_find_alpha_beta(static_cast(1), static_cast(-0.5), static_cast(0.5)); // x < 0 + + // No longer allow any parameter to be NaN or inf, so all these tests should throw. if (std::numeric_limits::has_quiet_NaN) { // Attempt to construct from non-finite should throw. @@ -545,6 +571,37 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(cdf(complement(w, +inf)), std::domain_error); // x = + inf BOOST_MATH_CHECK_THROW(quantile(w, +inf), std::domain_error); // p = + inf BOOST_MATH_CHECK_THROW(quantile(complement(w, +inf)), std::domain_error); // p = + inf + + test_find_alpha_beta(static_cast(1), static_cast(0.5), inf); // p = inf + test_find_alpha_beta(inf, static_cast(0.5)); // mean = inf + test_find_alpha_beta(static_cast(0.5), inf); // var = inf + test_find_alpha_beta(static_cast(1), inf, static_cast(0.5)); // x = inf + + // Check isnan policies + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef beta_distribution ignore_error_beta; + + BOOST_CHECK((boost::math::isnan)(ignore_error_beta::find_alpha(inf, static_cast(0.5)))); // mean = inf + BOOST_CHECK((boost::math::isnan)(ignore_error_beta::find_beta(inf, static_cast(0.5)))); // mean = inf + BOOST_CHECK((boost::math::isnan)(ignore_error_beta::find_alpha(inf, static_cast(0.5), static_cast(0.5)))); // beta = inf + BOOST_CHECK((boost::math::isnan)(ignore_error_beta::find_beta(inf, static_cast(0.5), static_cast(0.5)))); // alpha = inf + BOOST_CHECK((boost::math::isnan)(mode(ignore_error_beta(static_cast(2), static_cast(1))))); // beta = 1 + BOOST_CHECK((boost::math::isnan)(mode(ignore_error_beta(static_cast(1), static_cast(2))))); // alpha = 1 + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_beta(static_cast(1), static_cast(1)), 2))); // x > 1 + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_beta(static_cast(1), static_cast(1)), 2))); // x > 1 + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_beta(static_cast(1), static_cast(1)), 2)))); // x > 1 + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_beta(static_cast(1), static_cast(1)), 2))); // p > 1 + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_beta(static_cast(1), static_cast(1)), 2)))); // p > 1 } // has_infinity // Error handling checks: @@ -600,6 +657,8 @@ BOOST_AUTO_TEST_CASE( test_main ) BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.5), 0.5, 2 * std::numeric_limits::epsilon()); BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.9), 0.9, 2 * std::numeric_limits::epsilon()); BOOST_CHECK_EQUAL(cdf(mybeta11, 1), 1.); // Exact unity expected. + BOOST_CHECK_EQUAL(cdf(complement(mybeta11, 1)), 0.0); // Exact 0 expected. + BOOST_CHECK_EQUAL(cdf(complement(mybeta11, 0)), 1.0); // Exact unity expected. double tol = std::numeric_limits::epsilon() * 10; BOOST_CHECK_EQUAL(pdf(mybeta22, 1), 0); // is dome shape. @@ -627,6 +686,11 @@ BOOST_AUTO_TEST_CASE( test_main ) // quantile. BOOST_CHECK_CLOSE_FRACTION(quantile(mybeta22, 0.028), 0.1, tol); BOOST_CHECK_CLOSE_FRACTION(quantile(complement(mybeta22, 1 - 0.028)), 0.1, tol); + BOOST_CHECK_EQUAL(quantile(mybeta22, 0), 0.); + BOOST_CHECK_EQUAL(quantile(mybeta22, 1), 1.); + BOOST_CHECK_EQUAL(quantile(complement(mybeta22, 0)), 1.); + BOOST_CHECK_EQUAL(quantile(complement(mybeta22, 1)), 0.); + BOOST_CHECK_EQUAL(kurtosis(mybeta11), 3+ kurtosis_excess(mybeta11)); // Check kurtosis_excess = kurtosis - 3; BOOST_CHECK_CLOSE_FRACTION(variance(mybeta22), 0.05, tol); BOOST_CHECK_CLOSE_FRACTION(mean(mybeta22), 0.5, tol); diff --git a/test/test_binomial.cpp b/test/test_binomial.cpp index ef7f171723..b76c8e03e4 100644 --- a/test/test_binomial.cpp +++ b/test/test_binomial.cpp @@ -599,7 +599,111 @@ void test_spots(RealType T) binomial_distribution(static_cast(8), static_cast(1)), static_cast(7)), static_cast(0) ); + BOOST_CHECK_EQUAL( + cdf( + complement(binomial_distribution(static_cast(8), static_cast(0)), + static_cast(7))), static_cast(0) + ); + BOOST_CHECK_EQUAL( + cdf( + complement(binomial_distribution(static_cast(8), static_cast(1)), + static_cast(7))), static_cast(1) + ); + // Check Error handling and edge cases + BOOST_CHECK_EQUAL(binomial_distribution::find_lower_bound_on_p(static_cast(5), + static_cast(0), + static_cast(0.5)), + static_cast(0)); // success = 0 + BOOST_CHECK_THROW(binomial_distribution::find_lower_bound_on_p(static_cast(5), + static_cast(2), + static_cast(-0.5)), + std::domain_error); // probability < 0 + BOOST_CHECK_THROW(binomial_distribution::find_lower_bound_on_p(static_cast(5), + static_cast(2), + static_cast(1.5)), + std::domain_error); // probability > 1 + BOOST_CHECK_EQUAL(binomial_distribution::find_upper_bound_on_p(static_cast(5), + static_cast(5), + static_cast(0.5)), + static_cast(1)); // trials = successes + BOOST_CHECK_THROW(binomial_distribution::find_upper_bound_on_p(static_cast(5), + static_cast(2), + static_cast(-0.5)), + std::domain_error); // probability < 0 + BOOST_CHECK_THROW(binomial_distribution::find_upper_bound_on_p(static_cast(5), + static_cast(2), + static_cast(1.5)), + std::domain_error); // probability > 1 + BOOST_CHECK_THROW(binomial_distribution::find_minimum_number_of_trials(static_cast(10), + static_cast(1.5), + static_cast(0.05)), + std::domain_error); // probability > 1 + BOOST_CHECK_THROW(binomial_distribution::find_minimum_number_of_trials(static_cast(10), + static_cast(-1), + static_cast(0.05)), + std::domain_error); // probability < 0 + BOOST_CHECK_THROW(binomial_distribution::find_minimum_number_of_trials(static_cast(10), + static_cast(0.5), + static_cast(-0.05)), + std::domain_error); // alpha < 0 + BOOST_CHECK_THROW(binomial_distribution::find_minimum_number_of_trials(static_cast(10), + static_cast(0.5), + static_cast(1.314)), + std::domain_error); // alpha > 1 + BOOST_CHECK_THROW(binomial_distribution::find_maximum_number_of_trials(static_cast(10), + static_cast(1.5), + static_cast(0.05)), + std::domain_error); // probability > 1 + BOOST_CHECK_THROW(binomial_distribution::find_maximum_number_of_trials(static_cast(10), + static_cast(-1), + static_cast(0.05)), + std::domain_error); // probability < 0 + BOOST_CHECK_THROW(binomial_distribution::find_maximum_number_of_trials(static_cast(10), + static_cast(0.5), + static_cast(-0.05)), + std::domain_error); // alpha < 0 + BOOST_CHECK_THROW(binomial_distribution::find_maximum_number_of_trials(static_cast(10), + static_cast(0.5), + static_cast(1.5)), + std::domain_error); // alpha > 1 + + if (std::numeric_limits::has_infinity) + { + // Check isnan policies + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef binomial_distribution ignore_error_binomial; + + BOOST_CHECK((boost::math::isnan)(ignore_error_binomial::find_lower_bound_on_p(static_cast(5), + static_cast(2), + static_cast(-0.5)))); // probability < 0 + BOOST_CHECK((boost::math::isnan)(ignore_error_binomial::find_upper_bound_on_p(static_cast(5), + static_cast(2), + static_cast(-0.5)))); // probability < 0 + BOOST_CHECK((boost::math::isnan)(ignore_error_binomial::find_minimum_number_of_trials(static_cast(10), + static_cast(-1), + static_cast(0.05)))); // probability < 0 + BOOST_CHECK((boost::math::isnan)(ignore_error_binomial::find_maximum_number_of_trials(static_cast(10), + static_cast(-1), + static_cast(0.05)))); // probability < 0 + + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_binomial(static_cast(-8), static_cast(0.25)), static_cast(8)))); + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_binomial(static_cast(8), static_cast(-0.25)), static_cast(8)))); + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_binomial(static_cast(8), static_cast(0.25)), static_cast(9)))); + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_binomial(static_cast(8), static_cast(0.25)), static_cast(9)))); + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_binomial(static_cast(8), static_cast(0.25)), static_cast(9))))); + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_binomial(static_cast(8), static_cast(0.25)), static_cast(-1)))); + } #endif { diff --git a/test/test_cauchy.cpp b/test/test_cauchy.cpp index 8136999e37..4e8286c60e 100644 --- a/test/test_cauchy.cpp +++ b/test/test_cauchy.cpp @@ -744,7 +744,61 @@ void test_spots(RealType T) check_out_of_range >(0, 1); // (All) valid constructor parameter values. - + if (std::numeric_limits::has_infinity) + { + BOOST_CHECK_EQUAL( + boost::math::pdf( + cauchy_distribution(-2, 0.25), + std::numeric_limits::infinity()), // x + static_cast(0) // probability + ); + BOOST_CHECK_EQUAL( + boost::math::cdf( + cauchy_distribution(-2, 0.25), + std::numeric_limits::infinity()), // x + static_cast(1) // probability + ); + BOOST_CHECK_EQUAL( + boost::math::cdf( + cauchy_distribution(-2, 0.25), + -std::numeric_limits::infinity()), // x + static_cast(0) // probability + ); + BOOST_CHECK_EQUAL( + boost::math::quantile( + cauchy_distribution(-2, 0.25), + static_cast(0.5)), // p + static_cast(-2) // location + ); + + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef boost::math::cauchy_distribution ignore_error_cauchy; + + // PDF + BOOST_CHECK((boost::math::isnan)(boost::math::pdf(ignore_error_cauchy(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(boost::math::pdf(ignore_error_cauchy(std::numeric_limits::infinity(), static_cast(1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(boost::math::pdf(ignore_error_cauchy(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + // CDF + BOOST_CHECK((boost::math::isnan)(boost::math::cdf(ignore_error_cauchy(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(boost::math::cdf(ignore_error_cauchy(std::numeric_limits::infinity(), static_cast(1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(boost::math::cdf(ignore_error_cauchy(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + // Quantile + BOOST_CHECK((boost::math::isnan)(boost::math::quantile(ignore_error_cauchy(static_cast(0), static_cast(-1)), static_cast(0.25)))); + BOOST_CHECK((boost::math::isnan)(boost::math::quantile(ignore_error_cauchy(std::numeric_limits::infinity(), static_cast(1)), static_cast(0.25)))); + BOOST_CHECK((boost::math::isnan)(boost::math::quantile(ignore_error_cauchy(static_cast(0), static_cast(1)), static_cast(-0.25)))); + } } // template void test_spots(RealType) diff --git a/test/test_chi_squared.cpp b/test/test_chi_squared.cpp index a64926ca58..360f6947de 100644 --- a/test/test_chi_squared.cpp +++ b/test/test_chi_squared.cpp @@ -546,6 +546,9 @@ void test_spots(RealType T) chi_squared_distribution(static_cast(8)), static_cast(1.1))), std::domain_error ); + BOOST_MATH_CHECK_THROW( + mode(chi_squared_distribution(static_cast(1))), std::domain_error + ); // This first test value is taken from an example here: // http://www.itl.nist.gov/div898/handbook/prc/section2/prc232.htm @@ -566,6 +569,32 @@ void test_spots(RealType T) check_out_of_range >(1); // (All) valid constructor parameter values. + // NaN handling + if (std::numeric_limits::has_infinity) + { + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef boost::math::chi_squared_distribution ignore_error_chi_squared; + + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_chi_squared(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_chi_squared(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_chi_squared(static_cast(-1)), static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_chi_squared(static_cast(-1)), static_cast(0.5)))); + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_chi_squared(static_cast(-1)), static_cast(0.5))))); + BOOST_CHECK((boost::math::isnan)(mode(ignore_error_chi_squared(static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(ignore_error_chi_squared::find_degrees_of_freedom(10, -1, 0.01f, 100))); + BOOST_CHECK((boost::math::isnan)(ignore_error_chi_squared::find_degrees_of_freedom(10, 0.05f, -1, 100))); + + } } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main ) diff --git a/test/test_ellint_2.hpp b/test/test_ellint_2.hpp index 29a73c9961..3b565a493d 100644 --- a/test/test_ellint_2.hpp +++ b/test/test_ellint_2.hpp @@ -97,7 +97,7 @@ void test_spots(T, const char* type_name) // called a second time. // #ifndef TEST_UDT - static const std::array, 28> data1 = {{ + static const std::array, 31> data1 = {{ {{ SC_(0.0), SC_(0.0), SC_(0.0) }}, {{ SC_(-10.0), SC_(0.0), SC_(-10.0) }}, {{ SC_(-1.0), SC_(-1.0), SC_(-0.84147098480789650665250232163029899962256306079837) }}, @@ -127,6 +127,10 @@ void test_spots(T, const char* type_name) {{ -2 * boost::math::constants::pi(), SC_(1.0), SC_(-4.0) }}, {{ -(20 * boost::math::constants::pi()) / 21, SC_(1.0), SC_(-1.85095773382382555307064528472278244309033056100177750424044959888147151556780639619867349409562910808166870200808) }}, {{ -(20 * boost::math::constants::pi()) / 19, SC_(1.0), SC_(-2.16459459028073389414365205908793841951217248335965412335147127404907690285247170629778731438009399303864259295133) }}, + // Test at n * pi / 2, see https://github.com/boostorg/math/issues/1377 + {{ boost::math::constants::half_pi(), SC_(0.25), boost::math::ellint_2(static_cast(0.25)) }}, + {{ 3 * boost::math::constants::half_pi(), SC_(0.25), 3 * boost::math::ellint_2(static_cast(0.25)) }}, + {{ 6 * boost::math::constants::half_pi(), SC_(0.25), 6 * boost::math::ellint_2(static_cast(0.25)) }}, }}; do_test_ellint_e2(data1, type_name, "Elliptic Integral E: Mathworld Data"); diff --git a/test/test_exponential_dist.cpp b/test/test_exponential_dist.cpp index 124e2929ad..ecf5106d01 100644 --- a/test/test_exponential_dist.cpp +++ b/test/test_exponential_dist.cpp @@ -362,6 +362,48 @@ void test_spots(RealType T) BOOST_CHECK_EQUAL(pdf(exponential_distribution(2), inf), 0); BOOST_CHECK_EQUAL(cdf(exponential_distribution(2), inf), 1); BOOST_CHECK_EQUAL(cdf(complement(exponential_distribution(2), inf)), 0); + + // Test NaN policies + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef boost::math::exponential_distribution ignore_error_exponential; + + // PDF + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_exponential(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_exponential(static_cast(1)), static_cast(-1)))); + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_exponential(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_exponential(static_cast(1)), static_cast(-1)))); + + // CDF + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_exponential(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_exponential(static_cast(1)), static_cast(-1)))); + BOOST_CHECK((boost::math::isnan)(logcdf(ignore_error_exponential(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(logcdf(ignore_error_exponential(static_cast(1)), static_cast(-1)))); + + // Complement CDF + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_exponential(static_cast(-1)), static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_exponential(static_cast(1)), static_cast(-1))))); + BOOST_CHECK((boost::math::isnan)(logcdf(complement(ignore_error_exponential(static_cast(-1)), static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(logcdf(complement(ignore_error_exponential(static_cast(1)), static_cast(-1))))); + + // Quantile + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_exponential(static_cast(-1)), static_cast(0.5)))); + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_exponential(static_cast(1)), static_cast(-0.5)))); + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_exponential(static_cast(-1)), static_cast(0.5))))); + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_exponential(static_cast(1)), static_cast(-0.5))))); + + // Mean and Standard Deviation + BOOST_CHECK((boost::math::isnan)(mean(ignore_error_exponential(static_cast(-1))))); + BOOST_CHECK((boost::math::isnan)(standard_deviation(ignore_error_exponential(static_cast(-1))))); } } // template void test_spots(RealType) diff --git a/test/test_extreme_value.cpp b/test/test_extreme_value.cpp index 255777914f..dd0c32ae47 100644 --- a/test/test_extreme_value.cpp +++ b/test/test_extreme_value.cpp @@ -253,6 +253,65 @@ void test_spots(RealType) BOOST_CHECK_EQUAL(cdf(complement(extreme_value_distribution(), inf)), 0); BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution(), -inf), 0); BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution(), inf), 1); + + // Test NaN policy + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef boost::math::extreme_value_distribution ignore_error_extreme; + + // PDF + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_extreme(inf, static_cast(1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + // log(PDF) + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_extreme(inf, static_cast(1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + // CDF + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_extreme(inf, static_cast(1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0))))); + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_extreme(inf, static_cast(1)), static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN())))); + + // log(CDF) + BOOST_CHECK((boost::math::isnan)(logcdf(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(logcdf(ignore_error_extreme(inf, static_cast(1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(logcdf(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + BOOST_CHECK((boost::math::isnan)(logcdf(complement(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0))))); + BOOST_CHECK((boost::math::isnan)(logcdf(complement(ignore_error_extreme(inf, static_cast(1)), static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(logcdf(complement(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN())))); + + // Quantile + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0.5)))); + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_extreme(inf, static_cast(1)), static_cast(0.5)))); + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_extreme(static_cast(0), static_cast(1)), static_cast(-1)))); + + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0.5))))); + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_extreme(inf, static_cast(1)), static_cast(0.5))))); + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_extreme(static_cast(0), static_cast(1)), static_cast(-1))))); + + // Mean + BOOST_CHECK((boost::math::isnan)(mean(ignore_error_extreme(static_cast(0), static_cast(-1))))); + BOOST_CHECK((boost::math::isnan)(mean(ignore_error_extreme(inf, static_cast(1))))); + + // Standard deviation + BOOST_CHECK((boost::math::isnan)(standard_deviation(ignore_error_extreme(static_cast(0), static_cast(-1))))); + BOOST_CHECK((boost::math::isnan)(standard_deviation(ignore_error_extreme(inf, static_cast(1))))); } // // Bug reports: diff --git a/test/test_find_location.cpp b/test/test_find_location.cpp index 7e0fb543a5..4d376dc254 100644 --- a/test/test_find_location.cpp +++ b/test/test_find_location.cpp @@ -74,18 +74,18 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(find_location(0., 2., 0.), std::domain_error); // p above 0 to 1. BOOST_MATH_CHECK_THROW(find_location(numeric_limits::infinity(), 0.5, 0.), std::domain_error); // z not finite. - BOOST_MATH_CHECK_THROW(find_location(numeric_limits::quiet_NaN(), -1., 0.), + BOOST_MATH_CHECK_THROW(find_location(numeric_limits::quiet_NaN(), 0.5, 0.), std::domain_error); // z not finite - BOOST_MATH_CHECK_THROW(find_location(0., -1., numeric_limits::quiet_NaN()), + BOOST_MATH_CHECK_THROW(find_location(0., 0.5, numeric_limits::quiet_NaN()), std::domain_error); // scale not finite BOOST_MATH_CHECK_THROW(find_location(complement(0., -1., 0.)), std::domain_error); // p below 0 to 1. BOOST_MATH_CHECK_THROW(find_location(complement(0., 2., 0.)), std::domain_error); // p above 0 to 1. BOOST_MATH_CHECK_THROW(find_location(complement(numeric_limits::infinity(), 0.5, 0.)), std::domain_error); // z not finite. - BOOST_MATH_CHECK_THROW(find_location(complement(numeric_limits::quiet_NaN(), -1., 0.)), + BOOST_MATH_CHECK_THROW(find_location(complement(numeric_limits::quiet_NaN(), 0.5, 0.)), std::domain_error); // z not finite - BOOST_MATH_CHECK_THROW(find_location(complement(0., -1., numeric_limits::quiet_NaN())), + BOOST_MATH_CHECK_THROW(find_location(complement(0., 0.5, numeric_limits::quiet_NaN())), std::domain_error); // scale not finite //// Check for ab-use with unsuitable distribution(s) when concept check implemented. @@ -131,12 +131,40 @@ void test_spots(RealType) // Check that can use the complement version. RealType q = 1 - p; // complement. // cout << "find_location >(complement(z, q, sd)) = " << endl; + l = find_location >(complement(z, q, sd, policy<>())); l = find_location >(complement(z, q, sd)); normal_distribution np95pc(l, sd); // Same standard_deviation (scale) but with mean(location) shifted // cout << "Normal distribution with mean = " << l << " has " << "fraction <= " << z << " = " << cdf(np95pc, z) << endl; BOOST_CHECK_CLOSE_FRACTION(q, cdf(np95pc, z), tolerance); + // Test NaN handling + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + if (std::numeric_limits::has_infinity) + { + RealType inf = std::numeric_limits::infinity(); + BOOST_CHECK((boost::math::isnan)(find_location >(static_cast(0), static_cast(-1), static_cast(1), ignore_all_policy()))); + BOOST_CHECK((boost::math::isnan)(find_location >(inf, static_cast(0.5), static_cast(1), ignore_all_policy()))); + BOOST_CHECK((boost::math::isnan)(find_location >(static_cast(0), static_cast(0.5), inf, ignore_all_policy()))); + // Negative scale does not throw. Seems weird... waiting to see if this should be fixed + BOOST_CHECK((boost::math::isnan)(find_location >(static_cast(0), static_cast(0.5), static_cast(-1), ignore_all_policy()))); + + // Complement + BOOST_CHECK((boost::math::isnan)(find_location >(complement(static_cast(0), static_cast(-1), static_cast(1), ignore_all_policy())))); + BOOST_CHECK((boost::math::isnan)(find_location >(complement(inf, static_cast(0.5), static_cast(1), ignore_all_policy())))); + BOOST_CHECK((boost::math::isnan)(find_location >(complement(static_cast(0), static_cast(0.5), inf, ignore_all_policy())))); + } + + + } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main ) diff --git a/test/test_find_scale.cpp b/test/test_find_scale.cpp index f31e47969d..a3c62f39e8 100644 --- a/test/test_find_scale.cpp +++ b/test/test_find_scale.cpp @@ -84,9 +84,9 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(find_scale(complement(0., 2., 0.)), std::domain_error); // p above 0 to 1. BOOST_MATH_CHECK_THROW(find_scale(complement(numeric_limits::infinity(), 0.5, 0.)), std::domain_error); // z not finite. - BOOST_MATH_CHECK_THROW(find_scale(complement(numeric_limits::quiet_NaN(), -1., 0.)), + BOOST_MATH_CHECK_THROW(find_scale(complement(numeric_limits::quiet_NaN(), 0.5, 0.)), std::domain_error); // z not finite - BOOST_MATH_CHECK_THROW(find_scale(complement(0., -1., numeric_limits::quiet_NaN())), + BOOST_MATH_CHECK_THROW(find_scale(complement(0., 0.5, numeric_limits::quiet_NaN())), std::domain_error); // scale not finite BOOST_MATH_CHECK_THROW(find_scale(complement(0., -1., 0.)), std::domain_error); // p below 0 to 1. @@ -110,12 +110,36 @@ void test_spots(RealType) #ifndef BOOST_NO_EXCEPTIONS BOOST_CHECK_NO_THROW(find_scale(0, -1, 1, ignore_domain_policy())); // probability outside [0, 1] - BOOST_CHECK_NO_THROW(find_scale(numeric_limits::infinity(), -1, 1, + BOOST_CHECK_NO_THROW(find_scale(numeric_limits::infinity(), 0.5, 1, ignore_domain_policy())); // z not finite. BOOST_CHECK_NO_THROW(find_scale(complement(0, -1, 1, ignore_domain_policy()))); // probability outside [0, 1] - BOOST_CHECK_NO_THROW(find_scale(complement(numeric_limits::infinity(), -1, 1, + BOOST_CHECK_NO_THROW(find_scale(complement(numeric_limits::infinity(), 0.5, 1, ignore_domain_policy()))); // z not finite. #endif + + // Test NaN handling + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + if (std::numeric_limits::has_infinity) + { + RealType inf = std::numeric_limits::infinity(); + BOOST_CHECK((boost::math::isnan)(find_scale >(static_cast(0), static_cast(-1), static_cast(1), ignore_all_policy()))); + BOOST_CHECK((boost::math::isnan)(find_scale >(inf, static_cast(0.5), static_cast(1), ignore_all_policy()))); + BOOST_CHECK((boost::math::isnan)(find_scale >(static_cast(0), static_cast(0.5), inf, ignore_all_policy()))); + + // Complement + BOOST_CHECK((boost::math::isnan)(find_scale >(complement(static_cast(0), static_cast(-1), static_cast(1), ignore_all_policy())))); + BOOST_CHECK((boost::math::isnan)(find_scale >(complement(inf, static_cast(0.5), static_cast(1), ignore_all_policy())))); + BOOST_CHECK((boost::math::isnan)(find_scale >(complement(static_cast(0), static_cast(0.5), inf, ignore_all_policy())))); + } + RealType l = 0.; // standard normal distribution. RealType sd = static_cast(1); // normal default standard deviation = 1. normal_distribution n01(l, sd); // mean(location) = 0, standard_deviation (scale) = 1. diff --git a/test/test_ibeta.cpp b/test/test_ibeta.cpp index 563844a06c..0059b46619 100644 --- a/test/test_ibeta.cpp +++ b/test/test_ibeta.cpp @@ -380,18 +380,18 @@ BOOST_AUTO_TEST_CASE( test_main ) gsl_set_error_handler_off(); #endif #ifdef TEST_FLOAT - test_spots(0.0F); + test_spots(0.0F, "float"); #endif #ifdef TEST_DOUBLE - test_spots(0.0); + test_spots(0.0, "double"); #endif #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS #ifdef TEST_LDOUBLE - test_spots(0.0L); + test_spots(0.0L, "long double"); #endif #if !BOOST_WORKAROUND(BOOST_BORLANDC, BOOST_TESTED_AT(0x582)) #if defined(TEST_REAL_CONCEPT) && !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS) - test_spots(boost::math::concepts::real_concept(0.1)); + test_spots(boost::math::concepts::real_concept(0.1), "real_concept"); #endif #endif #endif diff --git a/test/test_ibeta.hpp b/test/test_ibeta.hpp index 2f0ea47107..13ff157f9d 100644 --- a/test/test_ibeta.hpp +++ b/test/test_ibeta.hpp @@ -152,12 +152,13 @@ void test_beta(T, const char* name) } template -void test_spots(T) +void test_spots(T, const char* name) { // // basic sanity checks, tolerance is 30 epsilon expressed as a percentage: // Spot values are from http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BetaRegularized // using precision of 50 decimal digits. + std::cout << "Testing spot values with type " << name << std::endl; T tolerance = boost::math::tools::epsilon() * 3000; if (boost::math::tools::digits() > 100) tolerance *= 2; @@ -487,5 +488,35 @@ void test_spots(T) } BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast(2), static_cast(1), static_cast(0)), 0); BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast(1), static_cast(2), static_cast(0)), 0); + + // Bug testing for large a,b and x close to a / (a+b). See PR 1363. + // The values for a,b are just too large for floats to handle. The tests here only + // check if ibeta is monotonically increasing for a,b fixed and x increasing. Both + // Mathematica and mpmath (in Python) are not able to evaluate ibeta for such large + // values of a,b so spot testing wasn't possible. The accuracy of the fix for + // these large values is very sensitive to the computer architecture. + // Macos with arm64 does the worst but linux and windows pass a larger range of tests. + if (!std::is_same::value) + { + // Larger values of a,b become more numerically unstable. Larger/smaller values + // of delta also become unstable. + T a_values[1] = {10000000272564224}; + T b_values[1] = {9965820922822656}; + T delta[7] = {0, 1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10}; + T a, b, x; + for (unsigned int i=0; i<1; i++){ + a = a_values[i]; + b = b_values[i]; + x = a / (a+b); // roughly the median of ibeta + + for (unsigned int j=0; j < 6; j++) + { + BOOST_CHECK_MESSAGE(boost::math::ibeta(a, b, x + delta[j+1]) > boost::math::ibeta(a, b, x + delta[j]), + "ibeta not monotonically increasing above a/(a+b) for ibeta(" << a << ", " << b << ", " << x << ") and delta=" << delta[j] << ": [" << boost::math::ibeta(a, b, x + delta[j+1]) << " >= " << boost::math::ibeta(a, b, x + delta[j]) << "]"); + BOOST_CHECK_MESSAGE(boost::math::ibeta(a, b, x - delta[j]) > boost::math::ibeta(a, b, x - delta[j+1]), + "ibeta not monotonically increasing below a/(a+b) for ibeta(" << a << ", " << b << ", " << x << ") and delta=" << delta[j] << ": [" << boost::math::ibeta(a, b, x - delta[j]) << " >= " << boost::math::ibeta(a, b, x - delta[j+1]) << "]"); + } + } + } } diff --git a/test/test_inverse_gaussian.cpp b/test/test_inverse_gaussian.cpp index 3825da9397..0fe64223dc 100644 --- a/test/test_inverse_gaussian.cpp +++ b/test/test_inverse_gaussian.cpp @@ -252,6 +252,13 @@ BOOST_AUTO_TEST_CASE( test_main ) // but better than R that gives up completely! // R Error in SuppDists::qinverse_gaussian(4.87914430108515e-219, 1, 1) : Infinite value in NewtonRoot() + inverse_gaussian w_big(66.99652081); + BOOST_CHECK_CLOSE_FRACTION( + quantile(w_big, 0.97969), 591.567880739988823, 15 * tolfeweps); + BOOST_CHECK_CLOSE_FRACTION( + quantile(complement(w_big, 1 - 0.97969)), 591.567880739988823, 10 * tolfeweps); + + BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 0.5), static_cast(0.87878257893544476), tolfeweps); // pdf BOOST_CHECK_CLOSE_FRACTION( diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 6b7c2347ab..7ad45c1ba3 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -138,13 +138,13 @@ void test_spot( BOOST_CHECK_CLOSE( cdf(complement(dist, x)), Q, tol); BOOST_CHECK_CLOSE( - quantile(dist, P), x, tol * 10); + quantile(dist, P), x, tol * 10); BOOST_CHECK_CLOSE( - quantile(complement(dist, Q)), x, tol * 10); + quantile(complement(dist, Q)), x, tol * 10); BOOST_CHECK_CLOSE( - dist.find_non_centrality(x, a, b, P), ncp, tol * 10); + dist.find_non_centrality(x, a, b, P), ncp, tol * 10); BOOST_CHECK_CLOSE( - dist.find_non_centrality(boost::math::complement(x, a, b, Q)), ncp, tol * 10); + dist.find_non_centrality(boost::math::complement(x, a, b, Q)), ncp, tol * 10); } if(boost::math::tools::digits() > 50) { @@ -369,6 +369,48 @@ void test_spots(RealType, const char* name = nullptr) } } } + + // Quick spot check for finding degrees of freedom. When checking for two degrees + // of freedom for real_concept types, the cdf at large/small v2 can be greater than 1 + // or less than 0. + if (!std::is_same::value){ + RealType v1 = 10; + RealType v2 = 5; + nc = 1; + x = 6; + boost::math::non_central_f_distribution ref(v1, v2, nc); + RealType P = cdf(ref, x); + BOOST_CHECK_CLOSE(ref.find_v2(x, v1, nc, P), v2, tolerance); + BOOST_CHECK_CLOSE(ref.find_v2(boost::math::complement(x, v1, nc, 1-P)), v2, tolerance); + BOOST_CHECK_CLOSE(ref.find_v1(x, v2, nc, P), v1, tolerance); + BOOST_CHECK_CLOSE(ref.find_v1(boost::math::complement(x, v2, nc, 1-P)), v1, tolerance); + } + + // Check case where two degrees of freedom solve the inversion problem + BOOST_MATH_CHECK_THROW(dist.find_v1(RealType(1.5), RealType(2.0), RealType(1.0), RealType(0.49845842011686358665786775091245664L)), boost::math::evaluation_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(RealType(3.51), RealType(5), RealType(0), RealType(0.85802971653663762108266155337333L)), boost::math::evaluation_error); + + // Check find_v1/v2 edge cases + // Case when P=1 or P=0 + nc = 2; + BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 1), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 0), std::domain_error); + // Case when Q=1 or Q=0 + BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 1)), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 0)), std::domain_error); + // Check very small values of x an evaluation error is thrown + x = boost::math::tools::epsilon() / 10; + BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 0.5)), boost::math::evaluation_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 0.5), boost::math::evaluation_error); + + BOOST_MATH_CHECK_THROW(dist.find_v2(x, b, nc, 1), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v2(x, b, nc, 0), std::domain_error); + // Case when Q=1 or Q=0 + BOOST_MATH_CHECK_THROW(dist.find_v2(boost::math::complement(x, b, nc, 1)), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v2(boost::math::complement(x, b, nc, 0)), std::domain_error); + // Check very small values of x an evaluation error is thrown + BOOST_MATH_CHECK_THROW(dist.find_v2(boost::math::complement(x, b, nc, 0.5)), boost::math::evaluation_error); + BOOST_MATH_CHECK_THROW(dist.find_v2(x, b, nc, 0.5), boost::math::evaluation_error); } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main ) diff --git a/test/test_pFq.hpp b/test/test_pFq.hpp index 0e7f3a0c54..df054caad3 100644 --- a/test/test_pFq.hpp +++ b/test/test_pFq.hpp @@ -298,6 +298,31 @@ void test_spots_2F2(T, const char*) } } +template +void test_special_cases(T, const char*) +{ + typedef boost::math::policies::policy, boost::math::policies::promote_double > nothrow_policy; + typedef boost::math::policies::policy, boost::math::policies::promote_double > throw_policy; + T error_rate; + BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ 2 }, { 3 }, static_cast(1e10), &error_rate, throw_policy()), std::overflow_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ static_cast(2) }, { static_cast(3) }, static_cast(1e10), &error_rate, throw_policy()), std::overflow_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ 2 }, { 3 }, boost::math::tools::max_value(), &error_rate, throw_policy()), std::overflow_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ static_cast(2) }, { static_cast(3) }, boost::math::tools::max_value(), &error_rate, throw_policy()), std::overflow_error); + BOOST_MATH_IF_CONSTEXPR(std::numeric_limits::is_specialized) + { + // Can't figure out how to make this work with real_concept (we get an evaluation_error): + BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ boost::math::tools::max_value() }, { static_cast(3) }, static_cast(0.5f), &error_rate, throw_policy()), std::overflow_error); + } + BOOST_MATH_IF_CONSTEXPR(std::numeric_limits::has_infinity) + { + BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ 2 }, { 3 }, static_cast(1e10), &error_rate, nothrow_policy()), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ static_cast(2) }, { static_cast(3) }, static_cast(1e10), &error_rate, nothrow_policy()), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ 2 }, { 3 }, boost::math::tools::max_value(), &error_rate, nothrow_policy()), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ static_cast(2) }, { static_cast(3) }, boost::math::tools::max_value(), &error_rate, nothrow_policy()), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ boost::math::tools::max_value() }, { static_cast(3) }, static_cast(0.5f), &error_rate, nothrow_policy()), std::numeric_limits::infinity()); + } +} + template void test_spots(T z, const char* type_name) { @@ -309,4 +334,5 @@ void test_spots(T z, const char* type_name) test_spots_1F2(z, type_name); test_spots_2F2(z, type_name); test_spots_2F1(z, type_name); + test_special_cases(z, type_name); } diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index a08cf75a01..3c63d3d639 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -548,6 +548,150 @@ void test_spots(RealType) static_cast(1.0))), 9); + // Tests for find_degrees_of_freedom(t, p) overload. + // Each case is derived from the CDF spot tests above: the exact df is known, + // so we verify that inverting CDF(x; df) = p recovers df to tight tolerance. + // float has ~7 significant digits; large-df inversion is ill-conditioned at that precision, + // so use a looser tolerance for single precision. + RealType tol_inv_df = std::is_same::value + ? static_cast(0.1) // float: 0.1% + : static_cast(0.01); // double+: 0.01% + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(-6.96455673428326), + static_cast(0.01)), + static_cast(2), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(-3.36492999890721), + static_cast(0.01)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(-0.559429644), + static_cast(0.3)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(1.475884049), + static_cast(0.9)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(-1.475884049), + static_cast(0.1)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(-5.2410429995425), + static_cast(0.00001)), + static_cast(25), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(6.96455673428326), + static_cast(0.99)), + static_cast(2), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(0.610822886098362)), + static_cast(0.1), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(0.777242554908434)), + static_cast(0.5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(0.822925875908677)), + static_cast(0.75), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(0.977114826753374)), + static_cast(1000), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(1), + static_cast(0.8413326478347855)), + static_cast(1e4), + tol_inv_df * static_cast(5)); // Looser tolerance for ill-conditioned problem with very large df + // Small df edge cases where approximation becomes problematic: would need different values for float + // For now just test double and long double to have test overage of these code paths + if (!std::is_same::value) + { + // Small df case 1: Edgeworth expansion breaks down for small degrees of freedom, use fallback + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(1e20), + static_cast(0.5244796002843015)), + static_cast(1e-3), + tol_inv_df); + // Small df case 2: Edgeworth expansion succeeds but is inaccurate, use fallback + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(2.0), + static_cast(0.500000000644961)), + static_cast(1e-10), + tol_inv_df); + } + // Analytical test: df=1 (Cauchy) case + { + boost::math::cauchy_distribution cauchy(0, 1); + RealType x = static_cast(1.0); + RealType p = cdf(cauchy, x); + RealType df_result = students_t_distribution::find_degrees_of_freedom(x, p); + BOOST_CHECK_EQUAL(df_result, static_cast(1.0)); + } + + + // Domain error: p outside (0,1) +#ifndef BOOST_NO_EXCEPTIONS + BOOST_MATH_CHECK_THROW( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(-0.1)), + std::domain_error); + BOOST_MATH_CHECK_THROW( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(1.1)), + std::domain_error); + // x == 0, p != 0.5: no df can satisfy CDF(0; df) == p -> domain error + BOOST_MATH_CHECK_THROW( + students_t_distribution::find_degrees_of_freedom( + static_cast(0), + static_cast(0.3)), + std::domain_error); + // x == 0, p == 0.5: CDF(0; df) == 0.5 for all df -> overflow (infinite solutions) + BOOST_MATH_CHECK_THROW( + students_t_distribution::find_degrees_of_freedom( + static_cast(0), + static_cast(0.5)), + std::overflow_error); + { + // Analytical test: df=infinity (Normal) case returns overflow error + boost::math::normal_distribution norm(0, 1); + RealType x = static_cast(1.0); + RealType p = cdf(norm, x); + BOOST_MATH_CHECK_THROW( + students_t_distribution::find_degrees_of_freedom(x, p), + std::overflow_error); + } +#endif + // Test for large degrees of freedom when should be same as normal. RealType inf = std::numeric_limits::infinity(); RealType nan = std::numeric_limits::quiet_NaN(); diff --git a/test/test_students_t_overflow.cpp b/test/test_students_t_overflow.cpp new file mode 100644 index 0000000000..d4020d3bf4 --- /dev/null +++ b/test/test_students_t_overflow.cpp @@ -0,0 +1,133 @@ +// Copyright Anton Leontev 2026. +// 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) +// +// Regression test for spurious intermediate overflow of x*x in the +// Student's t pdf() and cdf(): for |x| > sqrt(max_value) the functions +// must return the exact tail values (0 for the pdf, 0/1 for the cdf and +// its complement) without raising FE_OVERFLOW and without invoking any +// error handler. + +#define BOOST_TEST_MAIN +#include +#include +#include + +#include +#include +#include + +#ifndef BOOST_MATH_STANDALONE +#include // for real_concept +#include // for cpp_bin_float_50 +#endif + +#include +#include +#include + +#pragma STDC FENV_ACCESS ON + +template +RealType overflowing_deviate() +{ + BOOST_MATH_STD_USING // ADL sqrt for built-in and UDT types + using boost::math::tools::max_value; + return 2 * sqrt(max_value()); +} + +// For |x| > sqrt(max_value) the mathematically exact results have already +// saturated: pdf == 0, cdf == 0 or 1. These must be returned directly. +// A throwing policy is used so that any error-handler invocation would +// fail the test with an unexpected exception. +template +void test_tail_values(RealType) +{ + using namespace boost::math; + using namespace boost::math::policies; + + typedef policy > strict_policy; + typedef students_t_distribution dist_t; + + dist_t dist(static_cast(5)); + const RealType big = overflowing_deviate(); + + BOOST_CHECK((boost::math::isfinite)(big)); + + BOOST_CHECK_EQUAL(pdf(dist, big), static_cast(0)); + BOOST_CHECK_EQUAL(pdf(dist, -big), static_cast(0)); + + BOOST_CHECK_EQUAL(cdf(dist, big), static_cast(1)); + BOOST_CHECK_EQUAL(cdf(dist, -big), static_cast(0)); + + BOOST_CHECK_EQUAL(cdf(complement(dist, big)), static_cast(0)); + BOOST_CHECK_EQUAL(cdf(complement(dist, -big)), static_cast(1)); +} + +// Built-in types: the calls must not leave FE_OVERFLOW set -- avoiding the +// spurious intermediate overflow is the whole point of the guard. +template +void test_no_fp_exceptions(RealType) +{ + using namespace boost::math; + + students_t_distribution dist(static_cast(5)); + const RealType big = overflowing_deviate(); + + std::feclearexcept(FE_ALL_EXCEPT); + RealType p = pdf(dist, big); + RealType c1 = cdf(dist, big); + RealType c2 = cdf(dist, -big); + BOOST_CHECK_EQUAL(std::fetestexcept(FE_OVERFLOW), 0); + BOOST_CHECK_EQUAL(p, static_cast(0)); + BOOST_CHECK_EQUAL(c1, static_cast(1)); + BOOST_CHECK_EQUAL(c2, static_cast(0)); +} + +// Ordinary arguments must be unaffected by the guard. +template +void test_no_regression(RealType) +{ + using namespace boost::math; + + students_t_distribution dist(static_cast(10)); + const RealType tol = boost::math::tools::epsilon() * 100; + + BOOST_CHECK_CLOSE_FRACTION(cdf(dist, static_cast(0)), + static_cast(0.5), tol); + + RealType p = pdf(dist, static_cast(1.5)); + BOOST_CHECK((boost::math::isfinite)(p)); + BOOST_CHECK(p > 0); + + RealType big_but_safe = static_cast(1e6); + BOOST_CHECK_CLOSE_FRACTION(cdf(dist, big_but_safe), + static_cast(1), tol); +} + +BOOST_AUTO_TEST_CASE(students_t_overflow_test) +{ + test_tail_values(0.0F); + test_tail_values(0.0); + test_tail_values(0.0L); + + test_no_fp_exceptions(0.0F); + test_no_fp_exceptions(0.0); + test_no_fp_exceptions(0.0L); + + test_no_regression(0.0F); + test_no_regression(0.0); + test_no_regression(0.0L); + + // real_concept and multiprecision compute distribution constants via + // lexical_cast, which standalone mode disables, so both are skipped there. +#ifndef BOOST_MATH_STANDALONE + test_tail_values(boost::math::concepts::real_concept(0)); + test_no_regression(boost::math::concepts::real_concept(0)); + + test_tail_values(boost::multiprecision::cpp_bin_float_50(0)); + test_no_regression(boost::multiprecision::cpp_bin_float_50(0)); +#endif +}