# DoubleFloats.jl -\> Double128 | Quad64?
**URL:** https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049
**Category:** Numerics
**Created:** [December 1, 2020, 3:33pm UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049 "2020-12-01T15:33:07Z")
**Posts on this page:** 20
**Page:** 2
### Author:  [@abulak](https://discourse.julialang.org/u/abulak)
#### Post date: [December 3, 2020, 4:58pm UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/21 "2020-12-03T16:58:21Z")
Agreeable, `Arblib` is very, very young and you have the great opportunity 😃 to be the first serious user. One thing that is lacking is making it play nice with julia `AbstractFloat` protocol (? I know there is nothing formal like this, but maybe someone has hints ?), I opened an issue here: [https://github.com/kalmarek/Arblib.jl/issues/101](https://github.com/kalmarek/Arblib.jl/issues/101) please feel free to drop in and report missing bits and pieces 🙂
Eg., You need to define `Base.expm1(x::T) where T<:Arblib.ArbOrRef = Arblib.expm1!(T(prec=precision(x)), x)`. I’m sure there will be much more.
btw. instead of printing it’s much better to use `@debug` statements and enable/disable them setting `ENV["JULIA_DEBUG"] = "ThorinDistributions"`
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [December 3, 2020, 5:01pm UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/22 "2020-12-03T17:01:01Z")
> [@abulak](#):
>
> Agreeable, `Arblib` is very, very young and you have the great opportunity 😃 to be the first serious user. One thing that is lacking is making it play nice with julia `AbstractFloat` protocol (? I know there is nothing formal like this, but maybe someone has hints ?), I opened an issue here: [AbstractFloat interface functions · Issue #101 · kalmarek/Arblib.jl · GitHub](https://github.com/kalmarek/Arblib.jl/issues/101) please feel free to drop in and report missing bits and pieces 🙂
> Eg., You need to define `Base.expm1(x::T) where T<:Arblib.ArbOrRef = Arblib.expm1!(T(prec=precision(x)), x)` . I’m sure there will be much more.
>
> btw. instead of printing it’s much better to use `@debug` statements and enable/disable them setting `ENV["JULIA_DEBUG"] = "ThorinDistributions"`
Hey, you know what? I’ll try the upper example on Arblib and report issues right now.
---
### Author:  [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)
#### Post date: [December 3, 2020, 6:45pm UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/23 "2020-12-03T18:45:07Z")
I got these results running your example with the additional type `Float128` from Quadmath.jl.
256 bits, 40x40
setprecision(BigFloat, 256)
setworkingprecision(ArbFloat, 256); setextrabits(0)
julia\> display(Errors)
Dict{Type,Tuple{Float64,Float64,Float64}} with 6 entries:
Float64 =\> (1.01502e20, 0.277835, 1.37131e8)
ArbFloat =\> (-1.5526e-6, 1.94999, 2.22239e8)
Double64 =\> (0.000136315, 0.352428, 2.7385e8)
BigFloat =\> (0.000182493, 1.53269, 2.9704e8)
Float128 =\> (6.09923e-5, 0.510003, 2.95367e8)
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [December 3, 2020, 8:11pm UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/24 "2020-12-03T20:11:41Z")
> [@JeffreySarnoff](#):
>
> got these results running your example with the additional type `Float128` from Quadmath.jl.
>
> 256 bits, 40x40
> setprecision(BigFloat, 256)
> setworkingprecision(ArbFloat, 256); setextrabits(0)
>
> julia\> display(Errors)
> Dict{Type,Tuple{Float64,Float64,Float64}} with 6 entries:
> Float64 =\> (1.01502e20, 0.277835, 1.37131e8)
> ArbFloat =\> (-1.5526e-6, 1.94999, 2.22239e8)
> Double64 =\> (0.000136315, 0.352428, 2.7385e8)
> BigFloat =\> (0.000182493, 1.53269, 2.9704e8)
> Float128 =\> (6.09923e-5, 0.510003, 2.95367e8)
There was a mistake in the code I posted, the Radon.seed!() Calls needs to be inside the loop otherwise we are comparing it on different data.
---
### Author:  [@abulak](https://discourse.julialang.org/u/abulak)
#### Post date: [December 4, 2020, 1:49pm UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/25 "2020-12-04T13:49:22Z")
I glimpsed over `MultivariateGammaConvolution` and the first thing I see is that you use abstract fields:
[https://github.com/lrnv/ThorinDistributions.jl/blob/735be65ec76ca98968ac0e7f7bf38805e7a52f89/src/MultivariateGammaConvolution.jl#L21](https://github.com/lrnv/ThorinDistributions.jl/blob/735be65ec76ca98968ac0e7f7bf38805e7a52f89/src/MultivariateGammaConvolution.jl#L21)
it’d be much better to parametrize it concretely:
```julia
struct MultivariateGammaConvolution{T<:Real, V<:AbstractVector{T}, M<:AbstractMatrix{T}, XXX} <: Distributions.ContinuousMultivariateDistribution where T
α::V
θ::M
constants::XXX
end
```
then you can turn some functions into static ones:
```julia
#### eltype, length, support
eltype(d::MultivariateGammaConvolution{T, V}) where {T,V} = V
Base.length(d::MultivariateGammaConvolution) = size(d.θ,2)
Distributions.insupport(d::MultivariateGammaConvolution, x) = all(>(zero(eltype(x)), x)
# or all(>(zero(first(x))), x) ← this will probably be better when using Arbs that come with their own precision
#### Sampling
struct MGCSPL <: Distributions.Sampleable{Distributions.Multivariate,Distributions.Continuous}
Γs::Vector{MGSPL} # is MGSPL concrete?
end
Distributions.sampler(d::MultivariateGammaConvolution) = MGCSPL([Distributions.sampler(MultivariateGamma(d.α[i],d.θ[i,:])) for i in 1:length(d.α)])
```
and so on
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [December 4, 2020, 2:58pm UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/26 "2020-12-04T14:58:58Z")
Thanks for taking the time !
> [@abulak](#):
>
> it’d be much better to parametrize it concretely:
This is not the first time i hear that. Could you explain to me / point me to the doc where it is explained why ?
Following your advice, i’ll probably rewrite the whole package 😉
---
### Author:  [@Elrod](https://discourse.julialang.org/u/Elrod)
#### Post date: [December 4, 2020, 3:02pm UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/27 "2020-12-04T15:02:39Z")
> [@lrnv](#):
>
> This is not the first time i hear that. Could you explain to me / point me to the doc where it is explained why ?
[Performance Tips](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-fields-with-abstract-containers)
---
### Author:  [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)
#### Post date: [December 8, 2020, 1:19am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/28 "2020-12-08T01:19:49Z")
So, use `Float128` from `Quadmath.jl`.
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [December 8, 2020, 9:25am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/29 "2020-12-08T09:25:32Z")
> [@JeffreySarnoff](#):
>
> So, use `Float128` from `Quadmath.jl` .
Yeah no. As i said, 128 bits is not enough. There was a mistake in my MWE, we were comparing on differen data. This is fixed in this version :
> **MWE (click to show)**
>
> ```julia
> # You should ]add github.com/lrnv/ThorinDistributions.jl
> import Random, ThorinDistributions
>
> using DoubleFloats
> using ArbNumerics
> using Quadmath
> setprecision(256)
> setprecision(ArbFloat,256)
>
> types = (Float64,
> Double64,
> Float128,
> BigFloat,
> ArbFloat
> )
>
> n = 20
> d = 2
> m = Tuple(repeat([80],d))
> N = 1000
>
> Random.seed!(1234)
> α = rand(n)
> θ = reshape(10 .* rand(n*d),(n,d))
>
> init = (0.0, 0.0, 0.0)
> init_val = (zeros(m),zeros(m))
>
> Values = Dict(T => init_val for T in types)
> Errors = Dict(T => init for T in types)
> for T in types
> print(T,"\n")
>
> alpha = T.(α);
> scales = T.(θ);
> dist = ThorinDistributions.MultivariateGammaConvolution(alpha,scales);
> samples = zeros(T,(2,N));
> Random.seed!(123);
> Random.rand!(dist,samples);
>
> E = ThorinDistributions.empirical_coefs(samples,m); # empirical coefs
> coefs = ThorinDistributions.get_coefficients(alpha, scales,m); # theoreticla coefs
> aloc = @allocated ThorinDistributions.get_coefficients(alpha,scales,m);
> time = Base.@elapsed(ThorinDistributions.get_coefficients(alpha,scales,m));
>
>
> # All these should have the same magnitude, very small.
> Errors[T] = (Float64(sum((E - coefs)^2)),time,aloc);
> Values[T] = (Float64.(E), Float64.(coefs));
> print(Errors[T],"\n")
> end
>
> # Type => (Error, time1, allocations)
> display(Errors)
>
> ```
ArbFloat is not working proprerly at the moment (it’s probably my fault), but i do have the following results (after a few runs to get everything compiled):
```julia
Float128 => (4.51278e74, 1.77423, 467600.0)
BigFloat => (2.48274e-8, 11.7816, 8.16053e9)
Double64 => (2.10301e78, 0.470318, 467600.0)
Float64 => (2.68084e110, 0.121068, 286112.0)
```
11 seconds to evaluate the cost function is simply unnacceptble, I cannot afford this (Meaning: this is big enough that i can spend a few weeks trying to solve this issue 🤣 ). Therefore the problem remains. As you said earlyer, this may be fixable by an algorithm change, so that is what i’m trying to do right now (without a lot of success…); The source of overflow is in the `ThorinDistibutions.build_coefficients!` function if you wanna take a look.
---
### Author:  [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)
#### Post date: [December 8, 2020, 9:52am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/30 "2020-12-08T09:52:04Z")
One thing you could do is use `view` or `@view` (see online help) with your precomputed matricies and vectors. Currently, it appears that your code is doing a good deal of copying that can be avoided this way. Create the views where you precompute them, and use the names of those views rather than the names of the underlying arrays in your computational loops.
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [December 8, 2020, 10:27am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/31 "2020-12-08T10:27:41Z")
Okay, seems like a good idea. Anyway, the main bottneck of my code, when running in bigfloats, is `Base.MPFR.*`, taking approx 80% of my runtime.
I am trying to change the algorithm to both avoid overflow and reduce the number of multiplications, without success for the moment.
---
### Author:  [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)
#### Post date: [December 8, 2020, 11:01am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/32 "2020-12-08T11:01:02Z")
I don’t know how much of that overhead is from quietly copying BigFloats; if it is substantial, using views may help overall.
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [December 8, 2020, 11:02am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/33 "2020-12-08T11:02:57Z")
Looks like most of the overhead is in the bigfloat products, copying overhead happens _inside_ the bigfloat multiplication. Anayway a dirty `@view` fix reduced the runtime by 5-8%, which is still something.
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [December 8, 2020, 11:14am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/34 "2020-12-08T11:14:03Z")
Okay, i might have something. Using MultiFloats.jl, with a personal (dirty) fix for exp / log, i have the following 😛
> **Code that I ran**
>
> ```julia
> # ]add github.com/lrnv/ThorinDistributions.jl
> # ]add https://github.com/lrnv/MultiFloats.jl
> import Random, ThorinDistributions
>
> using DoubleFloats
> using ArbNumerics
> using Quadmath
> using MultiFloats
> setprecision(256)
> setprecision(ArbFloat,256)
>
> types = (Float64,
> Double64,
> Float128,
> BigFloat,
> Float64x1,
> Float64x2,
> Float64x3,
> Float64x4,
> Float64x5,
> Float64x6,
> Float64x7,
> Float64x8
> )
>
> n = 20
> d = 2
> m = Tuple(repeat([80],d))
> N = 1000
>
> Random.seed!(1234)
> α = rand(n)
> θ = reshape(10 .* rand(n*d),(n,d))
>
> init = (0.0, 0.0, 0.0)
> init_val = (zeros(m),zeros(m))
>
> Values = Dict(T => init_val for T in types)
> Errors = Dict(T => init for T in types)
> for T in types
> print(T,"\n")
>
> alpha = T.(α);
> scales = T.(θ);
> dist = ThorinDistributions.MultivariateGammaConvolution(alpha,scales);
> samples = zeros(T,(2,N));
> Random.seed!(123);
> Random.rand!(dist,samples);
>
> E = ThorinDistributions.empirical_coefs(samples,m); # empirical coefs
> coefs = ThorinDistributions.get_coefficients(alpha, scales,m); # theoreticla coefs
> aloc = @allocated ThorinDistributions.get_coefficients(alpha,scales,m);
> time = Base.@elapsed(ThorinDistributions.get_coefficients(alpha,scales,m));
>
>
> # All these should have the same magnitude, very small.
> Errors[T] = (Float64(sum((E - coefs)^2)),time,aloc);
> Values[T] = (Float64.(E), Float64.(coefs));
> print(Errors[T],"\n")
> end
>
> # Type => (Error, time1, allocations)
> display(Errors)
>
> ```
Which gives me the following result :
```julia
Float64 => (2.68084e110, 0.121672, 286112.0)
Float128 => (4.51278e74, 1.84456, 467600.0)
Double64 => (2.10301e78, 0.516504, 467600.0)
MultiFloat{Float64,1} => (4.1367e109, 0.124149, 296656.0)
MultiFloat{Float64,2} => (3.84925e78, 0.26139, 486208.0)
MultiFloat{Float64,3} => (-6.13933e44, 0.66715, 676464.0)
MultiFloat{Float64,4} => (3.1583e14, 1.33451, 864800.0)
MultiFloat{Float64,5} => (2.367e-8, 2.39511, 1.05234e6)
MultiFloat{Float64,6} => (2.36701e-8, 4.46531, 1.25013e6)
MultiFloat{Float64,7} => (2.36701e-8, 6.50383, 1.44102e6)
MultiFloat{Float64,8} => (2.36701e-8, 9.84249, 1.63157e6)
BigFloat => (2.48274e-8, 12.6092, 8.16053e9)
```
It seems like the `Float64x5` might be enough for this case. The dirty fix i used is that epx/log were not implemented, so i used a convertion to bigfloat back and forth to use the bigfloat implementation. @JeffreySarnoff Do you know if there is some cleaner algorithm that could be implemented for exp/log in this library ?
---
### Author:  [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)
#### Post date: [December 8, 2020, 11:17am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/35 "2020-12-08T11:17:31Z")
Here is a faster way to multiply BigFloats (it avoids the internal allocation).
```julia
function mul!(z::BigFloat, x::BigFloat, y::BigFloat)
ccall(("mpfr_mul",:libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode), z, x, y, Base.MPFR.ROUNDING_MODE[])
return z
end
mul! (generic function with 1 method)
# or
function mul!(z::BigFloat, x::BigFloat, y::BigFloat)
ccall(("mpfr_mul",:libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode), z, x, y, Base.MPFR.ROUNDING_MODE[])
return nothing
end
z=BigFloat(); x = BigFloat(pi); y = sqrt(BigFloat(2));
mul!(z, x, y)
z
```
(about 1.3x faster on a single multiply)
---
### Author:  [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)
#### Post date: [December 8, 2020, 11:25am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/36 "2020-12-08T11:25:02Z")
This is the first I have seen `MultiFloats.jl`. The algorithms for `exp` and `log` are interdependent and, for extended precision, generally difficult because the algebraic equivalences/relationships do not hold in floating point representations.
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [December 8, 2020, 11:29am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/37 "2020-12-08T11:29:01Z")
> [@JeffreySarnoff](#):
>
> Here is a faster way to multiply BigFloats
Seems great. Is there some equivalent hack for the `addmul` of 3 or 4 bigfloats at once ? The main overhead of my code is in two lines, of the form `x += z * y * t * u` and `x += z * y * t`.
---
### Author:  [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)
#### Post date: [December 8, 2020, 11:36am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/38 "2020-12-08T11:36:27Z")
Here is the equivalent for +:
```julia
function add!(z::BigFloat, x::BigFloat, y::BigFloat)
ccall(("mpfr_add",:libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode), z, x, y, Base.MPFR.ROUNDING_MODE[])
return z
end
```
---
### Author:  [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)
#### Post date: [December 8, 2020, 11:44am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/39 "2020-12-08T11:44:10Z")
This gives you z = fma(x,y,k)
```julia
function fma!(z::BigFloat, x::BigFloat, y::BigFloat, k::BigFloat)
ccall(("mpfr_fma",:libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode), z, x, y, k, Base.MPFR.ROUNDING_MODE[])
return z
end
```
This gives you z = x1_y1 + x2_y2 (works like fma does)
```julia
function fmma!(z::BigFloat, x1::BigFloat, y1::BigFloat, x2::BigFloat, y2::BigFloat)
ccall(("mpfr_fmma",:libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode), z, x1, y1, x2, y2, Base.MPFR.ROUNDING_MODE[])
return z
end
```
That’s all there is in the library that works directly.
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [December 8, 2020, 11:46am UTC](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049/40 "2020-12-08T11:46:07Z")
> [@JeffreySarnoff](#):
>
> That’s all there is in the library that works directly.
So for the rest i should compose these ones, right ?
[Previous page](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049.md?page=1)
[Next page](https://discourse.julialang.org/t/doublefloats-jl-double128-quad64/51049.md?page=3)