Some missing features about doublefloats.jl HOT 18 CLOSED

Thanks -- (also play with the stuff that appears to be working pls)

from doublefloats.jl.

oh -- there is no support for Complex as yet

from doublefloats.jl.

Sure? Simple complex operations work for more.

I have another thing, which is more like a wish. I currently use excessively muladd) to speed up calculcations. Currently Double just uses the fallback muladd(a, b, c) = a * b + c. Do you know whether somebody considered this special case so far?

from doublefloats.jl.

I do not understand your question about muladd. I require fma for everything in AccurateArithmetic. Are you suggesting we implement a Double aware fma and then use it in place of the default muladd?

we have the starter technology -- imagine this with hi, lo weaving

"""
    fma_(a, b, c) => (x, y, z)
Computes `x = fl(fma_(a, b, c))` and `y, z = fl(err(fma_(a, b, c)))`.
"""
function fma_(a::T, b::T, c::T) where {T<:AbstractFloat}
     x = fma(a, b, c)
     y, z = mul_(a, b)
     t, z = add_(c, z)
     t, u = add_(y, t)
     y = ((t - x) + u)
     y, z = add_hilo_(y, z)
     return x, y, z
end

from doublefloats.jl.

I thought whether there exists a special method for operations of the form a * b + c, similar to that we have square for a * a. But this is probably way too fancy...

from doublefloats.jl.

because * and + are error-free, there is usually a good outcome,
probably the way to do that would be to do Double( (a * b) + triple(c) ) -- that seems worthwhile for the Accurate stuff and I have the addition of a triple with a double already -- I will do that . Where do you think it helps most?

from doublefloats.jl.

or do you want exported that and call it fma?

from doublefloats.jl.

Cool! I think to overload Base.muladd is the best. One application of this is Horner's method to evaluate polynomials which basically boils down to nested muladds

julia> @macroexpand Base.Math.@horner a c1 c2 c3 c4 c5
quote
    #16031#t = a
    (Base.Math.muladd)(#16031#t, (Base.Math.muladd)(#16031#t, (Base.Math.muladd)(#16031#t, (Base.Math.muladd)(#16031#t, c5, c4), c3), c2), c1)
end

This could also be useful the evaluation of analytic functions.

from doublefloats.jl.

do that for Double and FastDouble?
or ...
do that for FastDouble and do triplemul tripleadd for Double
or ..
leave FastDouble alone

from doublefloats.jl.

julia> thi = Float64(tbig);tlo=Float64(tbig-thi);Double(thi,tlo)
Double(1.8674849396984466e14, 0.003872194160694882)

julia> muladd(d4,d2,d3)
Double(1.8674849396984466e14, 0.003872194160694882)

julia> d4*d2+d3
Double(1.8674849396984466e14, 0.0038721941606948826)

# fast double equivs

julia> fd4*fd2+fd3
FastDouble(1.8674849396984466e14, 0.0038721941606948826)

julia> muladd(fd4,fd2,fd3)
FastDouble(1.8674849396984466e14, 0.0038721941606948826)

from doublefloats.jl.

done ...

from doublefloats.jl.

now to look at your early notes re Inf NaN etc

from doublefloats.jl.

(I am on this today -- its the best time to question and ask)

from doublefloats.jl.

I noticed 3 more things:

• An implementation of rand would be nice. I think it doesn't need to be perfectly uniform for now. There is an implementation here, maybe this is useful.
• sin und cos currently allocate. I think this is because pio32, sin_pio32 and cos_pio32 are non-constant global variables. A const in front should fix this.
• I would also add an @inline in front of sin_circle and cos_circle such that the compiler can eliminate the redundant computations in calls of sincos.

Otherwise everything looks really nice :)

from doublefloats.jl.

If you can massage this, I will include the result:

import Random: rand
using Random

function Random.rand(rng::AbstractRNG, S::Type{Double{T,E}}) where {T<:AbstractFloat, E<:Emphasis}
    u = Random.rand(rng, UInt64)
    f = Float64(u)
    uf = UInt64(f)
    ur = uf > u ? uf - u : u - uf

    Double(E, T(5.421010862427522e-20 * f), T(5.421010862427522e-20 * Float64(ur)))
end


Random.rand(::Type{Double{T,E}}) where {T<:AbstractFloat, E<:Emphasis} = Random.rand(Base.Random.GLOBAL_RNG, Double{T,E})

from doublefloats.jl.

also, when you want the result of an external package's function to return FastDouble, it is a Good idea to call the func (e.g. rand) and then FastDouble(result) or (maybe check the result type eltype) and
map(FastDouble, resultvec)

from doublefloats.jl.

rand(Double), rand(FastDouble), rand(Double, n::Int), rand(FastDouble, n::Int)

current master is the pre-release

from doublefloats.jl.

v002 fixes some overlooked bugs and makes rand(Double, n) work better.
I need the FastDoubles to stay fast, as the Doubles now stay accurate.

from doublefloats.jl.