Re-implementation of Tang (1990) log procedures in pure float32

[sotlampr]

Base.Math.log1p calculation relies on two procedures , log_proc1 and log_proc2 from the literature Tang (1990). The current implementation upscales to Float64 internally and as such is incompatible with Metal.

This PR attempts to convert the two procedures to pure-float32.

See also relevant discussion #234

[sotlampr]

Julia Base uses Float64 tables

[sotlampr]

Similarly, logb in base returns a double

[sotlampr]

From Tang (1990):

(1) Define M := 12 for single precision, and define M := 24 for double precision.
(2) u1 := u rounded (or truncated) to M significant bits.
(3) fi := f rounded (or truncated) to M significant bits.

I'm not sure if this is what is meant

[sotlampr]

Non-upcasting alternative procedure from Tang (1990)

[maleadt]

cc @oscardssmith

[oscardssmith]

What's the error like for these implementations?

[sotlampr]

EDIT: CPU, running on metal has too high error
This log1p(using the modified procedures) compared with Base.Math.log1p (Float32) (ULP):

max 4.0 at x = -0.099609345
mean 0.12487844518189703

This log_proc2 compared with Base.Math.log_proc2 (Float32) (ULP):

log_proc2{Float32}, base=2
max 2.0 at x = 0.043909933
mean 0.1784229523275927

log_proc2{Float32}, base=ℯ
max 1.0 at x = 0.064494334
mean 0.004462545515167737

log_proc2{Float32}, base=10
max 2.0 at x = 0.03663287
mean 0.287801184326243

[sotlampr]

NOTE Running on Metal has very high error - something is not right

[oscardssmith]

Can you compare against Base.Math.log1p (Float64)? your current comparison isn't great since Base.Math.log1p (Float32) also has rounding error.

[sotlampr]

I assumed that we're interested in the difference between the upcasting float32 vs. the pure float32 implementation.

This log1p{Float32} v. Base.Math.log1p{Float64}

max 4.441346645355225 at x = -0.11138753
mean 0.17755624519091973

And for reference, Base.Math.log1p{Float32} v. Base.Math.log1p{Float64}

max 0.563554048538208 at x = 0.08200329
mean 0.11855523097385935

[sotlampr]

@oscardssmith do you think the error is acceptable? And most critically, running on Metal does not work - the error is way too high. Can you see any possible culprit in the code?

[oscardssmith]

the error is probably acceptable, but it is higher than I was expecting. I don't know why it's giving wrong answers on gpu. possibly the truncate is misbehaving?

[sotlampr]

The problems seems to be accessing t_log_Float32 - it returns 0

[oscardssmith]

now that you mention it, a table based method is probably not the best idea for a gpu.

[oscardssmith]

https://github.com/JuliaMath/openlibm/blob/master/src/s_log1pf.c might be a better approach.

[sotlampr]

Replaced by #239