Skip to content

add tests and fix bugs in modrem #250

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
JeffreySarnoff merged 1 commit into from
Mar 4, 2026
Merged

add tests and fix bugs in modrem #250

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
32 changes: 10 additions & 22 deletions src/math/arithmetic/modpi.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
## This work was supported in part by the U.S. Department of Energy, Office of Science, Office of Workforce Development for Teachers and Scientists (WDTS) under the Science Undergraduate Laboratory Internships Program (SULI) program at Oak Ridge National Laboratory, administered by the Oak Ridge Institute for Science and Education.

function mod2pi(x::DoubleFloat{T}) where {T<:IEEEFloat}
s = signbit(x)
s = x < 0
c = s ? -1 : 1
y = c * x
y < pi2o1(DoubleFloat{T}) && return s ? pi2o1(DoubleFloat{T}) - y : y
Expand Down Expand Up @@ -96,26 +96,6 @@ end

# ========================================== ^^^^^ =========================================

const pi_1o1_t64_hi = pi_1o1_t64[1]
const pi_1o1_t64_md = pi_1o1_t64[2]
const pi_1o1_t64_lo = pi_1o1_t64[3]

const pi_1o4_t64_hi = pi_1o4_t64[1]
const pi_1o4_t64_md = pi_1o4_t64[2]
const pi_1o4_t64_lo = pi_1o4_t64[3]

function mod1pipm(x::DoubleFloat{Float64})
abs(x) < onepi_d64 && return x
w1 = mul323(inv_pi_1o1_t64, HILO(x))
w2 = two_sum(w1[1] - trunc(w1[1]), w1[2], w1[3])
y = mul322(pi_1o1_t64, w2)
z = Double64(y)
return z
end

mod1pipm(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(mod1pipm(DoubleFloat{Float64}(x)))
mod1pipm(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(mod1pipm(DoubleFloat{Float64}(x)))

#=
rem2pi(x) = x - 2π*round(x/(2π),r)
rem2pi(x, RoundDown) == mod2pi(x)
Expand All @@ -127,9 +107,17 @@ mod1pipm(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(mod1pipm(DoubleFloat{Fl
• if r == RoundUp, then the result is in the interval [-2π, 0].
=#

function rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Nearest})
abs(x) < onepi && return x
w1 = mul323(inv_pi_2o1_t64, HILO(x))
w2 = two_sum(w1[1] - round(w1[1]), w1[2], w1[3])
y = mul322(pi_2o1_t64, w2)
z = Double64(y)
return z
end

rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Down}) = mod2pi(x)
rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Up}) = -rem2pi(-x, RoundDown)
rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Nearest}) = mod1pipm(x)
rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:ToZero}) =
signbit(x) ? rem2pi(x, RoundUp) : rem2pi(x, RoundDown)

Expand Down
11 changes: 11 additions & 0 deletions test/modrem.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,7 @@
@testset "modpi" begin

@test mod2pi(Double64(-13)) ≈ mod2pi(big(-13.)) atol = 1e-31

x=cbrt(41)*sqrt(Double64(pi))

@test DoubleFloats.mod2pi(x) == Double64(6.111805926475162, -1.667563077572613e-16)
Expand All @@ -13,4 +16,12 @@ end

@test DoubleFloats.rem2pi(x) == Double64(6.111805926475162, -1.667563077572613e-16)
@test DoubleFloats.rem2pi(-x) == Double64(-0.1713793807044248, 4.647966647701768e-18)

v = Double64.(-7:7)
bigv = big.(v)
for rounding in (RoundUp, RoundDown, RoundNearest, RoundToZero)
drem = rem2pi.(v, Ref(rounding))
bigrem = rem2pi.(bigv, Ref(rounding))
@test all(abs.(bigrem-drem) .< 1e-31)
end
end
Loading