[libc] Implement correctly rounded logf based on RLIBM library.

[lntue]

Summary

Implement correctly rounded logf based on RLIBM library: https://people.cs.rutgers.edu/~sn349/rlibm/.

[sivachandra]

Just a few cosmetic comments for now. Needs a rebase though to absorb a lot style changes.

libc/src/math/generic/logf.cpp

38

s/check out/refer to

39

s/ribm/rlibm

Also, seems like this is a yet to be published paper. So, we should update the link when it gets published.

47

For this and others, when they are already in an internal namespace which has a prefix of __, we shouldn't use __ prefix for globals. Also, can this be a constexpr?

[michaelrj]

libc/src/__support/FPUtil/x86_64/PolyEvalDouble.h

2

nit: formatting

libc/src/__support/FPUtil/x86_64/PolyEvalFloat.h

2

nit: formatting

[tschuett]

Do you believe this is a permalink? Maybe cite a paper of the algorithm?

[lntue]

Updating D115408: [libc] Implement correctly rounded logf based on RLIBM library.

[lntue]

[libc] Implement correctly rounded logf based on RLIBM library.

[lntue]

In D115408#3183323, @tschuett wrote:

Do you believe this is a permalink? Maybe cite a paper of the algorithm?

Thanks, I cited the full title of the paper instead.

libc/src/__support/FPUtil/x86_64/PolyEvalDouble.h

2

This file was moved to https://reviews.llvm.org/D115347

libc/src/__support/FPUtil/x86_64/PolyEvalFloat.h

2

This file was moved to https://reviews.llvm.org/D115347

libc/src/math/generic/logf.cpp

39

I added the full reference of the paper.

[lntue]

[libc] Implement correctly rounded logf based on RLIBM library.

[zimmermann6]

Do you know the cost in latency/throughput of the switch() with the 21 exceptional cases? Another way would be to perform a rounding test once the approximation of log has been computed, and go in the switch only if the test fails (which would happen very rarely).

[zimmermann6]

I confirm this version is correctly rounded for all binary32 inputs and rounding to nearest, by exhaustive testing. For other rounding modes I find 46 incorrect roundings for rounding towards zero, 46 for rounding towards +Inf, and 44 for rounding towards -Inf. I guess part of them are due to the hard-coded values for the 21 exceptional cases, which are on the wrong side with probability 1/2 each. Thus with little additional effort you could get a correctly rounded function for all rounding modes.

[santoshn]

One reason there are 21 exceptional inputs is because this patch uses a different polynomial evaluation than Horner's method that RLIBM uses. Let me generate a new polynomial using the polynomial evaluation in the patch. It will have significantly fewer exceptional inputs for a 5th-degree polynomial for log.

Addressing Paul Zimmerman's comment about other rounding modes, if you use the RLIBM polynomial generated with the round-to-odd for a 34-bit FP, it will produce correct results for all rounding modes. Isn't that a better solution than using the round-to-nearest polynomial and then tweaking it with exceptional inputs?

[santoshn]

If the goal is to generate 5th degree polynomial for log for the round-to-nearest with the specified polynomial evaluation in the patch.

Here is the best polynomial that I could generate as of now with 3 special case inputs;

Polynomial: y=0.0000000000000000000000000000000000000000000000000000000000000000000000e+00 x^(0) + 1.0000000000073561157165613622055388987064361572265625000000000000000000e+00 x^(1) + -5.0000000925149434838345996467978693544864654541015625000000000000000000e-01 x^(2) + 3.3333714042884438066849384085799101740121841430664062500000000000000000e-01 x^(3) + -2.5063110611886163514583358846721239387989044189453125000000000000000000e-01 x^(4) + 2.3503273830199439276000816789746750146150588989257812500000000000000000e-01 x^(5)

Special case inputs:

1.7966015845941743847569149750142969423905014991760253906250000000000000e-03
2.6190722430193863652647667805695164133794605731964111328125000000000000e-03
3.7279721707274009363797251381811292958445847034454345703125000000000000e-03

The constant term is 0.

[lntue]

Use a more accurate polynomial provided by @santoshn.

[lntue]

I've updated the patch using the polynomial provided by @santoshn . I've tried to benchmark various implementations, and it seems like the most efficient way to compute the quintic polynomial with zero free coeff that @santoshn provided is as follow:

• First compute it as a quartic polynomial r = a_5 x^4 + a4 x^3 + ... + a1 = polyeval(x, a1, a2, a3, a4, a5) = fma(polyeval(x, a2, a3, a4, a5), x, a1) // cubic polynomial is optimized in src/math/__support/FPUtil/x86_64/PolyEval.h
• Then use another fma to multiply r by x, and the add the extra part (m * log(2) + log(f)): result = fma(r, x, fma(m, log(2), log(f)))

With this computation, the exhaustive test found out that the number of exception cases reduced to 5: // Input decimal: 1.00639712810516357421875000000000000000000000000000

//   MPFR result: 0.00637675332837318907527892795527933097341434970404
// Input decimal: 9.47263622283935546875000000000000000000000000000000
//   MPFR result: 2.24840724468231192940413656173141803167922920169308
// Input decimal: 58037908.00000000000000000000000000000000000000000000000000
//   MPFR result: 17.87660694122314468798277511771537045730722298001695
// Input decimal: 127837836949849943048192.00000000000000000000000000000000000000000000000000
//   MPFR result: 53.20504951477050802794144423622214960099434498468168
// Input decimal: 54983060754563292101907316736.00000000000000000000000000000000000000000000000000
//   MPFR result: 66.17682266235352211194320433287158374936494916276426

The last 2 large exceptional value might come from evaluating m*log(2) + log(f) with fma.

@zimmermann6 : About latency: The early switch cases for exceptional values did slow down the whole function by 20-25%. I was playing around with it and found out that by moving the switch for exceptional cases close to the final fma call, the extra latency due to the switch is reduce to only ~5% of the whole function, which I think is acceptable.
About correct rounding for all the rounding modes, even though it would be a few extra switches, I'll need to update the testing infrastructure, especially the support for MPFR with all rounding modes, and the exhaustive tests for that. Also there is another design question needs to be considered is to support all rounding modes at compile-time and run-time. Since that works quite significant by itself, I would prefer to add the support for other rounding modes in a different patch.

[zimmermann6]

maybe I did something wrong, but with the latest version I get two failures for x=0x1.2f1fd6p+3 and x=0x1.bacb4ap+25.
If I disable the test for exceptional values I get five failures, those two and x=0x1.01a33ep+0, x=0x1.b121a6p+76 and 0x1.6351d8p+95.

With the following minimax polynomial found using Sollya I get only four failures (the value x=0x1.01a33ep+0 is CR):
0.19620128080245119450708557451434899121522903442383 * x^5 + -0.24996749678377125358785804110084427520632743835449 * x^4 + 0.33333320666203714033315463893814012408256530761719 * x^3 + -0.49999999978385922805301788685028441250324249267578 * x^2 + 0.9999999999998795408018281705153640359640121459961 * x

[lntue]

Fix double-rounding errors for exceptional cases.

[lntue]

In D115408#3188352, @zimmermann6 wrote:

maybe I did something wrong, but with the latest version I get two failures for x=0x1.2f1fd6p+3 and x=0x1.bacb4ap+25.
If I disable the test for exceptional values I get five failures, those two and x=0x1.01a33ep+0, x=0x1.b121a6p+76 and 0x1.6351d8p+95.

With the following minimax polynomial found using Sollya I get only four failures (the value x=0x1.01a33ep+0 is CR):
0.19620128080245119450708557451434899121522903442383 * x^5 + -0.24996749678377125358785804110084427520632743835449 * x^4 + 0.33333320666203714033315463893814012408256530761719 * x^3 + -0.49999999978385922805301788685028441250324249267578 * x^2 + 0.9999999999998795408018281705153640359640121459961 * x

Something is not right, since all of the failures that you found were already in the list of exceptional values:
(hexfloat, decimal, hex bits in float)
(0x0x1.01a33ep+0, 1.00639712810516357421875, 0x3f80d19f)
(0x1.2f1fd6p+3, 9.47263622283935546875, 0x41178feb)
(0x1.bacb4ap+25, 58037908.0, 0x4c5d65a5)
(0x1.b121a6p+76, 127837836949849943048192.0, 0x65d890d3)
(0x1.6351d8p+95, 54983060754563292101907316736.0, 0x6f31a8ec)

Nonetheless, after clearing the build cache and re-run the tests, I found 2 exceptional cases' results were off due to double-conversion when I tried to print their hexadecimal floats (inadvertently using intermediate doubles).
You might try clear the build cache and rerun the tests again.

[zimmermann6]

I confirm the new version is CR for all cases in rounding to nearest. A way to make the exceptional cases CR for directed rounding modes is the following:

case 0x3f80d19f:
  {
    volatile double h = 0x1.a1e82cp-8f, l = -0x1.fffd1ep-33f;
    return h + l;
  }

Or compile with -frounding-math to avoid the "volatile" trick.

[zimmermann6]

this revision is fine to me (for rounding to nearest), thanks!

[cqlauter]

libc/src/math/generic/logf.cpp

22

I think that this addition will burn a lot of accuracy in the case when the exponent m = -1 and 1.mant is close to 2. This is catastrophic cancellation.

There is a trick to avoid the cancellation, typically by bringing 1.mant into the range 0.75 <= 1.mant <= 1.5, which can be done by adding in a 1 into the IEEE754 representation at the right place.

[cqlauter]

libc/src/math/generic/logf.cpp

174

Overall, I think that there should be way more comments for math functions like this one. A line like the one I am highlighting is so hard to understand, even when you are working on that code. Imagine you open up the code in 5 years and need to reunderstand a line like this.

libc/config/linux/aarch64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.lrintf
     libc.src.math.lrintl
     libc.src.math.lround
     libc.src.math.lroundf
     libc.src.math.lroundl
     libc.src.math.ldexp
     libc.src.math.ldexpf
     libc.src.math.ldexpl
+    libc.src.math.logf
     libc.src.math.logb
     libc.src.math.logbf
     libc.src.math.logbl
     libc.src.math.modf
     libc.src.math.modff
     libc.src.math.modfl
     libc.src.math.nearbyint
     libc.src.math.nearbyintf

libc/config/linux/x86_64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.ldexpf
     libc.src.math.ldexpl
     libc.src.math.llrint
     libc.src.math.llrintf
     libc.src.math.llrintl
     libc.src.math.llround
     libc.src.math.llroundf
     libc.src.math.llroundl
+    libc.src.math.logf
     libc.src.math.logb
     libc.src.math.logbf
     libc.src.math.logbl
     libc.src.math.lrint
     libc.src.math.lrintf
     libc.src.math.lrintl
     libc.src.math.lround
     libc.src.math.lroundf

libc/config/windows/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.ilogbf
     libc.src.math.ilogbl
     libc.src.math.ldexp
     libc.src.math.ldexpf
     libc.src.math.ldexpl
     libc.src.math.llround
     libc.src.math.llroundf
     libc.src.math.llroundl
+    libc.src.math.logf
     libc.src.math.logb
     libc.src.math.logbf
     libc.src.math.logbl
     libc.src.math.lround
     libc.src.math.lroundf
     libc.src.math.lroundl
     libc.src.math.modf
     libc.src.math.modff

libc/spec/stdc.td

HeaderSpec Math = HeaderSpec<
           FunctionSpec<"ilogb", RetValSpec<IntType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"ilogbf", RetValSpec<IntType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"ilogbl", RetValSpec<IntType>, [ArgSpec<LongDoubleType>]>,
 
           FunctionSpec<"ldexp", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<IntType>]>,
           FunctionSpec<"ldexpf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<IntType>]>,
           FunctionSpec<"ldexpl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<IntType>]>,
 
+          FunctionSpec<"logf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+
           FunctionSpec<"logb", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"logbf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"logbl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>]>,
 
           FunctionSpec<"modf", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoublePtr>]>,
           FunctionSpec<"modff", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatPtr>]>,
           FunctionSpec<"modfl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoublePtr>]>,
 

libc/src/__support/FPUtil/PolyEval.h

 
 template <typename T> static inline T polyeval(T x, T a0) { return a0; }
 
 template <typename T, typename... Ts>
 static inline T polyeval(T x, T a0, Ts... a) {
   return fma(x, polyeval(x, a...), a0);
 }
 
+#ifdef LLVM_LIBC_ARCH_X86_64
+#include "immintrin.h"
+
+#include "x86_64/PolyEvalDouble.h"
+#include "x86_64/PolyEvalFloat.h"
+
+#endif // LLVM_LIBC_ARCH_X86_64
+
 } // namespace fputil
 } // namespace __llvm_libc
 
 #else
 
 namespace __llvm_libc {
 namespace fputil {
 

libc/src/__support/FPUtil/x86_64/FMA.h

 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_FMA_H
 
 #include "src/__support/architectures.h"
 
 #if !defined(LLVM_LIBC_ARCH_X86)
 #error "Invalid include"
 #endif
 
+#include "immintrin.h"
 #include "src/__support/CPP/TypeTraits.h"
 
 namespace __llvm_libc {
 namespace fputil {
 
 template <typename T>
 static inline cpp::EnableIfType<cpp::IsSame<T, float>::Value, T> fma(T x, T y,
                                                                      T z) {
-  float result = x;
-  __asm__ __volatile__("vfmadd213ss %x2, %x1, %x0"
-                       : "+x"(result)
-                       : "x"(y), "x"(z));
+  __m128 xmm = _mm_load_ss(&x);
+  __m128 ymm = _mm_load_ss(&y);
+  __m128 zmm = _mm_load_ss(&z);
+  __m128 r = _mm_fmadd_ss(xmm, ymm, zmm);
+  float result;
+  _mm_store_ss(&result, r);
   return result;
 }
 
 template <typename T>
 static inline cpp::EnableIfType<cpp::IsSame<T, double>::Value, T> fma(T x, T y,
                                                                       T z) {
-  double result = x;
-  __asm__ __volatile__("vfmadd213sd %x2, %x1, %x0"
-                       : "+x"(result)
-                       : "x"(y), "x"(z));
+  __m128d xmm = _mm_load_sd(&x);
+  __m128d ymm = _mm_load_sd(&y);
+  __m128d zmm = _mm_load_sd(&z);
+  __m128d r = _mm_fmadd_sd(xmm, ymm, zmm);
+  double result;
+  _mm_store_sd(&result, r);
   return result;
 }
 
 } // namespace fputil
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_FMA_H

libc/src/__support/FPUtil/x86_64/PolyEvalDouble.h

+//===-- Common header for PolyEval implementations for x86_64 double -- 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_SUPPORT_FPUTIL_X86_64_POLYEVALDOUBLE_H
+#define LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_POLYEVALDOUBLE_H
+
+template <>
+inline double polyeval(double x, double a0, double a1, double a2, double a3) {
+  __m256d xmm = _mm256_set1_pd(x);
+  __m256d a13 = _mm256_set_pd(0.0, x, a3, a1);
+  __m256d a02 = _mm256_set_pd(0.0, 0.0, a2, a0);
+  // r = (0, x^2, a3*x + a2, a1*x + a0)
+  __m256d r = _mm256_fmadd_pd(a13, xmm, a02);
+  return fma(r[2], r[1], r[0]);
+}
+
+template <>
+inline double polyeval(double x, double a0, double a1, double a2, double a3,
+                       double a4, double a5) {
+  __m256d xmm = _mm256_set1_pd(x);
+  __m256d a25 = _mm256_set_pd(0.0, x, a5, a2);
+  __m256d a14 = _mm256_set_pd(0.0, 0.0, a4, a1);
+  __m256d a03 = _mm256_set_pd(0.0, 0.0, a3, a0);
+  // r1 = (0, x^2, a5*x + a4, a2*x + a1)
+  __m256d r1 = _mm256_fmadd_pd(a25, xmm, a14);
+  // r2 = (0, x^3, (a5*x + a4)*x + a3, (a2*x + a1)*x + a0
+  __m256d r2 = _mm256_fmadd_pd(r1, xmm, a03);
+  return fma(r2[2], r2[1], r2[0]);
+}
+
+#endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_POLYEVALDOUBLE_H

libc/src/__support/FPUtil/x86_64/PolyEvalFloat.h

+//===-- Common header for PolyEval implementations for x86_64 float --- 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_SUPPORT_FPUTIL_X86_64_POLYEVALFLOAT_H
+#define LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_POLYEVALFLOAT_H
+
+template <>
+inline float polyeval(float x, float a0, float a1, float a2, float a3) {
+  __m128 xmm = _mm_set1_ps(x);
+  __m128 a13 = _mm_set_ps(0.0f, x, a3, a1);
+  __m128 a02 = _mm_set_ps(0.0f, 0.0f, a2, a0);
+  // r = (0, x^2, a3*x + a2, a1*x + a0)
+  __m128 r = _mm_fmadd_ps(a13, xmm, a02);
+  return fma(r[2], r[1], r[0]);
+}
+
+template <>
+inline float polyeval(float x, float a0, float a1, float a2, float a3, float a4,
+                      float a5) {
+  __m128 xmm = _mm_set1_ps(x);
+  __m128 a25 = _mm_set_ps(0.0f, x, a5, a2);
+  __m128 a14 = _mm_set_ps(0.0f, 0.0f, a4, a1);
+  __m128 a03 = _mm_set_ps(0.0f, 0.0f, a3, a0);
+  // r1 = (0, x^2, a5*x + a4, a2*x + a1)
+  __m128 r1 = _mm_fmadd_ps(a25, xmm, a14);
+  // r2 = (0, x^3, (a5*x + a4)*x + a3, (a2*x + a1)*x + a0
+  __m128 r2 = _mm_fmadd_ps(r1, xmm, a03);
+  return fma(r2[2], r2[1], r2[0]);
+}
+
+#endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_POLYEVALFLOAT_H

libc/src/math/CMakeLists.txt

 add_subdirectory(generic)
 if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/${LIBC_TARGET_ARCHITECTURE})
   add_subdirectory(${LIBC_TARGET_ARCHITECTURE})
 endif()
 
 function(add_math_entrypoint_object name)
   # We prefer machine specific implementation if available. Hence we check
   # that first and return early if we are able to add an alias target for the
   # machine specific implementation.
   get_fq_target_name("${LIBC_TARGET_ARCHITECTURE}.${name}" fq_machine_specific_target_name)
   if(TARGET ${fq_machine_specific_target_name})
     add_entrypoint_object(
       ${name}
       ALIAS
       DEPENDS
         .${LIBC_TARGET_ARCHITECTURE}.${name}
+      COMPILE_OPTIONS
+        -mfma
     )
     return()
   endif()
 
   get_fq_target_name("generic.${name}" fq_generic_target_name)
   if(TARGET ${fq_generic_target_name})
     add_entrypoint_object(
       ${name}
       ALIAS
       DEPENDS
         .generic.${name}
+      COMPILE_OPTIONS
+        -mfma
     )
     return()
   endif()
 
   # Add a dummy entrypoint object for missing implementations. They will be skipped
   # anyway as there will be no entry for them in the target entrypoints list.
   add_entrypoint_object(
     ${name}
add_entrypoint_object(
   SRCS
     fmaf.cpp
   HDRS
     fmaf.h
   DEPENDS
     libc.src.__support.FPUtil.fputil
   COMPILE_OPTIONS
     -O2
+    -mfma
 )
 
 add_entrypoint_object(
   fma
   SRCS
     fma.cpp
   HDRS
     fma.h
   DEPENDS
     libc.src.__support.FPUtil.fputil
   COMPILE_OPTIONS
     -O2
+    -mfma
 )
 
 add_math_entrypoint_object(ceil)
 add_math_entrypoint_object(ceilf)
 add_math_entrypoint_object(ceill)
 
 add_math_entrypoint_object(copysign)
 add_math_entrypoint_object(copysignf)

 add_math_entrypoint_object(ilogb)
 add_math_entrypoint_object(ilogbf)
 add_math_entrypoint_object(ilogbl)
 
 add_math_entrypoint_object(ldexp)
 add_math_entrypoint_object(ldexpf)
 add_math_entrypoint_object(ldexpl)
 
+add_math_entrypoint_object(logf)
+
 add_math_entrypoint_object(logb)
 add_math_entrypoint_object(logbf)
 add_math_entrypoint_object(logbl)
 
 add_math_entrypoint_object(llrint)
 add_math_entrypoint_object(llrintf)
 add_math_entrypoint_object(llrintl)
 

libc/src/math/generic/CMakeLists.txt

add_entrypoint_object(
   SRCS
     expm1f.cpp
   HDRS
     ../expm1f.h
   DEPENDS
     libc.include.math
     libc.src.math.expf
     libc.src.math.fabsf
+  COMPILE_OPTIONS
+    -mfma
 )
 
 add_entrypoint_object(
   copysign
   SRCS
     copysign.cpp
   HDRS
     ../copysign.h
HDRS
     ../ldexpl.h
   DEPENDS
     libc.src.__support.FPUtil.fputil
   COMPILE_OPTIONS
     -O2
 )
 
 add_entrypoint_object(
+  logf
+  SRCS
+    logf.cpp
+  HDRS
+    ../logf.h
+  DEPENDS
+    libc.src.__support.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O3
+    -mfma
+)
+
+add_entrypoint_object(
   logb
   SRCS
     logb.cpp
   HDRS
     ../logb.h
   DEPENDS
     libc.src.__support.FPUtil.fputil
   COMPILE_OPTIONS

libc/src/math/generic/logf.cpp

+//===-- Single-precision log(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/logf.h"
+#include "src/__support/FPUtil/BasicOperations.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 a correctly-rounded algorithm for log(x) in single precision with
+// round-to-nearest, tie-to-even mode 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 * log(2) + log(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:
+//   log(1.mant) = log(f) + log(1.mant / f)
+//               = log(f) + log(1 + d/f)
+//               ~ log(f) + P(d/f)
+// since d/f is sufficiently small.
+//   log(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, check out the following paper:
+// https://people.cs.rutgers.edu/~sn349/papers/ribmall-popl-2022.pdf
+
+namespace __llvm_libc {
+
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wc++17-extensions"
+
+// Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127.
+static const double __LOG_F[128] = {
+    0x0.0000000000000p+0, 0x1.fe02a6b106788p-8, 0x1.fc0a8b0fc03e3p-7,
+    0x1.7b91b07d5b11ap-6, 0x1.f829b0e783300p-6, 0x1.39e87b9febd5fp-5,
+    0x1.77458f632dcfcp-5, 0x1.b42dd711971bep-5, 0x1.f0a30c01162a6p-5,
+    0x1.16536eea37ae0p-4, 0x1.341d7961bd1d0p-4, 0x1.51b073f06183fp-4,
+    0x1.6f0d28ae56b4bp-4, 0x1.8c345d6319b20p-4, 0x1.a926d3a4ad563p-4,
+    0x1.c5e548f5bc743p-4, 0x1.e27076e2af2e5p-4, 0x1.fec9131dbeabap-4,
+    0x1.0d77e7cd08e59p-3, 0x1.1b72ad52f67a0p-3, 0x1.29552f81ff523p-3,
+    0x1.371fc201e8f74p-3, 0x1.44d2b6ccb7d1ep-3, 0x1.526e5e3a1b437p-3,
+    0x1.5ff3070a793d3p-3, 0x1.6d60fe719d21cp-3, 0x1.7ab890210d909p-3,
+    0x1.87fa06520c910p-3, 0x1.9525a9cf456b4p-3, 0x1.a23bc1fe2b563p-3,
+    0x1.af3c94e80bff2p-3, 0x1.bc286742d8cd6p-3, 0x1.c8ff7c79a9a21p-3,
+    0x1.d5c216b4fbb91p-3, 0x1.e27076e2af2e5p-3, 0x1.ef0adcbdc5936p-3,
+    0x1.fb9186d5e3e2ap-3, 0x1.0402594b4d040p-2, 0x1.0a324e27390e3p-2,
+    0x1.1058bf9ae4ad5p-2, 0x1.1675cababa60ep-2, 0x1.1c898c16999fap-2,
+    0x1.22941fbcf7965p-2, 0x1.2895a13de86a3p-2, 0x1.2e8e2bae11d30p-2,
+    0x1.347dd9a987d54p-2, 0x1.3a64c556945e9p-2, 0x1.404308686a7e3p-2,
+    0x1.4618bc21c5ec2p-2, 0x1.4be5f957778a0p-2, 0x1.51aad872df82dp-2,
+    0x1.5767717455a6cp-2, 0x1.5d1bdbf5809cap-2, 0x1.62c82f2b9c795p-2,
+    0x1.686c81e9b14aep-2, 0x1.6e08eaa2ba1e3p-2, 0x1.739d7f6bbd006p-2,
+    0x1.792a55fdd47a2p-2, 0x1.7eaf83b82afc3p-2, 0x1.842d1da1e8b17p-2,
+    0x1.89a3386c1425ap-2, 0x1.8f11e873662c7p-2, 0x1.947941c2116fap-2,
+    0x1.99d958117e08ap-2, 0x1.9f323ecbf984bp-2, 0x1.a484090e5bb0ap-2,
+    0x1.a9cec9a9a0849p-2, 0x1.af1293247786bp-2, 0x1.b44f77bcc8f62p-2,
+    0x1.b9858969310fbp-2, 0x1.beb4d9da71b7bp-2, 0x1.c3dd7a7cdad4dp-2,
+    0x1.c8ff7c79a9a21p-2, 0x1.ce1af0b85f3ebp-2, 0x1.d32fe7e00ebd5p-2,
+    0x1.d83e7258a2f3ep-2, 0x1.dd46a04c1c4a0p-2, 0x1.e24881a7c6c26p-2,
+    0x1.e744261d68787p-2, 0x1.ec399d2468cc0p-2, 0x1.f128f5faf06ecp-2,
+    0x1.f6123fa7028acp-2, 0x1.faf588f78f31ep-2, 0x1.ffd2e0857f498p-2,
+    0x1.02552a5a5d0fep-1, 0x1.04bdf9da926d2p-1, 0x1.0723e5c1cdf40p-1,
+    0x1.0986f4f573520p-1, 0x1.0be72e4252a82p-1, 0x1.0e44985d1cc8bp-1,
+    0x1.109f39e2d4c96p-1, 0x1.12f719593efbcp-1, 0x1.154c3d2f4d5e9p-1,
+    0x1.179eabbd899a0p-1, 0x1.19ee6b467c96ep-1, 0x1.1c3b81f713c24p-1,
+    0x1.1e85f5e7040d0p-1, 0x1.20cdcd192ab6dp-1, 0x1.23130d7bebf42p-1,
+    0x1.2555bce98f7cbp-1, 0x1.2795e1289b11ap-1, 0x1.29d37fec2b08ap-1,
+    0x1.2c0e9ed448e8bp-1, 0x1.2e47436e40268p-1, 0x1.307d7334f10bep-1,
+    0x1.32b1339121d71p-1, 0x1.34e289d9ce1d3p-1, 0x1.37117b54747b5p-1,
+    0x1.393e0d3562a19p-1, 0x1.3b68449fffc22p-1, 0x1.3d9026a7156fap-1,
+    0x1.3fb5b84d16f42p-1, 0x1.41d8fe84672aep-1, 0x1.43f9fe2f9ce67p-1,
+    0x1.4618bc21c5ec2p-1, 0x1.48353d1ea88dfp-1, 0x1.4a4f85db03ebbp-1,
+    0x1.4c679afccee39p-1, 0x1.4e7d811b75bb0p-1, 0x1.50913cc01686bp-1,
+    0x1.52a2d265bc5aap-1, 0x1.54b2467999497p-1, 0x1.56bf9d5b3f399p-1,
+    0x1.58cadb5cd7989p-1, 0x1.5ad404c359f2cp-1, 0x1.5cdb1dc6c1764p-1,
+    0x1.5ee02a9241675p-1, 0x1.60e32f44788d8p-1};
+
+// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127.
+static const double __ONE_OVER_F[128] = {
+    0x1.0000000000000p+0, 0x1.fc07f01fc07f0p-1, 0x1.f81f81f81f820p-1,
+    0x1.f44659e4a4271p-1, 0x1.f07c1f07c1f08p-1, 0x1.ecc07b301ecc0p-1,
+    0x1.e9131abf0b767p-1, 0x1.e573ac901e574p-1, 0x1.e1e1e1e1e1e1ep-1,
+    0x1.de5d6e3f8868ap-1, 0x1.dae6076b981dbp-1, 0x1.d77b654b82c34p-1,
+    0x1.d41d41d41d41dp-1, 0x1.d0cb58f6ec074p-1, 0x1.cd85689039b0bp-1,
+    0x1.ca4b3055ee191p-1, 0x1.c71c71c71c71cp-1, 0x1.c3f8f01c3f8f0p-1,
+    0x1.c0e070381c0e0p-1, 0x1.bdd2b899406f7p-1, 0x1.bacf914c1bad0p-1,
+    0x1.b7d6c3dda338bp-1, 0x1.b4e81b4e81b4fp-1, 0x1.b2036406c80d9p-1,
+    0x1.af286bca1af28p-1, 0x1.ac5701ac5701bp-1, 0x1.a98ef606a63bep-1,
+    0x1.a6d01a6d01a6dp-1, 0x1.a41a41a41a41ap-1, 0x1.a16d3f97a4b02p-1,
+    0x1.9ec8e951033d9p-1, 0x1.9c2d14ee4a102p-1, 0x1.999999999999ap-1,
+    0x1.970e4f80cb872p-1, 0x1.948b0fcd6e9e0p-1, 0x1.920fb49d0e229p-1,
+    0x1.8f9c18f9c18fap-1, 0x1.8d3018d3018d3p-1, 0x1.8acb90f6bf3aap-1,
+    0x1.886e5f0abb04ap-1, 0x1.8618618618618p-1, 0x1.83c977ab2beddp-1,
+    0x1.8181818181818p-1, 0x1.7f405fd017f40p-1, 0x1.7d05f417d05f4p-1,
+    0x1.7ad2208e0ecc3p-1, 0x1.78a4c8178a4c8p-1, 0x1.767dce434a9b1p-1,
+    0x1.745d1745d1746p-1, 0x1.724287f46debcp-1, 0x1.702e05c0b8170p-1,
+    0x1.6e1f76b4337c7p-1, 0x1.6c16c16c16c17p-1, 0x1.6a13cd1537290p-1,
+    0x1.6816816816817p-1, 0x1.661ec6a5122f9p-1, 0x1.642c8590b2164p-1,
+    0x1.623fa77016240p-1, 0x1.6058160581606p-1, 0x1.5e75bb8d015e7p-1,
+    0x1.5c9882b931057p-1, 0x1.5ac056b015ac0p-1, 0x1.58ed2308158edp-1,
+    0x1.571ed3c506b3ap-1, 0x1.5555555555555p-1, 0x1.5390948f40febp-1,
+    0x1.51d07eae2f815p-1, 0x1.5015015015015p-1, 0x1.4e5e0a72f0539p-1,
+    0x1.4cab88725af6ep-1, 0x1.4afd6a052bf5bp-1, 0x1.49539e3b2d067p-1,
+    0x1.47ae147ae147bp-1, 0x1.460cbc7f5cf9ap-1, 0x1.446f86562d9fbp-1,
+    0x1.42d6625d51f87p-1, 0x1.4141414141414p-1, 0x1.3fb013fb013fbp-1,
+    0x1.3e22cbce4a902p-1, 0x1.3c995a47babe7p-1, 0x1.3b13b13b13b14p-1,
+    0x1.3991c2c187f63p-1, 0x1.3813813813814p-1, 0x1.3698df3de0748p-1,
+    0x1.3521cfb2b78c1p-1, 0x1.33ae45b57bcb2p-1, 0x1.323e34a2b10bfp-1,
+    0x1.30d190130d190p-1, 0x1.2f684bda12f68p-1, 0x1.2e025c04b8097p-1,
+    0x1.2c9fb4d812ca0p-1, 0x1.2b404ad012b40p-1, 0x1.29e4129e4129ep-1,
+    0x1.288b01288b013p-1, 0x1.27350b8812735p-1, 0x1.25e22708092f1p-1,
+    0x1.2492492492492p-1, 0x1.23456789abcdfp-1, 0x1.21fb78121fb78p-1,
+    0x1.20b470c67c0d9p-1, 0x1.1f7047dc11f70p-1, 0x1.1e2ef3b3fb874p-1,
+    0x1.1cf06ada2811dp-1, 0x1.1bb4a4046ed29p-1, 0x1.1a7b9611a7b96p-1,
+    0x1.19453808ca29cp-1, 0x1.1811811811812p-1, 0x1.16e0689427379p-1,
+    0x1.15b1e5f75270dp-1, 0x1.1485f0e0acd3bp-1, 0x1.135c81135c811p-1,
+    0x1.12358e75d3033p-1, 0x1.1111111111111p-1, 0x1.0fef010fef011p-1,
+    0x1.0ecf56be69c90p-1, 0x1.0db20a88f4696p-1, 0x1.0c9714fbcda3bp-1,
+    0x1.0b7e6ec259dc8p-1, 0x1.0a6810a6810a7p-1, 0x1.0953f39010954p-1,
+    0x1.0842108421084p-1, 0x1.073260a47f7c6p-1, 0x1.0624dd2f1a9fcp-1,
+    0x1.05197f7d73404p-1, 0x1.0410410410410p-1, 0x1.03091b51f5e1ap-1,
+    0x1.0204081020408p-1, 0x1.0101010101010p-1};
+
+LLVM_LIBC_FUNCTION(float, logf, (float x)) {
+  const double __LOG_2 = 0x1.62e42fefa39efp-1;
+  using FPBits = typename fputil::FPBits<float>;
+  FPBits xbits(x);
+  int m = 0;
+
+  switch (xbits.uintval()) {
+  // Exceptional inputs:
+  case 0x3509dcf6U:
+    return -0x1.cf6b66p+3f;
+  case 0x3bf86ef0U:
+    return -0x1.38733cp+2f;
+  case 0x3ca1c99fU:
+    return -0x1.f65a4p+1f;
+  case 0x3d13e105U:
+    return -0x1.a922b8p+1f;
+  case 0x3f7ff1f2U:
+    return -0x1.c1cc58p-13f;
+  case 0x3f7fffffU:
+    return -0x1p-24f;
+  case 0x3f800006U:
+    return 0x1.7ffff8p-21f;
+  case 0x3f800014U:
+    return 0x1.3fffe8p-19f;
+  case 0x3f80001cU:
+    return 0x1.bfffdp-19f;
+  case 0x3f80c777U:
+    return 0x1.8db86cp-8f;
+  case 0x3f80ce72U:
+    return 0x1.9b986cp-8f;
+  case 0x3f80f7bfU:
+    return 0x1.eda0e4p-8f;
+  case 0x3f80fcfeU:
+    return 0x1.f80a86p-8f;
+  case 0x3f81feb4U:
+    return 0x1.fac3a6p-7f;
+  case 0x3f83d731U:
+    return 0x1.e45d6ap-6f;
+  case 0x3f90cb1dU:
+    return 0x1.f8f20cp-4f;
+  case 0x3fc55379U:
+    return 0x1.bb371p-2f;
+  case 0x3fd364d7U:
+    return 0x1.00dde6p-1f;
+  case 0x4e85f412U:
+    return 0x1.4d7022p+4f;
+  case 0x65d890d3U:
+    return 0x1.a9a3f2p+5f;
+  case 0x6f31a8ecU:
+    return 0x1.08b512p+6f;
+  }
+
+  if (xbits.uintval() < FPBits::minNormal ||
+      xbits.uintval() > FPBits::maxNormal) {
+    if (xbits.isZero()) {
+      return static_cast<float>(FPBits::negInf());
+    }
+    if (xbits.getSign() && !xbits.isNaN()) {
+      return FPBits::buildNaN(1 << (fputil::MantissaWidth<float>::value - 1));
+    }
+    if (xbits.isInfOrNaN()) {
+      return x;
+    }
+    // Normalize denormal inputs.
+    xbits.val *= 0x1.0p23f;
+    m = -23;
+  }
+
+  m += xbits.getExponent();
+  // Set bits to 1.m
+  xbits.setUnbiasedExponent(0x7F);
+  int fIndex = xbits.getMantissa() >> 16;
+
+  FPBits f(xbits.val);
+  f.bits &= ~0x0000'FFFF;
+
+  double d = static_cast<float>(xbits) - static_cast<float>(f);
+  d *= __ONE_OVER_F[fIndex];
+
+  double r = __llvm_libc::fputil::polyeval(
+      d, 0x1p-150, 0x1.ffffffffda0f7p-1, -0x1.fffffe19e207cp-2,
+      0x1.55517786e8551p-2, -0x1.f9e4d58f2259p-3);
+
+  r += __llvm_libc::fputil::fma(static_cast<double>(m), __LOG_2,
+                                __LOG_F[fIndex]);
+  return r;
+}
+
+#pragma clang diagnostic pop
+
+} // namespace __llvm_libc

libc/src/math/logf.h

+//===-- Implementation header for logf --------------------------*- 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_LOGF_H
+#define LLVM_LIBC_SRC_MATH_LOGF_H
+
+namespace __llvm_libc {
+
+float logf(float x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_LOGF_H

libc/test/src/math/CMakeLists.txt

SUITE
     libc-long-running-tests
   SRCS
     mod_k_pi_test.cpp
   DEPENDS
     libc.src.math.generic.dp_trig
     libc.src.__support.FPUtil.fputil
 )
 
+add_fp_unittest(
+  logf_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    logf_test.cpp
+  DEPENDS
+    libc.include.math
+    libc.src.math.logf
+    libc.src.__support.FPUtil.fputil
+)
+
 add_subdirectory(generic)
 add_subdirectory(exhaustive)
 add_subdirectory(differential_testing)

libc/test/src/math/differential_testing/CMakeLists.txt

add_diff_binary(
   DEPENDS
     .single_input_single_output_diff
     libc.src.math.floorf
   COMPILE_OPTIONS
     -fno-builtin
 )
 
 add_diff_binary(
+  logf_diff
+  SRCS
+    logf_diff.cpp
+  DEPENDS
+    .single_input_single_output_diff
+    libc.src.math.logf
+  COMPILE_OPTIONS
+    -fno-builtin
+)
+
+add_diff_binary(
+  logf_perf
+  SRCS
+    logf_perf.cpp
+  DEPENDS
+    .single_input_single_output_diff
+    libc.src.math.logf
+  COMPILE_OPTIONS
+    -fno-builtin
+)
+
+add_diff_binary(
   logbf_diff
   SRCS
     logbf_diff.cpp
   DEPENDS
     .single_input_single_output_diff
     libc.src.math.logbf
 )
 

libc/test/src/math/differential_testing/logf_diff.cpp

+//===-- Differential test for logf ----------------------------------------===//
+//
+// 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/logf.h"
+
+#include <math.h>
+
+SINGLE_INPUT_SINGLE_OUTPUT_DIFF(float, __llvm_libc::logf, ::logf,
+                                "logf_diff.log")

libc/test/src/math/differential_testing/logf_perf.cpp

+//===-- Differential test for logf ----------------------------------------===//
+//
+// 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/logf.h"
+
+#include <math.h>
+
+SINGLE_INPUT_SINGLE_OUTPUT_PERF(float, __llvm_libc::logf, ::logf,
+                                "logf_perf.log")

libc/test/src/math/exhaustive/CMakeLists.txt

SUITE
     libc_math_exhaustive_tests
   SRCS
   expm1f_test.cpp
   DEPENDS
     libc.include.math
     libc.src.math.expm1f
     libc.src.__support.FPUtil.fputil
 )
+
+add_fp_unittest(
+  logf_test
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    logf_test.cpp
+  DEPENDS
+    libc.include.math
+    libc.src.math.logf
+    libc.src.__support.FPUtil.fputil
+)

libc/test/src/math/exhaustive/logf_test.cpp

+//===-- Exhaustive test for logf ------------------------------------------===//
+//
+// 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/__support/FPUtil/FPBits.h"
+#include "src/math/logf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+TEST(LlvmLibcLogfExhaustiveTest, AllValues) {
+  // failed inputs: 3509dcf6, 3bf86ef0, 3ca1c99f, 3d13e105, 3f7ff1f2
+  uint32_t bits = 0x3f7ff1f3;
+  do {
+    FPBits xbits(bits);
+    float x = float(xbits);
+    EXPECT_MPFR_MATCH(mpfr::Operation::Log, x, __llvm_libc::logf(x), 0.5);
+  } while (bits++ < 0x7f7f'ffffU);
+}

libc/test/src/math/expm1f_test.cpp

for (uint32_t i = 0, v = 0; i <= count; ++i, v += step) {
     float result = __llvm_libc::expm1f(x);
 
     // If the computation resulted in an error or did not produce valid result
     // in the single-precision floating point range, then ignore comparing with
     // MPFR result as MPFR can still produce valid results because of its
     // wider precision.
     if (isnan(result) || isinf(result) || errno != 0)
       continue;
-    ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 1.5);
+    ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 3);
   }
 }

libc/test/src/math/logf_test.cpp

+//===-- Unittests for logf-----------------------------------------------===//
+//
+// 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/__support/FPUtil/FPBits.h"
+#include "src/math/logf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(float)
+
+TEST(LlvmLibcLogfTest, SpecialNumbers) {
+  EXPECT_FP_EQ(aNaN, __llvm_libc::logf(aNaN));
+  EXPECT_FP_EQ(inf, __llvm_libc::logf(inf));
+  EXPECT_TRUE(isnan(__llvm_libc::logf(negInf)));
+  EXPECT_FP_EQ(negInf, __llvm_libc::logf(0.0f));
+  EXPECT_FP_EQ(negInf, __llvm_libc::logf(-0.0f));
+  EXPECT_TRUE(isnan(__llvm_libc::logf(-1.0f)));
+  EXPECT_FP_EQ(zero, __llvm_libc::logf(1.0f));
+}
+
+TEST(LlvmLibcLogfTest, TrickyInputs) {
+  constexpr int N = 21;
+  constexpr uint32_t INPUTS[N] = {
+      0x3509dcf6U, 0x3bf86ef0U, 0x3ca1c99fU, 0x3d13e105U, 0x3f7ff1f2U,
+      0x3f7fffffU, 0x3f800006U, 0x3f800014U, 0x3f80001cU, 0x3f80c777U,
+      0x3f80ce72U, 0x3f80f7bfU, 0x3f80fcfeU, 0x3f81feb4U, 0x3f83d731U,
+      0x3f90cb1dU, 0x3fc55379U, 0x3fd364d7U, 0x4e85f412U, 0x65d890d3U,
+      0x6f31a8ecU,
+  };
+  for (int i = 0; i < N; ++i) {
+    float x = float(FPBits(INPUTS[i]));
+    ASSERT_MPFR_MATCH(mpfr::Operation::Log, x, __llvm_libc::logf(x), 0.5);
+  }
+}
+
+TEST(LlvmLibcLogfTest, InFloatRange) {
+  constexpr uint32_t count = 1000000;
+  constexpr uint32_t step = UINT32_MAX / count;
+  for (uint32_t i = 0, v = 0; i <= count; ++i, v += step) {
+    float x = float(FPBits(v));
+    if (isnan(x) || isinf(x))
+      continue;
+    errno = 0;
+    float result = __llvm_libc::logf(x);
+    // If the computation resulted in an error or did not produce valid result
+    // in the single-precision floating point range, then ignore comparing with
+    // MPFR result as MPFR can still produce valid results because of its
+    // wider precision.
+    if (isnan(result) || isinf(result) || errno != 0)
+      continue;
+    ASSERT_MPFR_MATCH(mpfr::Operation::Log, x, __llvm_libc::logf(x), 0.5);
+  }
+}

libc/utils/MPFRWrapper/MPFRUtils.h

enum class Operation : int {
   BeginUnaryOperationsSingleOutput,
   Abs,
   Ceil,
   Cos,
   Exp,
   Exp2,
   Expm1,
   Floor,
+  Log,
   Mod2PI,
   ModPIOver2,
   ModPIOver4,
   Round,
   Sin,
   Sqrt,
   Tan,
   Trunc,

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
   }
 
   MPFRNumber hypot(const MPFRNumber &b) {
     MPFRNumber result;
     mpfr_hypot(result.value, value, b.value, MPFR_RNDN);
     return result;
   }
 
+  MPFRNumber log() const {
+    MPFRNumber result;
+    mpfr_log(result.value, value, MPFR_RNDN);
+    return result;
+  }
+
   MPFRNumber remquo(const MPFRNumber &divisor, int &quotient) {
     MPFRNumber remainder;
     long q;
     mpfr_remquo(remainder.value, &q, value, divisor.value, MPFR_RNDN);
     quotient = q;
     return remainder;
   }
 
unaryOperation(Operation op, InputType input) {
   case Operation::Exp:
     return mpfrInput.exp();
   case Operation::Exp2:
     return mpfrInput.exp2();
   case Operation::Expm1:
     return mpfrInput.expm1();
   case Operation::Floor:
     return mpfrInput.floor();
+  case Operation::Log:
+    return mpfrInput.log();
   case Operation::Mod2PI:
     return mpfrInput.mod_2pi();
   case Operation::ModPIOver2:
     return mpfrInput.mod_pi_over_2();
   case Operation::ModPIOver4:
     return mpfrInput.mod_pi_over_4();
   case Operation::Round:
     return mpfrInput.round();