1 //===-- Utilities for trigonometric functions with FMA ----------*- 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_FMA_H
10 #define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H
11
12 #include "src/__support/FPUtil/FMA.h"
13 #include "src/__support/FPUtil/FPBits.h"
14 #include "src/__support/FPUtil/nearest_integer.h"
15 #include "src/__support/macros/config.h"
16
17 namespace LIBC_NAMESPACE_DECL {
18
19 namespace fma {
20
21 static constexpr uint32_t FAST_PASS_BOUND = 0x5600'0000U; // 2^45
22
23 // Digits of 32/pi, generated by Sollya with:
24 // > a0 = D(32/pi);
25 // > a1 = D(32/pi - a0);
26 // > a2 = D(32/pi - a0 - a1);
27 // > a3 = D(32/pi - a0 - a1 - a2);
28 static constexpr double THIRTYTWO_OVER_PI[5] = {
29 0x1.45f306dc9c883p+3, -0x1.6b01ec5417056p-51, -0x1.6447e493ad4cep-105,
30 0x1.e21c820ff28b2p-159, -0x1.508510ea79237p-214};
31
32 // Return k and y, where
33 // k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
small_range_reduction(double x,double & y)34 LIBC_INLINE int64_t small_range_reduction(double x, double &y) {
35 double kd = fputil::nearest_integer(x * THIRTYTWO_OVER_PI[0]);
36 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[0], -kd);
37 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[1], y);
38 return static_cast<int64_t>(kd);
39 }
40
41 // Return k and y, where
42 // k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
43 // This is used for sinf, cosf, sincosf.
44 LIBC_INLINE int64_t large_range_reduction(double x, int x_exp, double &y) {
45 // 2^45 <= |x| < 2^99
46 if (x_exp < 99) {
47 // - When x < 2^99, the full exact product of x * THIRTYTWO_OVER_PI[0]
48 // contains at least one integral bit <= 2^5.
49 // - When 2^45 <= |x| < 2^55, the lowest 6 unit bits are contained
50 // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[0]).
51 // - When |x| >= 2^55, the LSB of double(x * THIRTYTWO_OVER_PI[0]) is at
52 // least 2^6.
53 fputil::FPBits<double> prod_hi(x * THIRTYTWO_OVER_PI[0]);
54 prod_hi.set_uintval(prod_hi.uintval() &
55 ((x_exp < 55) ? (~0xfffULL) : (~0ULL))); // |x| < 2^55
56 double k_hi = fputil::nearest_integer(prod_hi.get_val());
57 double truncated_prod = fputil::fma<double>(x, THIRTYTWO_OVER_PI[0], -k_hi);
58 double prod_lo =
59 fputil::fma<double>(x, THIRTYTWO_OVER_PI[1], truncated_prod);
60 double k_lo = fputil::nearest_integer(prod_lo);
61 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[1], truncated_prod - k_lo);
62 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[2], y);
63 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[3], y);
64
65 return static_cast<int64_t>(k_lo);
66 }
67
68 // - When x >= 2^110, the full exact product of x * THIRTYTWO_OVER_PI[0] does
69 // not contain any of the lowest 6 unit bits, so we can ignore it completely.
70 // - When 2^99 <= |x| < 2^110, the lowest 6 unit bits are contained
71 // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[1]).
72 // - When |x| >= 2^110, the LSB of double(x * THIRTYTWO_OVER_PI[1]) is at
73 // least 64.
74 fputil::FPBits<double> prod_hi(x * THIRTYTWO_OVER_PI[1]);
75 prod_hi.set_uintval(prod_hi.uintval() &
76 ((x_exp < 110) ? (~0xfffULL) : (~0ULL))); // |x| < 2^110
77 double k_hi = fputil::nearest_integer(prod_hi.get_val());
78 double truncated_prod = fputil::fma<double>(x, THIRTYTWO_OVER_PI[1], -k_hi);
79 double prod_lo = fputil::fma<double>(x, THIRTYTWO_OVER_PI[2], truncated_prod);
80 double k_lo = fputil::nearest_integer(prod_lo);
81 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[2], truncated_prod - k_lo);
82 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[3], y);
83 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[4], y);
84
85 return static_cast<int64_t>(k_lo);
86 }
87
88 } // namespace fma
89
90 } // namespace LIBC_NAMESPACE_DECL
91
92 #endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H
93