Return type of entropy changes with the base

[sgaure]

I came across a somewhat unexpected property which I'm uncertain what I think of.

p = rand(Float32, 10)
p ./= sum(p)
typeof(entropy(p))    # Float32
typeof(entropy(p, 2))  # Float64
typeof(entropy(p, MathConstants.e)) # Float32

When specifying a base b, the result is "scaled by 1/log(b)" according to the docs, which is almost true. It is divided by log(b), which explains the last example. log(MathConstants.e) is an Int64, and dividing a Float32 by an Int64 yields a Float32.

This isn't a big problem, everything is type stable, but it's unexpected if one otherwise operates on Float32 (or Float16 or something else). It's of course easy to make it return Float32 in the second case by letting b be 2f0.

Anyway, should the return type of entropy depend on the type of the base b at all? Shouldn't it have the element type of p? Or the same type as the default base e entropy? I'm not suggesting a change right away, because it could possibly break existing software. But it should be considered.

[sgaure]

On second thought, it could be someone did something like this:

import ForwardDiff
ForwardDiff.derivative(x -> entropy(p, x), 2)

which would break if entropy(p,x) returns the type of entropy(p)