#!/usr/bin/gnuplot

# Aliasing a float32 to an int32 yields an approximation of its base-2
# logarithm.

# Parameters for IEEE-754 single-precision floats
L = 2**23   # mantissa scaling factor
B = 127     # exponent bias

# For m ∈ [0, 1), log₂(1 + m) ≈ 1 + σ
sigma = 0.0430357

# Assume x is a positive normal float32.
# Compute I(x): the int32 aliased to x.
log_of_2 = log(2)
log2(x) = log(x) / log_of_2
round(x) = floor(x + 0.5)
e(x) = floor(log2(x))               # exponent
m(x, e_x) = x / 2**e_x - 1          # mantissa, without the leading 1
M(x, e_x) = round(m(x, e_x) * L)    # scaled and rounded mantissa
I2(x, e_x) = (e_x + B) * L + M(x, e_x)
I(x) = I2(x, e(x))

# Best fit to a logarithm
scaled_log(x) = L * log2(x) + L * (B - sigma)

# Do the plot
set terminal svg size 640, 480 enhanced
set output 'Log_by_aliasing_to_int.svg'
set xrange [0:10]
set yrange [0x3d800000:]
set format y "0x%x"
set ytics L
set key top left Left reverse invert spacing 1.4
set samples 200
set style line 1 lw 3 lc rgb "gray80"
set style line 2 lw 1 lc rgb "blue"
set style increment user
set label 1 "L = 2^{23}"    at 6, 0x3fc00000
set label 2 "B = 127"       at 6, 0x3f600000
set label 3 "σ = 0.0430357" at 6, 0x3f000000
plot scaled_log(x) title "L log_2(x) + L (B − σ)", \
     I(x) title 'I_x' lt 3