llvm/libc/src/math/generic/log10.cpp

//===-- Double-precision log10(x) function --------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//

#include "src/math/log10.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/double_double.h"
#include "src/__support/FPUtil/dyadic_float.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/common.h"
#include "src/__support/integer_literals.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY

#include "common_constants.h"
#include "log_range_reduction.h"

namespace LIBC_NAMESPACE_DECL {

// 128-bit precision dyadic floating point numbers.
Float128;

operator""_u128;

namespace {

constexpr fputil::DoubleDouble LOG10_E = …;

// A simple upper bound for the error of e_x * log(2) - log(r).
constexpr double HI_ERR = …;

// Extra errors from P is from using x^2 to reduce evaluation latency.
constexpr double P_ERR = …;

// log10(2) with 128-bit precision generated by SageMath with:
// def format_hex(value):
//     l = hex(value)[2:]
//     n = 8
//     x = [l[i:i + n] for i in range(0, len(l), n)]
//     return "0x" + "'".join(x) + "_u128"
// (s, m, e) = RealField(128)(2).log10().sign_exponent_mantissa();
// print(format_hex(m));
constexpr Float128 LOG10_2(Sign::POS, /*exponent=*/-129, /*mantissa=*/
                           0x9a209a84'fbcff798'8f8959ac'0b7c9178_u128);

alignas(64) constexpr LogRR LOG10_TABLE = …;

// > P = fpminimax(log10(1 + x)/x, 3, [|128...|], [-0x1.0002143p-29 , 0x1p-29]);
// > P;
// > dirtyinfnorm(log10(1 + x)/x - P, [-0x1.0002143p-29 , 0x1p-29]);
// 0x1.64fb8...p-123
constexpr Float128 BIG_COEFFS[4]{ … };

// Reuse the output of the fast pass range reduction.
// -2^-8 <= m_x < 2^-7
double log10_accurate(int e_x, int index, double m_x) { … }

} // namespace

LLVM_LIBC_FUNCTION(double, log10, (double x)) { … }

} // namespace LIBC_NAMESPACE_DECL