93 lines · c
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-exception6//7//===----------------------------------------------------------------------===//8 9#ifndef LIBC_SRC___SUPPORT_MATH_RANGE_REDUCTION_FMA_H10#define LIBC_SRC___SUPPORT_MATH_RANGE_REDUCTION_FMA_H11 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 17namespace LIBC_NAMESPACE_DECL {18 19namespace fma {20 21static constexpr uint32_t FAST_PASS_BOUND = 0x5600'0000U; // 2^4522 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, where33// 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 * 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, where42// 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^9946 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 contained50 // 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 at52 // 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^5556 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] does69 // 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 contained71 // 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 at73 // 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^11077 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 fma89 90} // namespace LIBC_NAMESPACE_DECL91 92#endif // LIBC_SRC___SUPPORT_MATH_RANGE_REDUCTION_FMA_H93