Repository: JuliaLang/julia
File: base/math.jl (Julia)

Code before:
@constprop :aggressive function ^(x::Float64, y::Float64)
    xu = reinterpret(UInt64, x)
    xu == reinterpret(UInt64, 1.0) && return 1.0
    # Exponents greater than this will always overflow or underflow.
    # Note that NaN can pass through this, but that will end up fine.
    if !(abs(y)<0x1.8p62)
        isnan(y) && return y
        y = sign(y)*0x1.8p62
    end
    yint = unsafe_trunc(Int64, y) # This is actually safe since julia freezes the result
    yisint = y == yint
    if yisint
        yint == 0 && return 1.0
        use_power_by_squaring(yint) && return @noinline pow_body(x, yint)
    end
    2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x === +0.0 or -0.0 (Inf * false === 0.0)
    s = 1
    if x < 0
        !yisint && throw_exp_domainerror(x) # y isn't an integer
        s = ifelse(isodd(yint), -1, 1)
    end
    !isfinite(x) && return copysign(x,s)*(y>0 || isnan(x))           # x is inf or NaN
    return copysign(pow_body(abs(x), y), s)
end

@assume_effects :foldable @inline function pow_body(x::Float64, y::Float64)
    xu = reinterpret(UInt64, x)
    if xu < (UInt64(1)<<52) # x is subnormal
        xu = reinterpret(UInt64, x * 0x1p52) # normalize x
        xu &= ~sign_mask(Float64)
        xu -= UInt64(52) << 52 # mess with the exponent
    end
    logxhi,logxlo = _log_ext(xu)
    xyhi, xylo = two_mul(logxhi,y)
    xylo = muladd(logxlo, y, xylo)
    hi = xyhi+xylo
    return @inline Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
end

@constprop :aggressive function ^(x::T, y::T) where T <: Union{Float16, Float32}
    x == 1 && return one(T)
    # Exponents greater than this will always overflow or underflow.
    # Note that NaN can pass through this, but that will end up fine.
    max_exp = T == Float16 ? T(3<<14) : T(0x1.Ap30)
    if !(abs(y)<max_exp)
        isnan(y) && return y
        y = sign(y)*max_exp
    end
    yint = unsafe_trunc(Int32, y) # This is actually safe since julia freezes the result
    y == yint && return x^yint
    x < 0 && throw_exp_domainerror(x)

Review comment (suggestion):
I think we don't want this `@inline` but we do want the one on `exp_impl`. (because we want the `^`to inline so the compiler can const prop the `int` conversion etc.

Code after:
@constprop :aggressive function ^(x::Float64, y::Float64)
    xu = reinterpret(UInt64, x)
    xu == reinterpret(UInt64, 1.0) && return 1.0
    # Exponents greater than this will always overflow or underflow.
    # Note that NaN can pass through this, but that will end up fine.
    if !(abs(y)<0x1.8p62)
        isnan(y) && return y
        y = sign(y)*0x1.8p62
    end
    yint = unsafe_trunc(Int64, y) # This is actually safe since julia freezes the result
    yisint = y == yint
    if yisint
        yint == 0 && return 1.0
        use_power_by_squaring(yint) && return @noinline pow_body(x, yint)
    end
    2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x === +0.0 or -0.0 (Inf * false === 0.0)
    s = 1
    if x < 0
        !yisint && throw_exp_domainerror(x) # y isn't an integer
        s = ifelse(isodd(yint), -1, 1)
    end
    !isfinite(x) && return copysign(x,s)*(y>0 || isnan(x))           # x is inf or NaN
    return copysign(pow_body(abs(x), y), s)
end

@assume_effects :foldable @noinline function pow_body(x::Float64, y::Float64)
    xu = reinterpret(UInt64, x)
    if xu < (UInt64(1)<<52) # x is subnormal
        xu = reinterpret(UInt64, x * 0x1p52) # normalize x
        xu &= ~sign_mask(Float64)
        xu -= UInt64(52) << 52 # mess with the exponent
    end
    logxhi,logxlo = _log_ext(xu)
    xyhi, xylo = two_mul(logxhi,y)
    xylo = muladd(logxlo, y, xylo)
    hi = xyhi+xylo
    return @inline Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
end

@constprop :aggressive function ^(x::T, y::T) where T <: Union{Float16, Float32}
    x == 1 && return one(T)
    # Exponents greater than this will always overflow or underflow.
    # Note that NaN can pass through this, but that will end up fine.
    max_exp = T == Float16 ? T(3<<14) : T(0x1.Ap30)
    if !(abs(y)<max_exp)
        isnan(y) && return y
        y = sign(y)*max_exp
    end
    yint = unsafe_trunc(Int32, y) # This is actually safe since julia freezes the result
    y == yint && return x^yint
    x < 0 && throw_exp_domainerror(x)