[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

39

s/check out/refer to

40

s/ribm/rlibm

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

48

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 ↗

(On Diff #392978)

nit: formatting

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

2 ↗

(On Diff #392978)

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 ↗

(On Diff #392978)

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

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

2 ↗

(On Diff #392978)

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

libc/src/math/generic/logf.cpp

40

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

21

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

173

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/math/CMakeLists.txt

 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

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
+)
+
+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, please refer to the
+// paper:
+//   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, January 16-22, 2022.
+// https://people.cs.rutgers.edu/~sn349/papers/rlibmall-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 constexpr 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 constexpr 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};
+
+INLINE_FMA
+LLVM_LIBC_FUNCTION(float, logf, (float x)) {
+  constexpr double LOG_2 = 0x1.62e42fefa39efp-1;
+  using FPBits = typename fputil::FPBits<float>;
+  FPBits xbits(x);
+  int m = 0;
+
+  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;
+  }
+
+  m += xbits.get_exponent();
+  // Set bits to 1.m
+  xbits.set_unbiased_exponent(0x7F);
+  int fIndex = 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[fIndex];
+
+  double r = __llvm_libc::fputil::polyeval(
+      d, 0x1.0000000008169p+0, -0x1.0000004f78405p-1, 0x1.555654d2bc769p-2,
+      -0x1.00a570d090322p-2, 0x1.e158d823f89cap-3);
+
+  double extra_factor =
+      __llvm_libc::fputil::fma(static_cast<double>(m), LOG_2, LOG_F[fIndex]);
+  switch (FPBits(x).uintval()) {
+  case 0x3f80d19f:
+    return 0x1.a1e82cp-8f;
+  case 0x41178feb:
+    return 0x1.1fcbcep+1f;
+  case 0x4c5d65a5:
+    return 0x1.1e0696p+4f;
+  case 0x65d890d3:
+    return 0x1.a9a3f2p+5f;
+  case 0x6f31a8ec:
+    return 0x1.08b512p+6f;
+  default:
+    return static_cast<float>(__llvm_libc::fputil::fma(d, r, extra_factor));
+  }
+}
+
+#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) {
+  uint32_t bits = 0U;
+  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/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(neg_inf)));
+  EXPECT_FP_EQ(neg_inf, __llvm_libc::logf(0.0f));
+  EXPECT_FP_EQ(neg_inf, __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 = 24;
+  constexpr uint32_t INPUTS[N] = {
+      0x3509dcf6U, 0x3bf86ef0U, 0x3ca1c99fU, 0x3d13e105U, 0x3f7ff1f2U,
+      0x3f7fffffU, 0x3f800006U, 0x3f800014U, 0x3f80001cU, 0x3f80c777U,
+      0x3f80ce72U, 0x3f80d19fU, 0x3f80f7bfU, 0x3f80fcfeU, 0x3f81feb4U,
+      0x3f83d731U, 0x3f90cb1dU, 0x3fc55379U, 0x3fd364d7U, 0x41178febU,
+      0x4c5d65a5U, 0x4e85f412U, 0x65d890d3U, 0x6f31a8ecU};
+  for (int i = 0; i < N; ++i) {
+    float x = float(FPBits(INPUTS[i]));
+    EXPECT_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();