Correct gauss and gauss_kronrod integrators with strange-precision ty…

…pes.
boostorg · Feb 9, 2024 · 1d1a52c · 1d1a52c
1 parent 26e21c9
commit 1d1a52c
Show file tree

Hide file tree

Showing 3 changed files with 20 additions and 6 deletions.
diff --git a/include/boost/math/quadrature/gauss.hpp b/include/boost/math/quadrature/gauss.hpp
@@ -798,8 +798,8 @@ class gauss : public detail::gauss_detail<Real, N, detail::gauss_constant_catego
       Real L1 = abs(result);
       for (unsigned i = non_zero_start; i < base::abscissa().size(); ++i)
       {
-         K fp = f(base::abscissa()[i]);
-         K fm = f(-base::abscissa()[i]);
+         K fp = f(static_cast<Real>(base::abscissa()[i]));
+         K fm = f(static_cast<Real>(-base::abscissa()[i]));
          result += (fp + fm) * static_cast<Real>(base::weights()[i]);
          L1 += (abs(fp) + abs(fm)) * static_cast<Real>(base::weights()[i]);
       }

diff --git a/include/boost/math/quadrature/gauss_kronrod.hpp b/include/boost/math/quadrature/gauss_kronrod.hpp
@@ -1159,16 +1159,16 @@ class gauss_kronrod : public detail::gauss_kronrod_detail<Real, N, detail::gauss
       Real L1 = abs(kronrod_result);
       for (unsigned i = gauss_start; i < base::abscissa().size(); i += 2)
       {
-         fp = f(base::abscissa()[i]);
-         fm = f(-base::abscissa()[i]);
+         fp = f(static_cast<Real>(base::abscissa()[i]));
+         fm = f(static_cast<Real>(-base::abscissa()[i]));
          kronrod_result += (fp + fm) * static_cast<Real>(base::weights()[i]);
          L1 += (abs(fp) + abs(fm)) * static_cast<Real>(base::weights()[i]);
          gauss_result += (fp + fm) * static_cast<Real>(gauss<Real, (N - 1) / 2>::weights()[i / 2]);
       }
       for (unsigned i = kronrod_start; i < base::abscissa().size(); i += 2)
       {
-         fp = f(base::abscissa()[i]);
-         fm = f(-base::abscissa()[i]);
+         fp = f(static_cast<Real>(base::abscissa()[i]));
+         fm = f(static_cast<Real>(-base::abscissa()[i]));
          kronrod_result += (fp + fm) * static_cast<Real>(base::weights()[i]);
          L1 += (abs(fp) + abs(fm)) * static_cast<Real>(base::weights()[i]);
       }

diff --git a/test/git_issue_1075.cpp b/test/git_issue_1075.cpp
@@ -26,6 +26,17 @@ void test()
    std::cout << r << std::endl;
    r = boost::math::quadrature::gauss_kronrod<T, 61>::integrate(f, 0, 1, 0, 0, &error);
    std::cout << r << std::endl;
+
+   r = boost::math::quadrature::gauss<T, 7>::integrate(f, 0, 1, &error);
+   std::cout << r << std::endl;
+   r = boost::math::quadrature::gauss<T, 15>::integrate(f, 0, 1, &error);
+   std::cout << r << std::endl;
+   r = boost::math::quadrature::gauss<T, 20>::integrate(f, 0, 1, &error);
+   std::cout << r << std::endl;
+   r = boost::math::quadrature::gauss<T, 25>::integrate(f, 0, 1, &error);
+   std::cout << r << std::endl;
+   r = boost::math::quadrature::gauss<T, 30>::integrate(f, 0, 1, &error);
+   std::cout << r << std::endl;
 }
 
 
@@ -35,5 +46,8 @@ int main()
    test<boost::multiprecision::cpp_bin_float_double_extended>();
    test<boost::multiprecision::cpp_bin_float_quad>();
 
+   using strange = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<50, boost::multiprecision::backends::digit_base_2, void, std::int16_t, -16382, 16383>, boost::multiprecision::et_off>;
+   test<strange>();
+
    return 0;
 }