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___SUPPORT_MATH_RANGE_REDUCTION_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_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#include "src/__support/macros/config.h"
17
18namespace LIBC_NAMESPACE_DECL {
19
20namespace math {
21
22namespace trigonometric_func_utils_internal {
23
24LIBC_INLINE_VAR constexpr uint32_t FAST_PASS_BOUND = 0x4a80'0000U; // 2^22
25
26LIBC_INLINE_VAR constexpr int N_ENTRIES = 8;
27
28// We choose to split bits of 32/pi into 28-bit precision pieces, so that the
29// product of x * THIRTYTWO_OVER_PI_28[i] is exact.
30// These are generated by Sollya with:
31// > a1 = D(round(32/pi, 28, RN)); a1;
32// > a2 = D(round(32/pi - a1, 28, RN)); a2;
33// > a3 = D(round(32/pi - a1 - a2, 28, RN)); a3;
34// > a4 = D(round(32/pi - a1 - a2 - a3, 28, RN)); a4;
35// ...
36LIBC_INLINE_VAR constexpr double THIRTYTWO_OVER_PI_28[N_ENTRIES] = {
37 0x1.45f306ep+3, -0x1.b1bbeaep-28, 0x1.3f84ebp-57, -0x1.7056592p-87,
38 0x1.c0db62ap-116, -0x1.4cd8778p-145, -0x1.bef806cp-174, 0x1.63abdecp-204};
39
40// Exponents of the least significant bits of the corresponding entries in
41// THIRTYTWO_OVER_PI_28.
42LIBC_INLINE_VAR constexpr int THIRTYTWO_OVER_PI_28_LSB_EXP[N_ENTRIES] = {
43 -24, -55, -81, -114, -143, -170, -200, -230};
44
45// Return k and y, where
46// k = round(x * 16 / pi) and y = (x * 16 / pi) - k.
47LIBC_INLINE int64_t small_range_reduction(double x, double &y) {
48 double prod = x * THIRTYTWO_OVER_PI_28[0];
49 double kd = fputil::nearest_integer(x: prod);
50 y = prod - kd;
51 y = fputil::multiply_add(x, y: THIRTYTWO_OVER_PI_28[1], z: y);
52 y = fputil::multiply_add(x, y: THIRTYTWO_OVER_PI_28[2], z: y);
53 return static_cast<int64_t>(kd);
54}
55
56// Return k and y, where
57// k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
58// For large range, there are at most 2 parts of THIRTYTWO_OVER_PI_28
59// contributing to the lowest 6 binary digits (k & 63). If the least
60// significant bit of x * the least significant bit of THIRTYTWO_OVER_PI_28[i]
61// >= 64, we can completely ignore THIRTYTWO_OVER_PI_28[i].
62LIBC_INLINE int64_t large_range_reduction(double x, int x_exp, double &y) {
63 int idx = 0;
64 y = 0;
65 int x_lsb_exp_m4 = x_exp - fputil::FPBits<float>::FRACTION_LEN;
66
67 // Skipping the first parts of 32/pi such that:
68 // LSB of x * LSB of THIRTYTWO_OVER_PI_28[i] >= 32.
69 while (x_lsb_exp_m4 + THIRTYTWO_OVER_PI_28_LSB_EXP[idx] > 5)
70 ++idx;
71
72 double prod_hi = x * THIRTYTWO_OVER_PI_28[idx];
73 // Get the integral part of x * THIRTYTWO_OVER_PI_28[idx]
74 double k_hi = fputil::nearest_integer(x: prod_hi);
75 // Get the fractional part of x * THIRTYTWO_OVER_PI_28[idx]
76 double frac = prod_hi - k_hi;
77 double prod_lo = fputil::multiply_add(x, y: THIRTYTWO_OVER_PI_28[idx + 1], z: frac);
78 double k_lo = fputil::nearest_integer(x: prod_lo);
79
80 // Now y is the fractional parts.
81 y = prod_lo - k_lo;
82 y = fputil::multiply_add(x, y: THIRTYTWO_OVER_PI_28[idx + 2], z: y);
83 y = fputil::multiply_add(x, y: THIRTYTWO_OVER_PI_28[idx + 3], z: y);
84
85 return static_cast<int64_t>(k_hi) + static_cast<int64_t>(k_lo);
86}
87
88} // namespace trigonometric_func_utils_internal
89
90} // namespace math
91
92} // namespace LIBC_NAMESPACE_DECL
93
94#endif // LLVM_LIBC_SRC___SUPPORT_MATH_RANGE_REDUCTION_H
95