Reduce the argument of sin, cos and tan exactly and round once

[ibmibmibm]

See #1471.

Summary

On develop, sin, cos and tan lose digits from x = 1. Above about 1e9 they return values larger than 1, inf or NaN. This PR replaces the argument reduction and the kernels. The results are correctly rounded for the 32-bit and 64-bit types in all five rounding modes. All three functions are 1.8 to 32 times faster, and the code is less than half the size.

Problems on develop

• The reduction computes k = unsigned(2x / pi_v<T>). pi_v<T> has only the digits of the type, then the error of the remainder grows with k. The cast to unsigned overflows above x = 6.7e9.
• The decimal32 and decimal64 sine polynomials have a constant term and a bad fit. The decimal64 kernel has an error of 1.6e-10 at pi/4.
• The kernels are chains of fma calls. One decimal64 fma takes about 137 ns.
• tan calls sin and cos, then it reduces the argument three times. Near odd multiples of pi/2 it divides by zero.
• sin(±inf) and cos(±inf) return inf. IEEE 754 (9.2.1) and C Annex F require NaN.

Examples (decimal64_t unless noted). The correct values come from MPFR at 600 bits.

Call	develop	this PR	Correct value
sin(0.5)	0.4794255386073293	0.4794255386042030	0.4794255386042030
sin(1)	0.8414709848081057	0.8414709848078965	0.8414709848078965
tan(1)	1.557407724655289	1.557407724654902	1.557407724654902
sin(100)	-0.5063656411097701	-0.5063656411097588	-0.5063656411097588
sin(1e10)	-1.147951735515450e+28	-0.4875060250875107	-0.4875060250875107
cos(1e10)	7.135480319211579e+275	0.8731196226768560	0.8731196226768560
decimal32_t cos(100)	0.8623027	0.8623189	0.8623189
decimal32_t sin(1e10)	-1.147950e+28	-0.4875060	-0.4875060
decimal128_t sin(100)	-0.5063656411097587936565576104597925	-0.5063656411097587936565576104597854	-0.5063656411097587936565576104597854
sin(inf)	inf	nan	nan

The tests did not find these problems. test_sin_cos.cpp compares with float std::sin on [-2pi, 2pi] with an absolute tolerance of 4.8e-6. No test uses an argument above 22, a directed rounding mode, or a reference with more digits than the type.

Changes

• The new impl/trig_reduce.hpp does a Payne-Hanek reduction (an exact reduction with the digits of 2/pi) on integer words.
  • For |x| < 10^19, it uses 64-bit words of 2/pi and gives r in binary fixed point. The tables have 352 bytes. The guard words come from the worst cases below 2^64, with 27, 58 and 117 leading zero bits of r.
  • For larger x, it uses 2/pi in base-1e9 words. One table of 720 words covers the full range of all types. A binary table cannot do this, because the factor 5^e of 10^e is not a shift in base 2.
• The new impl/trig_fixed_point.hpp computes the kernels in binary fixed point. It uses one word for decimal32, one or two words for decimal64 and up to three words for decimal128. The result has one rounding, in the current rounding mode.
• The tests for 0 and for |x| < pi/4 use the integer significand, not decimal compares.
• tan does one reduction. The new impl/tan_impl.hpp has its kernel and its reciprocal.
• sin(±inf) and cos(±inf) return NaN. test_sin_cos.cpp and test_edges_and_behave.cpp now expect NaN.
• The new test test_trig_rounding.cpp compares with exact values from MPFR. It also has the worst cases below 10^19 and the inputs on each side of 10^19. It walks chains of adjacent inputs in three rounding modes. The results must not go against the direction of the exact function. The test skips the parts that need constant evaluation or that fast math turns off.

IEEE 754 (9.2) recommends these functions but does not require them. A version that conforms must be correctly rounded in the current rounding direction. C Annex F makes the rounding direction implementation-defined for sin, cos and tan. This PR meets the IEEE 754 level for the 32-bit and 64-bit types.

Accuracy

The sweep uses 4000 random full-precision inputs per range and per rounding mode, with 10 ranges per function from below epsilon to the largest value. The reference is MPFR.

	develop	this PR
Max error, round to nearest	10 to 3.3e5 ulp below 10, inf or NaN at large x	0.5 ulp in all ranges
Max error, directed modes	same	less than 1 ulp
Not correctly rounded, 32-bit and 64-bit types	almost all inputs above 1	0
Not correctly rounded, decimal128_t	almost all inputs above 1	12 of about 600,000

The 12 decimal128_t cases are near halfway points: 5 in sin on [1, 10), 2 in sin above 1e40, 3 in cos and 2 in tan near k*pi/2. The error is at most 1 ulp in these cases.

Monotonicity: for each function and type, the walk tests 147,900 pairs of adjacent inputs in three rounding modes.

	develop	this PR
Pairs that go against the exact direction	1,316 to 15,620 per function and type	0

Speed

Google Benchmark, median ns per call, one core, GCC -O3 -march=native. The fast types give similar numbers.

Function	Type	[0, pi/4)	[pi/4, 2pi)	[2pi, 1e3)	[1e3, 1e9)	[1e20, 1e40)
sin	decimal32_t	634 -> 66	899 -> 145	878 -> 148	653 -> 145	282 -> 157
cos	decimal32_t	813 -> 67	629 -> 141	822 -> 144	635 -> 138	535 -> 159
tan	decimal32_t	1771 -> 99	1949 -> 206	1724 -> 199	1174 -> 199	1148 -> 218
sin	decimal64_t	3018 -> 124	2936 -> 218	3325 -> 148	3381 -> 198	679 -> 238
cos	decimal64_t	2933 -> 135	3268 -> 175	3285 -> 150	2399 -> 197	1783 -> 240
tan	decimal64_t	6736 -> 208	6939 -> 233	6949 -> 286	6932 -> 293	2566 -> 336
sin	decimal128_t	3960 -> 341	3036 -> 452	4000 -> 470	3667 -> 430	1864 -> 423
cos	decimal128_t	3269 -> 338	3800 -> 477	2661 -> 463	3878 -> 364	2057 -> 423
tan	decimal128_t	8596 -> 607	6846 -> 738	8319 -> 740	8348 -> 605	5849 -> 663

Above 1e20, develop returns wrong values after a short computation. This PR does the full reduction there, then the gain is smaller.

Code size

Bytes that sin, cos and tan add to a program (.text and .rodata), GCC.

	develop	this PR
-O2	168,019	72,969
-Os	91,798	40,416

The tables add about 2 to 3 KB of .rodata. The removal of the fma chains makes the code smaller than this.

Tests

All single-file tests in test/Jamfile pass with GCC in C++17 and -Werror. test_format_fmtlib does not build with the local fmt 12, but this PR does not change it. The trig tests also pass with GCC and Clang in C++14, C++17 and C++20. Constant evaluation gives the same values in all compilers.

[codecov]

Codecov Report

❌ Patch coverage is 98.94180% with 8 lines in your changes missing coverage. Please review.
✅ Project coverage is 98.6%. Comparing base (1e058ab) to head (2ab9315).

Files with missing lines	Patch %	Lines
...de/boost/decimal/detail/cmath/impl/trig_reduce.hpp	97.2%	7 Missing ⚠️
...ost/decimal/detail/cmath/impl/trig_fixed_point.hpp	99.2%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff            @@
##           develop   #1475     +/-   ##
=========================================
- Coverage     98.6%   98.6%   -0.0%     
=========================================
  Files          308     312      +4     
  Lines        25727   26265    +538     
  Branches      2205    2264     +59     
=========================================
+ Hits         25362   25892    +530     
- Misses         365     373      +8     

Files with missing lines	Coverage Δ
include/boost/decimal/detail/cmath/cos.hpp	100.0% <100.0%> (ø)
...clude/boost/decimal/detail/cmath/impl/cos_impl.hpp	100.0% <100.0%> (ø)
...clude/boost/decimal/detail/cmath/impl/sin_impl.hpp	100.0% <100.0%> (ø)
...clude/boost/decimal/detail/cmath/impl/tan_impl.hpp	100.0% <100.0%> (ø)
include/boost/decimal/detail/cmath/sin.hpp	100.0% <100.0%> (+4.7%)	⬆️
include/boost/decimal/detail/cmath/tan.hpp	100.0% <100.0%> (+6.1%)	⬆️
test/test_edges_and_behave.cpp	100.0% <100.0%> (ø)
test/test_sin_cos.cpp	100.0% <100.0%> (ø)
test/test_trig_rounding.cpp	100.0% <100.0%> (ø)
...ost/decimal/detail/cmath/impl/trig_fixed_point.hpp	99.2% <99.2%> (ø)
... and 1 more

... and 3 files with indirect coverage changes

---

Continue to review full report in Codecov by Harness.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 1e058ab...2ab9315. Read the comment docs.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[mborland]

You have this note in your description: "test_format_fmtlib does not build with the local fmt 12, but this PR does not change it." Our {fmt} support should work from 11 so could you please investigate that when you have a chance? We install a number of different versions in the CI system for testing

[ibmibmibm]

You have this note in your description: "test_format_fmtlib does not build with the local fmt 12, but this PR does not change it." Our {fmt} support should work from 11 so could you please investigate that when you have a chance? We install a number of different versions in the CI system for testing

fmt_format.hpp:156:23: error: conversion to 'unsigned int' from 'int' may change the sign of the result [-Werror=sign-conversion]
fmt_format.hpp:156:28: error: conversion to 'int' from 'unsigned int' may change the sign of the result [-Werror=sign-conversion]
fmt_format.hpp:167:43: error: conversion to 'unsigned int' from 'int' may change the sign of the result [-Werror=sign-conversion]
fmt_format.hpp:167:48: error: conversion to 'int' from 'unsigned int' may change the sign of the result [-Werror=sign-conversion]

Looks like the *it returns a char32_tfor U"{}" format string in newest fmt, and sub with '0' makes an unsigned int