Fix the accuracy and monotonicity of asin, acos, atan and atan2 - #1478
Merged
Merged
Conversation
ibmibmibm
force-pushed
the
trig-asin-acos-atan
branch
from
September 29, 2026 11:07
969c056 to
9767201
Compare
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## develop #1478 +/- ##
=========================================
+ Coverage 98.6% 98.7% +0.1%
=========================================
Files 312 315 +3
Lines 26265 26454 +189
Branches 2264 2256 -8
=========================================
+ Hits 25892 26101 +209
+ Misses 373 353 -20
... and 8 files with indirect coverage changes Continue to review full report in Codecov by Harness.
🚀 New features to boost your workflow:
|
ibmibmibm
force-pushed
the
trig-asin-acos-atan
branch
from
September 30, 2026 03:50
9767201 to
6d9ee65
Compare
ibmibmibm
marked this pull request as ready for review
September 30, 2026 06:38
ibmibmibm
force-pushed
the
trig-asin-acos-atan
branch
from
September 30, 2026 07:41
6d9ee65 to
97aba40
Compare
Member
|
@ckormanyos do you want to take a look at this one as well? It seems reasonable to me |
Member
|
Hi Matt (@mborland) hi @ibmibmibm and thanks for looking into this. Actually I have gotten into the tradition of using Other than that, when changed, I think this is a good improvement. Hi @ibmibmibm if you agree when and where sensible, save just a few cycles in the functions we refactor by using |
- The asin polynomials for 64 and 128 bits were too short near 0.5, and a constant term broke small arguments. acos uses the same polynomials. - asin below cbrt(eps) rounded 1 + x^2/6 to the type, which lost digits. Add x^3/6 to x instead. - acos for x >= 0.5 subtracted pi/2 and added it back, which lost most digits near 1. - The 128-bit atan Pade kernel was too short above 0.3, and two atan constants were wrong. - pi_v<T> / 2 rounds twice, so asin(1), atan(inf) and atan2(y, 0) were one ulp off. atan2 took pi/4 and 3 pi/4 from numbers::pi_over_four_v, which has the wrong exponent for decimal128_t. - Use odd kernels x + x^3 P(x^2) from Chebyshev fits. Evaluate P in Q62 or Q126 fixed point and round x^2 P once, which also makes the functions faster. - Put pi/4, pi/2, 3 pi/4 and pi, split into a rounded part and the rest, in one table for asin, acos, atan and atan2. - asin and acos above 0.5 use 2 asin(sqrt((1 - x) / 2)), with a low part that corrects the rounding of the square root and of s + s. - atan uses the fdlibm reduction at 1/2, 1 and 3/2, and adds the rounding of the denominator and of the quotient to the low part. - atan2 rounded y / x before atan, and for a quotient q in [10^k, tan(10^k)) the grid of q is ten times coarser than that of the result. For these q only, give the remainder of the division to atan as a low part of its argument. - Write the integer constants in atan as T values, so that atan does not instantiate the operators for a decimal and an int. - Round each result only once at its own grid, so that consecutive inputs give results in order. - Three tests expected the old values of pi/2, pi/4, 3 pi/4 and asin(sqrt(eps)); update them. - Add a regression test with MPFR values, stored as a high and a low part, that requires an error below one ulp, and that walks consecutive inputs around each point where a formula changes. Refs boostorg#1472
ibmibmibm
force-pushed
the
trig-asin-acos-atan
branch
from
October 1, 2026 03:19
97aba40 to
5a8ee56
Compare
Contributor
Author
|
All five |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
See #1472.
Summary
On develop,
asin,acos,atanandatan2lose digits. The decimal64 functions have errors up to 9,230 ulp, and the decimal128 functions have errors up to 1.5e18 ulp. The results also go against the direction of the exact function for thousands of adjacent inputs. This PR replaces the kernels and the reductions. The error is now less than 1 ulp in round to nearest for all types, and the results are in order. All four functions are faster, by up to 20 times.Problems on develop
asinpolynomials are too short near 0.5, and they have a constant term. This constant term breaks small arguments.acosuses the same polynomials.asinrounds 1 + x^2/6 to the type, then it loses digits.acossubtracts pi/2 and adds it back. Near 1 this removes most digits.atanPade kernel is too short above 0.3. Two decimal128atanconstants are wrong.pi_v<T> / 2rounds two times, thenasin(1),atan(inf)andatan2(y, 0)are 1 ulp off.atan2takes pi/4 and 3 pi/4 fromnumbers::pi_over_four_v. For decimal128_t, this constant is 100 times too small.atan2computesatan(y / x). For q = y / x in [1, 1.557), q has a grid ten times coarser than the result. The rounding of q then adds up to 2.5 ulp.fmacalls on decimal values.Examples (decimal64_t unless noted). The correct values come from MPFR at 600 bits.
The tests did not find these problems.
test_asin.cpp,test_acos.cpp,test_atan.cppandtest_atan2.cppcompare withfloatstd::asin,std::acos,std::atanandstd::atan2, with a limit of 50 float ulps. Then they cannot see an error below the precision offloat. No test uses a directed rounding mode or a reference with more digits than the type.Changes
impl/fixed_point_series.hppcomputes P in binary fixed point: Q62 in anint64_tfor the 32-bit and 64-bit types, and Q126 in anint128_tfor the 128-bit types. The product x^2 P has one rounding.impl/split_pi.hpphas pi/4, pi/2, 3 pi/4 and pi as a rounded part and the rest.asin,acos,atanandatan2use this table. This PR does not changenumbers::pi_over_four_v.asinandacosuse 2 asin(sqrt((1 - x) / 2)). A low part corrects the rounding of the square root and of s + s.atanuses the fdlibm reduction at 1/2, 1 and 3/2. It adds the rounding of the denominator and of the quotient to the low part.atan2gives the remainder of y / x toatanas a low part of its argument. Other q do not need it, and they keep the speed of oneatancall.test_asin.cpp,test_atan.cppandtest_atan2.cppexpected the old values of pi/2, pi/4, 3 pi/4 and asin(sqrt(eps)). They now expect the correctly rounded values.test_inverse_trig_rounding.cppcompares with MPFR values, stored as a high and a low part. The error must be less than 1 ulp. The test also walks adjacent inputs around each point where a formula changes, and the results must be in order.Accuracy
The dense survey uses 20,000 random full-precision inputs per range and per rounding mode, on the ranges with the largest errors. The reference is MPFR.
Errors in ulp. The fast types give the same values.
For decimal32_t, a check of all inputs against MPFR gives these maximum errors in round to nearest: 0.89 ulp for
asin, 0.81 ulp foracosand 0.64 ulp foratan. No result has an error of 1 ulp or more. In the directed modes the maximum is 1.58 ulp.atan2survey: 5,000 inputs per band of |y/x|, in four quadrants and five rounding modes.In the directed modes, outside the ranges of q above, the rounding of q and the rounding of
atancan each add 1 ulp.Monotonicity: runs of adjacent inputs around each point where a formula changes, in five rounding modes. A plateau is a step where the result stays the same, but the exact value changes by 2 ulp or more.
atan, steps against the exact directionacos, plateausatan2, steps against the exact directionacos, steps against the exact directionThe one
acosstep is at x = -0.5 in upward mode. The result changes by 1 ulp in the wrong direction.Speed
Google Benchmark, median ns per call, one core, GCC
-O3. The fast types give similar numbers.The
atan2numbers are for x > 0. For x < 0 they are similar. Above 0.5,asinandacosspend most of the time insqrtand in the divisions.Code size
Bytes that the functions add to a program (.text and .rodata), for all six types.
asin,acosandatan, developatan2, developasin,acosandatanare smaller than on develop.atan2is larger. It has a secondatanpath with a low part, and it usesfmafor the remainder of y / x. For the 64-bit and 128-bit types,fmaneeds 256-bit arithmetic.Other functions
ellint_1andellint_2callatanin each step of the AGM. With this PR, their decimal128 error goes from up to 5.6e7 ulp to 11 to 26 ulp. They are also 10 to 20 percent faster. The decimal64 error does not change.Tests
All single-file tests in
test/Jamfilepass with GCC in C++17 and-Werror.test_format_fmtlibdoes not build with the local fmt 12, but this PR does not change it. The tests of these functions also pass with GCC and Clang in C++14 and C++20, in a 32-bit build and withBOOST_DECIMAL_FAST_MATH. The four functions also work in constant evaluation with GCC and Clang.