From 64286a442ece29eb50708db814d83a4cc240af8a Mon Sep 17 00:00:00 2001 From: Tim Esselens Date: Mon, 28 Sep 2026 18:18:57 +0200 Subject: [PATCH 1/2] Keep the sign in fast integer powers of negative Float32/Float16 numbers `@fastmath x^n` with an integer `n` outside -1:3 was computed as `__nv_fast_powf(x, Float32(n))`, that is exp2(n * log2(x)), which is NaN for x < 0: `(-0.5f0)^7` gave NaN where the CPU gives -0.0078125. This also affected literal exponents, which reach the same method through `pow_fast(x, ::Val)`. Float64 goes through `__nv_pow` and was not affected. Take the power of `abs(x)` and copy the sign of `x` for odd `n`, which is how LLVM's AMDGPU back-end expands a fast `pown`. The parity comes from the integer exponent, so an odd exponent that Float32 rounds to an even number (16_777_217) keeps its sign. The Float16 method now calls the Float32 one with the integer exponent, so it gets the same handling. Found while differentiating WaterLily.jl kernels with Enzyme on a GTX 1080. Co-Authored-By: Claude Opus 5.5 (1M context) --- CUDACore/src/device/intrinsics/math.jl | 8 ++++++-- test/core/device/intrinsics/math.jl | 15 +++++++++++++++ 2 files changed, 21 insertions(+), 2 deletions(-) diff --git a/CUDACore/src/device/intrinsics/math.jl b/CUDACore/src/device/intrinsics/math.jl index 7ec66b4d9b..1509bf571a 100644 --- a/CUDACore/src/device/intrinsics/math.jl +++ b/CUDACore/src/device/intrinsics/math.jl @@ -405,7 +405,11 @@ end y == 1 && return x y == 2 && return x*x y == 3 && return x*x*x - FastMath.pow_fast(x, Float32(y)) # uses __nv_fast_powf + # __nv_fast_powf computes exp2(y * log2(x)), which is NaN for x < 0, so take the power of + # abs(x) and copy the sign of x for odd y, like LLVM's AMDGPU expansion of fast `pown`. + # the parity comes from the integer y, since Float32(y) may round to an even number. + r = FastMath.pow_fast(abs(x), Float32(y)) # uses __nv_fast_powf + ifelse(isodd(y), copysign(r, x), r) end @device_override @assume_effects :foldable @inline function FastMath.pow_fast(x::Float16, y::Integer) y == -1 && return inv(x) @@ -413,7 +417,7 @@ end y == 1 && return x y == 2 && return x*x y == 3 && return x*x*x - Float16(FastMath.pow_fast(Float32(x), Float32(y))) + Float16(FastMath.pow_fast(Float32(x), y)) end # stay with Base's semantics and avoid the drift of libdevice's integer powers, diff --git a/test/core/device/intrinsics/math.jl b/test/core/device/intrinsics/math.jl index 1e5f1f28f0..ebaf7d327c 100644 --- a/test/core/device/intrinsics/math.jl +++ b/test/core/device/intrinsics/math.jl @@ -324,6 +324,21 @@ using SpecialFunctions @cuda threads=4 fastpow_kernel(A, 3) @test Array(A) == ones(T, 4) end + # larger exponents use `__nv_fast_powf`, which is NaN for negative bases + for T in (Float16, Float32, Float64), y in (4, 5, 7, -2, -3, Int32(5), Int32(-5)) + x = T[1.5, -0.5, 3, -2] + A = CuArray(x) + @cuda threads=4 fastpow_kernel(A, y) + @test all(isapprox.(Array(A), x .^ y; rtol=16*eps(T))) + end + # literal exponents take the same path, through `pow_fast(x, ::Val)` + fastpow7(x) = @fastmath x^7 + x = Float32[-0.5, -2] + @test all(isapprox.(Array(map(fastpow7, cu(x))), x .^ 7; rtol=16*eps(Float32))) + # the sign comes from the integer exponent: Float32(16_777_217) is even + A = CuArray(Float32[-1]) + @cuda threads=1 fastpow_kernel(A, 16_777_217) + @test Array(A) == Float32[-1] # Float16 hardware approximations: tanh.approx.f16 / ex2.approx.f16 on sm_75+ if capability(device()) >= v"7.5" From d28489203c4298b6aeda0c741afe49cbd1df0199 Mon Sep 17 00:00:00 2001 From: Tim Besard Date: Mon, 28 Sep 2026 18:18:57 +0200 Subject: [PATCH 2/2] Take the sign of precise integer powers from the integer exponent `^(::Float32, ::Integer)` converts the exponent and calls `__nv_powf`, which handles negative bases, but the conversion can round an odd exponent to an even one: `(-1f0)^16_777_217` gave 1. The same holds for Float64 beyond 2^53. Take the sign from the integer exponent, as the fast variant now does. Co-Authored-By: Claude Opus 5.5 (1M context) --- CUDACore/src/device/intrinsics/math.jl | 7 +++++-- test/core/device/intrinsics/math.jl | 7 +++++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/CUDACore/src/device/intrinsics/math.jl b/CUDACore/src/device/intrinsics/math.jl index 1509bf571a..0c6f143088 100644 --- a/CUDACore/src/device/intrinsics/math.jl +++ b/CUDACore/src/device/intrinsics/math.jl @@ -427,13 +427,15 @@ end # Base's `^(::Float, ::Integer)` calls `power_by_squaring`, whose # `trailing_zeros` loop emits `cttz_int`, so convert and defer to `__nv_pow`. +# the float exponent may round to an even number, so take the sign from y. @device_override @assume_effects :foldable @inline function Base.:(^)(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 - x ^ Float32(y) + r = abs(x) ^ Float32(y) + ifelse(isodd(y), copysign(r, x), r) end @device_override @assume_effects :foldable @inline function Base.:(^)(x::Float64, y::Integer) y == -1 && return inv(x) @@ -441,7 +443,8 @@ end y == 1 && return x y == 2 && return x*x y == 3 && return x*x*x - x ^ Float64(y) + r = abs(x) ^ Float64(y) + ifelse(isodd(y), copysign(r, x), r) end ## rounding and selection diff --git a/test/core/device/intrinsics/math.jl b/test/core/device/intrinsics/math.jl index ebaf7d327c..0038c0fcd7 100644 --- a/test/core/device/intrinsics/math.jl +++ b/test/core/device/intrinsics/math.jl @@ -19,9 +19,12 @@ using SpecialFunctions # ~1500 ulp on the cases below, while Int64 went via __nv_pow. @testset "integer exponent ($T, $I)" for T in (Float32, Float64), I in (Int32, Int64) - x = T[1.001, 0.9, 1.1, 1.0001, 2, 0.5] - n = I[500, 200, 30, 5000, 7, -3] + x = T[1.001, 0.9, 1.1, 1.0001, 2, 0.5, -2] + n = I[500, 200, 30, 5000, 7, -3, 5] @test Array(CuArray(x) .^ CuArray(n)) ≈ x .^ n rtol=8*eps(T) + # the sign comes from the integer exponent: Float32(16_777_217) is even + # (not compared with the CPU, which gets this wrong on Julia 1.13.0) + @test Array(CuArray(T[-1]) .^ CuArray(I[16_777_217])) == T[-1] end # Both widths have to agree with each other, not just with the CPU