Warning: This file is not a C or C++ file. It does not have highlighting.
1 | //===-- Utilities for trigonometric functions -------------------*- C++ -*-===// |
---|---|
2 | // |
3 | // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. |
4 | // See https://llvm.org/LICENSE.txt for license information. |
5 | // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception |
6 | // |
7 | //===----------------------------------------------------------------------===// |
8 | |
9 | #ifndef LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_H |
10 | #define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_H |
11 | |
12 | #include "src/__support/FPUtil/FPBits.h" |
13 | #include "src/__support/FPUtil/multiply_add.h" |
14 | #include "src/__support/FPUtil/nearest_integer.h" |
15 | #include "src/__support/common.h" |
16 | |
17 | namespace LIBC_NAMESPACE { |
18 | |
19 | namespace generic { |
20 | |
21 | static constexpr uint32_t FAST_PASS_BOUND = 0x4a80'0000U; // 2^22 |
22 | |
23 | static constexpr int N_ENTRIES = 8; |
24 | |
25 | // We choose to split bits of 32/pi into 28-bit precision pieces, so that the |
26 | // product of x * THIRTYTWO_OVER_PI_28[i] is exact. |
27 | // These are generated by Sollya with: |
28 | // > a1 = D(round(32/pi, 28, RN)); a1; |
29 | // > a2 = D(round(32/pi - a1, 28, RN)); a2; |
30 | // > a3 = D(round(32/pi - a1 - a2, 28, RN)); a3; |
31 | // > a4 = D(round(32/pi - a1 - a2 - a3, 28, RN)); a4; |
32 | // ... |
33 | static constexpr double THIRTYTWO_OVER_PI_28[N_ENTRIES] = { |
34 | 0x1.45f306ep+3, -0x1.b1bbeaep-28, 0x1.3f84ebp-57, -0x1.7056592p-87, |
35 | 0x1.c0db62ap-116, -0x1.4cd8778p-145, -0x1.bef806cp-174, 0x1.63abdecp-204}; |
36 | |
37 | // Exponents of the least significant bits of the corresponding entries in |
38 | // THIRTYTWO_OVER_PI_28. |
39 | static constexpr int THIRTYTWO_OVER_PI_28_LSB_EXP[N_ENTRIES] = { |
40 | -24, -55, -81, -114, -143, -170, -200, -230}; |
41 | |
42 | // Return k and y, where |
43 | // k = round(x * 16 / pi) and y = (x * 16 / pi) - k. |
44 | LIBC_INLINE int64_t small_range_reduction(double x, double &y) { |
45 | double prod = x * THIRTYTWO_OVER_PI_28[0]; |
46 | double kd = fputil::nearest_integer(prod); |
47 | y = prod - kd; |
48 | y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[1], y); |
49 | y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[2], y); |
50 | return static_cast<int64_t>(kd); |
51 | } |
52 | |
53 | // Return k and y, where |
54 | // k = round(x * 32 / pi) and y = (x * 32 / pi) - k. |
55 | // For large range, there are at most 2 parts of THIRTYTWO_OVER_PI_28 |
56 | // contributing to the lowest 6 binary digits (k & 63). If the least |
57 | // significant bit of x * the least significant bit of THIRTYTWO_OVER_PI_28[i] |
58 | // >= 64, we can completely ignore THIRTYTWO_OVER_PI_28[i]. |
59 | LIBC_INLINE int64_t large_range_reduction(double x, int x_exp, double &y) { |
60 | int idx = 0; |
61 | y = 0; |
62 | int x_lsb_exp_m4 = x_exp - fputil::FPBits<float>::FRACTION_LEN; |
63 | |
64 | // Skipping the first parts of 32/pi such that: |
65 | // LSB of x * LSB of THIRTYTWO_OVER_PI_28[i] >= 32. |
66 | while (x_lsb_exp_m4 + THIRTYTWO_OVER_PI_28_LSB_EXP[idx] > 5) |
67 | ++idx; |
68 | |
69 | double prod_hi = x * THIRTYTWO_OVER_PI_28[idx]; |
70 | // Get the integral part of x * THIRTYTWO_OVER_PI_28[idx] |
71 | double k_hi = fputil::nearest_integer(prod_hi); |
72 | // Get the fractional part of x * THIRTYTWO_OVER_PI_28[idx] |
73 | double frac = prod_hi - k_hi; |
74 | double prod_lo = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[idx + 1], frac); |
75 | double k_lo = fputil::nearest_integer(prod_lo); |
76 | |
77 | // Now y is the fractional parts. |
78 | y = prod_lo - k_lo; |
79 | y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[idx + 2], y); |
80 | y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[idx + 3], y); |
81 | |
82 | return static_cast<int64_t>(k_hi) + static_cast<int64_t>(k_lo); |
83 | } |
84 | |
85 | } // namespace generic |
86 | |
87 | } // namespace LIBC_NAMESPACE |
88 | |
89 | #endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_H |
90 |
Warning: This file is not a C or C++ file. It does not have highlighting.