Skip to content

Bugfix/Workaround sign error Asso. Legendre Function at z=-1 #76

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.

Sign up for GitHub

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

Merged
steppi merged 7 commits into from
Dec 5, 2025
Merged

Bugfix/Workaround sign error Asso. Legendre Function at z=-1 #76

Changes from 1 commit
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
f9251dc
Make the starting element in assoc_legendre_p_pm1 depend on the sign …
mberz Nov 18, 2025
64cdf17
BUG: adjust fix for z=-1, m=0 case of assoc_legendre_p
steppi Dec 5, 2025
860e9eb
BUG: apply previous fix to non-dual overload
steppi Dec 5, 2025
25a4ad7
Merge branch 'main' into bugfix/legendre_sign_error
steppi Dec 5, 2025
557710b
TST: Add test for assoc_legendre_p m=0 edge case
steppi Dec 5, 2025
d7d79e7
STY: apply clang-format
steppi Dec 5, 2025
1d0671e
STY: more clang-format
steppi Dec 5, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 5 additions & 1 deletion include/xsf/legendre.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -310,7 +310,11 @@ void assoc_legendre_p_pm1(NormPolicy norm, int n, int m, T z, int branch_cut, T
template <typename NormPolicy, typename T, size_t Order>
void assoc_legendre_p_pm1(NormPolicy norm, int n, int m, dual<T, Order> z, int branch_cut, dual<T, Order> &res) {
if (m == 0) {
res[0] = T(1);
if (real(z[0]) >= 0) {
res[0] = T(1);
} else {
res[0] = T(-1);
}
Comment on lines 317 to 321
Copy link
Collaborator

@steppi steppi Dec 5, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if (real(z[0]) >= 0) {
res[0] = T(1);
} else {
res[0] = T(-1);
}
if (real(z[0]) >= 0) {
res[0] = T(1);
} else {
res[0] = T(std::pow(-1, n));
}

The current patch breaks the even n case. You can get things to work again by adding dependence on n.
Also, it looks like the adjustment should also be added to the overload for non-dual types currently on lines 301-308:

xsf/include/xsf/legendre.h

Lines 301 to 308 in 98623fd

template <typename NormPolicy, typename T>
void assoc_legendre_p_pm1(NormPolicy norm, int n, int m, T z, int branch_cut, T &res) {
if (m == 0) {
res = T(1);
} else {
res = T(0);
}
}

} else {
res[0] = T(0);
}
Expand Down
Loading