Skip to content

Commit

Permalink
[libc] Implement log10f correctly rounded for all rounding modes.
Browse files Browse the repository at this point in the history
Based on RLIBM implementation similar to logf and log2f.  Most of the exceptional inputs are the exact powers of 10.

Reviewed By: sivachandra, zimmermann6, santoshn, jpl169

Differential Revision: https://reviews.llvm.org/D118093
  • Loading branch information
lntue committed Jan 25, 2022
1 parent a2505bd commit e581841
Show file tree
Hide file tree
Showing 15 changed files with 415 additions and 0 deletions.
1 change: 1 addition & 0 deletions libc/config/linux/aarch64/entrypoints.txt
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
libc.src.math.log10f
libc.src.math.log2f
libc.src.math.logf
libc.src.math.logb
Expand Down
1 change: 1 addition & 0 deletions libc/config/linux/x86_64/entrypoints.txt
Original file line number Diff line number Diff line change
Expand Up @@ -144,6 +144,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.llround
libc.src.math.llroundf
libc.src.math.llroundl
libc.src.math.log10f
libc.src.math.log2f
libc.src.math.logf
libc.src.math.logb
Expand Down
2 changes: 2 additions & 0 deletions libc/spec/stdc.td
Original file line number Diff line number Diff line change
Expand Up @@ -396,6 +396,8 @@ def StdC : StandardSpec<"stdc"> {
FunctionSpec<"ldexpf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<IntType>]>,
FunctionSpec<"ldexpl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<IntType>]>,

FunctionSpec<"log10f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,

FunctionSpec<"log2f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,

FunctionSpec<"logf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
Expand Down
2 changes: 2 additions & 0 deletions libc/src/math/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,8 @@ add_math_entrypoint_object(ldexp)
add_math_entrypoint_object(ldexpf)
add_math_entrypoint_object(ldexpl)

add_math_entrypoint_object(log10f)

add_math_entrypoint_object(log2f)

add_math_entrypoint_object(logf)
Expand Down
14 changes: 14 additions & 0 deletions libc/src/math/generic/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -666,6 +666,20 @@ add_object_library(
-Wno-c++17-extensions
)

add_entrypoint_object(
log10f
SRCS
log10f.cpp
HDRS
../log10f.h
DEPENDS
.common_constants
libc.src.__support.FPUtil.fputil
COMPILE_OPTIONS
-O3
-Wno-c++17-extensions
)

add_entrypoint_object(
log2f
SRCS
Expand Down
182 changes: 182 additions & 0 deletions libc/src/math/generic/log10f.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,182 @@
//===-- Single-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/log10f.h"
#include "common_constants.h" // Lookup table for (1/f)
#include "src/__support/FPUtil/BasicOperations.h"
#include "src/__support/FPUtil/FEnvUtils.h"
#include "src/__support/FPUtil/FMA.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/common.h"

// This is an algorithm for log10(x) in single precision which is
// correctly rounded for all rounding modes, based on the implementation of
// log10(x) from the RLIBM project at:
// https://people.cs.rutgers.edu/~sn349/rlibm

// Step 1 - Range reduction:
// For x = 2^m * 1.mant, log(x) = m * log10(2) + log10(1.m)
// If x is denormal, we normalize it by multiplying x by 2^23 and subtracting
// m by 23.

// Step 2 - Another range reduction:
// To compute log(1.mant), let f be the highest 8 bits including the hidden
// bit, and d be the difference (1.mant - f), i.e. the remaining 16 bits of the
// mantissa. Then we have the following approximation formula:
// log10(1.mant) = log10(f) + log10(1.mant / f)
// = log10(f) + log10(1 + d/f)
// ~ log10(f) + P(d/f)
// since d/f is sufficiently small.
// log10(f) and 1/f are then stored in two 2^7 = 128 entries look-up tables.

// Step 3 - Polynomial approximation:
// To compute P(d/f), we use a single degree-5 polynomial in double precision
// which provides correct rounding for all but few exception values.
// For more detail about how this polynomial is obtained, please refer to the
// papers:
// Lim, J. and Nagarakatte, S., "One Polynomial Approximation to Produce
// Correctly Rounded Results of an Elementary Function for Multiple
// Representations and Rounding Modes", Proceedings of the 49th ACM SIGPLAN
// Symposium on Principles of Programming Languages (POPL-2022), Philadelphia,
// USA, Jan. 16-22, 2022.
// https://people.cs.rutgers.edu/~sn349/papers/rlibmall-popl-2022.pdf
// Aanjaneya, M., Lim, J., and Nagarakatte, S., "RLibm-Prog: Progressive
// Polynomial Approximations for Fast Correctly Rounded Math Libraries",
// Dept. of Comp. Sci., Rutgets U., Technical Report DCS-TR-758, Nov. 2021.
// https://arxiv.org/pdf/2111.12852.pdf.

namespace __llvm_libc {

// Exact power of 10 in float:

// Lookup table for log10(f) = log10(1 + n*2^(-7)) where n = 0..127.
static constexpr double LOG10_F[128] = {
0x0.0000000000000p+0, 0x1.bafd47221ed26p-9, 0x1.b9476a4fcd10fp-8,
0x1.49b0851443684p-7, 0x1.b5e908eb13790p-7, 0x1.10a83a8446c78p-6,
0x1.45f4f5acb8be0p-6, 0x1.7adc3df3b1ff8p-6, 0x1.af5f92b00e610p-6,
0x1.e3806acbd058fp-6, 0x1.0ba01a8170000p-5, 0x1.25502c0fc314cp-5,
0x1.3ed1199a5e425p-5, 0x1.58238eeb353dap-5, 0x1.71483427d2a99p-5,
0x1.8a3fadeb847f4p-5, 0x1.a30a9d609efeap-5, 0x1.bba9a058dfd84p-5,
0x1.d41d5164facb4p-5, 0x1.ec6647eb58808p-5, 0x1.02428c1f08016p-4,
0x1.0e3d29d81165ep-4, 0x1.1a23445501816p-4, 0x1.25f5215eb594ap-4,
0x1.31b3055c47118p-4, 0x1.3d5d335c53179p-4, 0x1.48f3ed1df48fbp-4,
0x1.5477731973e85p-4, 0x1.5fe80488af4fdp-4, 0x1.6b45df6f3e2c9p-4,
0x1.769140a2526fdp-4, 0x1.81ca63d05a44ap-4, 0x1.8cf183886480dp-4,
0x1.9806d9414a209p-4, 0x1.a30a9d609efeap-4, 0x1.adfd07416be07p-4,
0x1.b8de4d3ab3d98p-4, 0x1.c3aea4a5c6effp-4, 0x1.ce6e41e463da5p-4,
0x1.d91d5866aa99cp-4, 0x1.e3bc1ab0e19fep-4, 0x1.ee4aba610f204p-4,
0x1.f8c9683468191p-4, 0x1.019c2a064b486p-3, 0x1.06cbd67a6c3b6p-3,
0x1.0bf3d0937c41cp-3, 0x1.11142f0811357p-3, 0x1.162d082ac9d10p-3,
0x1.1b3e71ec94f7bp-3, 0x1.204881dee8777p-3, 0x1.254b4d35e7d3cp-3,
0x1.2a46e8ca7ba2ap-3, 0x1.2f3b691c5a001p-3, 0x1.3428e2540096dp-3,
0x1.390f6844a0b83p-3, 0x1.3def0e6dfdf85p-3, 0x1.42c7e7fe3fc02p-3,
0x1.479a07d3b6411p-3, 0x1.4c65807e93338p-3, 0x1.512a644296c3dp-3,
0x1.55e8c518b10f8p-3, 0x1.5aa0b4b0988fap-3, 0x1.5f52447255c92p-3,
0x1.63fd857fc49bbp-3, 0x1.68a288b60b7fcp-3, 0x1.6d415eaf0906bp-3,
0x1.71da17c2b7e80p-3, 0x1.766cc40889e85p-3, 0x1.7af97358b9e04p-3,
0x1.7f80354d952a0p-3, 0x1.84011944bcb75p-3, 0x1.887c2e605e119p-3,
0x1.8cf183886480dp-3, 0x1.9161276ba2978p-3, 0x1.95cb2880f45bap-3,
0x1.9a2f95085a45cp-3, 0x1.9e8e7b0c0d4bep-3, 0x1.a2e7e8618c2d2p-3,
0x1.a73beaaaa22f4p-3, 0x1.ab8a8f56677fcp-3, 0x1.afd3e3a23b680p-3,
0x1.b417f49ab8807p-3, 0x1.b856cf1ca3105p-3, 0x1.bc907fd5d1c40p-3,
0x1.c0c5134610e26p-3, 0x1.c4f495c0002a2p-3, 0x1.c91f1369eb7cap-3,
0x1.cd44983e9e7bdp-3, 0x1.d165300e333f7p-3, 0x1.d580e67edc43dp-3,
0x1.d997c70da9b47p-3, 0x1.dda9dd0f4a329p-3, 0x1.e1b733b0c7381p-3,
0x1.e5bfd5f83d342p-3, 0x1.e9c3cec58f807p-3, 0x1.edc328d3184afp-3,
0x1.f1bdeeb654901p-3, 0x1.f5b42ae08c407p-3, 0x1.f9a5e79f76ac5p-3,
0x1.fd932f1ddb4d6p-3, 0x1.00be05b217844p-2, 0x1.02b0432c96ff0p-2,
0x1.04a054e139004p-2, 0x1.068e3fa282e3dp-2, 0x1.087a0832fa7acp-2,
0x1.0a63b3456c819p-2, 0x1.0c4b457d3193dp-2, 0x1.0e30c36e71a7fp-2,
0x1.1014319e661bdp-2, 0x1.11f594839a5bdp-2, 0x1.13d4f0862b2e1p-2,
0x1.15b24a0004a92p-2, 0x1.178da53d1ee01p-2, 0x1.1967067bb94b8p-2,
0x1.1b3e71ec94f7bp-2, 0x1.1d13ebb32d7f9p-2, 0x1.1ee777e5f0dc3p-2,
0x1.20b91a8e76105p-2, 0x1.2288d7a9b2b64p-2, 0x1.2456b3282f786p-2,
0x1.2622b0ee3b79dp-2, 0x1.27ecd4d41eb67p-2, 0x1.29b522a64b609p-2,
0x1.2b7b9e258e422p-2, 0x1.2d404b073e27ep-2, 0x1.2f032cf56a5bep-2,
0x1.30c4478f0835fp-2, 0x1.32839e681fc62p-2};

INLINE_FMA
LLVM_LIBC_FUNCTION(float, log10f, (float x)) {
constexpr double LOG10_2 = 0x1.34413509f79ffp-2;

using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
double m = 0.0;

// Exact powers of 10 and other hard-to-round cases.
switch (xbits.uintval()) {
case 0x4120'0000U: // x = 10
return 1.0f;
case 0x42c8'0000U: // x = 100
return 2.0f;
case 0x447a'0000U: // x = 1,000
return 3.0f;
case 0x461c'4000U: // x = 10,000
return 4.0f;
case 0x47c3'5000U: // x = 100,000
return 5.0f;
case 0x4974'2400U: // x = 1,000,000
return 6.0f;
case 0x4b18'9680U: // x = 10,000,000
return 7.0f;
case 0x4cbe'bc20U: // x = 100,000,000
return 8.0f;
case 0x4e6e'6b28U: // x = 1,000,000,000
return 9.0f;
case 0x5015'02f9U: // x = 10,000,000,000
return 10.0f;
case 0x4f13'4f83U: // x = 2471461632.0
if (fputil::get_round() == FE_UPWARD)
return 0x1.2c9314p+3f;
break;
case 0x7956'ba5eU: { // x = 69683218960000541503257137270226944.0
int round_mode = fputil::get_round();
if (round_mode == FE_DOWNWARD || round_mode == FE_TOWARDZERO)
return 0x1.16bebap+5f;
break;
}
}

if (xbits.uintval() < FPBits::MIN_NORMAL ||
xbits.uintval() > FPBits::MAX_NORMAL) {
if (xbits.is_zero()) {
return static_cast<float>(FPBits::neg_inf());
}
if (xbits.get_sign() && !xbits.is_nan()) {
return FPBits::build_nan(1 << (fputil::MantissaWidth<float>::VALUE - 1));
}
if (xbits.is_inf_or_nan()) {
return x;
}
// Normalize denormal inputs.
xbits.val *= 0x1.0p23f;
m -= 23.0;
}

m += static_cast<double>(xbits.get_exponent());
// Set bits to 1.m
xbits.set_unbiased_exponent(0x7F);
int f_index = xbits.get_mantissa() >> 16;

FPBits f(xbits.val);
f.bits &= ~0x0000'FFFF;

double d = static_cast<float>(xbits) - static_cast<float>(f);
d *= ONE_OVER_F[f_index];

double extra_factor = fputil::fma(m, LOG10_2, LOG10_F[f_index]);

double r = fputil::polyeval(d, extra_factor, 0x1.bcb7b1526e4c5p-2,
-0x1.bcb7b1518a5e9p-3, 0x1.287a72a6f716p-3,
-0x1.bcadb40b85565p-4, 0x1.5e0bc97f97e22p-4);

return static_cast<float>(r);
}

} // namespace __llvm_libc
18 changes: 18 additions & 0 deletions libc/src/math/log10f.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
//===-- Implementation header for log10f ------------------------*- C++ -*-===//
//
// 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
//
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC_MATH_LOG10F_H
#define LLVM_LIBC_SRC_MATH_LOG10F_H

namespace __llvm_libc {

float log10f(float x);

} // namespace __llvm_libc

#endif // LLVM_LIBC_SRC_MATH_LOG10F_H
13 changes: 13 additions & 0 deletions libc/test/src/math/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -1211,6 +1211,19 @@ add_fp_unittest(
libc.src.__support.FPUtil.fputil
)

add_fp_unittest(
log10f_test
NEED_MPFR
SUITE
libc_math_unittests
SRCS
log10f_test.cpp
DEPENDS
libc.include.math
libc.src.math.log10f
libc.src.__support.FPUtil.fputil
)

add_subdirectory(generic)
add_subdirectory(exhaustive)
add_subdirectory(differential_testing)
11 changes: 11 additions & 0 deletions libc/test/src/math/differential_testing/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -233,6 +233,17 @@ add_diff_binary(
-fno-builtin
)

add_diff_binary(
log10f_perf
SRCS
log10f_perf.cpp
DEPENDS
.single_input_single_output_diff
libc.src.math.log10f
COMPILE_OPTIONS
-fno-builtin
)

add_diff_binary(
log2f_diff
SRCS
Expand Down
16 changes: 16 additions & 0 deletions libc/test/src/math/differential_testing/log10f_perf.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
//===-- Differential test for log10f --------------------------------------===//
//
// 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 "SingleInputSingleOutputDiff.h"

#include "src/math/log10f.h"

#include <math.h>

SINGLE_INPUT_SINGLE_OUTPUT_PERF(float, __llvm_libc::log10f, ::log10f,
"log10f_perf.log")
17 changes: 17 additions & 0 deletions libc/test/src/math/exhaustive/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,23 @@ add_fp_unittest(
-lpthread
)

add_fp_unittest(
log10f_test
NO_RUN_POSTBUILD
NEED_MPFR
SUITE
libc_math_exhaustive_tests
SRCS
log10f_test.cpp
DEPENDS
.exhaustive_test
libc.include.math
libc.src.math.log10f
libc.src.__support.FPUtil.fputil
LINK_OPTIONS
-lpthread
)

add_fp_unittest(
log2f_test
NO_RUN_POSTBUILD
Expand Down
55 changes: 55 additions & 0 deletions libc/test/src/math/exhaustive/log10f_test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
//===-- Exhaustive test for log10f ----------------------------------------===//
//
// 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 "exhaustive_test.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/math/log10f.h"
#include "utils/MPFRWrapper/MPFRUtils.h"
#include "utils/UnitTest/FPMatcher.h"

using FPBits = __llvm_libc::fputil::FPBits<float>;

namespace mpfr = __llvm_libc::testing::mpfr;

struct LlvmLibcLog10fExhaustiveTest : public LlvmLibcExhaustiveTest<uint32_t> {
void check(uint32_t start, uint32_t stop,
mpfr::RoundingMode rounding) override {
mpfr::ForceRoundingMode r(rounding);
uint32_t bits = start;
do {
FPBits xbits(bits);
float x = float(xbits);
EXPECT_MPFR_MATCH(mpfr::Operation::Log10, x, __llvm_libc::log10f(x), 0.5,
rounding);
} while (bits++ < stop);
}
};

// Range: All non-negative;
static constexpr uint32_t START = 0x0000'0000U;
static constexpr uint32_t STOP = 0x7f80'0000U;
// Range: [1, 10];
// static constexpr uint32_t START = 0x3f80'0000U;
// static constexpr uint32_t STOP = 0x41c0'0000U;
static constexpr int NUM_THREADS = 16;

TEST_F(LlvmLibcLog10fExhaustiveTest, RoundNearestTieToEven) {
test_full_range(START, STOP, NUM_THREADS, mpfr::RoundingMode::Nearest);
}

TEST_F(LlvmLibcLog10fExhaustiveTest, RoundUp) {
test_full_range(START, STOP, NUM_THREADS, mpfr::RoundingMode::Upward);
}

TEST_F(LlvmLibcLog10fExhaustiveTest, RoundDown) {
test_full_range(START, STOP, NUM_THREADS, mpfr::RoundingMode::Downward);
}

TEST_F(LlvmLibcLog10fExhaustiveTest, RoundTowardZero) {
test_full_range(START, STOP, NUM_THREADS, mpfr::RoundingMode::TowardZero);
}

0 comments on commit e581841

Please sign in to comment.