Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
169 changes: 0 additions & 169 deletions src/doubletriple/double.jl
Original file line number Diff line number Diff line change
@@ -1,44 +1,3 @@
function double(::Type{T}, x::BigFloat) where {T<:IEEEFloat}
prec = precision(BigFloat)
setprecision(BigFloat, max(prec, 768))
hi = T(x)
lo = T(x - hi)
setprecision(BigFloat, prec)
return hi, lo
end

function double(::Type{T}, x::String) where {T<:IEEEFloat}
prec = precision(BigFloat)
setprecision(BigFloat, 768)
z = parse(BigFloat, x)
hi = T(z)
lo = T(z - hi)
setprecision(BigFloat, prec)
return hi, lo
end

double64(x::BigFloat) = double(Float64, x)
double32(x::BigFloat) = double(Float32, x)
double16(x::BigFloat) = double(Float16, x)

double64(x::String) = double(Float64, x)
double32(x::String) = double(Float32, x)
double16(x::String) = double(Float16, x)

function double_inv(::Type{T}, x::BigFloat) where {T<:IEEEFloat}
prec = precision(BigFloat)
setprecision(BigFloat, max(prec, 768))
x = inv(x)
hi = T(x)
lo = T(x - hi)
setprecision(BigFloat, prec)
return hi, lo
end

double64inv(x::BigFloat) = double_inv(Float64, x)
double32inv(x::BigFloat) = double_inv(Float32, x)
double16inv(x::BigFloat) = double_inv(Float16, x)

#=
algorithms from
Mioara Joldes, Jean-Michel Muller, Valentina Popescu.
Expand Down Expand Up @@ -88,16 +47,6 @@ function TwoSum(a::T, b::T) where {T<:AbstractFloat}
return s, t
end

function TwoDiff(a::T, b::T) where {T<:AbstractFloat}
s = a - b
a1 = s + b
b1 = s - a1
da = a - a1
db = b - b1
t = da + db
return s, t
end

# Algorithm 3 in ref: error-free transformation

@inline function Fast2Mult(a::T, b::T) where {T<:AbstractFloat}
Expand All @@ -115,47 +64,6 @@ end
return zₕᵢ, zₗₒ
end

@inline function DWMinusFP(xₕᵢ::T, xₗₒ::T, y::T) where {T<:AbstractFloat}
sₕᵢ, sₗₒ = TwoDiff(xₕᵢ, y)
v = xₗₒ + sₗₒ
zₕᵢ, zₗₒ = TwoSum(sₕᵢ, v)
return zₕᵢ, zₗₒ
end

# Algorithm 6 in ref: relerr 3u² + 13u³ [reltime 35]

function AccurateDWPlusDW(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
sₕᵢ, sₗₒ = TwoSum(xₕᵢ, yₕᵢ)
tₕᵢ, tₗₒ = TwoSum(xₗₒ, yₗₒ)
c = sₗₒ + tₕᵢ
vₕᵢ, vₗₒ = Fast2Sum(sₕᵢ, c)
w = tₗₒ + vₗₒ
zₕᵢ, zₗₒ = Fast2Sum(vₕᵢ, w)
return zₕᵢ, zₗₒ
end

function AccurateDWMinusDW(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
sₕᵢ, sₗₒ = TwoDiff(xₕᵢ, yₕᵢ)
tₕᵢ, tₗₒ = TwoDiff(xₗₒ, yₗₒ)
c = sₗₒ + tₕᵢ
vₕᵢ, vₗₒ = Fast2Sum(sₕᵢ, c)
w = tₗₒ + vₗₒ
zₕᵢ, zₗₒ = Fast2Sum(vₕᵢ, w)
return zₕᵢ, zₗₒ
end


# Algorithm 7 in ref: relerr (³/₂)u² + 4u³ [reltime 18]

@inline function DWTimesFP1(xₕᵢ::T, xₗₒ::T, y::T) where {T<:AbstractFloat}
cₕᵢ, c1 = Fast2Mult(xₕᵢ, y)
c2 = xₗₒ * y
tₕᵢ, t1 = Fast2Sum(cₕᵢ, c2)
t2 = t1 + c1
zₕᵢ, zₗₒ = Fast2Sum(tₕᵢ, t2)
return zₕᵢ, zₗₒ
end

# Algorithm 9 in ref: relerr 2u² [reltime 15]

@inline function DWTimesFP3(xₕᵢ::T, xₗₒ::T, y::T) where {T<:AbstractFloat}
Expand All @@ -165,86 +73,9 @@ end
return zₕᵢ, zₗₒ
end

# Algorithm 11 in ref: relerr 6u² [reltime 16]

function DWTimesDW2(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
cₕᵢ, c1 = Fast2Mult(xₕᵢ, yₕᵢ)
t0 = xₕᵢ * yₗₒ
c2 = fma(xₗₒ, yₕᵢ, t0)
c3 = c1 + c2
zₕᵢ, zₗₒ = Fast2Sum(cₕᵢ, c3)
return zₕᵢ, zₗₒ
end

# Algorithm 12 in ref: relerr 5u² [reltime 17]

function DWTimesDW3(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
cₕᵢ, c1 = Fast2Mult(xₕᵢ, yₕᵢ)
t0 = xₗₒ * yₗₒ
t1 = fma(xₕᵢ, yₗₒ, t0)
c2 = fma(xₗₒ, yₕᵢ, t1)
c3 = c1 + c2
zₕᵢ, zₗₒ = Fast2Sum(cₕᵢ, c3)
return zₕᵢ, zₗₒ
end

# Algorithm 15 in ref: relerr 3u² [reltime 27]

function DWDivFP3(xₕᵢ::T, xₗₒ::T, y::T) where {T<:AbstractFloat}
tₕᵢ = xₕᵢ / y
pₕᵢ, pₗₒ = Fast2Mult(tₕᵢ, y)
dₕᵢ = xₕᵢ - pₕᵢ
dₗₒ = dₕᵢ - pₗₒ
d = dₗₒ + xₗₒ
tₗₒ = d / y
zₕᵢ, zₗₒ = Fast2Sum(tₕᵢ, tₗₒ)
return zₕᵢ, zₗₒ
end

# Algorithm 17 in ref: relerr 15u² + 56u³ [reltime 50]

function DWDivDW2(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
tₕᵢ = xₕᵢ / yₕᵢ
rₕᵢ, rₗₒ = DWTimesFP1(yₕᵢ, yₗₒ, tₕᵢ)
dₕᵢ = xₕᵢ - rₕᵢ
dₗₒ = xₗₒ - rₗₒ
d = dₕᵢ + dₗₒ
tₗₒ = d / yₕᵢ
zₕᵢ, zₗₒ = Fast2Sum(tₕᵢ, tₗₒ)
return zₕᵢ, zₗₒ
end

# Algorithm 18 in ref: relerr < 10u² (6u² seen) [reltime 107]
# (note DWTimesDW3 replaces DWTimesDW2 per ref)

function DWDivDW3(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
tₕᵢ = inv(yₕᵢ)
rₕᵢ = fma(yₕᵢ, -tₕᵢ, one(T))
rₗₒ = -(yₗₒ * tₕᵢ)
eₕᵢ, eₗₒ = Fast2Sum(rₕᵢ, rₗₒ)
dₕᵢ, dₗₒ = DWTimesFP3(eₕᵢ, eₗₒ, tₕᵢ)
mₕᵢ, mₗₒ = DWPlusFP(dₕᵢ, dₗₒ, tₕᵢ)
zₕᵢ, zₗₒ = DWTimesDW3(xₕᵢ, xₗₒ, mₕᵢ, mₗₒ)
return zₕᵢ, zₗₒ
end

# inv(...) using Algorithms 17 and 18

# Algorithm 17 in ref: relerr 15u² + 56u³ [reltime 48]

function DWInvDW2(yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
tₕᵢ = one(T) / yₕᵢ
rₕᵢ, rₗₒ = DWTimesFP1(yₕᵢ, yₗₒ, tₕᵢ)
dₕᵢ = one(T) - rₕᵢ
dₗₒ = -rₗₒ
d = dₕᵢ + dₗₒ
tₗₒ = d / yₕᵢ
zₕᵢ, zₗₒ = Fast2Sum(tₕᵢ, tₗₒ)
return zₕᵢ, zₗₒ
end

# Algorithm 18 in ref: relerr < 10u² (6u² seen) [reltime 72]
# (note DWTimesDW3 replaces DWTimesDW2 per ref)

function DWInvDW3(yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
tₕᵢ = inv(yₕᵢ)
Expand Down
44 changes: 0 additions & 44 deletions src/doubletriple/double_consts.jl
Original file line number Diff line number Diff line change
@@ -1,41 +1,8 @@
const one_d64 = (1.0, 0.0)

const pi_4o1_d64 = (12.56637061435917200, 4.8985871965894130e-16)
const pi_2o1_d64 = ( 6.28318530717958600, 2.4492935982947064e-16)
const pi_1o1_d64 = ( 3.14159265358979300, 1.2246467991473532e-16)
const pi_1o2_d64 = ( 1.57079632679489660, 6.1232339957367660e-17)
const pi_15o32_d64 = ( 1.47262155637021560, -5.3616983752483473e-17)
const pi_7o16_d64 = ( 1.37444678594553450, 5.3578297462696700e-17)
const pi_1o3_d64 = ( 1.04719755119659790, -1.0720817664510910e-16)
const pi_1o4_d64 = ( 0.78539816339744830, 3.0616169978683830e-17)
const pi_1o6_d64 = ( 0.52359877559829890, -5.3604088322554550e-17)
const pi_1o8_d64 = ( 0.39269908169872414, 1.5308084989341915e-17)
const pi_1o12_d64 = ( 0.26179938779914946, -2.6802044161277275e-17)
const pi_1o16_d64 = ( 0.19634954084936207, 7.6540424946709580e-18)
const pi_1o24_d64 = ( 0.13089969389957473, -1.3401022080638637e-17)
const pi_1o32_d64 = ( 0.09817477042468103, 3.8270212473354790e-18)
const pi_1o64_d64 = ( 0.04908738521234052, 1.9135106236677394e-18)
const pi_1o128_d64 = ( 0.02454369260617026, 9.5675531183386970e-19)
const pi_1o256_d64 = ( 0.01227184630308513, 4.7837765591693480e-19)


const inv_pi_4o1_d64 = ( 0.07957747154594767, -4.9196691687956215e-18)
const inv_pi_2o1_d64 = ( 0.15915494309189535, -9.8393383375912430e-18)
const inv_pi_1o1_d64 = ( 0.31830988618379070, -1.9678676675182486e-17)
const inv_pi_1o2_d64 = ( 0.63661977236758140, -3.9357353350364970e-17)
const inv_pi_15o32_d64 = ( 0.67906109052542010, 2.4632204570453420e-17)
const inv_pi_7o16_d64 = ( 0.72756545413437870, 2.6011543692324530e-18)
const inv_pi_1o3_d64 = ( 0.95492965855137200, -3.5248787942896340e-18)
const inv_pi_1o4_d64 = ( 1.27323954473516280, -7.8714706700729940e-17)
const inv_pi_1o6_d64 = ( 1.90985931710274400, -7.0497575885792670e-18)
const inv_pi_1o8_d64 = ( 2.54647908947032550, -1.5742941340145989e-16)
const inv_pi_1o12_d64 = ( 3.81971863420548800, -1.4099515177158535e-17)
const inv_pi_1o16_d64 = ( 5.09295817894065100, -3.1485882680291977e-16)
const inv_pi_1o24_d64 = ( 7.63943726841097600, -2.8199030354317070e-17)
const inv_pi_1o32_d64 = ( 10.18591635788130200, -6.2971765360583950e-16)
const inv_pi_1o64_d64 = ( 20.37183271576260400, -1.2594353072116790e-15)
const inv_pi_1o128_d64 = ( 40.74366543152521000, -2.5188706144233580e-15)
const inv_pi_1o256_d64 = ( 81.48733086305042000, -5.0377412288467160e-15)

const pi_4o1_d32 = (12.566371f0, -3.496911f-7)
const pi_4o1_d16 = (Float16(12.56), Float16(0.00387))
Expand Down Expand Up @@ -71,14 +38,3 @@ pi1o2(::Type{Double16}) = Double16(pi_1o2_d16)
pi1o4(::Type{Double64}) = Double64(pi_1o4_d64)
pi1o4(::Type{Double32}) = Double32(pi_1o4_d32)
pi1o4(::Type{Double16}) = Double16(pi_1o4_d16)

# below is to be removed
const twopi_d64 = (6.283185307179586, 2.4492935982947064e-16)
const onepi_d64 = (3.141592653589793, 1.2246467991473532e-16)
const halfpi_d64 = (1.5707963267948966, 6.123233995736766e-17)
const qrtrpi_d64 = (0.7853981633974483, 3.061616997868383e-17)

const inv_twopi_d64 = (0.15915494309189535, -9.839338337591243e-18)
const inv_onepi_d64 = (0.3183098861837907, -1.9678676675182486e-17)
const inv_halfpi_d64 = (0.6366197723675814, -3.935735335036497e-17)
const inv_qrtrpi_d64 = (1.2732395447351628, -7.871470670072994e-17)
119 changes: 0 additions & 119 deletions src/doubletriple/quadword.jl

This file was deleted.

Loading
Loading