Skip to content
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.

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

improve concrete-foldability of core math functions #45613

Merged
merged 3 commits into from
Jun 11, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
28 changes: 14 additions & 14 deletions base/intfuncs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -246,20 +246,20 @@ end

# ^ for any x supporting *
to_power_type(x) = convert(Base._return_type(*, Tuple{typeof(x), typeof(x)}), x)
@noinline throw_domerr_powbysq(::Any, p) = throw(DomainError(p,
string("Cannot raise an integer x to a negative power ", p, '.',
"\nConvert input to float.")))
@noinline throw_domerr_powbysq(::Integer, p) = throw(DomainError(p,
string("Cannot raise an integer x to a negative power ", p, '.',
"\nMake x or $p a float by adding a zero decimal ",
"(e.g., 2.0^$p or 2^$(float(p)) instead of 2^$p), ",
"or write 1/x^$(-p), float(x)^$p, x^float($p) or (x//1)^$p")))
@noinline throw_domerr_powbysq(::AbstractMatrix, p) = throw(DomainError(p,
string("Cannot raise an integer matrix x to a negative power ", p, '.',
"\nMake x a float matrix by adding a zero decimal ",
"(e.g., [2.0 1.0;1.0 0.0]^$p instead ",
"of [2 1;1 0]^$p), or write float(x)^$p or Rational.(x)^$p")))
function power_by_squaring(x_, p::Integer)
@noinline throw_domerr_powbysq(::Any, p) = throw(DomainError(p, LazyString(
"Cannot raise an integer x to a negative power ", p, ".",
"\nConvert input to float.")))
@noinline throw_domerr_powbysq(::Integer, p) = throw(DomainError(p, LazyString(
"Cannot raise an integer x to a negative power ", p, ".",
"\nMake x or ", p, " a float by adding a zero decimal ",
"(e.g., 2.0^", p, " or 2^", float(p), " instead of 2^", p, ")",
"or write 1/x^", -p, ", float(x)^", p, ", x^float(", p, ") or (x//1)^", p, ".")))
@noinline throw_domerr_powbysq(::AbstractMatrix, p) = throw(DomainError(p, LazyString(
"Cannot raise an integer matrix x to a negative power ", p, ".",
"\nMake x a float matrix by adding a zero decimal ",
"(e.g., [2.0 1.0;1.0 0.0]^", p, " instead of [2 1;1 0]^", p, ")",
"or write float(x)^", p, " or Rational.(x)^", p, ".")))
@assume_effects :terminates_locally function power_by_squaring(x_, p::Integer)
x = to_power_type(x_)
if p == 1
return copy(x)
Expand Down
12 changes: 8 additions & 4 deletions base/special/exp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -175,7 +175,11 @@ const J_TABLE = (0x0000000000000000, 0xaac00b1afa5abcbe, 0x9b60163da9fb3335, 0xa
0xa66f0f9c1cb64129, 0x93af252b376bba97, 0xacdf3ac948dd7273, 0x99df50765b6e4540, 0x9faf6632798844f8,
0xa12f7bfdad9cbe13, 0xaeef91d802243c88, 0x874fa7c1819e90d8, 0xacdfbdba3692d513, 0x62efd3c22b8f71f1, 0x74afe9d96b2a23d9)

@inline function table_unpack(ind)
# XXX we want to mark :consistent-cy here so that this function can be concrete-folded,
# because the effect analysis currently can't prove it in the presence of `@inbounds` or
# `:boundscheck`, but still the access to `J_TABLE` is really safe here
Base.@assume_effects :consistent @inline function table_unpack(ind::Int32)
ind = ind & 255 + 1 # 255 == length(J_TABLE) - 1
j = @inbounds J_TABLE[ind]
jU = reinterpret(Float64, JU_CONST | (j&JU_MASK))
jL = reinterpret(Float64, JL_CONST | (j>>8))
Expand Down Expand Up @@ -211,7 +215,7 @@ end
r = muladd(N_float, LogBo256U(base, T), x)
r = muladd(N_float, LogBo256L(base, T), r)
k = N >> 8
jU, jL = table_unpack(N&255 + 1)
jU, jL = table_unpack(N)
aviatesk marked this conversation as resolved.
Show resolved Hide resolved
small_part = muladd(jU, expm1b_kernel(base, r), jL) + jU

if !(abs(x) <= SUBNORM_EXP(base, T))
Expand All @@ -236,7 +240,7 @@ end
r = muladd(N_float, LogBo256U(base, T), x)
r = muladd(N_float, LogBo256L(base, T), r)
k = N >> 8
jU, jL = table_unpack(N&255 + 1)
jU, jL = table_unpack(N)
very_small = muladd(jU, expm1b_kernel(base, r), jL)
small_part = muladd(jU,xlo,very_small) + jU
if !(abs(x) <= SUBNORM_EXP(base, T))
Expand Down Expand Up @@ -439,7 +443,7 @@ function expm1(x::Float64)
r = muladd(N_float, LogBo256U(Val(:ℯ), T), x)
r = muladd(N_float, LogBo256L(Val(:ℯ), T), r)
k = Int64(N >> 8)
jU, jL = table_unpack(N&255 +1)
jU, jL = table_unpack(N)
p = expm1b_kernel(Val(:ℯ), r)
twopk = reinterpret(Float64, (1023+k) << 52)
twopnk = reinterpret(Float64, (1023-k) << 52)
Expand Down
28 changes: 17 additions & 11 deletions base/special/log.jl
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,6 @@ const t_log_Float64 = ((0.0,0.0),(0.007782140442941454,-8.865052917267247e-13),
(0.6853040030982811,6.383161517064652e-13),(0.6892332812385575,2.5144230728376075e-13),
(0.6931471805601177,-1.7239444525614835e-13))


# Float32 lookup table
# to generate values:
# N=16
Expand Down Expand Up @@ -156,7 +155,12 @@ logbU(::Type{Float64},::Val{10}) = 0.4342944819032518
logbL(::Type{Float64},::Val{10}) = 1.098319650216765e-17

# Procedure 1
@inline function log_proc1(y::Float64,mf::Float64,F::Float64,f::Float64,jp::Int,base=Val(:ℯ))
# XXX we want to mark :consistent-cy here so that this function can be concrete-folded,
# because the effect analysis currently can't prove it in the presence of `@inbounds` or
# `:boundscheck`, but still the access to `t_log_Float64` is really safe here
Base.@assume_effects :consistent @inline function log_proc1(y::Float64,mf::Float64,F::Float64,f::Float64,base=Val(:ℯ))
jp = unsafe_trunc(Int,128.0*F)-127

## Steps 1 and 2
@inbounds hi,lo = t_log_Float64[jp]
l_hi = mf* 0.6931471805601177 + hi
Expand Down Expand Up @@ -211,8 +215,13 @@ end
return fma(m_hi, u, fma(m_lo, u, m_hi*fma(fma(-u,f,2(f-u)), g, q)))
end

# Procedure 1
# XXX we want to mark :consistent-cy here so that this function can be concrete-folded,
# because the effect analysis currently can't prove it in the presence of `@inbounds` or
# `:boundscheck`, but still the access to `t_log_Float32` is really safe here
Base.@assume_effects :consistent @inline function log_proc1(y::Float32,mf::Float32,F::Float32,f::Float32,base=Val(:ℯ))
jp = unsafe_trunc(Int,128.0f0*F)-127

@inline function log_proc1(y::Float32,mf::Float32,F::Float32,f::Float32,jp::Int,base=Val(:ℯ))
## Steps 1 and 2
@inbounds hi = t_log_Float32[jp]
l = mf*0.6931471805599453 + hi
Expand All @@ -232,6 +241,7 @@ end
Float32(logb(Float32, base)*(l + (u + q)))
end

# Procedure 2
@inline function log_proc2(f::Float32,base=Val(:ℯ))
## Step 1
# compute in higher precision
Expand Down Expand Up @@ -281,9 +291,8 @@ function _log(x::Float64, base, func)
mf = Float64(m)
F = (y + 3.5184372088832e13) - 3.5184372088832e13 # 0x1p-7*round(0x1p7*y)
f = y-F
jp = unsafe_trunc(Int,128.0*F)-127

return log_proc1(y,mf,F,f,jp,base)
return log_proc1(y,mf,F,f,base)
elseif x == 0.0
-Inf
elseif isnan(x)
Expand Down Expand Up @@ -317,9 +326,8 @@ function _log(x::Float32, base, func)
mf = Float32(m)
F = (y + 65536.0f0) - 65536.0f0 # 0x1p-7*round(0x1p7*y)
f = y-F
jp = unsafe_trunc(Int,128.0f0*F)-127

log_proc1(y,mf,F,f,jp,base)
log_proc1(y,mf,F,f,base)
elseif x == 0f0
-Inf32
elseif isnan(x)
Expand Down Expand Up @@ -352,9 +360,8 @@ function log1p(x::Float64)
mf = Float64(m)
F = (y + 3.5184372088832e13) - 3.5184372088832e13 # 0x1p-7*round(0x1p7*y)
f = (y - F) + c*s #2^m(F+f) = 1+x = z+c
jp = unsafe_trunc(Int,128.0*F)-127

log_proc1(y,mf,F,f,jp)
log_proc1(y,mf,F,f)
elseif x == -1.0
-Inf
elseif isnan(x)
Expand Down Expand Up @@ -385,9 +392,8 @@ function log1p(x::Float32)
mf = Float32(m)
F = (y + 65536.0f0) - 65536.0f0 # 0x1p-7*round(0x1p7*y)
f = (y - F) + s*c #2^m(F+f) = 1+x = z+c
jp = unsafe_trunc(Int,128.0*F)-127

log_proc1(y,mf,F,f,jp)
log_proc1(y,mf,F,f)
elseif x == -1f0
-Inf32
elseif isnan(x)
Expand Down
5 changes: 4 additions & 1 deletion base/special/rem_pio2.jl
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,10 @@ function fromfraction(f::Int128)
return (z1,z2)
end

function paynehanek(x::Float64)
# XXX we want to mark :consistent-cy here so that this function can be concrete-folded,
# because the effect analysis currently can't prove it in the presence of `@inbounds` or
# `:boundscheck`, but still the accesses to `INV_2PI` are really safe here
Base.@assume_effects :consistent function paynehanek(x::Float64)
# 1. Convert to form
#
# x = X * 2^k,
Expand Down
18 changes: 18 additions & 0 deletions test/math.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1449,3 +1449,21 @@ end
f44336()
@test (@allocated f44336()) == 0
end

# test constant-foldability
for fn in (:sin, :cos, :tan, :log, :log2, :log10, :log1p, :exponent, :sqrt, :cbrt,
# TODO :asin, :atan, :acos, :sinh, :cosh, :tanh, :asinh, :acosh, :atanh,
# TODO :exp, :exp2, :exp10, :expm1
)
for T in (Float32, Float64)
f = getfield(@__MODULE__, fn)
eff = Base.infer_effects(f, (T,))
if Core.Compiler.is_foldable(eff)
@test true
else
@error "bad effects found for $f(::$T)" eff
@test false
end
end
end
@test Core.Compiler.is_foldable(Base.infer_effects(^, (Float32,Int)))