-
-
Notifications
You must be signed in to change notification settings - Fork 28
BUG: amos.h: fix translation errors in asyi, unk1, and besi #158
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Changes from 3 commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -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) { | ||
|
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. |
||
| 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++) { | ||
|
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. |
||
| 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); | ||
|
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. |
||
| if (std::fabs(rs1) < elim) { | ||
| goto L120; | ||
| } | ||
|
|
||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -31,3 +31,33 @@ TEST_CASE("amos besj vectorized", "[amos][xsf_tests]") { | |
| REQUIRE(rel_error <= rtol); | ||
| } | ||
| } | ||
|
|
||
| TEST_CASE("amos besi vectorized", "[amos][xsf_tests]") { | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Can you please add a small regression test for
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Done. Your intuition was right, I did not find a way to hit the fixes in |
||
| // 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<complex<double>, 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<complex<double>> cy(n); | ||
| int ierr = 0; | ||
| int nz; | ||
|
|
||
| nz = xsf::amos::besi(z, fnu, kode, n, cy.data(), &ierr); | ||
|
|
||
| REQUIRE(ierr == 0); | ||
|
|
||
| complex<double> 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); | ||
| } | ||
| } | ||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
See https://github.com/scipy/scipy/blob/4edfcaa3ce8a387450b6efce968572def71be089/scipy/special/amos/zasyi.f#L153.