# This file is part of the IntervalArithmetic.jl package; MIT licensed

## Powers
for f in (:unsafe_pow_inner,)
    @eval begin
        ($f)(::EmptyInterval{T}, ::S) where {T, S <: Real} =
            emptyinterval(promote_type(T, S))
        ($f)(::T, ::EmptyInterval{S}) where {T <: Real, S} =
            emptyinterval(promote_type(T, S))
        ($f)(::ValidInterval{T}, ::ValidInterval{S}) where {T, S} =
            emptyinterval(promote_type(T, S))
    end
end

# Integer power:

# overwrite new behaviour for small integer powers from
# https://github.com/JuliaLang/julia/pull/24240:

Base.literal_pow(::typeof(^), x::Interval{T}, ::Val{p}) where {T,p} =
    ^(x, p)

^(x::EmptyInterval, ::Integer) = x

function ^(a::Interval{T}, n::Integer) where T <: AbstractFloat
    n < 0 && return a == zero(a) ? emptyinterval(a) : inv(a^(-n))
    n == 0 && return one(a)
    n == 1 && return a
    a == zero(a) && return a
    n == 2 && return sqr(a)
    # we put an arbitrary limit to recursion at n == 10
    if isodd(n)
        n > 10 && return a^unsafe_atomic(T, n)
        return a * a^(n - 1)
    else
        n > 10 && return nonnegatives(T) ∩ a^unsafe_atomic(T, n)
        return nonnegatives(T) ∩ (a * a^(n - 1))
    end
end

# input a must be positive
sqr(a::AbstractFloat, mode::RoundingMode) = *(a, a, mode)

sqr(a::Interval{T}) where T = zero_in(a) ?
    unsafe_interval(zero(T), @round_up(sqr(mag(a)))) :
    @round(sqr(mig(a)), sqr(mag(a)))

# Interval power of an interval:
function ^(a::ValidInterval{T}, x::ValidInterval{T}) where T
    domain = nonnegatives(T)
    return unsafe_pow_inner(a ∩ domain, x)
end

function unsafe_pow_inner(a::Interval{T}, x::Interval{T}) where T
    a == zero(a) && return !isempty(x ∩ nonnegatives(T)) ? a : emptyinterval(T)
    return exp(x * log(a))
end

function sqrt(a::ValidInterval{T}) where T
    domain = nonnegatives(T)
    return sqrt_inner(a ∩ domain)
end
sqrt_inner(a::Interval) = @round(sqrt(a.lo), sqrt(a.hi))
sqrt_inner(x::EmptyInterval)= x

cbrt(a::Interval{BigFloat}) = @round( cbrt(a.lo), cbrt(a.hi) )
cbrt(x::EmptyInterval)= x

hypot(a::Interval{T}, b::Interval{T}) where T = sqrt(sqr(a) + sqr(b))
