Merge pull request #45613 from JuliaLang/avi/concretemath

improve concrete-foldability of core math functions

Author: aviatesk

base/intfuncs.jl

@@ -248,20 +248,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)

base/special/exp.jl

@@ -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))
@@ -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)
     small_part =  muladd(jU, expm1b_kernel(base, r), jL) + jU
 
     if !(abs(x) <= SUBNORM_EXP(base, T))
@@ -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))
@@ -438,7 +442,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)

base/special/log.jl

@@ -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
@@ -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
@@ -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
@@ -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
@@ -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)
@@ -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)
@@ -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)
@@ -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)

base/special/rem_pio2.jl

@@ -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,

test/math.jl

@@ -1425,3 +1425,21 @@ end
 for T in (Float32, Float64)
     @test Core.Compiler.is_foldable(Base.infer_effects(^, (T,Int)))
 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)))