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

[lntue]

Summary

Implement log2f based on RLIBM library correctly rounded for all rounding modes.

[michaelrj]

libc/src/math/generic/common_constants.h

16

is this intended to be a long term solution? It feels like there should be an easier way to handle using hexadecimal float constants.

[lntue]

libc/src/math/generic/common_constants.h

16

I think the long term solution is to build libc with C++17 standard and get rid of these.

[sivachandra]

libc/src/math/generic/common_constants.h

16

Can you check if adding -std=c++17 to COMPILE_OPTIONS fixes the warning problem for you?

[lntue]

libc/src/math/generic/common_constants.h

16

I tried to add -std=c++17 to the CMakeList.txt files in src/math, src/math/generic, and test/src/math but the warnings still show up when the pragmas are removed.

[zimmermann6]

this patch does not apply to the current main branch (db5aceb):

$ patch -p1 -i /tmp/D115828.diff 
patching file libc/config/linux/aarch64/entrypoints.txt
Hunk #1 FAILED at 136.
1 out of 1 hunk FAILED -- saving rejects to file libc/config/linux/aarch64/entrypoints.txt.rej

It seems this patch was built on the branch which adds logf.

[sivachandra]

libc/src/math/generic/common_constants.h

16

I think we pass -std=c++14 explicitly also and that confuses the compiler. You should be able to add -Wno-c++17-extensions to COMPILE_OPTIONS.

[sivachandra]

libc/src/math/generic/common_constants.h

18

When you have large common global data, a better approach would be to put them in an object library like this: https://github.com/llvm/llvm-project/blob/main/libc/src/math/generic/CMakeLists.txt#L49.

And then use it as a dep for all users.

[zimmermann6]

a rebase is needed so that this patch can be applied on the 'main' branch

[lntue]

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

[lntue]

In D115828#3198002, @zimmermann6 wrote:

a rebase is needed so that this patch can be applied on the 'main' branch

I've rebased the patch. Can you check to see if it works now? You might need to patch this on top of the logf patch if possible.

[lntue]

Add -Wno-c++17-extensions compiler option and remove pragma usage.

[lntue]

Move common constants to an object library.

[sivachandra]

Accepting the mechanics of the change. Please wait for Paul and others for acceptance of the math parts.

[zimmermann6]

the new version applies cleanly to the main branch. I have tested it on x86_64 under Linux (haswell). I confirm it is CR for rounding to nearest, and I get 3 failures if I disable the 3 exceptional cases. For other rounding modes I get 8 failures for rounding towards zero (with the exceptional cases), 8 failures too for rounding towards -Inf, and 7 failures for rounding towards +Inf.

I tried with a polynomial generated by Sollya and with this polynomial we need no exceptional cases, and the routine is CR for all rounding modes (please can someone confirm?):

double r = __llvm_libc::fputil::polyeval(
    d, extra_factor, 0x1.71547652b7fefp+0, -0x1.715476500a42ep-1,
    0x1.ec70917f77152p-2, -0x1.71482b204ea69p-2, 0x1.21da0eb07c659p-2);

For the record this polynomial was obtained with the following input file (then run sollya log2f.sollya):

n = 5; /* polynomial degree */
P = 53; /* precision of the coefficients */

pretty = proc(u) {
  return ~(floor(u*1000)/1000);
};

d = [0, 1/2^7];
f = log2(1+x);
w = 1;
p =  remez(f, n, d, w);
pf = fpminimax(log2(1+x), [|1,2,3,4,5|], [|P...|], d, absolute, floating, 0, p)\
;
err_p = -log2(dirtyinfnorm(pf*w-f, d));
print (pf, pretty(err_p));

Sollya is available from https://www.sollya.org/. Would you consider using the Sollya polynomial?

[santoshn]

RLIBM's polynomial should also have 0 special case inputs and produce correctly rounded results for all inputs and for all rounding modes for log2f.

I think the current polynomial has special case inputs because the generated polynomial is generated with a different polynomial evaluation.

Here is a revised polynomial that uses the current polynomial evaluation and produces correct results for all rounding modes and all inputs with zero violated inputs:

Polynomial: y=-3.2945312494298684154217536821552091564491707891340621650044795387657359e-16 x^(0) + 1.4426950408890866217603843324468471109867095947265625000000000000000000e+00 x^(1) + -7.2134752022691861483849606884177774190902709960937500000000000000000000e-01 x^(2) + 4.8089833027421252653610395100258756428956985473632812500000000000000000e-01 x^(3) + -3.6069225263970772221711058591608889400959014892578125000000000000000000e-01 x^(4) + 2.8949201646411226729327381690382026135921478271484375000000000000000000e-01 x^(5)

[lntue]

In D115828#3199419, @zimmermann6 wrote:

the new version applies cleanly to the main branch. I have tested it on x86_64 under Linux (haswell). I confirm it is CR for rounding to nearest, and I get 3 failures if I disable the 3 exceptional cases. For other rounding modes I get 8 failures for rounding towards zero (with the exceptional cases), 8 failures too for rounding towards -Inf, and 7 failures for rounding towards +Inf.

I tried with a polynomial generated by Sollya and with this polynomial we need no exceptional cases, and the routine is CR for all rounding modes (please can someone confirm?):

double r = __llvm_libc::fputil::polyeval(
    d, extra_factor, 0x1.71547652b7fefp+0, -0x1.715476500a42ep-1,
    0x1.ec70917f77152p-2, -0x1.71482b204ea69p-2, 0x1.21da0eb07c659p-2);

For the record this polynomial was obtained with the following input file (then run sollya log2f.sollya):

n = 5; /* polynomial degree */
P = 53; /* precision of the coefficients */

pretty = proc(u) {
  return ~(floor(u*1000)/1000);
};

d = [0, 1/2^7];
f = log2(1+x);
w = 1;
p =  remez(f, n, d, w);
pf = fpminimax(log2(1+x), [|1,2,3,4,5|], [|P...|], d, absolute, floating, 0, p)\
;
err_p = -log2(dirtyinfnorm(pf*w-f, d));
print (pf, pretty(err_p));

Sollya is available from https://www.sollya.org/. Would you consider using the Sollya polynomial?

In D115828#3200306, @santoshn wrote:

RLIBM's polynomial should also have 0 special case inputs and produce correctly rounded results for all inputs and for all rounding modes for log2f.

I think the current polynomial has special case inputs because the generated polynomial is generated with a different polynomial evaluation.

Here is a revised polynomial that uses the current polynomial evaluation and produces correct results for all rounding modes and all inputs with zero violated inputs:

Polynomial: y=-3.2945312494298684154217536821552091564491707891340621650044795387657359e-16 x^(0) + 1.4426950408890866217603843324468471109867095947265625000000000000000000e+00 x^(1) + -7.2134752022691861483849606884177774190902709960937500000000000000000000e-01 x^(2) + 4.8089833027421252653610395100258756428956985473632812500000000000000000e-01 x^(3) + -3.6069225263970772221711058591608889400959014892578125000000000000000000e-01 x^(4) + 2.8949201646411226729327381690382026135921478271484375000000000000000000e-01 x^(5)

Thanks Paul and Santosh for looking into this! It looks like there are many polynomials that will make it correctly rounded for the entire range, which is a great news!

That also make me think that so actually the exception cases from logf (and maybe log10f in the near future) are coming from adding the extra factor: m*log(2) + log(f). So if we can make that addition more accurate, we shouldn't have any exceptional cases. Moreover, it's entirely possible that we can reduce the degree of our polynomials while maintaining the accuracy.

I've discussed with Christoph and see if he can get any update on those directions.

[zimmermann6]

Dear Santosh,

Here is a revised polynomial that uses the current polynomial evaluation and produces correct results for all rounding modes and all inputs with zero violated inputs:

Polynomial: y=-3.2945312494298684154217536821552091564491707891340621650044795387657359e-16 x^(0) + 1.4426950408890866217603843324468471109867095947265625000000000000000000e+00 x^(1) + -7.2134752022691861483849606884177774190902709960937500000000000000000000e-01 x^(2) + 4.8089833027421252653610395100258756428956985473632812500000000000000000e-01 x^(3) + -3.6069225263970772221711058591608889400959014892578125000000000000000000e-01 x^(4) + 2.8949201646411226729327381690382026135921478271484375000000000000000000e-01 x^(5)

this polynomial has a non-zero constant term, unlike the polynomial used in
this patch, and the one produced by Sollya. How can it fit the current
framework (see below)?

double r = __llvm_libc::fputil::polyeval(
    d, extra_factor, 0x1.71547652c2801p+0, -0x1.715476ec167eep-1,
    0x1.ec72eb4428a1ap-2, -0x1.72fd9daa7714fp-2, 0x1.8be682a823a9bp-2);

Best regards,
Paul

[santoshn]

Dear Paul,

Yes, it has a non-zero constant term because it generates the correctly rounded round-to-odd result for a 34-bit floating point. When this is rounded to any FP representation with the 8 exponent bits and for any rounding mode in the standard, it produces the correctly rounded result.
Here is our paper describing the result: https://people.cs.rutgers.edu/~sn349/papers/rlibmall-popl-2022.pdf

We plugged this polynomial locally in our infrastructure which uses the same range reduction and the same quintic polynomial evaluation, it produces correctly rounded results for all inputs and all inputs. Are you observing that it is not producing the correct result for some inputs?

Thanks,
Santosh

In D115828#3202691, @zimmermann6 wrote:

Dear Santosh,

this polynomial has a non-zero constant term, unlike the polynomial used in
this patch, and the one produced by Sollya. How can it fit the current
framework (see below)?

double r = __llvm_libc::fputil::polyeval(
    d, extra_factor, 0x1.71547652c2801p+0, -0x1.715476ec167eep-1,
    0x1.ec72eb4428a1ap-2, -0x1.72fd9daa7714fp-2, 0x1.8be682a823a9bp-2);

Best regards,
Paul

[michaelrj]

LGTM with nits from a formatting perspective.

libc/src/math/generic/log2f.cpp

122

this variable doesn't match the new formatting rules (should be f_index).

libc/test/src/math/log2f_test.cpp

44

These are constexpr and should be capitalized (so COUNT and STEP

[zimmermann6]

Dear Santosh,

Yes, it has a non-zero constant term because it generates the correctly rounded round-to-odd result for a 34-bit floating point. When this is rounded to any FP representation with the 8 exponent bits and for any rounding mode in the standard, it produces the correctly rounded result.
Here is our paper describing the result: https://people.cs.rutgers.edu/~sn349/papers/rlibmall-popl-2022.pdf

We plugged this polynomial locally in our infrastructure which uses the same range reduction and the same quintic polynomial evaluation, it produces correctly rounded results for all inputs and all inputs. Are you observing that it is not producing the correct result for some inputs?

sorry I was not clear enough: can you generate with RLIBM a degree-5
polynomial with double coefficients and a zero constant term that you can
plug into the llvm-libc infrastructure (lines 131-132 of the current
patch) and which produces correct-rounding for all rounding modes (like
the polynomial generated by Sollya which I gave above).

Best regards,
Paul

[santoshn]

Dear Paul,

Thanks, Here is a polynomial with zero constant term that I regenerated with RLIBM. It produces correctly rounded results for all rounding modes and inputs in our infrastructure. It should produce correct results for all representations that has 8 bits of exponents and precisions starting from 10 bits to 32-bits and for all rounding modes.

Polynomial: y=1.4426950408890186761112772728665731847286224365234375000000000000000000e+00 x^(1) + -7.2134752026802795299431636522058397531509399414062500000000000000000000e-01 x^(2) + 4.8089833837385143056053493637591600418090820312500000000000000000000000e-01 x^(3) + -3.6069150703943497759951242187526077032089233398437500000000000000000000e-01 x^(4) + 2.8934750971542422259830118491663597524166107177734375000000000000000000e-01 x^(5)

In D115828#3204271, @zimmermann6 wrote:

Dear Santosh,

Yes, it has a non-zero constant term because it generates the correctly rounded round-to-odd result for a 34-bit floating point. When this is rounded to any FP representation with the 8 exponent bits and for any rounding mode in the standard, it produces the correctly rounded result.
Here is our paper describing the result: https://people.cs.rutgers.edu/~sn349/papers/rlibmall-popl-2022.pdf

We plugged this polynomial locally in our infrastructure which uses the same range reduction and the same quintic polynomial evaluation, it produces correctly rounded results for all inputs and all inputs. Are you observing that it is not producing the correct result for some inputs?

sorry I was not clear enough: can you generate with RLIBM a degree-5
polynomial with double coefficients and a zero constant term that you can
plug into the llvm-libc infrastructure (lines 131-132 of the current
patch) and which produces correct-rounding for all rounding modes (like
the polynomial generated by Sollya which I gave above).

Best regards,
Paul

[zimmermann6]

Dear Santosh,

Thanks, Here is a polynomial with zero constant term that I regenerated with RLIBM. It produces correctly rounded results for all rounding modes and inputs in our infrastructure. It should produce correct results for all representations that has 8 bits of exponents and precisions starting from 10 bits to 32-bits and for all rounding modes.

Polynomial: y=1.4426950408890186761112772728665731847286224365234375000000000000000000e+00 x^(1) + -7.2134752026802795299431636522058397531509399414062500000000000000000000e-01 x^(2) + 4.8089833837385143056053493637591600418090820312500000000000000000000000e-01 x^(3) + -3.6069150703943497759951242187526077032089233398437500000000000000000000e-01 x^(4) + 2.8934750971542422259830118491663597524166107177734375000000000000000000e-01 x^(5)

if I plug this polynomial into the llvm-libc framework as follows:

double r = __llvm_libc::fputil::polyeval(
    d, extra_factor, 0x1.71547652b83f7p+0, -0x1.7154765134294p-1,
    0x1.ec709d3010d4p-2, -0x1.71591d4ab7a18p-2, 0x1.284ab6ada08d6p-2);

then I get 1 incorrect rounding for x=0x1.fc64e8p-1 and rounding to nearest,
and 1 incorrect rounding for x=0x1.197472p+0 and rounding down/towards zero
(for rounding up, all results are correctly rounded).

Best regards,
Paul

[santoshn]

Paul,

Thanks. Let me investigate and test it out with the patch. I have been testing the polynomials with our framework. There seems to be subtle differences between the patch and our public RLIBM framework.

Santosh

In D115828#3204839, @zimmermann6 wrote:

Dear Santosh,

Thanks, Here is a polynomial with zero constant term that I regenerated with RLIBM. It produces correctly rounded results for all rounding modes and inputs in our infrastructure. It should produce correct results for all representations that has 8 bits of exponents and precisions starting from 10 bits to 32-bits and for all rounding modes.

Polynomial: y=1.4426950408890186761112772728665731847286224365234375000000000000000000e+00 x^(1) + -7.2134752026802795299431636522058397531509399414062500000000000000000000e-01 x^(2) + 4.8089833837385143056053493637591600418090820312500000000000000000000000e-01 x^(3) + -3.6069150703943497759951242187526077032089233398437500000000000000000000e-01 x^(4) + 2.8934750971542422259830118491663597524166107177734375000000000000000000e-01 x^(5)

if I plug this polynomial into the llvm-libc framework as follows:

double r = __llvm_libc::fputil::polyeval(
    d, extra_factor, 0x1.71547652b83f7p+0, -0x1.7154765134294p-1,
    0x1.ec709d3010d4p-2, -0x1.71591d4ab7a18p-2, 0x1.284ab6ada08d6p-2);

then I get 1 incorrect rounding for x=0x1.fc64e8p-1 and rounding to nearest,
and 1 incorrect rounding for x=0x1.197472p+0 and rounding down/towards zero
(for rounding up, all results are correctly rounded).

Best regards,
Paul

[santoshn]

Here is a new polynomial that is generated using the exact output compensation in the patch. Tue suggested to use the FMA based poly eval as he was observing performance regressions with the SIMD instruction in x86-64.

Polynomial: y=1.4426950408936214387267682468518614768981933593750000000000000000000000e+00 x^(1) + -7.2134752892795794831926059487159363925457000732421875000000000000000000e-01 x^(2) + 4.8090233829603024062748772848863154649734497070312500000000000000000000e-01 x^(3) + -3.6137987525825709944626851211069151759147644042968750000000000000000000e-01 x^(4) + 3.2929554893140711158139311010017991065979003906250000000000000000000000e-01 x^(5)

Polynomial evaluation used is as follows:

double t1 = fma(x, a5, a4);
double t2 = fma(x, t1, a3);
double t3 = fma(x, t2, a2);
double t4 = fma(x, t3, a1);

final result = fma(d, t4, extra_factor)

Can you check if it produces correctly rounded results for all inputs and all rounding modes?

In D115828#3204858, @santoshn wrote:

Paul,

Thanks. Let me investigate and test it out with the patch. I have been testing the polynomials with our framework. There seems to be subtle differences between the patch and our public RLIBM framework.

Santosh

[zimmermann6]

Dear Santosh,

Here is a new polynomial that is generated using the exact output compensation in the patch. Tue suggested to use the FMA based poly eval as he was observing performance regressions with the SIMD instruction in x86-64.

Polynomial: y=1.4426950408936214387267682468518614768981933593750000000000000000000000e+00 x^(1) + -7.2134752892795794831926059487159363925457000732421875000000000000000000e-01 x^(2) + 4.8090233829603024062748772848863154649734497070312500000000000000000000e-01 x^(3) + -3.6137987525825709944626851211069151759147644042968750000000000000000000e-01 x^(4) + 3.2929554893140711158139311010017991065979003906250000000000000000000000e-01 x^(5)

Polynomial evaluation used is as follows:

double t1 = fma(x, a5, a4);
double t2 = fma(x, t1, a3);
double t3 = fma(x, t2, a2);
double t4 = fma(x, t3, a1);

final result = fma(d, t4, extra_factor)

Can you check if it produces correctly rounded results for all inputs and all rounding modes?

sure. If I converted the coefficients properly to hexadecimal values,
there is still one incorrectly rounded result for rounding towards zero
or down (same input x):

libm wrong by up to 1.01e+00 ulp(s) [1] for x=0x1.03a16ap+0
log2 gives 0x1.4cdc4ap-6
mpfr_log2 gives 0x1.4cdc4cp-6

Best regards,
Paul

[santoshn]

In D115828#3207780, @zimmermann6 wrote:

Dear Santosh,

Here is a new polynomial that is generated using the exact output compensation in the patch. Tue suggested to use the FMA based poly eval as he was observing performance regressions with the SIMD instruction in x86-64.

Polynomial: y=1.4426950408936214387267682468518614768981933593750000000000000000000000e+00 x^(1) + -7.2134752892795794831926059487159363925457000732421875000000000000000000e-01 x^(2) + 4.8090233829603024062748772848863154649734497070312500000000000000000000e-01 x^(3) + -3.6137987525825709944626851211069151759147644042968750000000000000000000e-01 x^(4) + 3.2929554893140711158139311010017991065979003906250000000000000000000000e-01 x^(5)

Polynomial evaluation used is as follows:

double t1 = fma(x, a5, a4);
double t2 = fma(x, t1, a3);
double t3 = fma(x, t2, a2);
double t4 = fma(x, t3, a1);

final result = fma(d, t4, extra_factor)

Can you check if it produces correctly rounded results for all inputs and all rounding modes?

sure. If I converted the coefficients properly to hexadecimal values,
there is still one incorrectly rounded result for rounding towards zero
or down (same input x):

libm wrong by up to 1.01e+00 ulp(s) [1] for x=0x1.03a16ap+0
log2 gives 0x1.4cdc4ap-6
mpfr_log2 gives 0x1.4cdc4cp-6

Best regards,
Paul

Dear Paul,

Thanks for testing it out.

I am seeing the exact oracle result in the local build for input 0x1.03a16ap+0. The result produced by the implementation is 0x1.4cdc4cp-6.

Have you commented out lines src/__support/FPUtil/PolyEval.h:38-42.

If you have not, it is most likely using using x86-64 SIMD extensions that Tue suggested us not to use. It could be one reason for the divergence we are seeing.

Thinking out loud: the round-to-nearest result for this input is exactly the same as round-to-zero result 0x1.4cdc4cp-6. Given that the implementation is producing the correct round-to-nearest, and disagreeing only for round-to-zero, is it possible that there is a bug in the test harness that checks round-to-zero results?

Thanks,
Santosh

[zimmermann6]

Dear Santosh,

Have you commented out lines src/__support/FPUtil/PolyEval.h:38-42.

yes:

• a/libc/src/__support/FPUtil/PolyEval.h

+++ b/libc/src/__support/FPUtil/PolyEval.h
@@ -35,7 +35,7 @@ INLINE_FMA static inline T polyeval(T x, T a0, Ts... a) {
} namespace fputil
} namespace __llvm_libc

-#ifdef LLVM_LIBC_ARCH_X86_64
+#ifdef LLVM_LIBC_ARCH_X86_64_XXX

#include "x86_64/PolyEval.h"

Please find attached the log2f.cpp file I am using.

Best regards,
Paul

• log2f.cpp7 KBDownload

[santoshn]

Dear Paul,

I was wondering how you are testing other rounding modes other than round-to-nearest-ties-to-even.

The coefficients that we are using are identical.

double r = __llvm_libc::fputil::polyeval(

d, extra_factor, 0x1.71547652bd4fp+0, -0x1.7154769b978c7p-1,                                                                                                           
0x1.ec71a99e349c8p-2, -0x1.720d90e6aac6cp-2, 0x1.5132da3583dap-2);

Now this double value "r" needs to be rounded according to the target rounding mode.

The cast on line 142 in the current patch is doing static_cast<float>(r). I assume it is just rounding it to round-to-nearest ties to even.

I am seeing that these coefficients are producing correctly rounded results for the round-to-nearest-ties-to-even.

Further when the double value is specifically rounded to the target rounding mode, it is producing the correctly rounded results for them.

If possible, can you tell me the double value (r) returned by the libc polynomial with the above coefficients for various rounding modes corresponding to the input (x=0x1.03a16ap+0)?

It should produce the same double value "r" for all rounding modes.

In my build, it produces r = 0x1.4cdc4c8p-6, which when rounded to any rounding mode produces the correctly rounded result.

Thanks,
Santosh

In D115828#3207913, @zimmermann6 wrote:

Dear Santosh,

Have you commented out lines src/__support/FPUtil/PolyEval.h:38-42.

yes:

• a/libc/src/__support/FPUtil/PolyEval.h

+++ b/libc/src/__support/FPUtil/PolyEval.h
@@ -35,7 +35,7 @@ INLINE_FMA static inline T polyeval(T x, T a0, Ts... a) {
} namespace fputil
} namespace __llvm_libc

-#ifdef LLVM_LIBC_ARCH_X86_64
+#ifdef LLVM_LIBC_ARCH_X86_64_XXX

#include "x86_64/PolyEval.h"

Please find attached the log2f.cpp file I am using.

Best regards,
Paul

• log2f.cpp7 KBDownload

[zimmermann6]

Dear Santosh,

Date: Thu, 23 Dec 2021 17:37:50 +0000 (UTC)
From: Santosh Nagarakatte via Phabricator <[e-mail]>

Dear Paul,

I was wondering how you are testing other rounding modes other than round-to-nearest-ties-to-even.

as follows:

fesetround(...);
y = log2f (x);

The coefficients that we are using are identical.

double r = __llvm_libc::fputil::polyeval(

d, extra_factor, 0x1.71547652bd4fp+0, -0x1.7154769b978c7p-1,                                                                                                           
0x1.ec71a99e349c8p-2, -0x1.720d90e6aac6cp-2, 0x1.5132da3583dap-2);

Now this double value "r" needs to be rounded according to the target rounding mode.

The cast on line 142 in the current patch is doing static_cast<float>(r). I assume it is just rounding it to round-to-nearest ties to even.

I am seeing that these coefficients are producing correctly rounded results for the round-to-nearest-ties-to-even.

Further when the double value is specifically rounded to the target rounding mode, it is producing the correctly rounded results for them.

If possible, can you tell me the double value (r) returned by the libc polynomial with the above coefficients for various rounding modes corresponding to the input (x=0x1.03a16ap+0)?

It should produce the same double value "r" for all rounding modes.

In my build, it produces r = 0x1.4cdc4c80p-6, which when rounded to any rounding mode produces the correctly rounded result.

here is what I get with the different rounding modes:

FE_TONEAREST:
r=0x1.4cdc4c0000001p-6
FE_TOWARDZERO:
r=0x1.4cdc4bfffffffp-6
FE_UPWARD:
r=0x1.4cdc4c0000001p-6
FE_DOWNWARD:
r=0x1.4cdc4bfffffffp-6

The double value r differs because the rounding mode is also used internally
(for example for the polynomial evaluation).

The code could set internally the rounding mode to FE_TONEAREST, and
restore it before the last rounding, this would solve the issue (at least
for that x-value), but I guess that would be slower than dealing with one
exceptional case.

Best regards,
Paul

[lntue]

In D115828#3209093, @zimmermann6 wrote:

Dear Santosh,

Date: Thu, 23 Dec 2021 17:37:50 +0000 (UTC)
From: Santosh Nagarakatte via Phabricator <[e-mail]>

Dear Paul,

I was wondering how you are testing other rounding modes other than round-to-nearest-ties-to-even.

as follows:

fesetround(...);
y = log2f (x);

The coefficients that we are using are identical.

double r = __llvm_libc::fputil::polyeval(

d, extra_factor, 0x1.71547652bd4fp+0, -0x1.7154769b978c7p-1,                                                                                                           
0x1.ec71a99e349c8p-2, -0x1.720d90e6aac6cp-2, 0x1.5132da3583dap-2);

Now this double value "r" needs to be rounded according to the target rounding mode.

The cast on line 142 in the current patch is doing static_cast<float>(r). I assume it is just rounding it to round-to-nearest ties to even.

I am seeing that these coefficients are producing correctly rounded results for the round-to-nearest-ties-to-even.

Further when the double value is specifically rounded to the target rounding mode, it is producing the correctly rounded results for them.

If possible, can you tell me the double value (r) returned by the libc polynomial with the above coefficients for various rounding modes corresponding to the input (x=0x1.03a16ap+0)?

It should produce the same double value "r" for all rounding modes.

In my build, it produces r = 0x1.4cdc4c80p-6, which when rounded to any rounding mode produces the correctly rounded result.

here is what I get with the different rounding modes:

FE_TONEAREST:
r=0x1.4cdc4c0000001p-6
FE_TOWARDZERO:
r=0x1.4cdc4bfffffffp-6
FE_UPWARD:
r=0x1.4cdc4c0000001p-6
FE_DOWNWARD:
r=0x1.4cdc4bfffffffp-6

The double value r differs because the rounding mode is also used internally
(for example for the polynomial evaluation).

The code could set internally the rounding mode to FE_TONEAREST, and
restore it before the last rounding, this would solve the issue (at least
for that x-value), but I guess that would be slower than dealing with one
exceptional case.

Best regards,
Paul

Thanks Paul and Santosh for getting to the bottom of this! I'm updating the testing infrastructure to fully support all rounding modes and once it completes, I'll update this patch accordingly.

[lntue]

Update the polynomial and make log2f correctly rounded for all rounding modes.

[lntue]

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

[lntue]

@zimmermann6, @santoshn: I've updated the implementation to be correctly rounded for all rounding modes. I also added the exhaustive tests for all rounding modes. Thanks!

[lntue]

Fix variable and constant style.

[zimmermann6]

Dear Tue,

@zimmermann6, @santoshn: I've updated the implementation to be correctly rounded for all rounding modes. I also added the exhaustive tests for all rounding modes. Thanks!

I confirm the new function is correctly rounded for all rounding modes.
Great job!
Paul

[santoshn]

looks good me me

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.log2f
     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/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.log2f
     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/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.log2f
     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/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<"log2f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+
           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>]>,

libc/src/__support/FPUtil/PolyEval.h

INLINE_FMA static inline T polyeval(T x, T a0, Ts... a) {
   return fma(x, polyeval(x, a...), a0);
 }
 
 } // namespace fputil
 } // namespace __llvm_libc
 
 #ifdef LLVM_LIBC_ARCH_X86_64
 
-#include "x86_64/PolyEval.h"
+// [DISABLED] There is a regression with using vectorized version for polyeval
+// compared to the naive Horner's scheme with fma.  Need further investigation
+// #include "x86_64/PolyEval.h"
 
 #endif // LLVM_LIBC_ARCH_X86_64
 
 #else
 
 namespace __llvm_libc {
 namespace fputil {
 

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(log2f)
+
 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)

libc/src/math/generic/CMakeLists.txt

add_entrypoint_object(
   HDRS
     ../ldexpl.h
   DEPENDS
     libc.src.__support.FPUtil.fputil
   COMPILE_OPTIONS
     -O2
 )
 
+add_object_library(
+  common_constants
+  HDRS
+    common_constants.h
+  SRCS
+    common_constants.cpp
+  COMPILE_OPTIONS
+    -Wno-c++17-extensions
+)
+
+add_entrypoint_object(
+  log2f
+  SRCS
+    log2f.cpp
+  HDRS
+    ../log2f.h
+  DEPENDS
+    .common_constants
+    libc.src.__support.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O3
+    -Wno-c++17-extensions
+)
+
 add_entrypoint_object(
   logf
   SRCS
     logf.cpp
   HDRS
     ../logf.h
   DEPENDS
+    .common_constants
     libc.src.__support.FPUtil.fputil
   COMPILE_OPTIONS
     -O3
 )
 
 add_entrypoint_object(
   logb
   SRCS

libc/src/math/generic/common_constants.h

+//===-- Common constants for math functions ---------------------*- 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_GENERIC_COMMON_CONSTANTS_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_COMMON_CONSTANTS_H
+
+namespace __llvm_libc {
+
+// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127.
+extern const double ONE_OVER_F[128];
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_GENERIC_COMMON_CONSTANTS_H

libc/src/math/generic/common_constants.cpp

+//===-- Common constants for math functions ---------------------*- 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "common_constants.h"
+
+namespace __llvm_libc {
+
+// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127.
+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};
+
+} // namespace __llvm_libc

libc/src/math/generic/log2f.cpp

+//===-- Single-precision log2(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/log2f.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 a correctly-rounded algorithm for log2(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, log2(x) = m + log2(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:
+//   log2(1.mant) = log2(f) + log2(1.mant / f)
+//                = log2(f) + log2(1 + d/f)
+//                ~ log2(f) + P(d/f)
+// since d/f is sufficiently small.
+//   log2(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 {
+
+// Lookup table for log2(f) = log2(1 + n*2^(-7)) where n = 0..127.
+static constexpr double LOG2_F[128] = {
+    0x0.0000000000000p+0, 0x1.6fe50b6ef0851p-7, 0x1.6e79685c2d22ap-6,
+    0x1.11cd1d5133413p-5, 0x1.6bad3758efd87p-5, 0x1.c4dfab90aab5fp-5,
+    0x1.0eb389fa29f9bp-4, 0x1.3aa2fdd27f1c3p-4, 0x1.663f6fac91316p-4,
+    0x1.918a16e46335bp-4, 0x1.bc84240adabbap-4, 0x1.e72ec117fa5b2p-4,
+    0x1.08c588cda79e4p-3, 0x1.1dcd197552b7bp-3, 0x1.32ae9e278ae1ap-3,
+    0x1.476a9f983f74dp-3, 0x1.5c01a39fbd688p-3, 0x1.70742d4ef027fp-3,
+    0x1.84c2bd02f03b3p-3, 0x1.98edd077e70dfp-3, 0x1.acf5e2db4ec94p-3,
+    0x1.c0db6cdd94deep-3, 0x1.d49ee4c325970p-3, 0x1.e840be74e6a4dp-3,
+    0x1.fbc16b902680ap-3, 0x1.0790adbb03009p-2, 0x1.11307dad30b76p-2,
+    0x1.1ac05b291f070p-2, 0x1.24407ab0e073ap-2, 0x1.2db10fc4d9aafp-2,
+    0x1.37124cea4cdedp-2, 0x1.406463b1b0449p-2, 0x1.49a784bcd1b8bp-2,
+    0x1.52dbdfc4c96b3p-2, 0x1.5c01a39fbd688p-2, 0x1.6518fe4677ba7p-2,
+    0x1.6e221cd9d0cdep-2, 0x1.771d2ba7efb3cp-2, 0x1.800a563161c54p-2,
+    0x1.88e9c72e0b226p-2, 0x1.91bba891f1709p-2, 0x1.9a802391e232fp-2,
+    0x1.a33760a7f6051p-2, 0x1.abe18797f1f49p-2, 0x1.b47ebf73882a1p-2,
+    0x1.bd0f2e9e79031p-2, 0x1.c592fad295b56p-2, 0x1.ce0a4923a587dp-2,
+    0x1.d6753e032ea0fp-2, 0x1.ded3fd442364cp-2, 0x1.e726aa1e754d2p-2,
+    0x1.ef6d67328e220p-2, 0x1.f7a8568cb06cfp-2, 0x1.ffd799a83ff9bp-2,
+    0x1.03fda8b97997fp-1, 0x1.0809cf27f703dp-1, 0x1.0c10500d63aa6p-1,
+    0x1.10113b153c8eap-1, 0x1.140c9faa1e544p-1, 0x1.18028cf72976ap-1,
+    0x1.1bf311e95d00ep-1, 0x1.1fde3d30e8126p-1, 0x1.23c41d42727c8p-1,
+    0x1.27a4c0585cbf8p-1, 0x1.2b803473f7ad1p-1, 0x1.2f56875eb3f26p-1,
+    0x1.3327c6ab49ca7p-1, 0x1.36f3ffb6d9162p-1, 0x1.3abb3faa02167p-1,
+    0x1.3e7d9379f7016p-1, 0x1.423b07e986aa9p-1, 0x1.45f3a98a20739p-1,
+    0x1.49a784bcd1b8bp-1, 0x1.4d56a5b33cec4p-1, 0x1.510118708a8f9p-1,
+    0x1.54a6e8ca5438ep-1, 0x1.5848226989d34p-1, 0x1.5be4d0cb51435p-1,
+    0x1.5f7cff41e09afp-1, 0x1.6310b8f553048p-1, 0x1.66a008e4788ccp-1,
+    0x1.6a2af9e5a0f0ap-1, 0x1.6db196a76194ap-1, 0x1.7133e9b156c7cp-1,
+    0x1.74b1fd64e0754p-1, 0x1.782bdbfdda657p-1, 0x1.7ba18f93502e4p-1,
+    0x1.7f1322182cf16p-1, 0x1.82809d5be7073p-1, 0x1.85ea0b0b27b26p-1,
+    0x1.894f74b06ef8bp-1, 0x1.8cb0e3b4b3bbep-1, 0x1.900e6160002cdp-1,
+    0x1.9367f6da0ab2fp-1, 0x1.96bdad2acb5f6p-1, 0x1.9a0f8d3b0e050p-1,
+    0x1.9d5d9fd5010b3p-1, 0x1.a0a7eda4c112dp-1, 0x1.a3ee7f38e181fp-1,
+    0x1.a7315d02f20c8p-1, 0x1.aa708f58014d3p-1, 0x1.adac1e711c833p-1,
+    0x1.b0e4126bcc86cp-1, 0x1.b418734a9008cp-1, 0x1.b74948f5532dap-1,
+    0x1.ba769b39e4964p-1, 0x1.bda071cc67e6ep-1, 0x1.c0c6d447c5dd3p-1,
+    0x1.c3e9ca2e1a055p-1, 0x1.c7095ae91e1c7p-1, 0x1.ca258dca93316p-1,
+    0x1.cd3e6a0ca8907p-1, 0x1.d053f6d260896p-1, 0x1.d3663b27f31d5p-1,
+    0x1.d6753e032ea0fp-1, 0x1.d9810643d6615p-1, 0x1.dc899ab3ff56cp-1,
+    0x1.df8f02086af2cp-1, 0x1.e29142e0e0140p-1, 0x1.e59063c8822cep-1,
+    0x1.e88c6b3626a73p-1, 0x1.eb855f8ca88fbp-1, 0x1.ee7b471b3a950p-1,
+    0x1.f16e281db7630p-1, 0x1.f45e08bcf0655p-1, 0x1.f74aef0efafaep-1,
+    0x1.fa34e1177c233p-1, 0x1.fd1be4c7f2af9p-1};
+
+INLINE_FMA
+LLVM_LIBC_FUNCTION(float, log2f, (float x)) {
+  using FPBits = typename fputil::FPBits<float>;
+  FPBits xbits(x);
+  int m = 0;
+
+  // Hard to round value(s).
+  if (FPBits(x).uintval() == 0x3f81d0b5U) {
+    int rounding_mode = fputil::get_round();
+    if (rounding_mode == FE_DOWNWARD || rounding_mode == FE_TOWARDZERO) {
+      return 0x1.4cdc4cp-6f;
+    }
+  }
+
+  // Exceptional inputs.
+  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);
+  // Get the 8 highest bits, use 7 bits (excluding the implicit hidden bit) for
+  // lookup tables.
+  int f_index = xbits.get_mantissa() >> 16;
+
+  FPBits f(xbits.val);
+  // Clear the lowest 16 bits.
+  f.bits &= ~0x0000'FFFF;
+
+  double d = static_cast<float>(xbits) - static_cast<float>(f);
+  d *= ONE_OVER_F[f_index];
+
+  double extra_factor = static_cast<double>(m) + LOG2_F[f_index];
+  double r = __llvm_libc::fputil::polyeval(
+      d, extra_factor, 0x1.71547652bd4fp+0, -0x1.7154769b978c7p-1,
+      0x1.ec71a99e349c8p-2, -0x1.720d90e6aac6cp-2, 0x1.5132da3583dap-2);
+
+  return static_cast<float>(r);
+}
+
+} // namespace __llvm_libc

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 "common_constants.h" // Lookup table for (1/f)
 #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:
static constexpr double LOG_F[128] = {
     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 ||

libc/src/math/log2f.h

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

libc/test/src/math/CMakeLists.txt

add_fp_unittest(
   SRCS
     logf_test.cpp
   DEPENDS
     libc.include.math
     libc.src.math.logf
     libc.src.__support.FPUtil.fputil
 )
 
+add_fp_unittest(
+  log2f_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    log2f_test.cpp
+  DEPENDS
+    libc.include.math
+    libc.src.math.log2f
+    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(
+  log2f_diff
+  SRCS
+    log2f_diff.cpp
+  DEPENDS
+    .single_input_single_output_diff
+    libc.src.math.log2f
+  COMPILE_OPTIONS
+    -fno-builtin
+)
+
+add_diff_binary(
+  log2f_perf
+  SRCS
+    log2f_perf.cpp
+  DEPENDS
+    .single_input_single_output_diff
+    libc.src.math.log2f
+  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

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

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

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

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

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

add_fp_unittest(
   DEPENDS
     .exhaustive_test
     libc.include.math
     libc.src.math.logf
     libc.src.__support.FPUtil.fputil
   LINK_OPTIONS
     -lpthread
 )
+
+add_fp_unittest(
+  log2f_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    log2f_test.cpp
+  DEPENDS
+    .exhaustive_test
+    libc.include.math
+    libc.src.math.log2f
+    libc.src.__support.FPUtil.fputil
+  LINK_OPTIONS
+    -lpthread
+)

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

+//===-- Exhaustive test for log2f -----------------------------------------===//
+//
+// 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/log2f.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 LlvmLibcLog2fExhaustiveTest : 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::Log2, x, __llvm_libc::log2f(x), 0.5,
+                        rounding);
+    } while (bits++ < stop);
+  }
+};
+
+TEST_F(LlvmLibcLog2fExhaustiveTest, RoundNearestTieToEven) {
+  test_full_range(/*start=*/0U, /*stop=*/0x7f80'0000U, /*nthreads=*/16,
+                  mpfr::RoundingMode::Nearest);
+}
+
+TEST_F(LlvmLibcLog2fExhaustiveTest, RoundUp) {
+  test_full_range(/*start=*/0U, /*stop=*/0x7f80'0000U, /*nthreads=*/16,
+                  mpfr::RoundingMode::Upward);
+}
+
+TEST_F(LlvmLibcLog2fExhaustiveTest, RoundDown) {
+  test_full_range(/*start=*/0U, /*stop=*/0x7f80'0000U, /*nthreads=*/16,
+                  mpfr::RoundingMode::Downward);
+}
+
+TEST_F(LlvmLibcLog2fExhaustiveTest, RoundTowardZero) {
+  test_full_range(/*start=*/0U, /*stop=*/0x7f80'0000U, /*nthreads=*/16,
+                  mpfr::RoundingMode::TowardZero);
+}

libc/test/src/math/log2f_test.cpp

+//===-- Unittests for log2f -----------------------------------------------===//
+//
+// 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/log2f.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(LlvmLibcLog2fTest, SpecialNumbers) {
+  EXPECT_FP_EQ(aNaN, __llvm_libc::log2f(aNaN));
+  EXPECT_FP_EQ(inf, __llvm_libc::log2f(inf));
+  EXPECT_TRUE(FPBits(__llvm_libc::log2f(neg_inf)).is_nan());
+  EXPECT_FP_EQ(neg_inf, __llvm_libc::log2f(0.0f));
+  EXPECT_FP_EQ(neg_inf, __llvm_libc::log2f(-0.0f));
+  EXPECT_TRUE(FPBits(__llvm_libc::log2f(-1.0f)).is_nan());
+  EXPECT_FP_EQ(zero, __llvm_libc::log2f(1.0f));
+}
+
+TEST(LlvmLibcLog2fTest, TrickyInputs) {
+  constexpr int N = 9;
+  constexpr uint32_t INPUTS[N] = {0x3f7d57f5U, 0x3f7ed848U, 0x3f7fd6ccU,
+                                  0x3f7fffffU, 0x3f80079bU, 0x3f81d0b5U,
+                                  0x3f82e602U, 0x3f83c98dU, 0x3f8cba39U};
+
+  for (int i = 0; i < N; ++i) {
+    float x = float(FPBits(INPUTS[i]));
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log2, x,
+                                   __llvm_libc::log2f(x), 0.5);
+  }
+}
+
+TEST(LlvmLibcLog2fTest, 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::log2f(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_ALL_ROUNDING(mpfr::Operation::Log2, x,
+                                   __llvm_libc::log2f(x), 0.5);
+  }
+}

libc/utils/MPFRWrapper/MPFRUtils.h

enum class Operation : int {
   Abs,
   Ceil,
   Cos,
   Exp,
   Exp2,
   Expm1,
   Floor,
   Log,
+  Log2,
   Mod2PI,
   ModPIOver2,
   ModPIOver4,
   Round,
   Sin,
   Sqrt,
   Tan,
   Trunc,
#define EXPECT_MPFR_MATCH_ROUNDING(op, input, match_value, ulp_tolerance,      \
   EXPECT_THAT(match_value, __llvm_libc::testing::mpfr::get_mpfr_matcher<op>(   \
                                input, match_value, ulp_tolerance, rounding))
 
 #define EXPECT_MPFR_MATCH(...)                                                 \
   GET_MPFR_MACRO(__VA_ARGS__, EXPECT_MPFR_MATCH_ROUNDING,                      \
                  EXPECT_MPFR_MATCH_DEFAULT, GET_MPFR_DUMMY_ARG)                \
   (__VA_ARGS__)
 
+#define EXPECT_MPFR_MATCH_ALL_ROUNDING(op, input, match_value, ulp_tolerance)  \
+  {                                                                            \
+    namespace mpfr = __llvm_libc::testing::mpfr;                               \
+    mpfr::ForceRoundingMode __r1(mpfr::RoundingMode::Nearest);                 \
+    EXPECT_MPFR_MATCH(op, input, match_value, ulp_tolerance,                   \
+                      mpfr::RoundingMode::Nearest);                            \
+    mpfr::ForceRoundingMode __r2(mpfr::RoundingMode::Upward);                  \
+    EXPECT_MPFR_MATCH(op, input, match_value, ulp_tolerance,                   \
+                      mpfr::RoundingMode::Upward);                             \
+    mpfr::ForceRoundingMode __r3(mpfr::RoundingMode::Downward);                \
+    EXPECT_MPFR_MATCH(op, input, match_value, ulp_tolerance,                   \
+                      mpfr::RoundingMode::Downward);                           \
+    mpfr::ForceRoundingMode __r4(mpfr::RoundingMode::TowardZero);              \
+    EXPECT_MPFR_MATCH(op, input, match_value, ulp_tolerance,                   \
+                      mpfr::RoundingMode::TowardZero);                         \
+  }
+
 #define ASSERT_MPFR_MATCH_DEFAULT(op, input, match_value, ulp_tolerance)       \
   ASSERT_THAT(match_value,                                                     \
               __llvm_libc::testing::mpfr::get_mpfr_matcher<op>(                \
                   input, match_value, ulp_tolerance,                           \
                   __llvm_libc::testing::mpfr::RoundingMode::Nearest))
 
 #define ASSERT_MPFR_MATCH_ROUNDING(op, input, match_value, ulp_tolerance,      \
                                    rounding)                                   \
   ASSERT_THAT(match_value, __llvm_libc::testing::mpfr::get_mpfr_matcher<op>(   \
                                input, match_value, ulp_tolerance, rounding))
 
 #define ASSERT_MPFR_MATCH(...)                                                 \
   GET_MPFR_MACRO(__VA_ARGS__, ASSERT_MPFR_MATCH_ROUNDING,                      \
                  ASSERT_MPFR_MATCH_DEFAULT, GET_MPFR_DUMMY_ARG)                \
   (__VA_ARGS__)
 
+#define ASSERT_MPFR_MATCH_ALL_ROUNDING(op, input, match_value, ulp_tolerance)  \
+  {                                                                            \
+    namespace mpfr = __llvm_libc::testing::mpfr;                               \
+    mpfr::ForceRoundingMode __r1(mpfr::RoundingMode::Nearest);                 \
+    ASSERT_MPFR_MATCH(op, input, match_value, ulp_tolerance,                   \
+                      mpfr::RoundingMode::Nearest);                            \
+    mpfr::ForceRoundingMode __r2(mpfr::RoundingMode::Upward);                  \
+    ASSERT_MPFR_MATCH(op, input, match_value, ulp_tolerance,                   \
+                      mpfr::RoundingMode::Upward);                             \
+    mpfr::ForceRoundingMode __r3(mpfr::RoundingMode::Downward);                \
+    ASSERT_MPFR_MATCH(op, input, match_value, ulp_tolerance,                   \
+                      mpfr::RoundingMode::Downward);                           \
+    mpfr::ForceRoundingMode __r4(mpfr::RoundingMode::TowardZero);              \
+    ASSERT_MPFR_MATCH(op, input, match_value, ulp_tolerance,                   \
+                      mpfr::RoundingMode::TowardZero);                         \
+  }
+
 #endif // LLVM_LIBC_UTILS_TESTUTILS_MPFRUTILS_H

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
   }
 
   MPFRNumber log() const {
     MPFRNumber result(*this);
     mpfr_log(result.value, value, mpfr_rounding);
     return result;
   }
 
+  MPFRNumber log2() const {
+    MPFRNumber result(*this);
+    mpfr_log2(result.value, value, mpfr_rounding);
+    return result;
+  }
+
   MPFRNumber remquo(const MPFRNumber &divisor, int &quotient) {
     MPFRNumber remainder(*this);
     long q;
     mpfr_remquo(remainder.value, &q, value, divisor.value, mpfr_rounding);
     quotient = q;
     return remainder;
   }
 
unary_operation(Operation op, InputType input, unsigned int precision,
   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::Log2:
+    return mpfrInput.log2();
   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();