Reduce the argument of sin, cos and tan exactly and round once (#1475) 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)`. `pi_v` 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.

Review and improve sin(), cos() and tan() #1471

Description

Activity

ckormanyos commented on Sep 29, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

Review and improve `sin()`, `cos()` and `tan()` #1471