diff --git a/include/xsf/amos/amos.h b/include/xsf/amos/amos.h index 49fdd2e59f..6c1a873283 100644 --- a/include/xsf/amos/amos.h +++ b/include/xsf/amos/amos.h @@ -1140,7 +1140,7 @@ namespace amos { return nz; } ck = std::exp(cz); - for (int i = 0; i < (nn + 1); i++) { + for (int i = 0; i < nn; i++) { y[i] *= ck; } /* 90 */ @@ -1770,7 +1770,7 @@ namespace amos { *ierr = 2; return 0; } - if (xx > 0.0) { + if (xx >= 0.0) { return nz; } // @@ -5622,7 +5622,7 @@ namespace amos { return -1; } nz = n; - for (i = 0; i < (n + 1); i++) { + for (i = 0; i < n; i++) { y[i] = 0.0; } return nz; @@ -6067,8 +6067,7 @@ namespace amos { // REFINE ESTIMATE AND TEST // aphi = std::abs(phid); - aarg = std::abs(argd); - rs1 += std::log(aphi) - 0.25 * std::log(aarg) - aic; + rs1 += std::log(aphi); if (std::fabs(rs1) < elim) { goto L120; } diff --git a/tests/xsf_tests/test_amos.cpp b/tests/xsf_tests/test_amos.cpp index e4bb069928..cb3d4afa4f 100644 --- a/tests/xsf_tests/test_amos.cpp +++ b/tests/xsf_tests/test_amos.cpp @@ -31,3 +31,60 @@ TEST_CASE("amos besj vectorized", "[amos][xsf_tests]") { REQUIRE(rel_error <= rtol); } } + +TEST_CASE("amos besi vectorized", "[amos][xsf_tests]") { + // tests the functionality of amos to return multiple consecutive orders for besi + // by comparing to the versions returning only a single order + using std::complex; + + using test_case = std::tuple, double, int, int, double>; + auto [z, fnu, kode, n, rtol] = GENERATE( + test_case{complex{0.0, -15.0}, 0.0, 1, 10, 1e-15} // gh-158 + ); + + std::vector> cy(n); + int ierr = 0; + int nz; + + nz = xsf::amos::besi(z, fnu, kode, n, cy.data(), &ierr); + + REQUIRE(ierr == 0); + + complex ref; + + for (int i = 0; i < n; ++i) { + nz = xsf::amos::besi(z, fnu + i, kode, 1, &ref, &ierr); + REQUIRE(ierr == 0); + + const auto rel_error = xsf::extended_relative_error(cy[i], ref); + CAPTURE(i, cy[i], ref, rel_error); + REQUIRE(rel_error <= rtol); + } +} + +TEST_CASE("amos asyi buffer overflow gh-158", "[amos][xsf_tests]") { + using std::complex; + + // parameters for asyi (computed as in besi) + const double tol = 2.220446049250313e-16; + const double rl = 24.6; + const double elim = 700.0; + const double alim = 664.0; + + const complex z{700.0, 10.0}; + const double fnu = 0.0; + const int n = 5; + const int kode = 1; + + // allocate n+1 elements, initialize extra to sentinel to detect overflow + const complex sentinel{12345.67, 98765.43}; + std::vector> y(n + 1, sentinel); + + int nz = xsf::amos::asyi(z, fnu, kode, n, y.data(), rl, tol, elim, alim); + + REQUIRE(nz == 0); + + // check if the extra element (index n) was touched + CAPTURE(y[n]); + CHECK(y[n] == sentinel); +}