Approximating log2(float) by separating significand and exponent

Question

Problem

I am trying to implement a fast float = log2(float). I wrote a simple C program to compare my results with log2 and I get a small error that I can't find a source for.

Background

I am taking the approach of identifying the floating point representation (ignoring sign bit) is 2^(exponent) * 1.significand. Using properties of logs I get log2(float) = exp + log2(1.significand). Eventually I'll truncate the significand for a table lookup, but for now I wanted to verify the correct result.

For further background reading on inspiration for this: http://www.icsi.berkeley.edu/pubs/techreports/TR-07-002.pdf

code

Herein lies the problem. This is a simple program that extracts floating point bits and adds the exponent to the log2(significand).

#include <math.h>
#include <stdint.h>
#include <stdio.h>

int main()
{
    typedef union {
      int32_t i;
      float f;
    } poly32_t;

    float x = 31415926535.8;
    poly32_t one;
    one.f = 1.0f;

    uint32_t ii;
    uint32_t num_iter = 15;
    for(ii=0; ii < num_iter; ++ii) {
        poly32_t poly_x;
        poly32_t poly_x_exponent;
        poly32_t poly_x_significand;

        // extract the exponent and significand
        poly_x.f = x;
        poly_x_significand.i = (0x007fffff & poly_x.i);
        poly_x_exponent.i = (0xff & (poly_x.i >> 23) ) - 127;

        // recover the hidden 1 of significand
        poly_x_significand.f = 1.0 + ((float)poly_x_significand.i)/10000000;
        // log2(2^exp * sig) = exponent + log2(significand)
        float log_sig = log2(poly_x_significand.f);
        float y_approx = (float)poly_x_exponent.i + log_sig;

        // Get the actual value
        float y_math = log2(x);
        printf("math::log2(%16.4f)=%8.4f  ;  approx=%8.4f  ;  diff=%.4f\n",
                                x, y_math,      y_approx, y_math-y_approx);
        x *= 0.1;
    }
    return 1;
}

Output

math::log2(31415926784.0000)= 34.8708  ;  approx= 34.7614  ;  diff=0.1094
math::log2( 3141592576.0000)= 31.5488  ;  approx= 31.4733  ;  diff=0.0755
math::log2(  314159264.0000)= 28.2269  ;  approx= 28.1927  ;  diff=0.0342
math::log2(   31415926.0000)= 24.9050  ;  approx= 24.7924  ;  diff=0.1126
math::log2(    3141592.5000)= 21.5831  ;  approx= 21.5036  ;  diff=0.0794
math::log2(     314159.2500)= 18.2611  ;  approx= 18.2221  ;  diff=0.0390
math::log2(      31415.9258)= 14.9392  ;  approx= 14.8235  ;  diff=0.1158
math::log2(       3141.5925)= 11.6173  ;  approx= 11.5340  ;  diff=0.0833
math::log2(        314.1592)=  8.2954  ;  approx=  8.2517  ;  diff=0.0437
math::log2(         31.4159)=  4.9734  ;  approx=  4.8546  ;  diff=0.1188
math::log2(          3.1416)=  1.6515  ;  approx=  1.5644  ;  diff=0.0871
math::log2(          0.3142)= -1.6704  ;  approx= -1.7187  ;  diff=0.0483
math::log2(          0.0314)= -4.9924  ;  approx= -4.9936  ;  diff=0.0012
math::log2(          0.0031)= -8.3143  ;  approx= -8.4050  ;  diff=0.0907
math::log2(          0.0003)=-11.6362  ;  approx=-11.6890  ;  diff=0.0528

The "approximation" should be exactly the same as the clib's log2 at this point; any help identifying my error would be greatly appreciated.

why don't you use the `frexp` function from the standard library — Cheers and hth. - Alf, Oct 23 '14 at 19:08
Two reasons. 1) I wasn't aware of it until now. 2) My end goal is a SIMD implementation where I can't use that function. Thanks for bringing it to my attention. — n-west, Oct 23 '14 at 19:34
Note: The `1.0` in `poly_x_significand.f = 1.0 + ((float)...` should be 0.0 when `(0xff & (poly_x.i >> 23) )` is 0. — chux - Reinstate Monica, Oct 23 '14 at 19:43
Yes, I know according to IEEE 754 denormalized numbers have a special interpretation. I happen to not care about denormalized numbers because my end-goal is targeting NEON which flushes denormalized numbers to 0. — n-west, Oct 23 '14 at 20:16

John Bollinger · Accepted Answer · 2014-10-23T19:13:34.567

5

Edited (2) to remove incorrect information about the number of significand bits

Furthermore, I think this is wrong:

    // recover the hidden 1 of significand
    poly_x_significand.f = 1.0 + ((float)poly_x_significand.i)/10000000;

The (explicit part of the) significand is a binary fraction, not a decimal one. You should be dividing by (float) (1 << 23).

Note, too, that it doesn't look like your implementation deals correctly with subnormal numbers.

edited Oct 23 '14 at 19:13

answered Oct 23 '14 at 18:56

John Bollinger

160,171
8
81
157

2

A 32-bit float according to IEEE 754 has a 23-bit significand, 8 bit exponent, and 1 sign bit. – n-west Oct 23 '14 at 18:59
You're right @n-west, I don't know what I was thinking. I edited that part out. – John Bollinger Oct 23 '14 at 19:05
1

Thanks @John, the binary fraction was the trick. However, because of the hidden 1 the fraction needs to be divided by `(float)(2<22)`. With that edit, you've got the solution. – n-west Oct 23 '14 at 19:11
1

Yes, 1 << 23, or (2 << 22). Evidently I am sloppy today. Fixed. – John Bollinger Oct 23 '14 at 19:14
@n-west [binary32](http://en.wikipedia.org/wiki/Single-precision_floating-point_format) has a 24-bit significand (23 explicitly stored). – chux - Reinstate Monica Oct 23 '14 at 19:46
1

@n-west Did you measure the speed compared to library log2(), I'm curious what effect would table cache misses have under different conditions. – 2501 Oct 23 '14 at 19:59
1

I, too, am curious about the prospects for improving on the library version of `log2()`. I would expect it already to be implemented pretty efficiently. – John Bollinger Oct 23 '14 at 20:15
@2501 : My use case doesn't requirean exact calculation, so it's fairly easy to get faster than a library that has a requirement of giving exact results. I used a very small table (6-bit lookup) to test out the approximation in C, then moved to NEON. Right now I'm at a 10x speedup, but I'll need to tweak it some more because there's high % error around the fringes (near 1 and near 2) of my significand polynomial. See https://github.com/n-west/gnuradio/blob/f75b9b3008551dac04255588f55d92b3116a6a59/volk/kernels/volk/volk_32f_binary_log_32f.h#L186 – n-west Oct 25 '14 at 00:55
@JohnBollinger, see previous comment for accuracy. Apparently SO only lets you tag one user/post. – n-west Oct 25 '14 at 00:56

Approximating log2(float) by separating significand and exponent

Problem

Background

code

Output

1 Answers1