From 0ac153fcd3cb1b12394ee0e11b4840c2dbf0d883 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Tue, 8 Mar 2022 16:28:06 -0500 Subject: [PATCH 01/10] don't screw up pow --- base/math.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/math.jl b/base/math.jl index 15bee7f1fccb4..3f46fb8b5c1e2 100644 --- a/base/math.jl +++ b/base/math.jl @@ -1011,7 +1011,7 @@ end x==0 && return abs(y)*Inf*(!(y>0)) logxhi,logxlo = Base.Math._log_ext(x) xyhi = logxhi*y - xylo = logxlo*y + xylo = muladd(logxlo, y, fma(logxhi, y, -xyhi)) hi = xyhi+xylo return Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ)) end From ca007ee2713a468f0d1e522f3fbbdb7c30bbcaa8 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Tue, 8 Mar 2022 16:31:28 -0500 Subject: [PATCH 02/10] better performance for cpus without fma --- base/math.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/base/math.jl b/base/math.jl index 3f46fb8b5c1e2..80c9c9b3a47dd 100644 --- a/base/math.jl +++ b/base/math.jl @@ -1010,8 +1010,8 @@ end !isfinite(x) && return x*(y>0 || isnan(x)) x==0 && return abs(y)*Inf*(!(y>0)) logxhi,logxlo = Base.Math._log_ext(x) - xyhi = logxhi*y - xylo = muladd(logxlo, y, fma(logxhi, y, -xyhi)) + xyhi, xylo = two_mul(logxhi,y) + xylo = muladd(logxlo, y, xylo) hi = xyhi+xylo return Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ)) end From 00e7522a873edb98f3d85385633a622e6448f678 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Wed, 9 Mar 2022 10:52:21 -0500 Subject: [PATCH 03/10] improve accuracy for x^-3 --- base/intfuncs.jl | 1 - base/math.jl | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index 3c2d9b4beec7b..44c7be0626126 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -321,7 +321,6 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}} @inline literal_pow(::typeof(^), x::HWNumber, ::Val{3}) = x*x*x @inline literal_pow(::typeof(^), x::HWNumber, ::Val{-1}) = inv(x) @inline literal_pow(::typeof(^), x::HWNumber, ::Val{-2}) = (i=inv(x); i*i) -@inline literal_pow(::typeof(^), x::HWNumber, ::Val{-3}) = (i=inv(x); i*i*i) # don't use the inv(x) transformation here since float^p is slightly more accurate @inline literal_pow(::typeof(^), x::AbstractFloat, ::Val{p}) where {p} = x^p diff --git a/base/math.jl b/base/math.jl index 80c9c9b3a47dd..d8cd3d7920116 100644 --- a/base/math.jl +++ b/base/math.jl @@ -1041,6 +1041,7 @@ end @assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer) y = 1.0 xnlo = ynlo = 0.0 + n == 3 && return x*x*x # keep compatibility with literal_pow if n < 0 rx = inv(x) n==-2 && return rx*rx #keep compatability with literal_pow @@ -1048,7 +1049,6 @@ end x = rx n = -n end - n == 3 && return x*x*x # keep compatibility with literal_pow while n > 1 if n&1 > 0 err = muladd(y, xnlo, x*ynlo) From c6d120556e4af54fb8179099b2aa2846691a33a7 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Wed, 9 Mar 2022 11:34:51 -0500 Subject: [PATCH 04/10] update doctest (similar accuracy) --- doc/src/manual/complex-and-rational-numbers.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/manual/complex-and-rational-numbers.md b/doc/src/manual/complex-and-rational-numbers.md index 94ad70982bbae..7edbc99932a78 100644 --- a/doc/src/manual/complex-and-rational-numbers.md +++ b/doc/src/manual/complex-and-rational-numbers.md @@ -36,7 +36,7 @@ julia> (-1 + 2im)^2 -3 - 4im julia> (-1 + 2im)^2.5 -2.7296244647840084 - 6.960664459571898im +2.729624464784009 - 6.9606644595719im julia> (-1 + 2im)^(1 + 1im) -0.27910381075826657 + 0.08708053414102428im From 1baab26dc5121122a2a4a54f34886029374adcff Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 9 Mar 2022 14:29:26 -0500 Subject: [PATCH 05/10] fix whitespace --- doc/src/manual/complex-and-rational-numbers.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/manual/complex-and-rational-numbers.md b/doc/src/manual/complex-and-rational-numbers.md index 7edbc99932a78..ac48e5b420f5e 100644 --- a/doc/src/manual/complex-and-rational-numbers.md +++ b/doc/src/manual/complex-and-rational-numbers.md @@ -36,7 +36,7 @@ julia> (-1 + 2im)^2 -3 - 4im julia> (-1 + 2im)^2.5 -2.729624464784009 - 6.9606644595719im +2.729624464784009 - 6.9606644595719im julia> (-1 + 2im)^(1 + 1im) -0.27910381075826657 + 0.08708053414102428im From ae4f0e3241e78e81e4acdbfd28b82572b99f8a7a Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 10 Mar 2022 01:37:00 -0500 Subject: [PATCH 06/10] add tests --- test/math.jl | 29 +++++++++++++++++++++++------ 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/test/math.jl b/test/math.jl index 34a21ca3335ad..d5edb3820200c 100644 --- a/test/math.jl +++ b/test/math.jl @@ -155,12 +155,6 @@ end @test x^y === T(big(x)^big(y)) @test x^1 === x @test x^yi === T(big(x)^yi) - # test (-x)^y for y larger than typemax(Int) - @test T(-1)^floatmax(T) === T(1) - @test prevfloat(T(-1))^floatmax(T) === T(Inf) - @test nextfloat(T(-1))^floatmax(T) === T(0.0) - # test for large negative exponent where error compensation matters - @test 0.9999999955206014^-1.0e8 == 1.565084574870928 @test (-x)^yi == x^yi @test (-x)^(yi+1) == -(x^(yi+1)) @test acos(x) ≈ acos(big(x)) @@ -1323,6 +1317,29 @@ end end end +@testset "pow" begin + for T in (Float16, Float32, Float64) + for x in (0.0, -0.0, 1.0, 10.0, 2.0, Inf, NaN, -Inf, -NaN) + for y in (0.0, -0.0, 1.0, -3.0,-10.0 , Inf, NaN, -Inf, -NaN) + @test T(x)^T(y) === T(big(x)^y) + end + end + for _ in 1:2^16 + x=rand(T)*100; y=rand(T)*200-100 + got, expected = x^y, widen(x)^y + if isfinite(eps(T(expected))) + @test abs(expected-got) <= 1.3*eps(T(expected)) + end + end + end + # test (-x)^y for y larger than typemax(Int) + @test T(-1)^floatmax(T) === T(1) + @test prevfloat(T(-1))^floatmax(T) === T(Inf) + @test nextfloat(T(-1))^floatmax(T) === T(0.0) + # test for large negative exponent where error compensation matters + @test 0.9999999955206014^-1.0e8 == 1.565084574870928 +end + # Test that sqrt behaves correctly and doesn't exhibit fp80 double rounding. # This happened on old glibc versions. # Test case from https://sourceware.org/bugzilla/show_bug.cgi?id=14032. From 1a8373104d51a56c8af7938ec93e49d377b5f41c Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 10 Mar 2022 12:30:37 -0500 Subject: [PATCH 07/10] fix tests --- test/math.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/test/math.jl b/test/math.jl index d5edb3820200c..828494511c7a9 100644 --- a/test/math.jl +++ b/test/math.jl @@ -1321,7 +1321,8 @@ end for T in (Float16, Float32, Float64) for x in (0.0, -0.0, 1.0, 10.0, 2.0, Inf, NaN, -Inf, -NaN) for y in (0.0, -0.0, 1.0, -3.0,-10.0 , Inf, NaN, -Inf, -NaN) - @test T(x)^T(y) === T(big(x)^y) + got, expected = T(x)^T(y), T(big(x))^T(y) + @test isnan_type(T, got) && isnan_type(T, expected) || (got === expected) end end for _ in 1:2^16 @@ -1331,11 +1332,11 @@ end @test abs(expected-got) <= 1.3*eps(T(expected)) end end + # test (-x)^y for y larger than typemax(Int) + @test T(-1)^floatmax(T) === T(1) + @test prevfloat(T(-1))^floatmax(T) === T(Inf) + @test nextfloat(T(-1))^floatmax(T) === T(0.0) end - # test (-x)^y for y larger than typemax(Int) - @test T(-1)^floatmax(T) === T(1) - @test prevfloat(T(-1))^floatmax(T) === T(Inf) - @test nextfloat(T(-1))^floatmax(T) === T(0.0) # test for large negative exponent where error compensation matters @test 0.9999999955206014^-1.0e8 == 1.565084574870928 end From 347d64725bb1afae8407c35e7218d5386af63195 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 10 Mar 2022 14:45:15 -0500 Subject: [PATCH 08/10] fix float32 for -n, test debug --- base/math.jl | 2 +- test/math.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/base/math.jl b/base/math.jl index d8cd3d7920116..6fd6a2f2a8be2 100644 --- a/base/math.jl +++ b/base/math.jl @@ -1065,8 +1065,8 @@ end end function ^(x::Float32, n::Integer) - n < 0 && return inv(x)^(-n) n == 3 && return x*x*x #keep compatibility with literal_pow + n < 0 && return Float32(Base.power_by_squaring(inv(Float64(x)),n)) Float32(Base.power_by_squaring(Float64(x),n)) end @inline ^(x::Float16, y::Integer) = Float16(Float32(x) ^ y) diff --git a/test/math.jl b/test/math.jl index 828494511c7a9..3ae8703c09100 100644 --- a/test/math.jl +++ b/test/math.jl @@ -1329,7 +1329,7 @@ end x=rand(T)*100; y=rand(T)*200-100 got, expected = x^y, widen(x)^y if isfinite(eps(T(expected))) - @test abs(expected-got) <= 1.3*eps(T(expected)) + @test abs(expected-got) <= 1.3*eps(T(expected)) || (x,y) end end # test (-x)^y for y larger than typemax(Int) From 0ac47534c52dcc6b198507712f5b9f9c97ee6c3f Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 11 Mar 2022 10:56:08 -0500 Subject: [PATCH 09/10] typo --- base/math.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/math.jl b/base/math.jl index 6fd6a2f2a8be2..40317d2ddb3a1 100644 --- a/base/math.jl +++ b/base/math.jl @@ -1066,7 +1066,7 @@ end function ^(x::Float32, n::Integer) n == 3 && return x*x*x #keep compatibility with literal_pow - n < 0 && return Float32(Base.power_by_squaring(inv(Float64(x)),n)) + n < 0 && return Float32(Base.power_by_squaring(inv(Float64(x)),-n)) Float32(Base.power_by_squaring(Float64(x),n)) end @inline ^(x::Float16, y::Integer) = Float16(Float32(x) ^ y) From 1a4f0557e0fff1d27f12ebcfa6f5e92fe2c7e1c3 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 11 Mar 2022 15:56:56 -0500 Subject: [PATCH 10/10] bugfix --- base/math.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/base/math.jl b/base/math.jl index 40317d2ddb3a1..0d20c7bf1cca5 100644 --- a/base/math.jl +++ b/base/math.jl @@ -1065,6 +1065,7 @@ end end function ^(x::Float32, n::Integer) + n == -2 && return (i=inv(x); i*i) n == 3 && return x*x*x #keep compatibility with literal_pow n < 0 && return Float32(Base.power_by_squaring(inv(Float64(x)),-n)) Float32(Base.power_by_squaring(Float64(x),n))