Skip to content

Commit

Permalink
Add mixed precition arithmetic functions
Browse files Browse the repository at this point in the history
  • Loading branch information
@eschnett
eschnett committed Nov 26, 2018
1 parent 38c91e3 commit 2c89e7f
Show file tree
Hide file tree
Showing 4 changed files with 142 additions and 1 deletion.
8 changes: 8 additions & 0 deletions src/QD128.jl
View file
Expand Up @@ -391,6 +391,14 @@ end



function Base.rand(::Type{Float128})
r = Ref{D2}()
ccall((:c_dd_rand, libqd), Cvoid, (Ref{D2},), r)
Float128(r[])
end



# Constants

const c128_0 = Float128(0.0)
Expand Down
63 changes: 62 additions & 1 deletion src/QD256.jl
View file
Expand Up @@ -29,7 +29,7 @@ function Float256(x::Float64)
end
function Float256(x::Float128)
r = Ref{D4}()
ccall((:c_qd_copy_qd, libqd), Cvoid, (Ref{D2}, Ref{D4}), x.d2, r)
ccall((:c_qd_copy_dd, libqd), Cvoid, (Ref{D2}, Ref{D4}), x.d2, r)
Float256(r[])
end
function Float256(x::Union{Int8, Int16, Int32, UInt8, UInt16, UInt32,
Expand Down Expand Up @@ -105,12 +105,24 @@ function Base. +(x::Float256, y::Float256)
ccall((:c_qd_add, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{D4}), x.d4, y.d4, r)
Float256(r[])
end
function Base. +(x::Float128, y::Float256)
r = Ref{D4}()
ccall((:c_qd_add_dd_qd, libqd), Cvoid, (Ref{D2}, Ref{D4}, Ref{D4}),
x.d2, y.d4, r)
Float256(r[])
end
function Base. +(x::Float64, y::Float256)
r = Ref{D4}()
ccall((:c_qd_add_d_qd, libqd), Cvoid, (Cdouble, Ref{D4}, Ref{D4}),
x, y.d4, r)
Float256(r[])
end
function Base. +(x::Float256, y::Float128)
r = Ref{D4}()
ccall((:c_qd_add_qd_dd, libqd), Cvoid, (Ref{D4}, Ref{D2}, Ref{D4}),
x.d4, y.d2, r)
Float256(r[])
end
function Base. +(x::Float256, y::Float64)
r = Ref{D4}()
ccall((:c_qd_add_qd_d, libqd), Cvoid, (Ref{D4}, Cdouble, Ref{D4}),
Expand All @@ -125,12 +137,24 @@ function Base. -(x::Float256, y::Float256)
ccall((:c_qd_sub, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{D4}), x.d4, y.d4, r)
Float256(r[])
end
function Base. -(x::Float128, y::Float256)
r = Ref{D4}()
ccall((:c_qd_sub_dd_qd, libqd), Cvoid, (Ref{D2}, Ref{D4}, Ref{D4}),
x.d2, y.d4, r)
Float256(r[])
end
function Base. -(x::Float64, y::Float256)
r = Ref{D4}()
ccall((:c_qd_sub_d_qd, libqd), Cvoid, (Cdouble, Ref{D4}, Ref{D4}),
x, y.d4, r)
Float256(r[])
end
function Base. -(x::Float256, y::Float128)
r = Ref{D4}()
ccall((:c_qd_sub_qd_dd, libqd), Cvoid, (Ref{D4}, Ref{D2}, Ref{D4}),
x.d4, y.d2, r)
Float256(r[])
end
function Base. -(x::Float256, y::Float64)
r = Ref{D4}()
ccall((:c_qd_sub_qd_d, libqd), Cvoid, (Ref{D4}, Cdouble, Ref{D4}),
Expand All @@ -145,12 +169,27 @@ function Base. *(x::Float256, y::Float256)
ccall((:c_qd_mul, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{D4}), x.d4, y.d4, r)
Float256(r[])
end
# This function as well as its counterpart below are much less
# accurate than expected. Their accuracy is that of Float128, not
# Float256. We therefore disable them.
# function Base. *(x::Float128, y::Float256)
# r = Ref{D4}()
# ccall((:c_qd_mul_dd_qd, libqd), Cvoid, (Ref{D2}, Ref{D4}, Ref{D4}),
# x.d2, y.d4, r)
# Float256(r[])
# end
function Base. *(x::Float64, y::Float256)
r = Ref{D4}()
ccall((:c_qd_mul_d_qd, libqd), Cvoid, (Cdouble, Ref{D4}, Ref{D4}),
x, y.d4, r)
Float256(r[])
end
# function Base. *(x::Float256, y::Float128)
# r = Ref{D4}()
# ccall((:c_qd_mul_qd_dd, libqd), Cvoid, (Ref{D4}, Ref{D2}, Ref{D4}),
# x.d4, y.d2, r)
# Float256(r[])
# end
function Base. *(x::Float256, y::Float64)
r = Ref{D4}()
ccall((:c_qd_mul_qd_d, libqd), Cvoid, (Ref{D4}, Cdouble, Ref{D4}),
Expand All @@ -165,12 +204,24 @@ function Base. /(x::Float256, y::Float256)
ccall((:c_qd_div, libqd), Cvoid, (Ref{D4}, Ref{D4}, Ref{D4}), x.d4, y.d4, r)
Float256(r[])
end
function Base. /(x::Float128, y::Float256)
r = Ref{D4}()
ccall((:c_qd_div_dd_qd, libqd), Cvoid, (Ref{D2}, Ref{D4}, Ref{D4}),
x.d2, y.d4, r)
Float256(r[])
end
function Base. /(x::Float64, y::Float256)
r = Ref{D4}()
ccall((:c_qd_div_d_qd, libqd), Cvoid, (Cdouble, Ref{D4}, Ref{D4}),
x, y.d4, r)
Float256(r[])
end
function Base. /(x::Float256, y::Float128)
r = Ref{D4}()
ccall((:c_qd_div_qd_dd, libqd), Cvoid, (Ref{D4}, Ref{D2}, Ref{D4}),
x.d4, y.d2, r)
Float256(r[])
end
function Base. /(x::Float256, y::Float64)
r = Ref{D4}()
ccall((:c_qd_div_qd_d, libqd), Cvoid, (Ref{D4}, Cdouble, Ref{D4}),
Expand All @@ -179,7 +230,9 @@ function Base. /(x::Float256, y::Float64)
end

Base. \(x::Float256, y::Float256) = y / x
Base. \(x::Float128, y::Float256) = y / x
Base. \(x::Float64, y::Float256) = y / x
Base. \(x::Float256, y::Float128) = y / x
Base. \(x::Float256, y::Float64) = y / x


Expand Down Expand Up @@ -403,6 +456,14 @@ end



function Base.rand(::Type{Float256})
r = Ref{D4}()
ccall((:c_qd_rand, libqd), Cvoid, (Ref{D4},), r)
Float256(r[])
end



# Constants

const c256_0 = Float256(0.0)
Expand Down
28 changes: 28 additions & 0 deletions test/test128.jl
View file
Expand Up @@ -48,6 +48,34 @@ for n in 1:1000
@test approxeq(Float128, 2, [x, y], big(x64 \ y), big(x64) \ big(y))
@test approxeq(Float128, 2, [x, y], big(x \ y64), big(x) \ big(y64))

@test isequal(x == x, big(x) == big(x))
@test isequal(x != x, big(x) != big(x))
@test isequal(x < x, big(x) < big(x))
@test isequal(x <= x, big(x) <= big(x))
@test isequal(x > x, big(x) > big(x))
@test isequal(x >= x, big(x) >= big(x))

@test isequal(x == y, big(x) == big(y))
@test isequal(x != y, big(x) != big(y))
@test isequal(x < y, big(x) < big(y))
@test isequal(x <= y, big(x) <= big(y))
@test isequal(x > y, big(x) > big(y))
@test isequal(x >= y, big(x) >= big(y))

@test isequal(x64 == y, big(x64) == big(y))
@test isequal(x64 != y, big(x64) != big(y))
@test isequal(x64 < y, big(x64) < big(y))
@test isequal(x64 <= y, big(x64) <= big(y))
@test isequal(x64 > y, big(x64) > big(y))
@test isequal(x64 >= y, big(x64) >= big(y))

@test isequal(x == y64, big(x) == big(y64))
@test isequal(x != y64, big(x) != big(y64))
@test isequal(x < y64, big(x) < big(y64))
@test isequal(x <= y64, big(x) <= big(y64))
@test isequal(x > y64, big(x) > big(y64))
@test isequal(x >= y64, big(x) >= big(y64))

@test isequal(big(floor(1.0e10 * x)), floor(big(1.0e10 * x)))
@test isequal(big(ceil(1.0e10 * x)), ceil(big(1.0e10 * x)))
@test isequal(big(round(1.0e10 * x)), round(big(1.0e10 * x)))
Expand Down
44 changes: 44 additions & 0 deletions test/test256.jl
View file
Expand Up @@ -21,8 +21,14 @@ for n in 1:1000

@test isequal(Float64(Float256(x64)), x64)

x128 = Float128(rand(BigFloat))
y128 = Float128(rand(BigFloat))

@test isequal(Float128(Float256(x128)), x128)

xbig = rand(BigFloat)
@test isequal(big(Float256(x64)), big(x64))
@test isequal(big(Float256(x128)), big(x128))
@test abs(big(Float256(xbig)) - xbig) <= eps(Float256)

x = Float256(rand(BigFloat))
Expand All @@ -34,25 +40,63 @@ for n in 1:1000
@test approxeq(Float256, 2, [x], big(inv(x)), inv(big(x)))

@test abs(big(x + y) - (big(x) + big(y))) <= eps(Float256)
@test abs(big(x128 + y) - (big(x128) + big(y))) <= eps(Float256)
@test abs(big(x64 + y) - (big(x64) + big(y))) <= eps(Float256)
@test abs(big(x + y128) - (big(x) + big(y128))) <= eps(Float256)
@test abs(big(x + y64) - (big(x) + big(y64))) <= eps(Float256)

@test abs(big(x - y) - (big(x) - big(y))) <= eps(Float256)
@test abs(big(x128 - y) - (big(x128) - big(y))) <= eps(Float256)
@test abs(big(x64 - y) - (big(x64) - big(y))) <= eps(Float256)
@test abs(big(x - y128) - (big(x) - big(y128))) <= eps(Float256)
@test abs(big(x - y64) - (big(x) - big(y64))) <= eps(Float256)

@test abs(big(x * y) - (big(x) * big(y))) <= eps(Float256)
@test abs(big(x128 * y) - (big(x128) * big(y))) <= eps(Float256)
@test abs(big(x64 * y) - (big(x64) * big(y))) <= eps(Float256)
@test abs(big(x * y128) - (big(x) * big(y128))) <= eps(Float256)
@test abs(big(x * y64) - (big(x) * big(y64))) <= eps(Float256)

@test approxeq(Float256, 2, [x, y], big(x / y), big(x) / big(y))
@test approxeq(Float256, 2, [x, y], big(x128 / y), big(x128) / big(y))
@test approxeq(Float256, 2, [x, y], big(x64 / y), big(x64) / big(y))
@test approxeq(Float256, 2, [x, y], big(x / y128), big(x) / big(y128))
@test approxeq(Float256, 2, [x, y], big(x / y64), big(x) / big(y64))

@test approxeq(Float256, 2, [x, y], big(x \ y), big(x) \ big(y))
@test approxeq(Float256, 2, [x, y], big(x128 \ y), big(x128) \ big(y))
@test approxeq(Float256, 2, [x, y], big(x64 \ y), big(x64) \ big(y))
@test approxeq(Float256, 2, [x, y], big(x \ y128), big(x) \ big(y128))
@test approxeq(Float256, 2, [x, y], big(x \ y64), big(x) \ big(y64))

@test isequal(x == x, big(x) == big(x))
@test isequal(x != x, big(x) != big(x))
@test isequal(x < x, big(x) < big(x))
@test isequal(x <= x, big(x) <= big(x))
@test isequal(x > x, big(x) > big(x))
@test isequal(x >= x, big(x) >= big(x))

@test isequal(x == y, big(x) == big(y))
@test isequal(x != y, big(x) != big(y))
@test isequal(x < y, big(x) < big(y))
@test isequal(x <= y, big(x) <= big(y))
@test isequal(x > y, big(x) > big(y))
@test isequal(x >= y, big(x) >= big(y))

@test isequal(x64 == y, big(x64) == big(y))
@test isequal(x64 != y, big(x64) != big(y))
@test isequal(x64 < y, big(x64) < big(y))
@test isequal(x64 <= y, big(x64) <= big(y))
@test isequal(x64 > y, big(x64) > big(y))
@test isequal(x64 >= y, big(x64) >= big(y))

@test isequal(x == y64, big(x) == big(y64))
@test isequal(x != y64, big(x) != big(y64))
@test isequal(x < y64, big(x) < big(y64))
@test isequal(x <= y64, big(x) <= big(y64))
@test isequal(x > y64, big(x) > big(y64))
@test isequal(x >= y64, big(x) >= big(y64))

@test isequal(big(floor(1.0e10 * x)), floor(big(1.0e10 * x)))
@test isequal(big(ceil(1.0e10 * x)), ceil(big(1.0e10 * x)))
@test isequal(big(round(1.0e10 * x)), round(big(1.0e10 * x)))
Expand Down

0 comments on commit 2c89e7f

Please sign in to comment.