Skip to content

Commit

Permalink
RoundNearest argument for rem
Browse files Browse the repository at this point in the history
  • Loading branch information
simonbyrne committed Apr 23, 2015
1 parent 7e63e6a commit 1b84e72
Show file tree
Hide file tree
Showing 5 changed files with 47 additions and 6 deletions.
1 change: 1 addition & 0 deletions base/exports.jl
Original file line number Diff line number Diff line change
Expand Up @@ -407,6 +407,7 @@ export
reinterpret,
rem,
rem1,
rem2pi,
round,
sec,
secd,
Expand Down
2 changes: 1 addition & 1 deletion base/linalg/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ det(D::Diagonal) = prod(D.diag)
logdet{T<:Real}(D::Diagonal{T}) = sum(log(D.diag))
function logdet{T<:Complex}(D::Diagonal{T}) #Make sure branch cut is correct
x = sum(log(D.diag))
-pi<imag(x)<pi ? x : real(x)+(mod2pi(imag(x)+pi)-pi)*im
complex(real(x),rem2pi(imag(x),RoundNearest))
end
# identity matrices via eye(Diagonal{type},n)
eye{T}(::Type{Diagonal{T}}, n::Int) = Diagonal(ones(T,n))
Expand Down
4 changes: 1 addition & 3 deletions base/linalg/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -175,9 +175,7 @@ function logdet{T<:Complex,S}(A::LU{T,S})
if isodd(sum(A.ipiv .!= 1:n))
s = Complex(real(s), imag(s)+π)
end
r, a = reim(s)
a = π-mod2pi-a) #Take principal branch with argument (-pi,pi]
complex(r,a)
complex(real(s),rem2pi(imag(s),RoundNearest))
end

inv!{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}) = @assertnonsingular LAPACK.getri!(A.factors, A.ipiv) A.info
Expand Down
39 changes: 37 additions & 2 deletions base/math.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan,
cbrt, sqrt, erf, erfc, erfcx, erfi, dawson,
significand,
lgamma, hypot, gamma, lfact, max, min, minmax, ldexp, frexp,
clamp, modf, ^, mod2pi,
clamp, modf, ^, mod2pi, rem2pi,
airy, airyai, airyprime, airyaiprime, airybi, airybiprime, airyx,
besselj0, besselj1, besselj, besseljx,
bessely0, bessely1, bessely, besselyx,
Expand All @@ -22,7 +22,7 @@ export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan,

import Base: log, exp, sin, cos, tan, sinh, cosh, tanh, asin,
acos, atan, asinh, acosh, atanh, sqrt, log2, log10,
max, min, minmax, ^, exp2,
max, min, minmax, ^, exp2, rem,
exp10, expm1, log1p,
sign_mask, exponent_mask, exponent_one, exponent_half,
significand_mask, significand_bits, exponent_bits, exponent_bias
Expand Down Expand Up @@ -247,6 +247,12 @@ function frexp{T<:FloatingPoint}(A::Array{T})
return (f, e)
end


rem(x::Float64, y::Float64, ::RoundingMode{:Nearest}) =
ccall((:remainder, libm),Float64,(Float64,Float64),x,y)
rem(x::Float32, y::Float32, ::RoundingMode{:Nearest}) =
ccall((:remainderf, libm),Float32,(Float32,Float32),x,y)

modf(x) = rem(x,one(x)), trunc(x)

const _modff_temp = Float32[0]
Expand Down Expand Up @@ -324,6 +330,35 @@ const pi3o2_l = 1.8369701987210297e-16 # convert(Float64, pi * BigFloat(3/2) -
const pi4o2_h = 6.283185307179586 # convert(Float64, pi * BigFloat(2))
const pi4o2_l = 2.4492935982947064e-16 # convert(Float64, pi * BigFloat(2) - pi4o2_h)

function rem2pi(x::Float64, ::RoundingMode{:Nearest})
abs(x) < pi && return x

(n,y) = ieee754_rem_pio2(x)

if iseven(n)
if n&2 == 2 # add pi
return add22condh(y[1],y[2],pi2o2_h,pi2o2_l)
else # add 0
return y[1]
end
else
if n&2 == 2 # subtract pi/2
return add22condh(y[1],y[2],-pi1o2_h,-pi1o2_l)
else # add pi/2
return add22condh(y[1],y[2],pi1o2_h,pi1o2_l)
end
end
end

rem2pi(x::Float32, ::RoundingMode{:Nearest}) = Float32(rem2pi(Float64(x),RoundNearest))
rem2pi(x::Int32, ::RoundingMode{:Nearest}) = rem2pi(Float64(x),RoundNearest)
function rem2pi(x::Int64, ::RoundingMode{:Nearest})
fx = Float64(x)
fx == x || throw(ArgumentError("Int64 argument to rem2pi is too large: $x"))
rem2pi(fx,RoundNearest)
end


function mod2pi(x::Float64) # or modtau(x)
# with r = mod2pi(x)
# a) 0 <= r < 2π (note: boundary open or closed - a bit fuzzy, due to rem_pio2 implementation)
Expand Down
7 changes: 7 additions & 0 deletions base/mpfr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -541,6 +541,13 @@ function rem(x::BigFloat, y::BigFloat)
return z
end

function rem(x::BigFloat, y::BigFloat, ::RoundingMode{:Nearest})
z = BigFloat()
ccall((:mpfr_remainder, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, &y, ROUNDING_MODE[end])
return z
end


function sum(arr::AbstractArray{BigFloat})
z = BigFloat(0)
for i in arr
Expand Down

0 comments on commit 1b84e72

Please sign in to comment.