95 lines · c
1//===-- Extra range reduction steps for accurate pass of logarithms -------===//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 LLVM_LIBC_SRC_MATH_GENERIC_LOG_RANGE_REDUCTION_H10#define LLVM_LIBC_SRC_MATH_GENERIC_LOG_RANGE_REDUCTION_H11 12#include "src/__support/FPUtil/dyadic_float.h"13#include "src/__support/macros/config.h"14#include "src/__support/math/common_constants.h"15#include "src/__support/uint128.h"16 17namespace LIBC_NAMESPACE_DECL {18 19// Struct to store -log*(r) for 4 range reduction steps.20struct LogRR {21 fputil::DyadicFloat<128> step_1[128];22 fputil::DyadicFloat<128> step_2[193];23 fputil::DyadicFloat<128> step_3[161];24 fputil::DyadicFloat<128> step_4[130];25};26 27// Perform logarithm range reduction steps 2-4.28// Inputs from the first step of range reduction:29// m_x : the reduced argument after the first step of range reduction30// satisfying -2^-8 <= m_x < 2^-7 and ulp(m_x) >= 2^-60.31// idx1: index of the -log(r1) table from the first step.32// Outputs of the extra range reduction steps:33// sum: adding -log(r1) - log(r2) - log(r3) - log(r4) to the resulted sum.34// return value: the reduced argument v satisfying:35// -0x1.0002143p-29 <= v < 0x1p-29, and ulp(v) >= 2^(-125).36LIBC_INLINE fputil::DyadicFloat<128>37log_range_reduction(double m_x, const LogRR &log_table,38 fputil::DyadicFloat<128> &sum) {39 using namespace common_constants_internal;40 using Float128 = typename fputil::DyadicFloat<128>;41 using MType = typename Float128::MantissaType;42 43 int64_t v = static_cast<int64_t>(m_x * 0x1.0p60); // ulp = 2^-6044 45 // Range reduction - Step 246 // Output range: vv2 in [-0x1.3ffcp-15, 0x1.3e3dp-15].47 // idx2 = trunc(2^14 * (v + 2^-8 + 2^-15))48 size_t idx2 = static_cast<size_t>((v + 0x10'2000'0000'0000) >> 46);49 sum = fputil::quick_add(sum, log_table.step_2[idx2]);50 51 int64_t s2 = static_cast<int64_t>(S2[idx2]); // |s| <= 2^-7, ulp = 2^-1652 int64_t sv2 = s2 * v; // |s*v| < 2^-14, ulp = 2^(-60-16) = 2^-7653 int64_t spv2 = (s2 << 44) + v; // |s + v| < 2^-14, ulp = 2^-6054 int64_t vv2 = (spv2 << 16) + sv2; // |vv2| < 2^-14, ulp = 2^-7655 56 // Range reduction - Step 357 // Output range: vv3 in [-0x1.01928p-22 , 0x1p-22]58 // idx3 = trunc(2^21 * (v + 80*2^-21 + 2^-22))59 size_t idx3 = static_cast<size_t>((vv2 + 0x2840'0000'0000'0000) >> 55);60 sum = fputil::quick_add(sum, log_table.step_3[idx3]);61 62 int64_t s3 = static_cast<int64_t>(S3[idx3]); // |s| < 2^-13, ulp = 2^-2163 int64_t spv3 = (s3 << 55) + vv2; // |s + v| < 2^-21, ulp = 2^-7664 // |s*v| < 2^-27, ulp = 2^(-76-21) = 2^-9765 Int128 sv3 = static_cast<Int128>(s3) * static_cast<Int128>(vv2);66 // |vv3| < 2^-21, ulp = 2^-9767 Int128 vv3 = (static_cast<Int128>(spv3) << 21) + sv3;68 69 // Range reduction - Step 470 // Output range: vv4 in [-0x1.0002143p-29 , 0x1p-29]71 // idx4 = trunc(2^21 * (v + 65*2^-28 + 2^-29))72 size_t idx4 = static_cast<size_t>((static_cast<int>(vv3 >> 68) + 131) >> 1);73 74 sum = fputil::quick_add(sum, log_table.step_4[idx4]);75 76 Int128 s4 = static_cast<Int128>(S4[idx4]); // |s| < 2^-21, ulp = 2^-2877 // |s + v| < 2^-28, ulp = 2^-9778 Int128 spv4 = (s4 << 69) + vv3;79 // |s*v| < 2^-42, ulp = 2^(-97-28) = 2^-12580 Int128 sv4 = s4 * vv3;81 // |vv4| < 2^-28, ulp = 2^-12582 Int128 vv4 = (spv4 << 28) + sv4;83 84 return (vv4 < 0) ? Float128(Sign::NEG, -125,85 MType({static_cast<uint64_t>(-vv4),86 static_cast<uint64_t>((-vv4) >> 64)}))87 : Float128(Sign::POS, -125,88 MType({static_cast<uint64_t>(vv4),89 static_cast<uint64_t>(vv4 >> 64)}));90}91 92} // namespace LIBC_NAMESPACE_DECL93 94#endif // LLVM_LIBC_SRC_MATH_GENERIC_LOG_RANGE_REDUCTION_H95