From 4f0550ab32472136bfc4b447b3c00e42d2f76faa Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 27 Aug 2021 09:37:47 -0500 Subject: [PATCH 01/12] ^(::Float, ::Integer) uses compensated power by squaring for Float64. For Float32, this just using power by squaring with Float64. .5 ULP accuracy for both types. --- base/math.jl | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/base/math.jl b/base/math.jl index 3857b1b1e8c10..e4954b19ef0ca 100644 --- a/base/math.jl +++ b/base/math.jl @@ -926,21 +926,28 @@ end end @inline ^(x::Float16, y::Float16) = Float16(Float32(x)^Float32(y)) # TODO: optimize +# compensated power by squaring @inline function ^(x::Float64, y::Integer) - y == -1 && return inv(x) - y == 0 && return one(x) - y == 1 && return x - y == 2 && return x*x - y == 3 && return x*x*x - ccall("llvm.pow.f64", llvmcall, Float64, (Float64, Float64), x, Float64(y)) + n < 0 && return inv(myintpow(x,-n)) + y = 1.0 + xnlo = ynlo = 0.0 + while n > 1 + if n&1 > 0 + yn = x*y + ynlo = fma(x, y , -yn) + muladd(y, xnlo, x*ynlo) + y = yn + end + xn = x * x + xnlo = muladd(x, 2*xnlo, fma(x, x, -xn)) + x = xn + n >>>= 1 + end + !isfinite(x) && return x*y + return muladd(x, y, muladd(y, xnlo, x*ynlo)) end @inline function ^(x::Float32, y::Integer) - y == -1 && return inv(x) - y == 0 && return one(x) - y == 1 && return x - y == 2 && return x*x - y == 3 && return x*x*x - ccall("llvm.pow.f32", llvmcall, Float32, (Float32, Float32), x, Float32(y)) + y < 0 && return inv(myintpow(x,-y)) + Float32(Base.power_by_squaring(Float64(x),y)) end @inline ^(x::Float16, y::Integer) = Float16(Float32(x) ^ y) @inline literal_pow(::typeof(^), x::Float16, ::Val{p}) where {p} = Float16(literal_pow(^,Float32(x),Val(p))) From 7b87679cb56d4606cff3a54bc4181bceab99844c Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 27 Aug 2021 11:35:26 -0500 Subject: [PATCH 02/12] fix bug --- base/math.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/base/math.jl b/base/math.jl index e4954b19ef0ca..86f082d53afe9 100644 --- a/base/math.jl +++ b/base/math.jl @@ -928,7 +928,7 @@ end # compensated power by squaring @inline function ^(x::Float64, y::Integer) - n < 0 && return inv(myintpow(x,-n)) + n < 0 && return inv(x^(-n)) y = 1.0 xnlo = ynlo = 0.0 while n > 1 @@ -946,7 +946,7 @@ end return muladd(x, y, muladd(y, xnlo, x*ynlo)) end @inline function ^(x::Float32, y::Integer) - y < 0 && return inv(myintpow(x,-y)) + y < 0 && return inv(x^(-y)) Float32(Base.power_by_squaring(Float64(x),y)) end @inline ^(x::Float16, y::Integer) = Float16(Float32(x) ^ y) From c4798cf99531e043f5ec6e2d8a7770933089f196 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 27 Aug 2021 12:51:21 -0500 Subject: [PATCH 03/12] Update base/math.jl Co-authored-by: Kristoffer Carlsson --- base/math.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/math.jl b/base/math.jl index 86f082d53afe9..83d046004c632 100644 --- a/base/math.jl +++ b/base/math.jl @@ -927,7 +927,7 @@ end @inline ^(x::Float16, y::Float16) = Float16(Float32(x)^Float32(y)) # TODO: optimize # compensated power by squaring -@inline function ^(x::Float64, y::Integer) +@inline function ^(x::Float64, n::Integer) n < 0 && return inv(x^(-n)) y = 1.0 xnlo = ynlo = 0.0 From eb75934a74e3e5f2b166d6c434a325afdfa8a042 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 17 Sep 2021 07:52:31 -0500 Subject: [PATCH 04/12] fix negative numbers, and tests --- base/math.jl | 39 +++++++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/base/math.jl b/base/math.jl index 83d046004c632..cd8297dcd9776 100644 --- a/base/math.jl +++ b/base/math.jl @@ -911,6 +911,8 @@ function modf(x::Float64) end @inline function ^(x::Float64, y::Float64) + yint = round(Int, y) + y == yint && return x^yint z = ccall("llvm.pow.f64", llvmcall, Float64, (Float64, Float64), x, y) if isnan(z) & !isnan(x+y) throw_exp_domainerror(x) @@ -918,6 +920,8 @@ end z end @inline function ^(x::Float32, y::Float32) + yint = round(Int, y) + y == yint && return x^yint z = ccall("llvm.pow.f32", llvmcall, Float32, (Float32, Float32), x, y) if isnan(z) & !isnan(x+y) throw_exp_domainerror(x) @@ -928,22 +932,25 @@ end # compensated power by squaring @inline function ^(x::Float64, n::Integer) - n < 0 && return inv(x^(-n)) - y = 1.0 - xnlo = ynlo = 0.0 - while n > 1 - if n&1 > 0 - yn = x*y - ynlo = fma(x, y , -yn) + muladd(y, xnlo, x*ynlo) - y = yn - end - xn = x * x - xnlo = muladd(x, 2*xnlo, fma(x, x, -xn)) - x = xn - n >>>= 1 - end - !isfinite(x) && return x*y - return muladd(x, y, muladd(y, xnlo, x*ynlo)) + if n > 0 + y = 1.0 + xnlo = ynlo = 0.0 + while n > 1 + if n&1 > 0 + yn = x*y + ynlo = fma(x, y , -yn) + muladd(y, xnlo, x*ynlo) + y = yn + end + xn = x * x + xnlo = muladd(x, 2*xnlo, fma(x, x, -xn)) + x = xn + n >>>= 1 + end + !isfinite(x) && return x*y + return muladd(x, y, muladd(y, xnlo, x*ynlo)) + end + n < 0 && return inv(x^(-n)) + return one(n) # n == 0 end @inline function ^(x::Float32, y::Integer) y < 0 && return inv(x^(-y)) From b00e3b0b80a9815ed4386b7894d251f2c349590b Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 17 Sep 2021 09:09:35 -0500 Subject: [PATCH 05/12] fix subnormals --- base/math.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/math.jl b/base/math.jl index cd8297dcd9776..fdd29e446b2db 100644 --- a/base/math.jl +++ b/base/math.jl @@ -949,7 +949,7 @@ end !isfinite(x) && return x*y return muladd(x, y, muladd(y, xnlo, x*ynlo)) end - n < 0 && return inv(x^(-n)) + n < 0 && return inv(x)^(-n) return one(n) # n == 0 end @inline function ^(x::Float32, y::Integer) From 002188984a2844f582425c81f9a4f495e5f081e9 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 17 Sep 2021 09:35:07 -0500 Subject: [PATCH 06/12] fix subnormals --- base/math.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/math.jl b/base/math.jl index fdd29e446b2db..e1a788b2f3d95 100644 --- a/base/math.jl +++ b/base/math.jl @@ -950,7 +950,7 @@ end return muladd(x, y, muladd(y, xnlo, x*ynlo)) end n < 0 && return inv(x)^(-n) - return one(n) # n == 0 + return one(x) # n == 0 end @inline function ^(x::Float32, y::Integer) y < 0 && return inv(x^(-y)) From 1730a3226af2ef95d1e91de4ce34fcbff8daf8ff Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 17 Sep 2021 09:36:19 -0500 Subject: [PATCH 07/12] fix whitespace --- base/math.jl | 46 +++++++++++++++++++++++----------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/base/math.jl b/base/math.jl index e1a788b2f3d95..79a6967e436d2 100644 --- a/base/math.jl +++ b/base/math.jl @@ -911,8 +911,8 @@ function modf(x::Float64) end @inline function ^(x::Float64, y::Float64) - yint = round(Int, y) - y == yint && return x^yint + yint = round(Int, y) + y == yint && return x^yint z = ccall("llvm.pow.f64", llvmcall, Float64, (Float64, Float64), x, y) if isnan(z) & !isnan(x+y) throw_exp_domainerror(x) @@ -920,8 +920,8 @@ end z end @inline function ^(x::Float32, y::Float32) - yint = round(Int, y) - y == yint && return x^yint + yint = round(Int, y) + y == yint && return x^yint z = ccall("llvm.pow.f32", llvmcall, Float32, (Float32, Float32), x, y) if isnan(z) & !isnan(x+y) throw_exp_domainerror(x) @@ -932,25 +932,25 @@ end # compensated power by squaring @inline function ^(x::Float64, n::Integer) - if n > 0 - y = 1.0 - xnlo = ynlo = 0.0 - while n > 1 - if n&1 > 0 - yn = x*y - ynlo = fma(x, y , -yn) + muladd(y, xnlo, x*ynlo) - y = yn - end - xn = x * x - xnlo = muladd(x, 2*xnlo, fma(x, x, -xn)) - x = xn - n >>>= 1 - end - !isfinite(x) && return x*y - return muladd(x, y, muladd(y, xnlo, x*ynlo)) - end - n < 0 && return inv(x)^(-n) - return one(x) # n == 0 + if n > 0 + y = 1.0 + xnlo = ynlo = 0.0 + while n > 1 + if n&1 > 0 + yn = x*y + ynlo = fma(x, y , -yn) + muladd(y, xnlo, x*ynlo) + y = yn + end + xn = x * x + xnlo = muladd(x, 2*xnlo, fma(x, x, -xn)) + x = xn + n >>>= 1 + end + !isfinite(x) && return x*y + return muladd(x, y, muladd(y, xnlo, x*ynlo)) + end + n < 0 && return inv(x)^(-n) + return one(x) # n == 0 end @inline function ^(x::Float32, y::Integer) y < 0 && return inv(x^(-y)) From ec107090994ebc44b5102ed25047f8cbb9d73565 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 17 Sep 2021 10:06:40 -0500 Subject: [PATCH 08/12] fix infinity --- base/math.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/base/math.jl b/base/math.jl index 79a6967e436d2..9362e66414019 100644 --- a/base/math.jl +++ b/base/math.jl @@ -911,8 +911,10 @@ function modf(x::Float64) end @inline function ^(x::Float64, y::Float64) - yint = round(Int, y) - y == yint && return x^yint + if isfinite(x) + yint = round(Int, y) + y == yint && return x^yint + end z = ccall("llvm.pow.f64", llvmcall, Float64, (Float64, Float64), x, y) if isnan(z) & !isnan(x+y) throw_exp_domainerror(x) @@ -920,8 +922,10 @@ end z end @inline function ^(x::Float32, y::Float32) - yint = round(Int, y) - y == yint && return x^yint + if isfinite(x) + yint = round(Int, y) + y == yint && return x^yint + end z = ccall("llvm.pow.f32", llvmcall, Float32, (Float32, Float32), x, y) if isnan(z) & !isnan(x+y) throw_exp_domainerror(x) From 804b71ccbcb230301bcfd3313dded144062e8fca Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 17 Sep 2021 10:36:42 -0500 Subject: [PATCH 09/12] fix big powers --- base/math.jl | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/base/math.jl b/base/math.jl index 9362e66414019..d881775d7c347 100644 --- a/base/math.jl +++ b/base/math.jl @@ -911,10 +911,8 @@ function modf(x::Float64) end @inline function ^(x::Float64, y::Float64) - if isfinite(x) - yint = round(Int, y) - y == yint && return x^yint - end + yint = unsafe_trunc(Int, y) + y == yint && return x^yint z = ccall("llvm.pow.f64", llvmcall, Float64, (Float64, Float64), x, y) if isnan(z) & !isnan(x+y) throw_exp_domainerror(x) @@ -922,10 +920,8 @@ end z end @inline function ^(x::Float32, y::Float32) - if isfinite(x) - yint = round(Int, y) - y == yint && return x^yint - end + yint = unsafe_trunc(Int, y) + y == yint && return x^yint z = ccall("llvm.pow.f32", llvmcall, Float32, (Float32, Float32), x, y) if isnan(z) & !isnan(x+y) throw_exp_domainerror(x) From a6ca2f4ac8cc7b06c048bf522067322ae0c36714 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 30 Sep 2021 01:11:08 -0500 Subject: [PATCH 10/12] fixes --- base/math.jl | 43 ++++++++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/base/math.jl b/base/math.jl index d881775d7c347..3e71f7744050f 100644 --- a/base/math.jl +++ b/base/math.jl @@ -911,7 +911,7 @@ function modf(x::Float64) end @inline function ^(x::Float64, y::Float64) - yint = unsafe_trunc(Int, y) + yint = unsafe_trunc(Int, y) # Note, this is actually safe since julia freezes the result y == yint && return x^yint z = ccall("llvm.pow.f64", llvmcall, Float64, (Float64, Float64), x, y) if isnan(z) & !isnan(x+y) @@ -920,7 +920,7 @@ end z end @inline function ^(x::Float32, y::Float32) - yint = unsafe_trunc(Int, y) + yint = unsafe_trunc(Int, y) # Note, this is actually safe since julia freezes the result y == yint && return x^yint z = ccall("llvm.pow.f32", llvmcall, Float32, (Float32, Float32), x, y) if isnan(z) & !isnan(x+y) @@ -932,25 +932,30 @@ end # compensated power by squaring @inline function ^(x::Float64, n::Integer) - if n > 0 - y = 1.0 - xnlo = ynlo = 0.0 - while n > 1 - if n&1 > 0 - yn = x*y - ynlo = fma(x, y , -yn) + muladd(y, xnlo, x*ynlo) - y = yn - end - xn = x * x - xnlo = muladd(x, 2*xnlo, fma(x, x, -xn)) - x = xn - n >>>= 1 + n == 0 && return one(x) + y = 1.0 + xnlo = ynlo = 0.0 + + if n < 0 + rx = inv(x) + x = rx + xnlo = fma(x, rx, -1.) * rx + n = -n + end + n==3 && return x*x*x #keep compatability with literal_pow + while n > 1 + if n&1 > 0 + yn = x*y + ynlo = fma(x, y , -yn) + muladd(y, xnlo, x*ynlo) + y = yn end - !isfinite(x) && return x*y - return muladd(x, y, muladd(y, xnlo, x*ynlo)) + xn = x * x + xnlo = muladd(x, 2*xnlo, fma(x, x, -xn)) + x = xn + n >>>= 1 end - n < 0 && return inv(x)^(-n) - return one(x) # n == 0 + !isfinite(x) && return x*y + return muladd(x, y, muladd(y, xnlo, x*ynlo)) end @inline function ^(x::Float32, y::Integer) y < 0 && return inv(x^(-y)) From 125330acfcc0c481b97db1b32dd6a35b728486fe Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 30 Sep 2021 23:13:44 -0500 Subject: [PATCH 11/12] fix failing tests --- base/math.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/base/math.jl b/base/math.jl index 3e71f7744050f..0b6adc996bc81 100644 --- a/base/math.jl +++ b/base/math.jl @@ -935,11 +935,10 @@ end n == 0 && return one(x) y = 1.0 xnlo = ynlo = 0.0 - if n < 0 rx = inv(x) - x = rx xnlo = fma(x, rx, -1.) * rx + x = rx n = -n end n==3 && return x*x*x #keep compatability with literal_pow @@ -957,9 +956,10 @@ end !isfinite(x) && return x*y return muladd(x, y, muladd(y, xnlo, x*ynlo)) end -@inline function ^(x::Float32, y::Integer) - y < 0 && return inv(x^(-y)) - Float32(Base.power_by_squaring(Float64(x),y)) +@inline function ^(x::Float32, n::Integer) + n < 0 && return inv(x)^(-n) + n==3 && return x*x*x #keep compatability with literal_pow + Float32(Base.power_by_squaring(Float64(x),n)) end @inline ^(x::Float16, y::Integer) = Float16(Float32(x) ^ y) @inline literal_pow(::typeof(^), x::Float16, ::Val{p}) where {p} = Float16(literal_pow(^,Float32(x),Val(p))) From 965a2c3eb5bcfd2b3edcfe345040f4111f00bb52 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 1 Oct 2021 01:15:43 -0500 Subject: [PATCH 12/12] fix Inf^-1 --- base/math.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/math.jl b/base/math.jl index 0b6adc996bc81..4a88dd4e3a33c 100644 --- a/base/math.jl +++ b/base/math.jl @@ -937,7 +937,7 @@ end xnlo = ynlo = 0.0 if n < 0 rx = inv(x) - xnlo = fma(x, rx, -1.) * rx + isfinite(x) && (xnlo = fma(x, rx, -1.) * rx) x = rx n = -n end