Support for QuadGK

[giordano]

It'd cool to be able to use this package with QuadGK, which enables numerical integration with custom Julia types.

However, this currently doesn't work because the numeric type must be an AbstractFloat or a type that can be converted to an AbstractFloat. This initial issue would be solved by #14.

After that, the first issue I've faced is that there is no eps method for posit numbers:

julia> quadgk(x -> x ^ 2, Posit16(0), Posit16(1))
ERROR: MethodError: no method matching eps(::Type{Posit16})

[milankl]

Okay, I think this is already reason enough to make AbstractPosit <: AbstractFloat, i'll change that. Fingers crossed, that we don't come across another good reason to change it back. From your error message it seems that there is only the eps(::Type) functionality needed and not the eps(::AbstractFloat) which seems to be based on nextfloat and prevfloat. I'm happy to implement that, Posits number are also sorted by bits, that means a bit-wise plus one should do the job for nextfloat(::Posit), and vice versa for prevfloat(::Posit).

[giordano]

I think you need both eps(::Type) and eps(::Posit), because eps for an instance gives the precision around that number, while eps(T) gives the precision around 1

[giordano]

Good, now that there is eps we get another error:

julia> quadgk(x -> x ^ 2, Posit16(0), Posit16(1))
ERROR: MethodError: no method matching exponent(::Posit16)

[milankl]

Any idea whether the exponent function is directly needed or only inside the eps(::Posit) function? Still wondering whether the best way to extract the regime and exponent bits and calculate the base 2 exponent with integer arithmetics or to implement a nextfloat, prevfloat variant for Posits (which should be doable by reinterpreting it as UInt and adding 1) and then to calculate eps based on these.

[giordano]

I think it's only inside eps(::Posit), this is the function that is being called:

eps(x::AbstractFloat) = isfinite(x) ? abs(x) >= floatmin(x) ? ldexp(eps(typeof(x)), exponent(x)) : nextfloat(zero(x)) : oftype(x, NaN)

[milankl]

isfinite, nextfloat, prevfloat, sign, signbit are now implemented. However, following the underlying idea of posit arithmetic (look at the posit circle), we should have a wrap-around behaviour for posits similar to integers

julia> Float64(nextfloat(nextfloat(Posit8(64))))
-64.0

(64 is floatmax(Posit8)). Do you think this causes any issues? Alternatively, we could also aim for a saturation behaviour such that notareal(Posit8) == nextfloat(notareal(Posit8)).

[giordano]

Do you think this causes any issues?

No idea at all! But this could be a good occasion to check if it makes sense or not.

Anyway, a comprehensive set of tests is really needed to make sure that everything is behaving as expected at any time.

[milankl]

I plan to implement the eps(::AbstractPosit) function simply as

eps(p::AbstractPosit) == max(p-prevfloat(p),nextfloat(p)-p)

which avoids the exponent function that you mentioned earlier raises and error as non-existent. Yes, we could also implement the exponent function directly, however, for posits that means we would need to count the regime bits, which is certainly done in the C library somewhere, but I think that here is an easier way for the beginning.

[giordano]

That might work as an interim solution!

[milankl]

Btw, this works now

julia> quadgk(x -> x ^ 2, Posit16(0), Posit16(1))
(Posit16(0x255b), Posit16(0x0080))
julia> Float64.(ans)
(0.33367919921875, 6.103515625e-5)

[giordano]

This is very nice! Comparison with Float16:

julia> using SoftPosit, QuadGK

julia> Float64.(quadgk(x -> x ^ 2, Posit16(0), Posit16(1)))
(0.33367919921875, 6.103515625e-5)

julia> Float64.(quadgk(x -> x ^ 2, Float16(0), Float16(1)))
(0.3330078125, 0.00048828125)