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/common.h"
16
17namespace LIBC_NAMESPACE {
18
19namespace fma {
20
21static 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);
28static 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.
34LIBC_INLINE int64_t small_range_reduction(double x, double &y) {
35 double kd = fputil::nearest_integer(x: x * THIRTYTWO_OVER_PI[0]);
36 y = fputil::fma(x, y: THIRTYTWO_OVER_PI[0], z: -kd);
37 y = fputil::fma(x, y: THIRTYTWO_OVER_PI[1], z: 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.
44LIBC_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(x: prod_hi.get_val());
57 double truncated_prod = fputil::fma(x, y: THIRTYTWO_OVER_PI[0], z: -k_hi);
58 double prod_lo = fputil::fma(x, y: THIRTYTWO_OVER_PI[1], z: truncated_prod);
59 double k_lo = fputil::nearest_integer(x: prod_lo);
60 y = fputil::fma(x, y: THIRTYTWO_OVER_PI[1], z: truncated_prod - k_lo);
61 y = fputil::fma(x, y: THIRTYTWO_OVER_PI[2], z: y);
62 y = fputil::fma(x, y: THIRTYTWO_OVER_PI[3], z: y);
63
64 return static_cast<int64_t>(k_lo);
65 }
66
67 // - When x >= 2^110, the full exact product of x * THIRTYTWO_OVER_PI[0] does
68 // not contain any of the lowest 6 unit bits, so we can ignore it completely.
69 // - When 2^99 <= |x| < 2^110, the lowest 6 unit bits are contained
70 // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[1]).
71 // - When |x| >= 2^110, the LSB of double(x * THIRTYTWO_OVER_PI[1]) is at
72 // least 64.
73 fputil::FPBits<double> prod_hi(x * THIRTYTWO_OVER_PI[1]);
74 prod_hi.set_uintval(prod_hi.uintval() &
75 ((x_exp < 110) ? (~0xfffULL) : (~0ULL))); // |x| < 2^110
76 double k_hi = fputil::nearest_integer(x: prod_hi.get_val());
77 double truncated_prod = fputil::fma(x, y: THIRTYTWO_OVER_PI[1], z: -k_hi);
78 double prod_lo = fputil::fma(x, y: THIRTYTWO_OVER_PI[2], z: truncated_prod);
79 double k_lo = fputil::nearest_integer(x: prod_lo);
80 y = fputil::fma(x, y: THIRTYTWO_OVER_PI[2], z: truncated_prod - k_lo);
81 y = fputil::fma(x, y: THIRTYTWO_OVER_PI[3], z: y);
82 y = fputil::fma(x, y: THIRTYTWO_OVER_PI[4], z: y);
83
84 return static_cast<int64_t>(k_lo);
85}
86
87} // namespace fma
88
89} // namespace LIBC_NAMESPACE
90
91#endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H
92

source code of libc/src/math/generic/range_reduction_fma.h