Skip to content

Fix complex asin, acos, tan, tanh and log1p mismatches with std - #1456

Merged
DiamonDinoia merged 2 commits into
xtensor-stack:masterfrom
emrys53:fix/complex-math-std-mismatch-1432
Nov 13, 2026
Merged

DiamonDinoia merged 2 commits into
xtensor-stack:masterfrom
emrys53:fix/complex-math-std-mismatch-1432

Conversation

@emrys53

@emrys53 emrys53 commented Nov 11, 2026 •

Copy link
Copy Markdown
Contributor

Fixes #1432.
Supersedes #1434, #1435, #1436, #1437, #1438, #1439, #1440.

The implementations are mostly taken from FreeBSD implementation (scalar, converted to simd version) and some from Julia. Main idea is that we always try to sum up variables that have the same sign to avoid cancellations and in some cases we check for some special fixed bounds so that the term won't become Inf. I also added tests, that tests for the inputs mentioned in #1432.

Performance: ns per element, master vs this PR

Lower is faster; the last column of each architecture is PR / master (below 1 means the PR is faster).

function type SSE2 master SSE2 PR SSE2 PR/master AVX2+FMA master AVX2+FMA PR AVX2+FMA PR/master AVX512F master AVX512F PR AVX512F PR/master
asin float 11.33 8.53 0.76 5.52 4.33 0.78 3.11 2.82 0.90
asin double 28.99 25.79 0.90 14.84 13.33 0.90 12.78 11.46 0.90
acos float 11.25 8.51 0.76 5.54 4.35 0.79 3.14 2.76 0.87
acos double 28.95 26.20 0.90 14.81 13.36 0.91 12.67 11.54 0.91
asinh float 11.26 8.36 0.74 5.60 4.35 0.78 3.12 2.75 0.88
asinh double 29.13 25.71 0.90 14.83 13.42 0.91 12.76 11.49 0.90
acosh float 11.18 8.47 0.78 5.54 4.25 0.77 3.14 2.79 0.88
acosh double 28.81 26.12 0.90 14.81 13.43 0.91 12.79 11.53 0.90
tan float 14.99 11.51 0.76 9.11 6.45 0.71 9.81 6.29 0.64
tan double 23.41 20.39 0.88 7.76 8.00 1.04 6.23 6.64 1.07
tanh float 14.46 10.97 0.78 9.26 6.30 0.68 10.14 6.49 0.63
tanh double 22.87 20.18 0.87 7.47 8.11 1.07 5.70 6.80 1.20
log1p float 5.08 5.42 1.07 2.39 2.67 1.05 1.48 1.49 1.04
log1p double 13.81 14.51 1.05 7.00 7.48 1.04 5.06 5.13 1.01
atan2 float 1.75 1.76 1.00 0.77 0.81 1.05 0.55 0.57 1.04
atan2 double 4.55 4.98 1.10 2.31 2.24 1.03 2.35 2.30 0.98
  • Machine: Intel Xeon w5-3435X. Compiler: GCC 14.2, -O3 -DNDEBUG.

Fix xtensor-stack#1432: asin/acos/asinh/acosh use Kahan's formulas, tan/tanh use
FreeBSD's ctanh, log1p keeps the sign of a -0 imaginary part, and
atan2(+-0, +-0) no longer returns NaN. Add tests with every input
from the issue.

Assisted-By: Claude Opus 5.5 <noreply@anthropic.com>
@emrys53

emrys53 commented Nov 11, 2026

Copy link
Copy Markdown
Contributor Author

All failing CI jobs except one use libc++, whose trigonometric functions use different method and are less accurate than the current xsimd for large |z|. I checked some of the values with "exact" results and xsimd is a lot closer to "exact" results than what libc++ provides. I don't know how to fix this issue. Maybe instead of testing against std:: implementations, we can test it against precomputed values for each case so it doesn't change between different environments.

The one exception is clang-cl x64 /fp:fast. The test hits log1p(-1), which should be -inf. With fast-math on, xsimd's log doesn't handle infinities, so it returns a finite value instead. I think with -ffast-math we can assume there should be no NaN/Inf for arguments/results so in my opinion xsimd's approach seems better than what clang does. Maybe we can disable these -inf cases for this specific envrionment.

@DiamonDinoia

DiamonDinoia commented Nov 11, 2026 •

Copy link
Copy Markdown
Contributor

We can change the tests to be accurate to machine precision multiplied by the conditioning of the problem.
I would say that we need to be backward stable (threfeten & bau). If we account for the conditioning of the problem we lost roughly one digit for every power of 10 of the magnitude of the argument.

Once we do that we can just measure the relative error point-wise to that accuracy.

fast-math is best effort user that enable it know that they cannot distinguish inf from a large value.

@DiamonDinoia

Copy link
Copy Markdown
Contributor

I do not like testing against generated data from python.
Nor I like a 3rd party dependency for those.
How hard is to hardcode the corner cases we care? There should be very few of those.

@DiamonDinoia

Copy link
Copy Markdown
Contributor

Also, can you check the license of the code you took? I do not want to violate any of them. We can just implement the textbook algorithm for these things and fit the coefficients with Remez.

@DiamonDinoia DiamonDinoia left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Over engineered tests :)

@emrys53

emrys53 commented Nov 12, 2026

Copy link
Copy Markdown
Contributor Author

Also, can you check the license of the code you took? I do not want to violate any of them. We can just implement the textbook algorithm for these things and fit the coefficients with Remez.

Licences for FreeBSD are BSD-2 and for Julia it is MIT. I will add the Licences and credit on top of each function in the next commit.

@emrys53
emrys53 force-pushed the fix/complex-math-std-mismatch-1432 branch from 27bf679 to 1208d75 Compare November 12, 2026 21:05
Check asin, acos, asinh, acosh, tan, tanh and log1p on the inputs of
issue 1432 that gave a wrong value or a wrong sign of zero before the
fix: signed zeros, inputs where the sign of a zero imaginary part
decides the sign of the result, points near the poles of tan and
tanh, and large |z|. The expected values are the exact results
rounded to double, written in the tests. A result passes if it is
within 4 * cond * eps of the expected value, where cond is the
condition number |z f'(z) / f(z)|. Infinite values and the sign of
zeros must match; with fast-math, infinite values are skipped.

Cite the papers whose formulas complex asin and tanh use (Kahan 1987;
Hull, Fairgrieve and Tang 1997) instead of the Julia and FreeBSD code.

Assisted-By: Claude Opus 5.5 <noreply@anthropic.com>
@emrys53
emrys53 force-pushed the fix/complex-math-std-mismatch-1432 branch from 1208d75 to ac774c8 Compare November 12, 2026 21:26
Comment thread test/test_utils.hpp Resolved

@DiamonDinoia DiamonDinoia left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

since there are citations I do not need to restate the equations. Also I can read the equations from the code itself.

@emrys53

emrys53 commented Nov 12, 2026

Copy link
Copy Markdown
Contributor Author

since there are citations I do not need to restate the equations. Also I can read the equations from the code itself.

I tried to keep the comments to minimum and only explain math and why some steps are necessary and why is it done this way to make it easier to review (and also myself to keep track of what is going on). I tried to keep them to minimum but if you think they are too much i can remove them. But in my opinion it makes the code more clear and it is not too much comment.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Complex asin, acos, asinh, acosh, sqrt, log1p, tan, tanh return very different results from std::

2 participants