From e99f547a01ea1d7c93f3e491c6f28d6840ccce88 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sat, 28 May 2022 15:38:43 -0500 Subject: [PATCH 1/5] Fixes in power with rationals and nthroot --- src/intervals/arithmetic/power.jl | 93 ++++++++++++------------------- 1 file changed, 37 insertions(+), 56 deletions(-) diff --git a/src/intervals/arithmetic/power.jl b/src/intervals/arithmetic/power.jl index 5ba05e048..f40d76563 100644 --- a/src/intervals/arithmetic/power.jl +++ b/src/intervals/arithmetic/power.jl @@ -130,6 +130,8 @@ function ^(a::F, x::Rational{R}) where {F<:Interval, R<:Integer} isempty(a) && return emptyinterval(a) iszero(x) && return one(a) + # x < 0 && return inv(a^(-x)) + x < 0 && return F( inv(bigequiv(a^(-x)) ) ) if isthinzero(a) x > zero(x) && return zero(a) @@ -140,44 +142,26 @@ function ^(a::F, x::Rational{R}) where {F<:Interval, R<:Integer} x == (1//2) && return sqrt(a) - alo = inf(a) - ahi = sup(a) - - if x >= 0 - if alo ≥ 0 - abig = @biginterval(a) - ui = convert(Culong, q) - low = BigFloat() - high = BigFloat() - ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , inf(abig) , ui, MPFRRoundDown) - ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , sup(abig) , ui, MPFRRoundUp) - b = interval(low, high) - b = convert(F, b) - return b^p - end - - if alo < 0 && ahi ≥ 0 - a = a ∩ F(0, Inf) - abig = @biginterval(a) - ui = convert(Culong, q) - low = BigFloat() - high = BigFloat() - ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , inf(abig) , ui, MPFRRoundDown) - ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , sup(abig) , ui, MPFRRoundUp) - b = interval(low, high) - b = convert(F, b) - return b^p - end - - if ahi < 0 - return emptyinterval(a) - end + alo, ahi = bounds(a) + if ahi < 0 + return emptyinterval(a) end - if x < 0 - return inv(a^(-x)) + if alo < 0 && ahi ≥ 0 + a = a ∩ F(0, Inf) end + + abig = bigequiv(a) + abiglo, abighi = bounds(abig) + ui = convert(Culong, q) + low = BigFloat() + high = BigFloat() + ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , abiglo , ui, MPFRRoundDown) + ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , abighi , ui, MPFRRoundUp) + b = interval(low, high) ^ p + return F(b) + end # Interval power of an interval: @@ -297,35 +281,32 @@ end Compute the real n-th root of Interval. """ function nthroot(a::Interval{BigFloat}, n::Integer) + isempty(a) && return a n == 1 && return a n == 2 && return sqrt(a) - n < 0 && isthinzero(a) && return emptyinterval(a) - isempty(a) && return a - if n > 0 - alo = inf(a) - ahi = sup(a) - ahi < 0 && iseven(n) && return emptyinterval(BigFloat) - if alo < 0 && ahi >= 0 && iseven(n) - a = a ∩ Interval{BigFloat}(0, Inf) - end - ui = convert(Culong, n) - low = BigFloat() - high = BigFloat() - ccall((:mpfr_rootn_ui, :libmpfr), Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , a.lo , ui, MPFRRoundDown) - ccall((:mpfr_rootn_ui, :libmpfr), Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , a.hi , ui, MPFRRoundUp) - b = interval(low , high) - return b - elseif n < 0 - return inv(nthroot(a, -n)) - elseif n == 0 - return emptyinterval(a) + n == 0 && return emptyinterval(a) + # n < 0 && isthinzero(a) && return emptyinterval(a) + n < 0 && return inv(nthroot(a, -n)) + + alo, ahi = bounds(a) + ahi < 0 && iseven(n) && return emptyinterval(BigFloat) + if alo < 0 && ahi >= 0 && iseven(n) + a = a ∩ Interval{BigFloat}(0, Inf) + alo, ahi = bounds(a) end + ui = convert(Culong, n) + low = BigFloat() + high = BigFloat() + ccall((:mpfr_rootn_ui, :libmpfr), Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , alo , ui, MPFRRoundDown) + ccall((:mpfr_rootn_ui, :libmpfr), Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , ahi , ui, MPFRRoundUp) + return interval(low , high) end -function nthroot(a::Interval{T}, n::Integer) where T +function nthroot(a::F, n::Integer) where {T, F<:Interval{T}} n == 1 && return a n == 2 && return sqrt(a) + n < 0 && return inv(nthroot(a, -n)) abig = @biginterval(a) b = nthroot(abig, n) - return convert(Interval{T}, b) + return F(b) end From 33055dacb27be68738c3756e767303aa1368992f Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sun, 29 May 2022 14:12:58 -0500 Subject: [PATCH 2/5] Fix nthroot for subnormals --- src/intervals/arithmetic/power.jl | 34 +++++++++++++++---------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/intervals/arithmetic/power.jl b/src/intervals/arithmetic/power.jl index f40d76563..002e91ef7 100644 --- a/src/intervals/arithmetic/power.jl +++ b/src/intervals/arithmetic/power.jl @@ -10,7 +10,7 @@ # Write explicitly like this to avoid ambiguity warnings: for T in (:Integer, :Float64, :BigFloat, :Interval) - @eval ^(a::Interval{Float64}, x::$T) = Interval{Float64}(bigequiv(a)^x) + @eval ^(a::Interval{Float64}, x::$T) = Interval{Float64}((@biginterval(a))^x) end @@ -29,11 +29,11 @@ Base.literal_pow(::typeof(^), x::Interval{T}, ::Val{p}) where {T,p} = x^p Implement the `pow` function of the IEEE Std 1788-2015 (Table 9.1). """ -^(a::F, b::F) where {F<:Interval} = F(bigequiv(a)^b) +^(a::F, b::F) where {F<:Interval} = F((@biginterval(a))^b) ^(a::F, x::AbstractFloat) where {F<:Interval{BigFloat}} = a^big(x) for T in (:AbstractFloat, :Integer) - @eval ^(a::F, x::$T) where {F<:Interval} = F(bigequiv(a)^x) + @eval ^(a::F, x::$T) where {F<:Interval} = F((@biginterval(a))^x) end function ^(a::F, n::Integer) where {F<:Interval{BigFloat}} @@ -131,7 +131,7 @@ function ^(a::F, x::Rational{R}) where {F<:Interval, R<:Integer} isempty(a) && return emptyinterval(a) iszero(x) && return one(a) # x < 0 && return inv(a^(-x)) - x < 0 && return F( inv(bigequiv(a^(-x)) ) ) + x < 0 && return F( inv( (@biginterval(a))^(-x) ) ) if isthinzero(a) x > zero(x) && return zero(a) @@ -152,16 +152,11 @@ function ^(a::F, x::Rational{R}) where {F<:Interval, R<:Integer} a = a ∩ F(0, Inf) end - abig = bigequiv(a) - abiglo, abighi = bounds(abig) - ui = convert(Culong, q) - low = BigFloat() - high = BigFloat() - ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , low , abiglo , ui, MPFRRoundDown) - ccall((:mpfr_rootn_ui, :libmpfr) , Int32 , (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode) , high , abighi , ui, MPFRRoundUp) - b = interval(low, high) ^ p - return F(b) + b = nthroot( @biginterval(a), q) + p == 1 && return F(b) + + return F(b^p) end # Interval power of an interval: @@ -276,11 +271,11 @@ function log1p(a::F) where {T, F<:Interval{T}} end """ - nthroot(a::Interval{BigFloat}, n::Integer) + nthroot(a::Interval, n::Integer) Compute the real n-th root of Interval. """ -function nthroot(a::Interval{BigFloat}, n::Integer) +function nthroot(a::F, n::Integer) where {F<:Interval{BigFloat}} isempty(a) && return a n == 1 && return a n == 2 && return sqrt(a) @@ -291,7 +286,7 @@ function nthroot(a::Interval{BigFloat}, n::Integer) alo, ahi = bounds(a) ahi < 0 && iseven(n) && return emptyinterval(BigFloat) if alo < 0 && ahi >= 0 && iseven(n) - a = a ∩ Interval{BigFloat}(0, Inf) + a = a ∩ F(0, Inf) alo, ahi = bounds(a) end ui = convert(Culong, n) @@ -305,8 +300,13 @@ end function nthroot(a::F, n::Integer) where {T, F<:Interval{T}} n == 1 && return a n == 2 && return sqrt(a) - n < 0 && return inv(nthroot(a, -n)) + abig = @biginterval(a) + if n < 0 + issubnormal(mag(a)) && return inv(nthroot(a, -n)) + return F( inv(nthroot(abig, -n)) ) + end + b = nthroot(abig, n) return F(b) end From e3c21f0029f2d353895e19dd6169b62efd3b7764 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sun, 29 May 2022 14:16:37 -0500 Subject: [PATCH 3/5] Remove '@biginterval' in favor of lbigequiv' --- src/intervals/arithmetic/power.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/intervals/arithmetic/power.jl b/src/intervals/arithmetic/power.jl index 002e91ef7..855d90412 100644 --- a/src/intervals/arithmetic/power.jl +++ b/src/intervals/arithmetic/power.jl @@ -10,7 +10,7 @@ # Write explicitly like this to avoid ambiguity warnings: for T in (:Integer, :Float64, :BigFloat, :Interval) - @eval ^(a::Interval{Float64}, x::$T) = Interval{Float64}((@biginterval(a))^x) + @eval ^(a::Interval{Float64}, x::$T) = Interval{Float64}((bigequiv(a))^x) end @@ -29,11 +29,11 @@ Base.literal_pow(::typeof(^), x::Interval{T}, ::Val{p}) where {T,p} = x^p Implement the `pow` function of the IEEE Std 1788-2015 (Table 9.1). """ -^(a::F, b::F) where {F<:Interval} = F((@biginterval(a))^b) +^(a::F, b::F) where {F<:Interval} = F((bigequiv(a))^b) ^(a::F, x::AbstractFloat) where {F<:Interval{BigFloat}} = a^big(x) for T in (:AbstractFloat, :Integer) - @eval ^(a::F, x::$T) where {F<:Interval} = F((@biginterval(a))^x) + @eval ^(a::F, x::$T) where {F<:Interval} = F((bigequiv(a))^x) end function ^(a::F, n::Integer) where {F<:Interval{BigFloat}} @@ -131,7 +131,7 @@ function ^(a::F, x::Rational{R}) where {F<:Interval, R<:Integer} isempty(a) && return emptyinterval(a) iszero(x) && return one(a) # x < 0 && return inv(a^(-x)) - x < 0 && return F( inv( (@biginterval(a))^(-x) ) ) + x < 0 && return F( inv( (bigequiv(a))^(-x) ) ) if isthinzero(a) x > zero(x) && return zero(a) @@ -152,7 +152,7 @@ function ^(a::F, x::Rational{R}) where {F<:Interval, R<:Integer} a = a ∩ F(0, Inf) end - b = nthroot( @biginterval(a), q) + b = nthroot( bigequiv(a), q) p == 1 && return F(b) @@ -301,7 +301,7 @@ function nthroot(a::F, n::Integer) where {T, F<:Interval{T}} n == 1 && return a n == 2 && return sqrt(a) - abig = @biginterval(a) + abig = bigequiv(a) if n < 0 issubnormal(mag(a)) && return inv(nthroot(a, -n)) return F( inv(nthroot(abig, -n)) ) From ab117e3664eb9953b47891152e5aef4787dfd24a Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Wed, 1 Jun 2022 12:41:31 -0500 Subject: [PATCH 4/5] Small improvement in nthroot --- src/intervals/arithmetic/power.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/intervals/arithmetic/power.jl b/src/intervals/arithmetic/power.jl index 855d90412..7c0627c8b 100644 --- a/src/intervals/arithmetic/power.jl +++ b/src/intervals/arithmetic/power.jl @@ -301,12 +301,11 @@ function nthroot(a::F, n::Integer) where {T, F<:Interval{T}} n == 1 && return a n == 2 && return sqrt(a) - abig = bigequiv(a) if n < 0 issubnormal(mag(a)) && return inv(nthroot(a, -n)) return F( inv(nthroot(abig, -n)) ) end - b = nthroot(abig, n) + b = nthroot(bigequiv(a), n) return F(b) end From d9ffe0e1a0b8f76af7dbd05787fbd880e559f16f Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Wed, 1 Jun 2022 12:51:20 -0500 Subject: [PATCH 5/5] Revert the last commit :-) --- src/intervals/arithmetic/power.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/intervals/arithmetic/power.jl b/src/intervals/arithmetic/power.jl index 7c0627c8b..855d90412 100644 --- a/src/intervals/arithmetic/power.jl +++ b/src/intervals/arithmetic/power.jl @@ -301,11 +301,12 @@ function nthroot(a::F, n::Integer) where {T, F<:Interval{T}} n == 1 && return a n == 2 && return sqrt(a) + abig = bigequiv(a) if n < 0 issubnormal(mag(a)) && return inv(nthroot(a, -n)) return F( inv(nthroot(abig, -n)) ) end - b = nthroot(bigequiv(a), n) + b = nthroot(abig, n) return F(b) end