diff --git a/LICENSES/berkeley-softfloat-BSD-3-Clause.txt b/LICENSES/berkeley-softfloat-BSD-3-Clause.txt new file mode 100644 index 000000000..7ab5488e6 --- /dev/null +++ b/LICENSES/berkeley-softfloat-BSD-3-Clause.txt @@ -0,0 +1,42 @@ +The Berkeley SoftFloat lookup tables (APPROX_RECIP_SQRT_K0S/K1S, +APPROX_RECIP_K0S/K1S) and the recip / recipSqrt approximation +functions in dist/softfloat64.metal and shaders/softfloat.metal are +derived from John R. Hauser's Berkeley SoftFloat-3e: + + Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017 + The Regents of the University of California. + All Rights Reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions, and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above + copyright notice, this list of conditions, and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + 3. Neither the name of the University nor the names of its + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, + THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + PARTICULAR PURPOSE, ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF + USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT + OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + SUCH DAMAGE. + +The vendored Berkeley SoftFloat-3e and TestFloat-3e sources under +vendor/ are also covered by this license. See +vendor/SoftFloat-3e/COPYING.txt and vendor/TestFloat-3e/COPYING.txt +for the unmodified upstream license texts. diff --git a/LICENSES/julia-MIT.txt b/LICENSES/julia-MIT.txt new file mode 100644 index 000000000..9823029ed --- /dev/null +++ b/LICENSES/julia-MIT.txt @@ -0,0 +1,22 @@ +MIT License + +Copyright (c) 2009-2024: Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors: https://github.com/JuliaLang/julia/contributors + +Permission is hereby granted, free of charge, to any person obtaining +a copy of this software and associated documentation files (the +"Software"), to deal in the Software without restriction, including +without limitation the rights to use, copy, modify, merge, publish, +distribute, sublicense, and/or sell copies of the Software, and to +permit persons to whom the Software is furnished to do so, subject to +the following conditions: + +The above copyright notice and this permission notice shall be +included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE +LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/LICENSES/metal-softfloat-MIT.txt b/LICENSES/metal-softfloat-MIT.txt new file mode 100644 index 000000000..a8ac0c7a3 --- /dev/null +++ b/LICENSES/metal-softfloat-MIT.txt @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 metal-softfloat contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md index d3fde856d..66489d3a2 100644 --- a/README.md +++ b/README.md @@ -53,6 +53,11 @@ identifies cached kernels. Compiled kernels are not invalidated when a library c Julia library methods are recompiled when redefined, but cached kernels keep using the previous version, and bitcode files are read once per session. +Software binary64 support is an optional library; see +[`src/softfloat/README.md`](src/softfloat/README.md). Its legalization runs at this same +boundary for both the caller and Julia libraries to keep late-imported runtime ABIs +consistent. + ### Testing GPUCompiler.jl's own test suite runs on GitHub Actions. In addition, Buildkite runs the test diff --git a/src/GPUCompiler.jl b/src/GPUCompiler.jl index 8a5903978..03ca66ea7 100644 --- a/src/GPUCompiler.jl +++ b/src/GPUCompiler.jl @@ -84,6 +84,7 @@ include("devicelib.jl") include("mcgen.jl") include("debug.jl") include("driver.jl") +include("softfloat/SoftFloat.jl") # other reusable functionality include("execution.jl") diff --git a/src/softfloat/README.md b/src/softfloat/README.md new file mode 100644 index 000000000..20bdaac61 --- /dev/null +++ b/src/softfloat/README.md @@ -0,0 +1,72 @@ +# Software floating point + +This directory implements IEEE 754 binary64 (`Float64`) support for GPU targets without +native double-precision arithmetic, using only 64-bit integer operations. Back-ends opt in +by returning `softfloat_libraries()...` from `GPUCompiler.device_libraries(job)` and calling +`GPUCompiler.SoftFloat.legalize!(job, mod)` in `lower_intrinsics!`. +The same hook legalizes Julia libraries imported after optimization, so imported runtime functions +and their callers agree on integer ABIs. External bitcode is an ABI boundary: unsupported +binary64 calls are rejected rather than silently changing their calling convention. +The optimized integer routines are linked at the late library stage, from two libraries: +small routines (sign operations, comparisons, rounding, most conversions) are forced inline, +while the arithmetic and narrowing conversions remain out of line to avoid duplicating +their implementation at every operation (see `OUTLINED` in `provider.jl`). + +The emulation operates on ordinary `Float64` code in the final LLVM module, including +Base's pure-Julia math functions. The usual GPU restrictions still apply: external +floating-point library calls need backend overrides, and unsupported LLVM operations +are rejected. `legalize.jl` first outlines every floating-point operation on +`double` values into a call to one of the integer-only routines, and then rebuilds the +module with `double` replaced by `i64` everywhere (signatures, aggregates, constants, +globals, attributes). The routines are linked in afterwards, through the device-library +machinery in `../devicelib.jl`, from the table in `provider.jl`. + +The routines themselves are the methods of `SoftFloat64`, a `primitive type` with the bit +representation of `Float64` that can be used as an ordinary number type on the CPU as well +(which is how `test/softfloat.jl` validates it, bit for bit against native arithmetic). +`binary64/types.jl` defines the type, using Base's format traits (`significand_bits`, +`exponent_mask`, ...), and `unpack`, which yields the sign, exponent and normalized +significand of a value. Every operation in `arithmetic.jl` and `conversions.jl` then +computes an exact or sticky-extended significand for its result and hands it to +`round_pack` in `rounding.jl`, which is where rounding happens: a single `round_up` decision +implements all `RoundingMode`s, for `round_pack`, the narrowing conversions to `Float32` and +`Float16`, and `round` to integral values. The operators default to round-to-nearest-even +but accept a trailing rounding mode, e.g. `+(a, b, RoundUp)`, which specializes at compile +time. `uint128.jl` provides the two-word `U128` that significand products need, since many +back-ends cannot lower `UInt128`. + +The large routines (`+`, `*`, `/`, `sqrt`, `fma`) are `@noinline` so that they are shared +between call sites; the others are meant to be inlined. + +Semantics, as used by the legalizer: round-to-nearest-even, gradual underflow, signed zeros, +infinities, quiet comparisons, and canonical quiet NaN arithmetic results +(`0x7ff8000000000000`). Sign operations preserve NaN payloads. Not supported: exception +flags, `Float64` atomics and `frem` +(Julia's `rem` does not use it). `Float64`→`Float16` rounds directly rather than through +`Float32`, which would double-round. Widening `Float16` and `Float32` uses their bits +directly, preserving subnormals even when native arithmetic flushes them to zero. + +`paynehanek.jl` is unrelated to the emulation itself: it is a version of Base's Payne-Hanek +argument reduction (used by `sin`, `cos`, and `tan` for large arguments) that avoids the +`UInt128` arithmetic many GPU back-ends cannot lower. Back-ends should override +`Base.Math.paynehanek(::Float64)` with it. + +## Package boundary + +The numerical core is `uint128.jl` and `binary64/`; it depends only on Base and can move +into a standalone package. `provider.jl` and `legalize.jl` are GPUCompiler integration, +while `paynehanek.jl` is a backend workaround for Base's implementation. + +`SoftFloat64` currently supplies the arithmetic and conversions needed by that integration, +plus basic numeric traits, promotion and hashing. A standalone package should complete +and test the remaining number interface (checked integer conversions, adjacent values, +parsing and broader promotion) before presenting it as a general replacement for `Float64`. + +## Provenance + +The algorithms in `binary64/` are those of +[Berkeley SoftFloat 3e](https://github.com/ucb-bar/berkeley-softfloat-3) (BSD-3-Clause), +including its reciprocal tables, by way of +[metal-softfloat](https://github.com/guyfischman/metal-softfloat) +(commit `8b6c592e2e383040fe2778bed8dda7904df284b1`, MIT). `paynehanek.jl` is adapted from +Julia's `base/special/rem_pio2.jl` (MIT). The notices are in `../../LICENSES/`. diff --git a/src/softfloat/SoftFloat.jl b/src/softfloat/SoftFloat.jl new file mode 100644 index 000000000..e13d22295 --- /dev/null +++ b/src/softfloat/SoftFloat.jl @@ -0,0 +1,22 @@ +# software floating-point support for targets without native double precision + +module SoftFloat + +using LLVM, LLVM.IR, LLVM.Build, LLVM.Passes +import ..GPUCompiler +using Base: sign_mask, exponent_mask, significand_mask, exponent_one, significand_bits, + exponent_bits, exponent_bias, exponent_raw_max, uinttype + +include("uint128.jl") +include("binary64/types.jl") +include("binary64/rounding.jl") +include("binary64/arithmetic.jl") +include("binary64/comparisons.jl") +include("binary64/conversions.jl") +include("paynehanek.jl") +include("provider.jl") +include("legalize.jl") + +export SoftFloat64, softfloat_libraries + +end diff --git a/src/softfloat/binary64/arithmetic.jl b/src/softfloat/binary64/arithmetic.jl new file mode 100644 index 000000000..ab95a69e8 --- /dev/null +++ b/src/softfloat/binary64/arithmetic.jl @@ -0,0 +1,283 @@ +# Binary64 arithmetic: addition, multiplication, division, square root and fused +# multiply-add, following Berkeley SoftFloat 3e (BSD-3-Clause; see +# LICENSES/berkeley-softfloat-BSD-3-Clause.txt), by way of its port to Metal in +# metal-softfloat (MIT; see LICENSES/metal-softfloat-MIT.txt). +# +# Every operation first handles the special values, then computes a significand for the +# result that is either exact or extended with a sticky bit, and hands it to `round_pack`. +# The operators without a rounding-mode argument are the device library's entry points: +# they are `@noinline` so that the generated device code shares them between call sites. + +Base.:-(x::SoftFloat64) = reinterpret(SoftFloat64, bits(x) ⊻ sign_mask(SoftFloat64)) +Base.abs(x::SoftFloat64) = reinterpret(SoftFloat64, bits(x) & ~sign_mask(SoftFloat64)) +Base.copysign(x::SoftFloat64, y::SoftFloat64) = + reinterpret(SoftFloat64, (bits(x) & ~sign_mask(SoftFloat64)) | (bits(y) & sign_mask(SoftFloat64))) + + +## addition + +Base.@noinline Base.:+(a::SoftFloat64, b::SoftFloat64) = +(a, b, RoundNearest) +Base.:-(a::SoftFloat64, b::SoftFloat64) = a + -b +Base.:-(a::SoftFloat64, b::SoftFloat64, mode::RoundingMode) = +(a, -b, mode) + +@inline function Base.:+(a::SoftFloat64, b::SoftFloat64, mode::RoundingMode) + if !isfinite(a) | !isfinite(b) + (isnan(a) | isnan(b)) && return quiet_nan(SoftFloat64) + (isinf(a) & isinf(b) & (signbit(a) != signbit(b))) && return quiet_nan(SoftFloat64) + return isinf(a) ? a : b + end + if iszero(a) | iszero(b) + (iszero(a) & iszero(b)) && return signbit(a) == signbit(b) ? a : cancellation_zero(mode) + return iszero(a) ? b : a + end + + sign, exp, sig = unpack(a) + sign_b, exp_b, sig_b = unpack(b) + # order the operands by magnitude, so that only the smaller one needs to be aligned + if (exp, sig) < (exp_b, sig_b) + sign, exp, sig, sign_b, exp_b, sig_b = sign_b, exp_b, sig_b, sign, exp, sig + end + # ten bits below the significands hold the rounding bits, and the sticky bit from aligning + sig <<= 10 + sig_b = shift_right_jam(sig_b << 10, exp - exp_b) + if sign == sign_b + sum = sig + sig_b + if sum >= one(UInt64) << 63 + sum = shift_right_jam(sum, 1) + exp += 1 + end + return round_pack(SoftFloat64, sign, exp, sum, mode) + else + diff = sig - sig_b + iszero(diff) && return cancellation_zero(mode) + # the leading bit of the difference can be anywhere; its lowest bit has weight + # 2^(exp - 62), the significand's leading bit having been at 62 + return norm_round_pack(SoftFloat64, sign, exp - 62, diff, mode) + end +end + + +## multiplication + +Base.@noinline Base.:*(a::SoftFloat64, b::SoftFloat64) = *(a, b, RoundNearest) + +@inline function Base.:*(a::SoftFloat64, b::SoftFloat64, mode::RoundingMode) + sign = signbit(a) ⊻ signbit(b) + if !isfinite(a) | !isfinite(b) + (isnan(a) | isnan(b) | iszero(a) | iszero(b)) && return quiet_nan(SoftFloat64) + return infinity(SoftFloat64, sign) + end + (iszero(a) | iszero(b)) && return signed_zero(SoftFloat64, sign) + + _, exp_a, sig_a = unpack(a) + _, exp_b, sig_b = unpack(b) + # the top word of the 106-bit product has its leading bit at 62 or 61, at weight + # 2^(exp_a + exp_b + 1); the low word only contributes to the sticky bit + sig = sticky_high(widemul(sig_a << 10, sig_b << 11)) + exp = exp_a + exp_b + 1 + if sig < one(UInt64) << 62 + sig <<= 1 + exp -= 1 + end + round_pack(SoftFloat64, sign, exp, sig, mode) +end + + +## division + +# Approximate reciprocals from the tables of Berkeley SoftFloat 3e: `a` is a fixed-point +# number with its leading bit set, with one integer bit (`1 <= a < 2`), or two for the +# square root when `two_integer_bits` is set (`2 <= a < 4`). The result is a pure fraction +# that never exceeds the true reciprocal, to within about 2 ulps. + +const APPROX_RECIP_K0S = (0xffc4, 0xf0be, 0xe363, 0xd76f, 0xccad, 0xc2f0, 0xba16, 0xb201, + 0xaa97, 0xa3c6, 0x9d7a, 0x97a6, 0x923c, 0x8d32, 0x887e, 0x8417) +const APPROX_RECIP_K1S = (0xf0f1, 0xd62c, 0xbfa1, 0xac77, 0x9c0a, 0x8ddb, 0x8185, 0x76ba, + 0x6d3b, 0x64d4, 0x5d5c, 0x56b1, 0x50b6, 0x4b55, 0x4679, 0x4211) + +@inline function approx_recip32(a::UInt32) + index = (a >> 27) & 0xf + eps = (a >> 11) % UInt16 + @inbounds r0 = APPROX_RECIP_K0S[index + 1] - + ((UInt32(APPROX_RECIP_K1S[index + 1]) * eps) >> 20) % UInt16 + sigma0 = ~((widemul(UInt32(r0), a) >> 7) % UInt32) + r = (UInt32(r0) << 16) + (widemul(UInt32(r0), sigma0) >> 24) % UInt32 + sqr_sigma0 = (widemul(sigma0, sigma0) >> 32) % UInt32 + r + (widemul(r, sqr_sigma0) >> 48) % UInt32 +end + +const APPROX_RECIP_SQRT_K0S = (0xb4c9, 0xffab, 0xaa7d, 0xf11c, 0xa1c5, 0xe4c7, 0x9a43, 0xda29, + 0x93b5, 0xd0e5, 0x8ded, 0xc8b7, 0x88c6, 0xc16d, 0x8424, 0xbae1) +const APPROX_RECIP_SQRT_K1S = (0xa5a5, 0xea42, 0x8c21, 0xc62d, 0x788f, 0xaa7f, 0x6928, 0x94b6, + 0x5cc7, 0x8335, 0x52a6, 0x74e2, 0x4a3e, 0x68fe, 0x432b, 0x5efd) + +@inline function approx_recip_sqrt32(a::UInt32, two_integer_bits::Bool) + index = ((a >> 27) & 0xe) + !two_integer_bits + eps = (a >> 12) % UInt16 + @inbounds r0 = APPROX_RECIP_SQRT_K0S[index + 1] - + ((UInt32(APPROX_RECIP_SQRT_K1S[index + 1]) * eps) >> 20) % UInt16 + e_sqr_r0 = UInt32(r0) * UInt32(r0) + two_integer_bits && (e_sqr_r0 <<= 1) + sigma0 = ~((widemul(e_sqr_r0, a) >> 23) % UInt32) + r = (UInt32(r0) << 16) + (widemul(UInt32(r0), sigma0) >> 25) % UInt32 + sqr_sigma0 = (widemul(sigma0, sigma0) >> 32) % UInt32 + r += (widemul((r >> 1) + (r >> 3) - (UInt32(r0) << 14), sqr_sigma0) >> 48) % UInt32 + r & 0x8000_0000 == 0 ? UInt32(0x8000_0000) : r +end + +Base.@noinline Base.:/(a::SoftFloat64, b::SoftFloat64) = /(a, b, RoundNearest) + +@inline function Base.:/(a::SoftFloat64, b::SoftFloat64, mode::RoundingMode) + sign = signbit(a) ⊻ signbit(b) + if !isfinite(a) | !isfinite(b) | iszero(a) | iszero(b) + (isnan(a) | isnan(b)) && return quiet_nan(SoftFloat64) + ((isinf(a) & isinf(b)) | (iszero(a) & iszero(b))) && return quiet_nan(SoftFloat64) + (isinf(a) | iszero(b)) && return infinity(SoftFloat64, sign) + return signed_zero(SoftFloat64, sign) + end + + _, exp_a, sig_a = unpack(a) + _, exp_b, sig_b = unpack(b) + # scale the dividend such that the quotient's leading bit lands at 62 + exp = exp_a - exp_b + if sig_a < sig_b + exp -= 1 + sig_a <<= 11 + else + sig_a <<= 10 + end + sig_b <<= 11 + + # the upper half of the quotient from an approximate reciprocal of the divisor, then the + # lower half from the remainder, which is finally used to correct the last few bits and + # to set the sticky bit + recip = approx_recip32((sig_b >> 32) % UInt32) - 0x2 + sig_b_hi = (sig_b >> 32) % UInt32 + sig_b_lo = (sig_b % UInt32) >> 4 + q_hi = (widemul((sig_a >> 32) % UInt32, recip) >> 32) % UInt32 + double_term = q_hi << 1 + rem = ((sig_a - widemul(double_term, sig_b_hi)) << 28) - widemul(double_term, sig_b_lo) + q_lo = (widemul((rem >> 32) % UInt32, recip) >> 32) % UInt32 + 0x4 + sig = (UInt64(q_hi) << 32) + (UInt64(q_lo) << 4) + if sig & 0x1ff < 0x40 + q_lo &= ~UInt32(7) + sig &= ~UInt64(0x7f) + double_term = q_lo << 1 + rem = ((rem - widemul(double_term, sig_b_hi)) << 28) - widemul(double_term, sig_b_lo) + if rem & sign_mask(SoftFloat64) != 0 + sig -= 0x80 + elseif rem != 0 + sig |= 1 + end + end + round_pack(SoftFloat64, sign, exp, sig, mode) +end + + +## square root + +Base.@noinline Base.sqrt(x::SoftFloat64) = sqrt(x, RoundNearest) + +@inline function Base.sqrt(x::SoftFloat64, mode::RoundingMode) + if !isfinite(x) | iszero(x) | signbit(x) + (isnan(x) | (signbit(x) & !iszero(x))) && return quiet_nan(SoftFloat64) + return x # ±0 and +Inf + end + + _, exp_a, sig_a = unpack(x) + # halving an odd exponent leaves a factor of two in the significand + exp = exp_a >> 1 + odd_exp = isodd(exp_a) + + # the upper half of the root from an approximate reciprocal square root, then the lower + # half from the remainder, which is finally used to correct the last few bits and to + # set the sticky bit + sig32_a = (sig_a >> 21) % UInt32 + recip = approx_recip_sqrt32(sig32_a, odd_exp) + sig32_z = (widemul(sig32_a, recip) >> 32) % UInt32 + if odd_exp + sig_a <<= 9 + else + sig_a <<= 8 + sig32_z >>= 1 + end + rem = sig_a - widemul(sig32_z, sig32_z) + q = (widemul((rem >> 2) % UInt32, recip) >> 32) % UInt32 + sig = ((UInt64(sig32_z) << 32) | (UInt64(1) << 5)) + (UInt64(q) << 3) + if sig & 0x1ff < 0x22 + sig &= ~UInt64(0x3f) + shifted = sig >> 6 + rem = (sig_a << 52) - shifted * shifted + if rem & sign_mask(SoftFloat64) != 0 + sig -= 1 + elseif rem != 0 + sig |= 1 + end + end + round_pack(SoftFloat64, false, exp, sig, mode) +end + + +## fused multiply-add + +Base.@noinline Base.fma(a::SoftFloat64, b::SoftFloat64, c::SoftFloat64) = + fma(a, b, c, RoundNearest) + +@inline function Base.fma(a::SoftFloat64, b::SoftFloat64, c::SoftFloat64, mode::RoundingMode) + sign = signbit(a) ⊻ signbit(b) + if !isfinite(a) | !isfinite(b) | !isfinite(c) + (isnan(a) | isnan(b) | isnan(c)) && return quiet_nan(SoftFloat64) + if isinf(a) | isinf(b) + (iszero(a) | iszero(b)) && return quiet_nan(SoftFloat64) + (isinf(c) & (signbit(c) != sign)) && return quiet_nan(SoftFloat64) + return infinity(SoftFloat64, sign) + end + return c + end + if iszero(a) | iszero(b) + # the product is an exact zero + (iszero(c) & (signbit(c) != sign)) && return cancellation_zero(mode) + return c + end + + _, exp_a, sig_a = unpack(a) + _, exp_b, sig_b = unpack(b) + # the exact product, normalized to have its leading bit at 125, at weight 2^exp + product = widemul(sig_a << 10, sig_b << 10) + exp = exp_a + exp_b + 1 + if product.hi < one(UInt64) << 61 + product <<= 1 + exp -= 1 + end + iszero(c) && return round_pack(SoftFloat64, sign, exp, sticky_high(product) << 1, mode) + + # align the addend's significand to the product's, shifting whichever is smaller + sign_c, exp_c, sig_c = unpack(c) + addend = U128(sig_c << 9, 0) + if exp >= exp_c + addend = shift_right_jam(addend, exp - exp_c) + else + product = shift_right_jam(product, exp_c - exp) + exp = exp_c + end + + if sign == sign_c + sig = sticky_high(product + addend) # leading bit at 62 or 61 + if sig >= one(UInt64) << 62 + exp += 1 + else + sig <<= 1 + end + return round_pack(SoftFloat64, sign, exp, sig, mode) + else + # the larger magnitude determines the sign; cancellation can move the leading bit + # of the difference anywhere, so normalize from its lowest bit, at weight 2^(exp - 125) + if product < addend + sign = sign_c + product, addend = addend, product + end + diff = product - addend + iszero(diff) && return cancellation_zero(mode) + return norm_round_pack(SoftFloat64, sign, exp - 125, diff, mode) + end +end diff --git a/src/softfloat/binary64/comparisons.jl b/src/softfloat/binary64/comparisons.jl new file mode 100644 index 000000000..9e6fda1ee --- /dev/null +++ b/src/softfloat/binary64/comparisons.jl @@ -0,0 +1,46 @@ +# Binary64 comparisons and extrema, on the bit patterns. + +Base.:(==)(a::SoftFloat64, b::SoftFloat64) = + !(isnan(a) | isnan(b)) & ((bits(a) == bits(b)) | (iszero(a) & iszero(b))) + +function Base.:<(a::SoftFloat64, b::SoftFloat64) + (isnan(a) | isnan(b)) && return false + sign_a, sign_b = signbit(a), signbit(b) + if sign_a != sign_b + sign_a & !(iszero(a) & iszero(b)) + else + # the magnitudes order like their bit patterns, in reverse for negative values + (bits(a) != bits(b)) & (sign_a ⊻ (bits(a) < bits(b))) + end +end + +function Base.:<=(a::SoftFloat64, b::SoftFloat64) + (isnan(a) | isnan(b)) && return false + sign_a, sign_b = signbit(a), signbit(b) + if sign_a != sign_b + sign_a | (iszero(a) & iszero(b)) + else + (bits(a) == bits(b)) | (sign_a ⊻ (bits(a) < bits(b))) + end +end + +unordered(a::SoftFloat64, b::SoftFloat64) = isnan(a) | isnan(b) + +# `min` and `max` propagate NaNs and order -0.0 below +0.0, like Base's for Float64 and +# LLVM's `minimum` and `maximum` intrinsics +function Base.min(a::SoftFloat64, b::SoftFloat64) + unordered(a, b) && return quiet_nan(SoftFloat64) + (iszero(a) & iszero(b)) && return reinterpret(SoftFloat64, bits(a) | bits(b)) + a < b ? a : b +end +function Base.max(a::SoftFloat64, b::SoftFloat64) + unordered(a, b) && return quiet_nan(SoftFloat64) + (iszero(a) & iszero(b)) && return reinterpret(SoftFloat64, bits(a) & bits(b)) + a < b ? b : a +end + +# LLVM's `minnum` and `maxnum` intrinsics return the number when only one operand is NaN +minnum(a::SoftFloat64, b::SoftFloat64) = isnan(a) ? (isnan(b) ? quiet_nan(SoftFloat64) : b) : + isnan(b) ? a : min(a, b) +maxnum(a::SoftFloat64, b::SoftFloat64) = isnan(a) ? (isnan(b) ? quiet_nan(SoftFloat64) : b) : + isnan(b) ? a : max(a, b) diff --git a/src/softfloat/binary64/conversions.jl b/src/softfloat/binary64/conversions.jl new file mode 100644 index 000000000..e41456bfc --- /dev/null +++ b/src/softfloat/binary64/conversions.jl @@ -0,0 +1,100 @@ +# Conversions between binary64 and integers, integral values, and the narrower floating-point +# formats, which reuse `round_pack` for their rounding. + + +## integers + +@inline SoftFloat64(x::Int64, mode::RoundingMode=RoundNearest) = + # `abs(typemin(Int64))` wraps, but its reinterpretation is the right magnitude + iszero(x) ? zero(SoftFloat64) : + norm_round_pack(SoftFloat64, x < 0, 0, reinterpret(UInt64, abs(x)), mode) +@inline SoftFloat64(x::UInt64, mode::RoundingMode=RoundNearest) = + iszero(x) ? zero(SoftFloat64) : norm_round_pack(SoftFloat64, false, 0, x, mode) +SoftFloat64(x::Integer, mode::RoundingMode=RoundNearest) = + SoftFloat64((x isa Unsigned ? UInt64 : Int64)(x), mode) + +# The truncating conversions implement LLVM's `fptosi` and `fptoui`, which are only +# defined for values that fit; other values saturate here, but that must not be relied on. +function Base.unsafe_trunc(::Type{Int64}, x::SoftFloat64) + (isnan(x) | iszero(x)) && return Int64(0) + isinf(x) && return signbit(x) ? typemin(Int64) : typemax(Int64) + sign, exp, sig = unpack(x) + exp < 0 && return Int64(0) + exp >= 63 && return sign ? typemin(Int64) : typemax(Int64) + magnitude = sig << (exp - significand_bits(SoftFloat64)) # a negative shift is to the right + sign ? -(magnitude % Int64) : magnitude % Int64 +end +function Base.unsafe_trunc(::Type{UInt64}, x::SoftFloat64) + (isnan(x) | iszero(x) | signbit(x)) && return UInt64(0) + isinf(x) && return typemax(UInt64) + _, exp, sig = unpack(x) + exp < 0 && return UInt64(0) + exp >= 64 && return typemax(UInt64) + sig << (exp - significand_bits(SoftFloat64)) +end + + +## rounding to an integral value + +# (per mode, as Base has fallbacks for AbstractFloat with some of them, which a method +# for all modes would be ambiguous with) +for mode in (:Nearest, :NearestTiesAway, :NearestTiesUp, :ToZero, :Up, :Down, :FromZero) + @eval Base.round(x::SoftFloat64, mode::RoundingMode{$(QuoteNode(mode))}) = + round_integral(x, mode) +end +Base.round(x::SoftFloat64) = round(x, RoundNearest) + +function round_integral(x::SoftFloat64, mode::RoundingMode) + isnan(x) && return quiet_nan(SoftFloat64) + iszero(x) && return x + u = bits(x) + sign = signbit(x) + exp = ((u & exponent_mask(SoftFloat64)) >> significand_bits(SoftFloat64)) % Int - + exponent_bias(SoftFloat64) + exp >= significand_bits(SoftFloat64) && return x # integral already, or infinite + if exp < 0 + # a magnitude below one rounds to zero or one; its leading bit is the half bit when + # the exponent is -1, and any lower bit is a rest bit + up = round_up(mode, sign, false, exp == -1, + (exp < -1) | (u & significand_mask(SoftFloat64) != 0)) + return up ? copysign(one(SoftFloat64), x) : signed_zero(SoftFloat64, sign) + end + # the weight of one in the bit pattern: adding it to the integral part increments the + # value, carrying into the exponent when required + unit = one(UInt64) << (significand_bits(SoftFloat64) - exp) + fraction = u & (unit - 1) + iszero(fraction) && return x + up = round_up(mode, sign, u & unit != 0, fraction & (unit >> 1) != 0, + fraction & ((unit >> 1) - 1) != 0) + reinterpret(SoftFloat64, (u - fraction) + up * unit) +end + + +## narrower floating-point formats + +@inline function SoftFloat64(x::T) where {T<:Union{Float16,Float32}} + # Native comparisons and intermediate floating-point conversions may flush subnormals + # to zero on the target. Classify and widen the original bits instead. + magnitude = bits(x) & ~sign_mask(T) + sign = bits(x) & sign_mask(T) != 0 + magnitude > exponent_mask(T) && return quiet_nan(SoftFloat64) + magnitude == exponent_mask(T) && return infinity(SoftFloat64, sign) + iszero(magnitude) && return signed_zero(SoftFloat64, sign) + _, exp, sig = unpack(x) + # exact, so the rounding bits are zero + round_pack(SoftFloat64, sign, exp, UInt64(sig) << (62 - significand_bits(T)), + RoundNearest) +end + +@inline function narrow(::Type{T}, x::SoftFloat64, mode::RoundingMode) where {T<:Union{Float16,Float32}} + isnan(x) && return quiet_nan(T) + isinf(x) && return infinity(T, signbit(x)) + iszero(x) && return signed_zero(T, signbit(x)) + sign, exp, sig = unpack(x) + # move the leading bit to where `round_pack` expects it, jamming the bits that fall off + shift = significand_bits(SoftFloat64) - (8sizeof(T) - 2) + round_pack(T, sign, exp, shift_right_jam(sig, shift) % uinttype(T), mode) +end +@inline Base.Float32(x::SoftFloat64, mode::RoundingMode=RoundNearest) = narrow(Float32, x, mode) +# Float16 is rounded directly: rounding through Float32 could round twice +@inline Base.Float16(x::SoftFloat64, mode::RoundingMode=RoundNearest) = narrow(Float16, x, mode) diff --git a/src/softfloat/binary64/rounding.jl b/src/softfloat/binary64/rounding.jl new file mode 100644 index 000000000..634c3d54a --- /dev/null +++ b/src/softfloat/binary64/rounding.jl @@ -0,0 +1,84 @@ +# Rounding: the decision shared by all operations, parameterized on the rounding mode, and +# the rounding and packing of results. + +# rounding-mode predicates (Base.Rounding has similar ones from Julia 1.11) +rounds_to_nearest(::RoundingMode) = false +rounds_to_nearest(::RoundingMode{:Nearest}) = true +rounds_to_nearest(::RoundingMode{:NearestTiesAway}) = true +rounds_to_nearest(::RoundingMode{:NearestTiesUp}) = true +round_tie(::RoundingMode{:Nearest}, sign::Bool, odd::Bool) = odd +round_tie(::RoundingMode{:NearestTiesAway}, sign::Bool, odd::Bool) = true +round_tie(::RoundingMode{:NearestTiesUp}, sign::Bool, odd::Bool) = !sign +rounds_away_from_zero(::RoundingMode{:ToZero}, sign::Bool) = false +rounds_away_from_zero(::RoundingMode{:FromZero}, sign::Bool) = true +rounds_away_from_zero(::RoundingMode{:Up}, sign::Bool) = !sign +rounds_away_from_zero(::RoundingMode{:Down}, sign::Bool) = sign +overflows_to_infinity(mode::RoundingMode, sign::Bool) = + rounds_to_nearest(mode) || rounds_away_from_zero(mode, sign) + +# the zero produced by an exact cancellation, whose sign depends on the rounding mode +cancellation_zero(::RoundingMode) = signed_zero(SoftFloat64, false) +cancellation_zero(::RoundingMode{:Down}) = signed_zero(SoftFloat64, true) + +""" + round_up(mode, sign, odd, half, rest) -> Bool + +Whether rounding a truncated magnitude must increment its lowest kept bit, given that bit +(`odd`), the highest discarded bit (`half`), and whether any lower discarded bit is set +(`rest`). +""" +@inline function round_up(mode::RoundingMode, sign::Bool, odd::Bool, half::Bool, rest::Bool) + if rounds_to_nearest(mode) + half & (rest | round_tie(mode, sign, odd)) + else + rounds_away_from_zero(mode, sign) & (half | rest) + end +end + +""" + round_pack(T, sign, exp, sig, mode) -> T + +Round `±sig * 2^(exp - (8sizeof(T) - 2))` to format `T`: `sig::uinttype(T)` has its leading +bit at position `8sizeof(T) - 2` (so `exp` is the binary exponent of the value), followed by +the significand, and then `exponent_bits(T) - 1` rounding bits, the lowest of which must be +sticky, i.e. set whenever any lower-order bit of the exact result was dropped. +""" +@inline function round_pack(::Type{T}, sign::Bool, exp::Int, sig::U, + mode::RoundingMode) where {T, U<:Unsigned} + nbits = 8sizeof(U) + nround = nbits - 2 - significand_bits(T) + biased = exp + exponent_bias(T) + if biased < 1 + # subnormal: scale to the subnormal grid before rounding, to avoid rounding twice + sig = shift_right_jam(sig, 1 - biased) + biased = 1 + end + odd = sig & (one(U) << nround) != 0 + half = sig & (one(U) << (nround - 1)) != 0 + rest = sig & ((one(U) << (nround - 1)) - one(U)) != 0 + sig = (sig >> nround) + round_up(mode, sign, odd, half, rest) + + # The rounded significand still contains its leading bit, which adds one to the exponent + # field when packed. A subnormal result, whose leading bit is lower, thus leaves the field + # at zero, and a carry out of the rounding increments the field, as required. + field = biased - 1 + (sig >> significand_bits(T)) % Int + if field >= exponent_raw_max(T) + return overflows_to_infinity(mode, sign) ? infinity(T, sign) : largest_finite(T, sign) + end + reinterpret(T, (U(sign) << (nbits - 1)) + (((biased - 1) % U) << significand_bits(T)) + sig) +end + +""" + norm_round_pack(T, sign, exp, sig, mode) -> T + +Like `round_pack`, for a nonzero `sig` (a `UInt64` or `U128`) whose leading bit can be +anywhere, and whose lowest bit has weight `2^exp`. +""" +@inline function norm_round_pack(::Type{T}, sign::Bool, exp::Int, sig, + mode::RoundingMode) where {T} + shift = leading_zeros(sig) - 1 + sig = shift >= 0 ? sig << shift : shift_right_jam(sig, -shift) + round_pack(T, sign, exp + 8sizeof(sig) - 2 - shift, narrow(sig), mode) +end +@inline narrow(sig::UInt64) = sig +@inline narrow(sig::U128) = sticky_high(sig) diff --git a/src/softfloat/binary64/types.jl b/src/softfloat/binary64/types.jl new file mode 100644 index 000000000..e9e30f87e --- /dev/null +++ b/src/softfloat/binary64/types.jl @@ -0,0 +1,93 @@ +# The SoftFloat64 type, and the representation shared by all operations on it. + +""" + SoftFloat64 <: AbstractFloat + +An IEEE 754 binary64 value, represented like `Float64`, whose arithmetic, comparisons and +conversions are implemented with 64-bit integer operations only, for targets without native +double-precision support. `reinterpret` converts to and from `Float64` for free. + +Operations round to nearest, ties to even, unless a trailing `RoundingMode` is passed: the +arithmetic operators, `sqrt`, `fma`, `round`, and the conversions from integers and to +narrower floating-point types accept one. Arithmetic and rounding canonicalize NaN results +(`0x7ff8000000000000`). +Bit-preserving operations (`reinterpret`, negation, `abs`, and `copysign`) retain NaN +payloads. There are no exception flags. +""" +primitive type SoftFloat64 <: AbstractFloat 64 end + +SoftFloat64(x::Float64) = reinterpret(SoftFloat64, x) +Base.Float64(x::SoftFloat64) = reinterpret(Float64, x) +Base.show(io::IO, x::SoftFloat64) = print(io, "SoftFloat64(", Float64(x), ")") + +# the format traits of Float64 apply +for f in (:sign_mask, :exponent_mask, :significand_mask, :exponent_one, :exponent_half, + :significand_bits, :exponent_bits, :exponent_bias, :exponent_max, :exponent_raw_max) + @eval Base.$f(::Type{SoftFloat64}) = Base.$f(Float64) +end +Base.uinttype(::Type{SoftFloat64}) = UInt64 +Base.precision(::Type{SoftFloat64}; base::Integer=2) = precision(Float64; base) +Base.eps(::Type{SoftFloat64}) = SoftFloat64(eps(Float64)) +Base.eps(x::SoftFloat64) = abs(x - reinterpret(SoftFloat64, bits(x) ⊻ UInt64(1))) + +# Keep mixed arithmetic in software. Larger integer and arbitrary-precision types need +# their own conversion and promotion policies before they can participate. +Base.promote_rule(::Type{SoftFloat64}, ::Type{T}) where + {T<:Union{Bool,Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Float16,Float32,Float64}} = SoftFloat64 + +# The helpers below apply to any IEEE format `T` these traits describe, i.e. SoftFloat64 +# and the native Float16 and Float32, so that conversions between them can share code. + +bits(x::T) where {T} = reinterpret(uinttype(T), x) + +Base.signbit(x::SoftFloat64) = reinterpret(Int64, x) < 0 +Base.isnan(x::SoftFloat64) = bits(x) & ~sign_mask(SoftFloat64) > exponent_mask(SoftFloat64) +Base.isinf(x::SoftFloat64) = bits(x) & ~sign_mask(SoftFloat64) == exponent_mask(SoftFloat64) +Base.isfinite(x::SoftFloat64) = bits(x) & exponent_mask(SoftFloat64) != exponent_mask(SoftFloat64) +Base.iszero(x::SoftFloat64) = bits(x) & ~sign_mask(SoftFloat64) == 0 +Base.issubnormal(x::SoftFloat64) = (bits(x) & exponent_mask(SoftFloat64) == 0) & !iszero(x) + +# special values of format `T` +signed_zero(::Type{T}, sign::Bool) where {T} = + reinterpret(T, uinttype(T)(sign) << (8sizeof(T) - 1)) +infinity(::Type{T}, sign::Bool) where {T} = + reinterpret(T, bits(signed_zero(T, sign)) | exponent_mask(T)) +largest_finite(::Type{T}, sign::Bool) where {T} = + reinterpret(T, bits(infinity(T, sign)) - one(uinttype(T))) +quiet_nan(::Type{T}) where {T} = + reinterpret(T, exponent_mask(T) | (one(uinttype(T)) << (significand_bits(T) - 1))) + +Base.zero(::Type{SoftFloat64}) = signed_zero(SoftFloat64, false) +Base.one(::Type{SoftFloat64}) = reinterpret(SoftFloat64, exponent_one(SoftFloat64)) +Base.typemin(::Type{SoftFloat64}) = infinity(SoftFloat64, true) +Base.typemax(::Type{SoftFloat64}) = infinity(SoftFloat64, false) +Base.floatmax(::Type{SoftFloat64}) = largest_finite(SoftFloat64, false) +Base.floatmin(::Type{SoftFloat64}) = reinterpret(SoftFloat64, significand_mask(SoftFloat64) + 1) + +""" + unpack(x) -> (sign, exp, sig) + +Decompose a finite, nonzero value into its sign, binary exponent and integral significand, +such that `x == ±sig * 2^(exp - significand_bits(T))`: the significand's implicit leading +bit is made explicit, and subnormals are normalized. +""" +@inline function unpack(x::T) where {T} + u = bits(x) + sign = u & sign_mask(T) != 0 + exp = ((u & exponent_mask(T)) >> significand_bits(T)) % Int + sig = u & significand_mask(T) + if exp == 0 # subnormal + shift = leading_zeros(sig) - exponent_bits(T) + return sign, 1 - exponent_bias(T) - shift, sig << shift + end + return sign, exp - exponent_bias(T), sig | (significand_mask(T) + one(u)) +end + +# Base's numeric hashing uses an exact numerator * 2^exponent / denominator decomposition. +function Base.decompose(x::SoftFloat64) + isnan(x) && return Int64(0), 0, 0 + isinf(x) && return Int64(signbit(x) ? -1 : 1), 0, 0 + iszero(x) && return Int64(0), 0, signbit(x) ? -1 : 1 + sign, exp, sig = unpack(x) + return Int64(sig), exp - significand_bits(SoftFloat64), sign ? -1 : 1 +end diff --git a/src/softfloat/legalize.jl b/src/softfloat/legalize.jl new file mode 100644 index 000000000..74e026d8b --- /dev/null +++ b/src/softfloat/legalize.jl @@ -0,0 +1,679 @@ +# LLVM legalization of binary64 operations and values +# +# This happens in two steps. First, every floating-point operation on `double` values is +# outlined into a call to a placeholder function that still has the original `double` +# signature, e.g. `fadd double %a, %b` becomes `call double @gpu_softfloat_add64(double %a, +# double %b)`. Then, the module is rebuilt with every `double` type replaced by `i64` +# (including in aggregates, globals, constants, signatures and attributes), which turns the +# placeholders into declarations of the integer-only routines that are linked afterwards. +# +# Splitting the semantic and representational changes keeps each step simple: outlining +# only looks at individual instructions, and the rebuild is a type-preserving clone that +# does not need to understand any operation. + +mutable struct Float64TypeRemapper + cache::Dict{LLVM.LLVMType,LLVM.LLVMType} + suffix::String +end +Float64TypeRemapper(; suffix=".softfloat64") = + Float64TypeRemapper(Dict{LLVM.LLVMType,LLVM.LLVMType}(), suffix) + +function remap_type(remapper::Float64TypeRemapper, typ::LLVM.LLVMType) + haskey(remapper.cache, typ) && return remapper.cache[typ] + if typ == LLVM.DoubleType() + mapped = LLVM.Int64Type() + elseif typ isa LLVM.ArrayType + mapped = LLVM.ArrayType(remap_type(remapper, typ.element_type), typ.length) + elseif typ isa LLVM.VectorType + mapped = LLVM.VectorType(remap_type(remapper, typ.element_type), typ.length) + elseif typ isa LLVM.FunctionType + mapped = LLVM.FunctionType(remap_type(remapper, LLVM.return_type(typ)), + LLVM.LLVMType[remap_type(remapper, t) for t in LLVM.parameters(typ)]; + vararg=LLVM.isvararg(typ)) + elseif typ isa LLVM.PointerType + mapped = LLVM.isopaque(typ) ? typ : + LLVM.PointerType(remap_type(remapper, typ.element_type), LLVM.addrspace(typ)) + elseif typ isa LLVM.StructType + if LLVM.isopaque(typ) + mapped = typ + elseif LLVM.name(typ) === nothing + mapped = LLVM.StructType( + LLVM.LLVMType[remap_type(remapper, t) for t in LLVM.elements(typ)]; + packed=LLVM.ispacked(typ)) + else + mapped = LLVM.StructType(LLVM.name(typ) * remapper.suffix) + remapper.cache[typ] = mapped + LLVM.elements!(mapped, + LLVM.LLVMType[remap_type(remapper, t) for t in LLVM.elements(typ)]; + packed=LLVM.ispacked(typ)) + end + else + mapped = typ + end + remapper.cache[typ] = mapped + mapped +end + +function contains_double(typ::LLVM.LLVMType, seen=Set{LLVM.LLVMType}()) + typ == LLVM.DoubleType() && return true + typ in seen && return false + push!(seen, typ) + if typ isa Union{LLVM.ArrayType,LLVM.VectorType} + return contains_double(typ.element_type, seen) + elseif typ isa LLVM.FunctionType + return contains_double(LLVM.return_type(typ), seen) || + any(t -> contains_double(t, seen), LLVM.parameters(typ)) + elseif typ isa LLVM.PointerType + return !LLVM.isopaque(typ) && contains_double(typ.element_type, seen) + elseif typ isa LLVM.StructType && !LLVM.isopaque(typ) + return any(t -> contains_double(t, seen), LLVM.elements(typ)) + end + false +end + +function placeholder!(mod::LLVM.Module, symbol::String, return_type::LLVM.LLVMType, + argument_types::Vector{<:LLVM.LLVMType}) + ft = LLVM.FunctionType(return_type, argument_types) + if haskey(LLVM.functions(mod), symbol) + fn = LLVM.functions(mod)[symbol] + LLVM.function_type(fn) == ft || error("soft-float placeholder ABI mismatch for $symbol") + return fn + end + fun = LLVM.Function(mod, symbol, ft) + fun.memory_effects = LLVM.MemoryEffects() + push!(fun.function_attributes, LLVM.EnumAttribute("nounwind")) + push!(fun.function_attributes, LLVM.EnumAttribute("willreturn")) + fun +end + +function copy_semantic_metadata!(new::LLVM.Instruction, old::LLVM.Instruction) + # LLVM.jl cannot enumerate instruction metadata keys, so copy the standardized kinds + # relevant to arithmetic and memory semantics explicitly. Debug location is installed by + # the IRBuilder before construction. + # `fpmath` describes accuracy of a floating-point result and is invalid after + # the call's result becomes an integer bit pattern. + for key in ("prof", "range", "tbaa", "tbaa.struct", "alias.scope", + "noalias", "nontemporal", "invariant.group") + kind = LLVM.MDKind(key) + haskey(LLVM.metadata(old), kind) || continue + LLVM.metadata(new)[kind] = LLVM.metadata(old)[kind] + end + new +end + +# Julia does not attach a location to every instruction (e.g. some `fcmp`s), but a call +# needs one when both functions have debug info: inlining the routine through a call +# without one leaves its locations in the routine's subprogram, which the verifier rejects. +function call_location(inst::LLVM.Instruction) + loc = inst.debug_location + loc === nothing || return loc + subprogram = LLVM.parent(LLVM.parent(inst)).subprogram + subprogram === nothing ? nothing : LLVM.DILocation(0, 0, subprogram) +end + +function outlined_call!(builder::LLVM.IRBuilder, inst::LLVM.Instruction, symbol::String, + args::Vector{<:LLVM.Value}, return_type::LLVM.LLVMType) + mod = LLVM.parent(LLVM.parent(LLVM.parent(inst))) + fn = placeholder!(mod, symbol, return_type, LLVM.LLVMType[LLVM.value_type(x) for x in args]) + LLVM.position!(builder, LLVM.before(inst)) + builder.debug_location = call_location(inst) + call = LLVM.call!(builder, LLVM.function_type(fn), fn, args, LLVM.name(inst)) + copy_semantic_metadata!(call, inst) + call +end + +# apply `f` per lane when `result_type` is a vector, extracting the lanes of vector `args` +function scalarize!(f, builder, result_type, args...) + result_type isa LLVM.VectorType || return f(args...) + result = LLVM.UndefValue(result_type) + for lane in 0:result_type.length-1 + index = LLVM.ConstantInt(Int32(lane)) + lane_args = map(args) do arg + LLVM.value_type(arg) isa LLVM.VectorType ? + LLVM.extract_element!(builder, arg, index) : arg + end + result = LLVM.insert_element!(builder, result, f(lane_args...), index) + end + result +end + +scalar_type(typ) = typ isa LLVM.VectorType ? typ.element_type : typ +is_binary64_type(typ) = scalar_type(typ) == LLVM.DoubleType() + +function comparison_call!(builder, inst, f, a, b) + # Julia's device ABI represents `Bool` returns as zero-extended i8, while LLVM `fcmp` + # produces i1. Make that boundary explicit; calling an i8 helper through an i1 prototype + # happens to verify with opaque pointers but returns incorrect predicates on Metal. + raw = outlined_call!(builder, inst, helper_name(f), LLVM.Value[a, b], LLVM.Int8Type()) + LLVM.icmp!(builder, LLVM.API.LLVMIntNE, raw, LLVM.ConstantInt(UInt8(0))) +end + +# lower an `fcmp` predicate to the eq/lt/le/unordered helpers +function comparison_value!(builder, inst, pred, a, b) + P = LLVM.API + eq() = comparison_call!(builder, inst, :eq64, a, b) + lt() = comparison_call!(builder, inst, :lt64, a, b) + le() = comparison_call!(builder, inst, :le64, a, b) + gt() = comparison_call!(builder, inst, :lt64, b, a) + ge() = comparison_call!(builder, inst, :le64, b, a) + uno() = comparison_call!(builder, inst, :unordered64, a, b) + if pred == P.LLVMRealPredicateFalse + LLVM.ConstantInt(false) + elseif pred == P.LLVMRealPredicateTrue + LLVM.ConstantInt(true) + elseif pred == P.LLVMRealOEQ + eq() + elseif pred == P.LLVMRealOGT + gt() + elseif pred == P.LLVMRealOGE + ge() + elseif pred == P.LLVMRealOLT + lt() + elseif pred == P.LLVMRealOLE + le() + elseif pred == P.LLVMRealONE + LLVM.and!(builder, LLVM.not!(builder, eq()), LLVM.not!(builder, uno())) + elseif pred == P.LLVMRealORD + LLVM.not!(builder, uno()) + elseif pred == P.LLVMRealUNO + uno() + elseif pred == P.LLVMRealUEQ + LLVM.or!(builder, eq(), uno()) + elseif pred == P.LLVMRealUGT + LLVM.or!(builder, gt(), uno()) + elseif pred == P.LLVMRealUGE + LLVM.or!(builder, ge(), uno()) + elseif pred == P.LLVMRealULT + LLVM.or!(builder, lt(), uno()) + elseif pred == P.LLVMRealULE + LLVM.or!(builder, le(), uno()) + elseif pred == P.LLVMRealUNE + LLVM.not!(builder, eq()) + else + error("unsupported binary64 fcmp predicate $pred") + end +end + +function extend_integer!(builder, value::LLVM.Value, signed::Bool) + typ = LLVM.value_type(value) + typ isa LLVM.IntegerType || error("binary64 conversion from non-integer $typ") + width = LLVM.width(typ) + width > 64 && error("binary64 conversion from i$width is unsupported") + width == 64 && return value + signed ? LLVM.sext!(builder, value, LLVM.Int64Type()) : + LLVM.zext!(builder, value, LLVM.Int64Type()) +end + +# LLVM intrinsics on double values, by name prefix (covering both scalar and vector forms) +const INTRINSICS = [ + "llvm.sqrt." => :sqrt64, + "llvm.fma." => :fma64, + "llvm.fmuladd." => :fma64, + "llvm.fabs." => :abs64, + "llvm.copysign." => :copysign64, + "llvm.trunc." => :trunc64, + "llvm.rint." => :rint64, + "llvm.nearbyint." => :rint64, + "llvm.roundeven." => :rint64, + "llvm.round." => :round64, + "llvm.floor." => :floor64, + "llvm.ceil." => :ceil64, + "llvm.minimum." => :minimum64, + "llvm.maximum." => :maximum64, + "llvm.minnum." => :minnum64, + "llvm.maxnum." => :maxnum64, +] + +# replace a binary64 operation with calls to the emulation routines, returning whether +# the instruction was outlined +function outline_instruction!(builder, inst::LLVM.Instruction) + double = LLVM.DoubleType() + result_type = LLVM.value_type(inst) + operands = collect(LLVM.operands(inst)) + input_type = isempty(operands) ? nothing : LLVM.value_type(first(operands)) + LLVM.position!(builder, LLVM.before(inst)) + builder.debug_location = inst.debug_location + call(f, args...) = outlined_call!(builder, inst, helper_name(f), LLVM.Value[args...], double) + + replacement = if inst isa LLVM.FAddInst && is_binary64_type(result_type) + scalarize!((a, b) -> call(:add64, a, b), builder, result_type, operands...) + elseif inst isa LLVM.FSubInst && is_binary64_type(result_type) + # subtraction is addition of the negated operand (a sign flip, after inlining) + scalarize!((a, b) -> call(:add64, a, call(:neg64, b)), builder, result_type, operands...) + elseif inst isa LLVM.FMulInst && is_binary64_type(result_type) + scalarize!((a, b) -> call(:mul64, a, b), builder, result_type, operands...) + elseif inst isa LLVM.FDivInst && is_binary64_type(result_type) + scalarize!((a, b) -> call(:div64, a, b), builder, result_type, operands...) + elseif inst isa LLVM.FNegInst && is_binary64_type(result_type) + scalarize!(a -> call(:neg64, a), builder, result_type, operands...) + elseif inst isa LLVM.FRemInst && is_binary64_type(result_type) + error("binary64 remainder (frem) is not supported") + elseif inst isa LLVM.FCmpInst && is_binary64_type(input_type) + pred = LLVM.predicate(inst) + scalarize!((a, b) -> comparison_value!(builder, inst, pred, a, b), + builder, result_type, operands...) + elseif inst isa Union{LLVM.SIToFPInst,LLVM.UIToFPInst} && is_binary64_type(result_type) + signed = inst isa LLVM.SIToFPInst + scalarize!(builder, result_type, operands...) do x + call(signed ? :i64_to_f64 : :u64_to_f64, extend_integer!(builder, x, signed)) + end + elseif inst isa Union{LLVM.FPToSIInst,LLVM.FPToUIInst} && is_binary64_type(input_type) + T = scalar_type(result_type) + T isa LLVM.IntegerType && LLVM.width(T) <= 64 || + error("unsupported conversion from Float64 to $T") + f = inst isa LLVM.FPToSIInst ? :f64_to_i64 : :f64_to_u64 + scalarize!(builder, result_type, operands...) do x + wide = outlined_call!(builder, inst, helper_name(f), LLVM.Value[x], LLVM.Int64Type()) + LLVM.width(T) == 64 ? wide : LLVM.trunc!(builder, wide, T) + end + elseif inst isa LLVM.FPExtInst && is_binary64_type(result_type) + T = scalar_type(input_type) + T in (LLVM.HalfType(), LLVM.FloatType()) || + error("unsupported conversion from $T to Float64") + scalarize!(builder, result_type, operands...) do x + # Widen the bits directly; native fpext may flush a half subnormal to zero. + f, bits_type = T == LLVM.HalfType() ? (:f16_to_f64, LLVM.Int16Type()) : + (:f32_to_f64, LLVM.Int32Type()) + call(f, LLVM.bitcast!(builder, x, bits_type)) + end + elseif inst isa LLVM.FPTruncInst && is_binary64_type(input_type) + T = scalar_type(result_type) + T in (LLVM.HalfType(), LLVM.FloatType()) || + error("unsupported conversion from Float64 to $T") + f, bits_type = T == LLVM.HalfType() ? (:f64_to_f16, LLVM.Int16Type()) : + (:f64_to_f32, LLVM.Int32Type()) + scalarize!(builder, result_type, operands...) do x + bits = outlined_call!(builder, inst, helper_name(f), LLVM.Value[x], bits_type) + LLVM.bitcast!(builder, bits, T) + end + elseif inst isa LLVM.CallInst && LLVM.called_operand(inst) isa LLVM.Function && + is_binary64_type(result_type) + callee = LLVM.name(LLVM.called_operand(inst)) + i = findfirst(((prefix, _),) -> startswith(callee, prefix), INTRINSICS) + i === nothing && return false + f = INTRINSICS[i][2] + scalarize!((args...) -> call(f, args...), builder, result_type, + collect(LLVM.arguments(inst))...) + else + return false + end + LLVM.replace_uses!(inst, replacement) + LLVM.erase!(inst) + true +end + +function outline_float64!(mod::LLVM.Module) + fpclass = tryparse(LLVM.Intrinsic, "llvm.is.fpclass") + for fun in collect(mod.functions) + GPUCompiler.expand_vector_reductions!(fun) + GPUCompiler.scalarize_vector_calls!((call, intrinsic) -> intrinsic == fpclass, fun) + GPUCompiler.lower_intrinsic_calls!(fun) do builder, call, intrinsic + intrinsic == fpclass || return nothing + first(call.arguments).value_type == LLVM.DoubleType() || return nothing + GPUCompiler.lower_is_fpclass!(builder, call) + end + end + @dispose builder=LLVM.IRBuilder() begin + worklist = LLVM.Instruction[] + for fn in LLVM.functions(mod), bb in LLVM.blocks(fn), inst in LLVM.instructions(bb) + push!(worklist, inst) + end + for inst in worklist + outline_instruction!(builder, inst) + end + end + mod +end + +function remap_attribute(remapper, attr::LLVM.Attribute) + attr isa LLVM.TypeAttribute || return attr + LLVM.TypeAttribute(attr.kind, remap_type(remapper, attr.value)) +end + +function copy_global_properties!(new::LLVM.GlobalValue, old::LLVM.GlobalValue) + new.linkage = old.linkage + LLVM.copy_attributes!(new, old) + new +end + +function copy_function_properties!(new::LLVM.Function, old::LLVM.Function) + copy_global_properties!(new, old) + LLVM.callconv!(new, LLVM.callconv(old)) + isempty(LLVM.gc(old)) || LLVM.gc!(new, LLVM.gc(old)) + for (old_arg, new_arg) in zip(LLVM.parameters(old), LLVM.parameters(new)) + LLVM.name!(new_arg, LLVM.name(old_arg)) + end + new +end + +# `nofpclass` constrains floating-point values, not their integer representations. +# Keep it for native floating-point arguments and results, but drop it for mapped ones. +function copy_attributes!(dest, source, remapper, typ=nothing) + attrs = collect(source) # source and destination may be the same call-site set + for attr in collect(dest) + delete!(dest, attr) + end + for attr in attrs + if typ !== nothing && scalar_type(typ) isa LLVM.IntegerType && + attr isa LLVM.EnumAttribute && + attr.kind === :nofpclass + continue + end + push!(dest, remap_attribute(remapper, attr)) + end +end + +function copy_function_attributes!(new::LLVM.Function, old::LLVM.Function, remapper) + # CloneFunctionInto copies attributes verbatim, including their embedded types. + copy_attributes!(LLVM.function_attributes(new), LLVM.function_attributes(old), remapper) + copy_attributes!(LLVM.return_attributes(new), LLVM.return_attributes(old), remapper, + LLVM.return_type(LLVM.function_type(new))) + for (i, arg) in enumerate(LLVM.parameters(new)) + copy_attributes!(new.parameter_attributes[i], old.parameter_attributes[i], + remapper, LLVM.value_type(arg)) + end + new +end + +function map_constant(remapper, value_map::Dict{LLVM.Value,LLVM.Value}, value::LLVM.Value) + haskey(value_map, value) && return value_map[value] + value isa LLVM.Constant || error("cannot materialize non-constant LLVM value $value") + typ = remap_type(remapper, LLVM.value_type(value)) + if value isa LLVM.ConstantFP && LLVM.value_type(value) == LLVM.DoubleType() + return LLVM.const_bitcast(value, LLVM.Int64Type()) + elseif value isa LLVM.ConstantAggregateZero + return LLVM.null(typ) + elseif value isa LLVM.PointerNull + return LLVM.PointerNull(typ) + elseif value isa LLVM.UndefValue + return LLVM.UndefValue(typ) + elseif value isa LLVM.PoisonValue + return LLVM.PoisonValue(typ) + elseif value isa LLVM.ConstantArray + # LLVM.jl's array interface descends to scalar leaves. Preserve one LLVM + # aggregate level here, including zero/undef subarrays, and recurse ourselves. + vals = LLVM.Constant[map_constant(remapper, value_map, x) for x in LLVM.operands(value)] + return LLVM.ConstantArray(typ.element_type, vals) + elseif value isa LLVM.ConstantDataArray + vals = LLVM.Constant[map_constant(remapper, value_map, x) for x in value.elements] + return LLVM.ConstantArray(typ.element_type, vals) + elseif value isa Union{LLVM.ConstantVector,LLVM.ConstantDataVector} + vals = LLVM.Constant[map_constant(remapper, value_map, x) for x in value.elements] + return LLVM.Value(LLVM.API.LLVMConstVector(vals, length(vals))) + elseif value isa LLVM.ConstantStruct + vals = LLVM.Constant[map_constant(remapper, value_map, x) for x in LLVM.operands(value)] + return LLVM.name(typ) === nothing ? LLVM.ConstantStruct(vals; packed=LLVM.ispacked(typ)) : + LLVM.ConstantStruct(typ, vals) + elseif value isa LLVM.ConstantExpr + ops = LLVM.Constant[map_constant(remapper, value_map, x) for x in LLVM.operands(value)] + opcode = LLVM.opcode(value) + if opcode == LLVM.API.LLVMBitCast + return LLVM.value_type(first(ops)) == typ ? first(ops) : LLVM.const_bitcast(first(ops), typ) + elseif opcode == LLVM.API.LLVMAddrSpaceCast + return LLVM.const_addrspacecast(first(ops), typ) + elseif opcode == LLVM.API.LLVMPtrToInt + return LLVM.const_ptrtoint(first(ops), typ) + elseif opcode == LLVM.API.LLVMIntToPtr + return LLVM.const_inttoptr(first(ops), typ) + elseif opcode == LLVM.API.LLVMGetElementPtr + source_type = remap_type(remapper, + LLVM.LLVMType(LLVM.API.LLVMGetGEPSourceElementType(value))) + return LLVM.const_gep(source_type, first(ops), ops[2:end]) + end + # e.g. `icmp` constant expressions (LLVM < 19) that do not involve binary64 + typ == LLVM.value_type(value) && ops == collect(LLVM.operands(value)) && return value + error("unsupported constant expression during binary64 legalization: opcode $opcode") + elseif remap_type(remapper, LLVM.value_type(value)) != LLVM.value_type(value) + error("unsupported changed constant $(typeof(value)) during binary64 legalization") + end + value +end + +function rebuild_module!(mod::LLVM.Module) + LLVM.API.LLVMGetFirstGlobalAlias(mod) == C_NULL || + error("binary64 legalization does not yet support LLVM global aliases") + LLVM.API.LLVMGetFirstGlobalIFunc(mod) == C_NULL || + error("binary64 legalization does not yet support LLVM ifuncs") + remapper = Float64TypeRemapper() + value_map = Dict{LLVM.Value,LLVM.Value}() + + # Only globals whose type changes, and functions with a body, need to be rebuilt. + # Other declarations are kept as they are: recreating an intrinsic declaration, in + # particular, is not safe with typed pointers. This does require that any remaining + # double-typed intrinsic declaration is unused (see `check_float64_calls`). + for fn in collect(LLVM.functions(mod)) + if LLVM.isdeclaration(fn) && startswith(LLVM.name(fn), "llvm.") && + contains_double(LLVM.function_type(fn)) + isempty(LLVM.uses(fn)) || error("unsupported binary64 intrinsic: $(LLVM.name(fn))") + LLVM.erase!(fn) + end + end + old_functions = LLVM.Function[] + for fn in LLVM.functions(mod) + if LLVM.isdeclaration(fn) && !contains_double(LLVM.function_type(fn)) + value_map[fn] = fn + else + push!(old_functions, fn) + end + end + old_globals = LLVM.GlobalVariable[] + kept_globals = LLVM.GlobalVariable[] + for gv in LLVM.globals(mod) + if contains_double(LLVM.global_value_type(gv)) + push!(old_globals, gv) + else + value_map[gv] = gv + push!(kept_globals, gv) + end + end + + for old in old_functions + LLVM.name!(old, LLVM.name(old) * ".softfloat64.old") + end + for old in old_globals + LLVM.name!(old, LLVM.name(old) * ".softfloat64.old") + end + + for old in old_globals + final_name = replace(LLVM.name(old), r"\.softfloat64\.old$" => "") + typ = remap_type(remapper, LLVM.global_value_type(old)) + new = LLVM.GlobalVariable(mod, typ, final_name, LLVM.addrspace(LLVM.value_type(old))) + copy_global_properties!(new, old) + new.constant = old.constant + value_map[old] = new + end + for old in old_functions + final_name = replace(LLVM.name(old), r"\.softfloat64\.old$" => "") + new = LLVM.Function(mod, final_name, remap_type(remapper, LLVM.function_type(old))) + copy_function_properties!(new, old) + value_map[old] = new + for (old_arg, new_arg) in zip(LLVM.parameters(old), LLVM.parameters(new)) + value_map[old_arg] = new_arg + end + end + + materializer(value) = map_constant(remapper, value_map, value) + for old in old_functions + LLVM.isdeclaration(old) && continue + LLVM.clone_into!(value_map[old], old; value_map, + changes=LLVM.API.LLVMCloneFunctionChangeTypeGlobalChanges, + type_mapper=t -> remap_type(remapper, t), materializer) + end + for old in old_functions + copy_function_attributes!(value_map[old], old, remapper) + for bb in LLVM.blocks(value_map[old]), inst in LLVM.instructions(bb) + inst isa LLVM.CallInst || continue + attrs = LLVM.function_attributes(inst) + copy_attributes!(attrs, attrs, remapper) + attrs = LLVM.return_attributes(inst) + copy_attributes!(attrs, attrs, remapper, LLVM.value_type(inst)) + for (i, arg) in enumerate(LLVM.arguments(inst)) + attrs = inst.argument_attributes[i] + copy_attributes!(attrs, attrs, remapper, LLVM.value_type(arg)) + end + end + end + for old in old_globals + init = LLVM.initializer(old) + init === nothing || LLVM.initializer!(value_map[old], materializer(init)) + end + # initializers of unchanged globals may still refer to rebuilt functions or globals + for gv in kept_globals + init = LLVM.initializer(gv) + init === nothing && continue + new_init = materializer(init) + new_init == init || LLVM.initializer!(gv, new_init) + end + + # Rewrite kernel annotations and other named metadata while both generations exist. + for old in old_functions + LLVM.replace_metadata_uses!(old, value_map[old]) + end + for old in old_globals + LLVM.replace_metadata_uses!(old, value_map[old]) + end + + # Remove all old bodies/initializers first so cross-references among old definitions no + # longer prevent safe erasure. + for old in old_functions + LLVM.isdeclaration(old) || empty!(old) + end + for old in old_globals + LLVM.initializer!(old, nothing) + end + foreach(LLVM.erase!, old_functions) + foreach(LLVM.erase!, old_globals) + return +end + +# Opaque pointers hide the allocated/indexed type from result and operand types, and calls +# carry a function type independently of their callee pointer. +function embedded_type(inst::LLVM.Instruction) + if inst isa LLVM.AllocaInst + LLVM.LLVMType(LLVM.API.LLVMGetAllocatedType(inst)) + elseif inst isa LLVM.GetElementPtrInst + LLVM.LLVMType(LLVM.API.LLVMGetGEPSourceElementType(inst)) + elseif inst isa LLVM.CallInst + LLVM.called_type(inst) + else + nothing + end +end + +typed_attributes(attrs) = [LLVM.value(attr) for attr in collect(attrs) if attr isa LLVM.TypeAttribute] + +function attribute_types(value::Union{LLVM.Function,LLVM.CallInst}) + types = LLVM.LLVMType[] + append!(types, typed_attributes(LLVM.function_attributes(value))) + append!(types, typed_attributes(LLVM.return_attributes(value))) + if value isa LLVM.Function + for i in eachindex(LLVM.parameters(value)) + append!(types, typed_attributes(value.parameter_attributes[i])) + end + else + for i in eachindex(LLVM.arguments(value)) + append!(types, typed_attributes(value.argument_attributes[i])) + end + end + types +end + +# call `f(typ, describe)` for every type in the module, where `describe()` returns a +# description of where the type occurs +function foreach_type(f, mod::LLVM.Module) + for gv in LLVM.globals(mod) + f(LLVM.global_value_type(gv), () -> "global $(LLVM.name(gv))") + end + for fn in LLVM.functions(mod) + f(LLVM.function_type(fn), () -> "function $(LLVM.name(fn))") + for typ in attribute_types(fn) + f(typ, () -> "attribute of function $(LLVM.name(fn))") + end + for bb in LLVM.blocks(fn), inst in LLVM.instructions(bb) + describe() = "instruction in $(LLVM.name(fn)): $inst" + f(LLVM.value_type(inst), describe) + typ = embedded_type(inst) + typ === nothing || f(typ, describe) + for op in LLVM.operands(inst) + f(LLVM.value_type(op), describe) + end + inst isa LLVM.CallInst || continue + for typ in attribute_types(inst) + f(typ, describe) + end + end + end +end + +function uses_float64(mod::LLVM.Module) + found = false + foreach_type(mod) do typ, _ + found |= contains_double(typ) + end + found +end + +# check that no trace of binary64 remains, returning a list of problems +function validate_module(mod::LLVM.Module) + problems = String[] + foreach_type(mod) do typ, describe + contains_double(typ) && push!(problems, describe() * " retains double") + end + for fn in LLVM.functions(mod) + if LLVM.isdeclaration(fn) && startswith(LLVM.name(fn), "llvm.") && + occursin(".f64", LLVM.name(fn)) && !isempty(LLVM.uses(fn)) + push!(problems, "f64 intrinsic survived: $(LLVM.name(fn))") + end + end + problems +end + +# reject calls with a binary64 ABI to functions we do not control +function check_float64_calls(mod::LLVM.Module, imports) + for fn in LLVM.functions(mod), bb in LLVM.blocks(fn), inst in LLVM.instructions(bb) + inst isa LLVM.CallInst || continue + callee = LLVM.called_operand(inst) + # byval/sret and other typed attributes are part of the ABI even when every + # explicit argument is an opaque pointer. Check the declaration as well, since + # not all attributes have to be repeated at the call site. + (contains_double(LLVM.called_type(inst)) || + any(contains_double, attribute_types(inst)) || + (callee isa LLVM.Function && any(contains_double, attribute_types(callee)))) || continue + # internal functions are rebuilt with their callers, but external functions, + # indirect calls and inline assembly have no such ABI agreement + if callee isa LLVM.Function + !LLVM.isdeclaration(callee) && continue + name = LLVM.name(callee) + any(m -> m.llvm_name == name, METHODS) && continue + name in imports && continue + end + error("unsupported binary64 call during software legalization: $inst") + end +end + +# replace all binary64 operations and values in `mod`. functions and globals keep their +# names, but are replaced by new ones when their type changes. +function legalize_module!(mod::LLVM.Module; imports=Set{String}()) + uses_float64(mod) || return + outline_float64!(mod) + check_float64_calls(mod, imports) + rebuild_module!(mod) + LLVM.verify(mod) + problems = validate_module(mod) + isempty(problems) || error("incomplete binary64 legalization:\n" * join(problems, "\n")) + return +end + +""" + legalize!(job, mod) + +Legalize binary64 in a target's `lower_intrinsics!` hook. Julia libraries selected by the +same target undergo this hook before they are imported late, so their exported ABIs are +rewritten together with the caller. External bitcode and unregistered functions remain +strict ABI boundaries and cannot be silently rewritten. +""" +function legalize!(job, mod::LLVM.Module) + imports = Set(method.llvm_name for library in GPUCompiler.selected_device_libraries(job) + if library.bitcode === nothing for method in library.methods) + legalize_module!(mod; imports) +end diff --git a/src/softfloat/paynehanek.jl b/src/softfloat/paynehanek.jl new file mode 100644 index 000000000..88d22eecc --- /dev/null +++ b/src/softfloat/paynehanek.jl @@ -0,0 +1,76 @@ +# Payne-Hanek argument reduction without 128-bit integers. +# +# Adapted from Julia's base/special/rem_pio2.jl (MIT, see LICENSES/julia-MIT.txt), which +# uses `UInt128` arithmetic that many GPU back-ends cannot lower. This version replaces it +# by the two-word `U128`, and is otherwise identical, so that back-ends can override +# `Base.Math.paynehanek(::Float64)` with it to make Base's trigonometric functions compile +# for large arguments. The results are bit-for-bit equal to Base's. + +# convert a signed 128-bit fraction to two Float64s, the first truncated to 26 bits +function fromfraction(f::U128) + iszero(f) && return (0.0, 0.0) + + s = f.hi & sign_mask(Float64) # sign bit + x = s == 0 ? f : -f # magnitude + + # 1. get leading term truncated to 26 bits + n1 = 128 - leading_zeros(x) + m1 = (x >> (n1 - 26)).lo << 27 + d1 = ((n1 - 128 + 1021) % UInt64) << 52 + z1 = reinterpret(Float64, s | (d1 + m1)) + + # 2. compute remaining term + x2 = x - (U128(m1) << (n1 - 53)) + iszero(x2) && return (z1, 0.0) + n2 = 128 - leading_zeros(x2) + m2 = (x2 >> (n2 - 53)).lo + d2 = ((n2 - 128 + 1021) % UInt64) << 52 + z2 = reinterpret(Float64, s | (d2 + m2)) + return (z1, z2) +end + +# see `Base.Math.paynehanek` for the derivation; only the arithmetic differs +function paynehanek(x::Float64) + # 1. convert to form x = X * 2^k, where X is a 53-bit integer + u = reinterpret(UInt64, x) + X = (u & significand_mask(Float64)) | (one(UInt64) << significand_bits(Float64)) + raw_exponent = ((u & exponent_mask(Float64)) >> significand_bits(Float64)) % Int + k = raw_exponent - exponent_bias(Float64) - significand_bits(Float64) + + # 2. extract the relevant 192-bit window of 1/2π (the caller guarantees |x| >= 2^20·π/2, + # so the table indices are in bounds, as in Base) + idx = k >> 6 + shift = k - (idx << 6) + INV_2PI = Base.Math.INV_2PI + @inbounds if shift == 0 + a1 = INV_2PI[idx+1] + a2 = INV_2PI[idx+2] + a3 = INV_2PI[idx+3] + else + a1 = (idx < 0 ? zero(UInt64) : INV_2PI[idx+1] << shift) | (INV_2PI[idx+2] >> (64 - shift)) + a2 = (INV_2PI[idx+2] << shift) | (INV_2PI[idx+3] >> (64 - shift)) + a3 = (INV_2PI[idx+3] << shift) | (INV_2PI[idx+4] >> (64 - shift)) + end + + # 3. multiply, keeping the fractional 128 bits of the quotient after division by 2π + w1 = U128(X * a1, 0) # overflow becomes integer + w2 = widemul(X, a2) + w3 = widemul(X, a3) >> 64 + w = w1 + w2 + w3 + signbit(x) && (w = -w) + + # 4. convert to quadrant, and the quotient fraction after division by π/2 + q = ((w >> 125).lo % Int + 1) >> 1 + f = w << 2 + + # 5. convert quotient fraction to split precision Float64 + z_hi, z_lo = fromfraction(f) + + # 6. multiply by π/2 + pio2 = 1.5707963267948966 + pio2_hi = 1.5707963407039642 + pio2_lo = -1.3909067614167116e-8 + y_hi = (z_hi + z_lo) * pio2 + y_lo = (((z_hi * pio2_hi - y_hi) + z_hi * pio2_lo) + z_lo * pio2_hi) + z_lo * pio2_lo + return q, Base.Math.DoubleFloat64(y_hi, y_lo) +end diff --git a/src/softfloat/provider.jl b/src/softfloat/provider.jl new file mode 100644 index 000000000..b7c50a655 --- /dev/null +++ b/src/softfloat/provider.jl @@ -0,0 +1,89 @@ +# device-library provider for binary64 emulation + +""" + softfloat_libraries() + +Device libraries implementing IEEE 754 binary64 (`Float64`) arithmetic using only +64-bit integer operations, for targets without native double-precision support. + +The selecting target calls `legalize!` from `lower_intrinsics!`, including when building +its Julia libraries. All `double` values in the LLVM module are replaced by `i64` bit +patterns, and floating-point operations on them by calls to the [`SoftFloat64`](@ref) +routines (see `legalize_module!`). This is transparent to Julia code: ordinary `Float64` +code, including the pure-Julia math functions in Base, compiles as usual. + +Supported semantics: round-to-nearest-even, gradual underflow, signed zeros, infinities, +quiet comparisons, and canonical quiet NaN results (`0x7ff8000000000000`). Not supported: +floating-point exception flags, other rounding modes, and `Float64` atomics. +""" +softfloat_libraries() = + (GPUCompiler.DeviceLibrary(:softfloat64; + methods=Tuple(m for m in METHODS if !(m.name in OUTLINED)), + optimize=true, always_inline=true), + GPUCompiler.DeviceLibrary(:softfloat64_outlined; + methods=Tuple(m for m in METHODS if m.name in OUTLINED), + optimize=true, always_inline=false)) + +# Routines kept out of line. GPUCompiler has no cost-based inliner, so the others are +# forced inline: they are at most ~45 instructions once optimized (sign operations 2-4, +# comparisons 7-33, rounding 19-40, other conversions 27-43), and are often called with +# constant operands that simplify further. These are 127 (`mul64`) to 444 (`fma64`) +# instructions, or 79 for the narrowing conversions, and large kernels (e.g. Base's `sin`) +# call `add64` and `mul64` over a hundred times each. +const OUTLINED = (:add64, :mul64, :div64, :sqrt64, :fma64, :f64_to_f32, :f64_to_f16) + +# name of the LLVM function implementing one of the library's routines +helper_name(name::Symbol) = "gpu_softfloat_$name" + +# The library's ABI exchanges binary64 values as `i64`, other floating-point values as +# their bits too, and `Bool` as `i8` like any Julia function. These shims adapt the +# SoftFloat64 methods where their signature differs. +rint(x::SoftFloat64) = round(x, RoundNearest) +round_ties_away(x::SoftFloat64) = round(x, RoundNearestTiesAway) +i64_to_f64(x::Int64) = SoftFloat64(x) +u64_to_f64(x::UInt64) = SoftFloat64(x) +f64_to_i64(x::SoftFloat64) = unsafe_trunc(Int64, x) +f64_to_u64(x::SoftFloat64) = unsafe_trunc(UInt64, x) +f16_to_f64(x::UInt16) = SoftFloat64(reinterpret(Float16, x)) +f32_to_f64(x::UInt32) = SoftFloat64(reinterpret(Float32, x)) +f64_to_f32(x::SoftFloat64) = reinterpret(UInt32, Float32(x)) +f64_to_f16(x::SoftFloat64) = reinterpret(UInt16, Float16(x)) + +const METHODS = map([ + # arithmetic + (:add64, +, SoftFloat64, (SoftFloat64, SoftFloat64)), + (:mul64, *, SoftFloat64, (SoftFloat64, SoftFloat64)), + (:div64, /, SoftFloat64, (SoftFloat64, SoftFloat64)), + (:sqrt64, sqrt, SoftFloat64, (SoftFloat64,)), + (:fma64, fma, SoftFloat64, (SoftFloat64, SoftFloat64, SoftFloat64)), + (:neg64, -, SoftFloat64, (SoftFloat64,)), + (:abs64, abs, SoftFloat64, (SoftFloat64,)), + (:copysign64, copysign, SoftFloat64, (SoftFloat64, SoftFloat64)), + # comparisons + (:eq64, ==, Bool, (SoftFloat64, SoftFloat64)), + (:lt64, <, Bool, (SoftFloat64, SoftFloat64)), + (:le64, <=, Bool, (SoftFloat64, SoftFloat64)), + (:unordered64, unordered, Bool, (SoftFloat64, SoftFloat64)), + (:minimum64, min, SoftFloat64, (SoftFloat64, SoftFloat64)), + (:maximum64, max, SoftFloat64, (SoftFloat64, SoftFloat64)), + (:minnum64, minnum, SoftFloat64, (SoftFloat64, SoftFloat64)), + (:maxnum64, maxnum, SoftFloat64, (SoftFloat64, SoftFloat64)), + # rounding to integral values + (:trunc64, trunc, SoftFloat64, (SoftFloat64,)), + (:rint64, rint, SoftFloat64, (SoftFloat64,)), + (:round64, round_ties_away, SoftFloat64, (SoftFloat64,)), + (:floor64, floor, SoftFloat64, (SoftFloat64,)), + (:ceil64, ceil, SoftFloat64, (SoftFloat64,)), + # conversions + (:i64_to_f64, i64_to_f64, SoftFloat64, (Int64,)), + (:u64_to_f64, u64_to_f64, SoftFloat64, (UInt64,)), + (:f64_to_i64, f64_to_i64, Int64, (SoftFloat64,)), + (:f64_to_u64, f64_to_u64, UInt64, (SoftFloat64,)), + (:f16_to_f64, f16_to_f64, SoftFloat64, (UInt16,)), + (:f32_to_f64, f32_to_f64, SoftFloat64, (UInt32,)), + (:f64_to_f32, f64_to_f32, UInt32, (SoftFloat64,)), + (:f64_to_f16, f64_to_f16, UInt16, (SoftFloat64,)), + ]) do (name, f, return_type, argument_types) + GPUCompiler.DeviceLibraryMethod(f, return_type, argument_types; + name, llvm_name=helper_name(name)) +end diff --git a/src/softfloat/uint128.jl b/src/softfloat/uint128.jl new file mode 100644 index 000000000..b681b1d25 --- /dev/null +++ b/src/softfloat/uint128.jl @@ -0,0 +1,74 @@ +# A two-word unsigned 128-bit integer, for the significand products of binary64 emulation. +# +# Many GPU back-ends cannot lower `UInt128` arithmetic, so the few operations the emulation +# needs are defined here on a pair of `UInt64`s, with the semantics of the corresponding +# `UInt128` operations: shifts by 128 bits or more yield zero, and negative shift counts +# shift the other way, as for Base's integers. + +struct U128 + hi::UInt64 + lo::UInt64 +end +U128(x::UInt64) = U128(0, x) + +Base.zero(::Type{U128}) = U128(0, 0) +Base.iszero(a::U128) = (a.hi | a.lo) == 0 +Base.:(==)(a::U128, b::U128) = (a.hi == b.hi) & (a.lo == b.lo) +Base.:<(a::U128, b::U128) = (a.hi < b.hi) | ((a.hi == b.hi) & (a.lo < b.lo)) +Base.:<=(a::U128, b::U128) = !(b < a) + +function Base.:+(a::U128, b::U128) + lo = a.lo + b.lo + U128(a.hi + b.hi + (lo < a.lo), lo) +end +function Base.:-(a::U128, b::U128) + lo = a.lo - b.lo + U128(a.hi - b.hi - (a.lo < b.lo), lo) +end +Base.:-(a::U128) = zero(U128) - a + +# (the operators are not defined in terms of each other, which would make them recursive +# and thus prevent them from being inlined) +shift_left(a::U128, n::Int) = + n < 64 ? U128((a.hi << n) | (a.lo >> (64 - n)), a.lo << n) : U128(a.lo << (n - 64), 0) +shift_right(a::U128, n::Int) = + n < 64 ? U128(a.hi >> n, (a.lo >> n) | (a.hi << (64 - n))) : U128(0, a.hi >> (n - 64)) +Base.:<<(a::U128, n::Int) = n >= 0 ? shift_left(a, n) : shift_right(a, -n) +Base.:>>(a::U128, n::Int) = n >= 0 ? shift_right(a, n) : shift_left(a, -n) + +Base.leading_zeros(a::U128) = a.hi == 0 ? 64 + leading_zeros(a.lo) : leading_zeros(a.hi) + +""" + widemul(a, b) + +The full product of two unsigned integers, like `Base.widemul`, except that the product of +two `UInt64`s is a `U128`. +""" +widemul(a::UInt32, b::UInt32) = UInt64(a) * UInt64(b) +function widemul(a::UInt64, b::UInt64) + # schoolbook multiplication of 32-bit halves, as in Base's 32-bit `widemul` fallback + a0, a1 = a & 0xffff_ffff, a >> 32 + b0, b1 = b & 0xffff_ffff, b >> 32 + p00 = a0 * b0 + mid = a1 * b0 + (p00 >> 32) + p01 = a0 * b1 + (mid & 0xffff_ffff) + U128(a1 * b1 + (mid >> 32) + (p01 >> 32), (p01 << 32) | (p00 & 0xffff_ffff)) +end + +""" + shift_right_jam(a, n) + +Shift `a` right by `n` bits, "jamming" any nonzero bits shifted out into the lowest bit of +the result, which thereby becomes a sticky bit for rounding. `n` may exceed the width of +`a`, in which case the result is `0` or `1`. +""" +@inline shift_right_jam(a::T, n::Int) where {T<:Unsigned} = + n < 8sizeof(T) ? (a >> n) | T((a << (8sizeof(T) - n)) != 0) : T(a != 0) +@inline function shift_right_jam(a::U128, n::Int) + n >= 128 && return U128(UInt64(!iszero(a))) + shifted = a >> n + U128(shifted.hi, shifted.lo | !iszero(a << (128 - n))) +end + +# the top word of a 128-bit significand, with the rest jammed into its lowest bit +sticky_high(a::U128) = a.hi | (a.lo != 0) diff --git a/test/helpers/native.jl b/test/helpers/native.jl index d0f83c3b3..8ff133dab 100644 --- a/test/helpers/native.jl +++ b/test/helpers/native.jl @@ -11,10 +11,11 @@ struct CompilerParams <: AbstractCompilerParams entry_safepoint::Bool method_table relocations::Symbol + softfloat::Bool CompilerParams(entry_safepoint::Bool=false, method_table=test_method_table, - relocations::Symbol=:bake) = - new(entry_safepoint, method_table, relocations) + relocations::Symbol=:bake, softfloat::Bool=false) = + new(entry_safepoint, method_table, relocations, softfloat) end module Runtime end @@ -77,11 +78,12 @@ end function create_job(@nospecialize(func), @nospecialize(types); entry_safepoint::Bool=false, method_table=test_method_table, - relocations::Symbol=:bake, jlruntime::Bool=true, kwargs...) + relocations::Symbol=:bake, jlruntime::Bool=true, + softfloat::Bool=false, kwargs...) config_kwargs, kwargs = split_kwargs(kwargs, GPUCompiler.CONFIG_KWARGS) source = methodinstance(typeof(func), Base.to_tuple_type(types), Base.get_world_counter()) target = NativeCompilerTarget(;jlruntime) - params = CompilerParams(entry_safepoint, method_table, relocations) + params = CompilerParams(entry_safepoint, method_table, relocations, softfloat) config = CompilerConfig(target, params; kernel=false, config_kwargs...) CompilerJob(source, config), kwargs end diff --git a/test/helpers/spirv.jl b/test/helpers/spirv.jl index adaf7942a..f045a097f 100644 --- a/test/helpers/spirv.jl +++ b/test/helpers/spirv.jl @@ -4,18 +4,29 @@ using ..GPUCompiler import LLVM import ..TestRuntime -struct CompilerParams <: AbstractCompilerParams end +struct CompilerParams <: AbstractCompilerParams + emulate_fp64::Bool + CompilerParams(emulate_fp64::Bool=false) = new(emulate_fp64) +end GPUCompiler.runtime_module(::CompilerJob{<:Any,CompilerParams}) = TestRuntime +GPUCompiler.device_libraries(job::CompilerJob{<:Any,CompilerParams}) = + job.config.params.emulate_fp64 ? + (GPUCompiler.runtime_library(TestRuntime), GPUCompiler.SoftFloat.softfloat_libraries()...) : + (GPUCompiler.runtime_library(TestRuntime),) +function GPUCompiler.lower_intrinsics!(job::CompilerJob{<:Any,CompilerParams}, mod::LLVM.Module) + job.config.params.emulate_fp64 && GPUCompiler.SoftFloat.legalize!(job, mod) + return +end function create_job(@nospecialize(func), @nospecialize(types); supports_fp16=true, supports_fp64=true, supports_bfloat16=false, backend::Symbol, driver::Symbol=:generic, extensions::String="", - atomics::SPIRVAtomics=SPIRVAtomics(), kwargs...) + atomics::SPIRVAtomics=SPIRVAtomics(), emulate_fp64=false, kwargs...) config_kwargs, kwargs = split_kwargs(kwargs, GPUCompiler.CONFIG_KWARGS) source = methodinstance(typeof(func), Base.to_tuple_type(types), Base.get_world_counter()) target = SPIRVCompilerTarget(; backend, validate=true, optimize=true, extensions, atomics, supports_fp16, supports_fp64, supports_bfloat16, driver) - params = CompilerParams() + params = CompilerParams(emulate_fp64) config = CompilerConfig(target, params; kernel=false, config_kwargs...) CompilerJob(source, config), kwargs end diff --git a/test/softfloat.jl b/test/softfloat.jl new file mode 100644 index 000000000..09937196b --- /dev/null +++ b/test/softfloat.jl @@ -0,0 +1,306 @@ +# numerical tests of the SoftFloat64 routines on the CPU, against native Float64 arithmetic + +using GPUCompiler.SoftFloat: SoftFloat64, unordered, minnum, maxnum, paynehanek +using Random, Test + +soft(x::Float64) = SoftFloat64(x) +soft(x::UInt64) = reinterpret(SoftFloat64, x) +# results are compared bit for bit, except that NaNs are canonicalized +canonical(x::Float64) = isnan(x) ? NaN : x +canonical(x::SoftFloat64) = canonical(Float64(x)) + +const EDGE_BITS = UInt64[ + 0x0000_0000_0000_0000, 0x8000_0000_0000_0000, # ±0 + 0x0000_0000_0000_0001, 0x000f_ffff_ffff_ffff, # smallest and largest subnormal + 0x0010_0000_0000_0000, 0x3ff0_0000_0000_0000, # smallest normal, 1 + 0x3ff0_0000_0000_0001, 0x7fef_ffff_ffff_ffff, # nextfloat(1), floatmax + 0x7ff0_0000_0000_0000, 0xfff0_0000_0000_0000, # ±Inf + 0x7ff0_0000_0000_0001, 0x7ff8_0000_0000_0042, # signaling and payloaded NaN +] + +rng = Xoshiro(0x5f64) +corpus = reinterpret(Float64, vcat(EDGE_BITS, rand(rng, UInt64, 25_000))) +edges = reinterpret(Float64, EDGE_BITS) +pairs = [(corpus[i], corpus[mod1(7919i, length(corpus))]) for i in eachindex(corpus)] + +@testset "arithmetic" begin + for (a, b) in pairs + @test canonical(soft(a) + soft(b)) === canonical(a + b) + @test canonical(soft(a) - soft(b)) === canonical(a - b) + @test canonical(soft(a) * soft(b)) === canonical(a * b) + @test canonical(soft(a) / soft(b)) === canonical(a / b) + end + + # uniformly random bits rarely exercise cancellation or the subnormal boundary + for i in 1:10_000 + a = rand(rng, UInt64) + b = (a & 0x7ff0_0000_0000_0000) | (rand(rng, UInt64) & 0x800f_ffff_ffff_ffff) + for (x, y) in ((a, b), (a & 0x801f_ffff_ffff_ffff, b & 0x801f_ffff_ffff_ffff)) + @test canonical(soft(x) + soft(y)) === canonical(reinterpret(Float64, x) + reinterpret(Float64, y)) + @test canonical(soft(x) - soft(y)) === canonical(reinterpret(Float64, x) - reinterpret(Float64, y)) + end + end + for a in edges, b in edges + @test canonical(soft(a) + soft(b)) === canonical(a + b) + @test canonical(soft(a) * soft(b)) === canonical(a * b) + @test canonical(soft(a) / soft(b)) === canonical(a / b) + end + + for a in corpus + @test canonical(sqrt(soft(a))) === canonical(a < 0 ? NaN : sqrt(a)) + @test Float64(abs(soft(a))) === abs(a) + @test Float64(-soft(a)) === -a + end + for (a, b) in pairs + @test Float64(copysign(soft(a), soft(b))) === copysign(a, b) + end + + for i in eachindex(corpus) + a = corpus[i] + b = corpus[mod1(3571i, length(corpus))] + c = corpus[mod1(12347i, length(corpus))] + @test canonical(fma(soft(a), soft(b), soft(c))) === canonical(fma(a, b, c)) + end + for a in edges, b in edges, c in edges + @test canonical(fma(soft(a), soft(b), soft(c))) === canonical(fma(a, b, c)) + end +end + +@testset "rounding modes" begin + # Directed rounding is checked against a BigFloat computation of the exact result (the + # precision covers the exact product and sum in an FMA) that is rounded once. The + # computation itself also happens under the mode, for the sign of exact zero results. + reference(mode, f, xs...) = setrounding(BigFloat, mode) do + setprecision(BigFloat, 4300) do + Float64(f(map(BigFloat, xs)...), mode) + end + end + modes = (RoundNearest, RoundToZero, RoundDown, RoundUp) + finite = filter(isfinite, corpus) + # add values near the overflow and underflow boundaries + boundaries = [floatmax(Float64), -floatmax(Float64), prevfloat(floatmax(Float64)), + floatmin(Float64), nextfloat(0.0), prevfloat(floatmin(Float64)), 0.0, -0.0, + 1.0, -1.0, 2.0^-1074 * 3, 1e-160, -1e-160, 0.5, 3.0] + inputs = vcat(finite[1:5000], boundaries) + for mode in modes, i in eachindex(inputs) + a = inputs[i] + b = inputs[mod1(7919i, length(inputs))] + c = inputs[mod1(3571i, length(inputs))] + @test Float64(+(soft(a), soft(b), mode)) === reference(mode, +, a, b) + @test Float64(-(soft(a), soft(b), mode)) === reference(mode, -, a, b) + @test Float64(*(soft(a), soft(b), mode)) === reference(mode, *, a, b) + iszero(b) || @test Float64(/(soft(a), soft(b), mode)) === reference(mode, /, a, b) + a < 0 || @test Float64(sqrt(soft(a), mode)) === reference(mode, sqrt, a) + @test Float64(fma(soft(a), soft(b), soft(c), mode)) === reference(mode, (x, y, z) -> x * y + z, a, b, c) + end + for a in boundaries, b in boundaries + for mode in modes + @test Float64(+(soft(a), soft(b), mode)) === reference(mode, +, a, b) + @test Float64(*(soft(a), soft(b), mode)) === reference(mode, *, a, b) + end + # exact cancellation gives -0.0 only when rounding down + @test Float64(+(soft(a), soft(-a), RoundDown)) === -0.0 + @test Float64(+(soft(a), soft(-a), RoundUp)) === 0.0 + end + for (a, b, c) in ((floatmax(Float64), 2.0, -floatmax(Float64)), + (floatmin(Float64), 0.5, -prevfloat(floatmin(Float64))), + (prevfloat(1.0), nextfloat(1.0), -1.0)) + for mode in modes + @test Float64(fma(soft(a), soft(b), soft(c), mode)) === + reference(mode, (x, y, z) -> x * y + z, a, b, c) + end + end + # Random triples almost never cancel a product. Use the rounded product as the + # opposing addend, including cases where the residual is subnormal. + for i in 1:2000 + a = ldexp(1 + rand(rng), rand(rng, -500:500)) + b = ldexp(1 + rand(rng), rand(rng, -500:500)) + c = -(a * b) + for mode in modes + @test Float64(fma(soft(a), soft(b), soft(c), mode)) === + reference(mode, (x, y, z) -> x * y + z, a, b, c) + end + end + # overflow rounds to the largest finite value when rounding towards zero + @test Float64(*(soft(floatmax(Float64)), soft(2.0), RoundToZero)) === floatmax(Float64) + @test Float64(*(soft(floatmax(Float64)), soft(2.0), RoundDown)) === floatmax(Float64) + @test Float64(*(soft(floatmax(Float64)), soft(2.0), RoundUp)) === Inf + @test Float64(*(soft(-floatmax(Float64)), soft(2.0), RoundUp)) === -floatmax(Float64) + @test Float64(*(soft(-floatmax(Float64)), soft(2.0), RoundDown)) === -Inf +end + +@testset "comparisons" begin + for (a, b) in vcat(pairs, [(a, b) for a in edges, b in edges][:]) + x, y = soft(a), soft(b) + @test (x == y) == (a == b) + @test (x < y) == (a < b) + @test (x <= y) == (a <= b) + @test unordered(x, y) == (isnan(a) || isnan(b)) + @test canonical(min(x, y)) === canonical(min(a, b)) + @test canonical(max(x, y)) === canonical(max(a, b)) + # minnum/maxnum return the number when only one operand is NaN + @test canonical(minnum(x, y)) === (isnan(a) ? canonical(b) : isnan(b) ? a : canonical(min(a, b))) + @test canonical(maxnum(x, y)) === (isnan(a) ? canonical(b) : isnan(b) ? a : canonical(max(a, b))) + end + for a in corpus + x = soft(a) + @test isnan(x) == isnan(a) + @test isinf(x) == isinf(a) + @test isfinite(x) == isfinite(a) + @test iszero(x) == iszero(a) + @test issubnormal(x) == issubnormal(a) + @test signbit(x) == signbit(a) + end +end + +@testset "rounding" begin + # random bits rarely have a fractional part, so also exercise moderate magnitudes + inputs = vcat(corpus, (rand(rng, 10_000) .- 0.5) .* 2.0 .^ rand(rng, -5:60, 10_000), + -4.0:0.25:4.0, [2.0^52 + 0.5, 2.0^52 - 0.5, -2.0^53 + 1]) + for a in inputs + x = soft(a) + @test canonical(trunc(x)) === canonical(trunc(a)) + @test canonical(round(x)) === canonical(round(a)) + @test canonical(round(x, RoundNearestTiesAway)) === canonical(round(a, RoundNearestTiesAway)) + @test canonical(floor(x)) === canonical(floor(a)) + @test canonical(ceil(x)) === canonical(ceil(a)) + end +end + +@testset "conversions" begin + for x in vcat(rand(rng, Int64, 25_000), typemin(Int64), typemax(Int64), 0, -1, 1) + @test Float64(SoftFloat64(x)) === Float64(x) + for mode in (RoundToZero, RoundDown, RoundUp) + @test Float64(SoftFloat64(x, mode)) === Float64(BigFloat(x), mode) + end + end + for x in vcat(rand(rng, UInt64, 25_000), typemax(UInt64), 0, 1, UInt64(1) << 63) + @test Float64(SoftFloat64(x)) === Float64(x) + for mode in (RoundToZero, RoundDown, RoundUp) + @test Float64(SoftFloat64(x, mode)) === Float64(BigFloat(x), mode) + end + end + @test Float64(SoftFloat64(Int32(-7))) === -7.0 + @test Float64(SoftFloat64(UInt8(255))) === 255.0 + # float-to-integer routines implement fptosi/fptoui, i.e. truncation of in-range values + for x in vcat([(rand(rng) - 0.5) * 2.0^63 for _ in 1:10_000], [rand(rng) * 100 for _ in 1:10_000], + -1.5, 0.5, -0.0, -2.0^63, prevfloat(2.0^63)) + @test unsafe_trunc(Int64, soft(x)) === unsafe_trunc(Int64, x) + end + for x in vcat([rand(rng) * 2.0^64 for _ in 1:10_000], [rand(rng) * 100 for _ in 1:10_000], + 0.5, 0.0, prevfloat(2.0^64)) + @test unsafe_trunc(UInt64, soft(x)) === unsafe_trunc(UInt64, x) + end + + for x in reinterpret(Float32, rand(rng, UInt32, 25_000)) + @test canonical(SoftFloat64(x)) === canonical(Float64(x)) + end + for h in reinterpret(Float16, UInt16(0):typemax(UInt16)) + @test canonical(SoftFloat64(h)) === canonical(Float64(h)) + end + for x in corpus + @test Float32(soft(x)) === (isnan(x) ? NaN32 : Float32(x)) + for mode in (RoundToZero, RoundDown, RoundUp) + isnan(x) || @test Float32(soft(x), mode) === Float32(BigFloat(x), mode) + end + end + # values above the half-way point in the 2^-150 bin round to the least Float32 subnormal + @test Float32(soft(0x369c_6d9b_5f38_1bd0)) === reinterpret(Float32, 0x0000_0001) + + # every Float16 round-trips + for h in reinterpret(Float16, UInt16(0):typemax(UInt16)) + x = Float64(h) + @test Float16(soft(x)) === (isnan(x) ? NaN16 : h) + end + # both sides of every positive finite Float16 midpoint, including the normal/subnormal + # boundary; an intermediate Float32 rounding would lose the side + for h in UInt16(0):UInt16(0x7bfe) + x = Float64(reinterpret(Float16, h)) + y = Float64(reinterpret(Float16, h + UInt16(1))) + midpoint = (x + y) / 2 + for z in (prevfloat(midpoint), midpoint, nextfloat(midpoint)) + @test Float16(soft(z)) === Float16(z) + @test Float16(soft(-z)) === Float16(-z) + for mode in (RoundToZero, RoundDown, RoundUp) + @test Float16(soft(z), mode) === Float16(BigFloat(z), mode) + @test Float16(soft(-z), mode) === Float16(BigFloat(-z), mode) + end + end + end +end + +@testset "paynehanek" begin + # bit-for-bit identical to Base's UInt128-based implementation + for x in vcat([2.0^20 * pi/2, 1e20, 1e100, floatmax(Float64), 2.0^20, nextfloat(2.0^20)], + [rand(rng) * 2.0^rand(rng, 21:1023) for _ in 1:10_000]) + for y in (x, -x) + @test paynehanek(y) === Base.Math.paynehanek(y) + end + end +end + +@testset "number interface" begin + for x in corpus + sx = soft(x) + @test isequal(eps(sx), soft(eps(x))) + @test hash(sx) == hash(x) + @test isequal(sx, x) + end + @test precision(SoftFloat64) == 53 + @test precision(soft(1.0); base=10) == precision(1.0; base=10) + @test eps(SoftFloat64) === soft(eps(Float64)) + for T in (Bool, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, + Float16, Float32, Float64) + @test soft(1.5) + one(T) === soft(2.5) + @test one(T) + soft(1.5) === soft(2.5) + @test promote_type(SoftFloat64, T) == SoftFloat64 + end + @test length(Set([soft(-0.0), soft(0.0), soft(NaN), soft(1.5)])) == 4 +end + +@testset "rounding ties and exact results" begin + for x in -4.0:0.25:4.0 + @test Float64(round(soft(x), RoundNearestTiesUp)) === round(x, RoundNearestTiesUp) + @test Float64(round(soft(x), RoundFromZero)) === round(x, RoundFromZero) + end + for mode in (RoundNearest, RoundToZero, RoundDown, RoundUp, RoundFromZero, + RoundNearestTiesAway, RoundNearestTiesUp) + for i in 1:1000 + # At most 26 significant bits, so the square is exactly representable. + x = Float64(rand(rng, UInt32) & 0x03ffffff) + y = ldexp(1.0, rand(rng, -900:900)) + @test Float64(/(soft(x*y), soft(y), mode)) === x + @test Float64(sqrt(soft(x*x), mode)) === x + end + end + for mode in (RoundNearestTiesAway, RoundNearestTiesUp) + @test Float64(+(soft(1.0), soft(2.0^-53), mode)) === nextfloat(1.0) + @test Float64(+(soft(-1.0), soft(-2.0^-53), mode)) === + (mode === RoundNearestTiesUp ? -1.0 : prevfloat(-1.0)) + @test Float64(*(soft(nextfloat(0.0)), soft(0.5), mode)) === nextfloat(0.0) + @test Float64(*(soft(-nextfloat(0.0)), soft(0.5), mode)) === + (mode === RoundNearestTiesUp ? -0.0 : -nextfloat(0.0)) + end +end + +@testset "widening with flush-to-zero enabled" begin + # Pass bits through a call boundary so constant folding cannot hide native FP operations. + @noinline widen_bits(x::T) where {T<:Union{UInt16,UInt32}} = + reinterpret(UInt64, SoftFloat64(reinterpret(T === UInt16 ? Float16 : Float32, x))) + inputs = (UInt16[0, 1, 0x03ff, 0x0400, 0x8001, 0x83ff, 0x7c00, 0x7e01], + UInt32[0, 1, 0x007fffff, 0x00800000, 0x80000001, 0x807fffff, 0x7f800000, 0x7fc00001]) + old = get_zero_subnormals() + try + set_zero_subnormals(false) + expected = map(xs -> map(widen_bits, xs), inputs) + if set_zero_subnormals(true) + for (xs, ys) in zip(inputs, expected), (x, y) in zip(xs, ys) + @test widen_bits(x) == y + end + else + @test_skip "flush-to-zero is unavailable on this CPU" + end + finally + set_zero_subnormals(old) + end +end diff --git a/test/softfloat/legalize.jl b/test/softfloat/legalize.jl new file mode 100644 index 000000000..25e9792f4 --- /dev/null +++ b/test/softfloat/legalize.jl @@ -0,0 +1,307 @@ +# LLVM-level tests of the binary64 legalization pass + +using GPUCompiler.SoftFloat: legalize_module!, validate_module +using LLVM, Test + +# The IR below is written with opaque pointers, which LLVM 15/16 contexts do not accept +# unless asked; the typed-pointer path of the pass is covered by the Metal tests. + +@testset "module legalization" begin + LLVM.Context(; opaque_pointers=true) do ctx + ir = """ + %node = type { double, ptr } + @values = constant [2 x double] [double -0.0, double 0x7FF8000000000042] + @vector = constant <2 x double> + + declare double @llvm.sqrt.f64(double) + declare double @llvm.fma.f64(double, double, double) + declare double @llvm.fabs.f64(double) + + define double @transport(double %a, double %b, ptr byval(double) %ignored) { + entry: + %slot = alloca %node + %field = getelementptr %node, ptr %slot, i32 0, i32 0 + store double %a, ptr %field + %loaded = load double, ptr %field + %sum = fadd double %loaded, %b + %product = fmul double %sum, %b + %root = call double @llvm.sqrt.f64(double %product) + %fused = call double @llvm.fma.f64(double %root, double %a, double 1.5) + %negative = fneg double %fused + %absolute = call double @llvm.fabs.f64(double %negative) + %cmp = fcmp ule double %absolute, %a + %selected = select i1 %cmp, double %absolute, double %a + br i1 %cmp, label %left, label %right + left: + br label %exit + right: + br label %exit + exit: + %phi = phi double [ %selected, %left ], [ 0x8000000000000000, %right ] + ret double %phi + } + + define <2 x double> @vectors(<2 x double> %a, <2 x double> %b) { + %x = fdiv <2 x double> %a, %b + ret <2 x double> %x + } + + !air.kernel = !{!0} + !0 = !{ptr @transport} + """ + mod = parse(LLVM.Module, ir) + legalize_module!(mod) + entry = LLVM.functions(mod)["transport"] + output = string(mod) + @test isempty(validate_module(mod)) + @test LLVM.verify(mod) === nothing + @test LLVM.function_type(entry) == LLVM.FunctionType( + LLVM.Int64Type(), [LLVM.Int64Type(), LLVM.Int64Type(), LLVM.PointerType()]) + @test occursin("@values = constant [2 x i64] [i64 -9223372036854775808, i64 9221120237041090626]", output) + @test occursin("%node.softfloat64 = type { i64, ptr }", output) + @test occursin("byval(i64)", output) + @test occursin("!air.kernel", output) + @test !occursin(r"\bdouble\b", output) + @test !occursin(r"\bi128\b", output) + dispose(mod) + end +end + +@testset "nested constant arrays" begin + LLVM.Context(; opaque_pointers=true) do ctx + mod = parse(LLVM.Module, """ + @table = constant [3 x [2 x double]] [ + [2 x double] zeroinitializer, + [2 x double] [double 1.0, double -0.0], + [2 x double] undef] + @cube = constant [1 x [2 x [2 x double]]] [[2 x [2 x double]] [ + [2 x double] [double 2.0, double 3.0], + [2 x double] [double 4.0, double 5.0]]] + define double @lookup(i64 %i, i64 %j) { + %ptr = getelementptr [3 x [2 x double]], ptr @table, i64 0, i64 %i, i64 %j + %x = load double, ptr %ptr + ret double %x + } + """) + legalize_module!(mod) + @test LLVM.verify(mod) === nothing + @test isempty(validate_module(mod)) + table = LLVM.initializer(LLVM.globals(mod)["table"]) + @test LLVM.value_type(table) == LLVM.ArrayType(LLVM.ArrayType(LLVM.Int64Type(), 2), 3) + elements = collect(LLVM.operands(table)) + @test elements[1] isa LLVM.ConstantAggregateZero + @test string(elements[2]) == "[2 x i64] [i64 4607182418800017408, i64 -9223372036854775808]" + @test elements[3] isa LLVM.UndefValue + cube = LLVM.initializer(LLVM.globals(mod)["cube"]) + @test LLVM.value_type(cube) == LLVM.ArrayType(LLVM.ArrayType(LLVM.ArrayType(LLVM.Int64Type(), 2), 2), 1) + @test occursin("getelementptr [3 x [2 x i64]]", string(mod)) + dispose(mod) + end +end + +@testset "call-site types and metadata" begin + LLVM.Context(; opaque_pointers=true) do ctx + mod = parse(LLVM.Module, """ + define void @callee(ptr byval(double) %x) { ret void } + define double @caller(double %x, ptr %p) { + call void @callee(ptr byval(double) %p) + %y = fadd double %x, 1.0, !fpmath !0 + ret double %y + } + !0 = !{float 2.5} + """) + legalize_module!(mod) + @test LLVM.verify(mod) === nothing + @test occursin("call void @callee(ptr byval(i64)", string(mod)) + @test !occursin("!fpmath", string(mod)) + dispose(mod) + end +end + +@testset "outlined calls have a debug location" begin + LLVM.Context(; opaque_pointers=true) do ctx + mod = parse(LLVM.Module, """ + define i1 @caller(double %x) !dbg !3 { + %y = fadd double %x, 1.0, !dbg !4 + %z = fcmp une double %y, 0.0 + ret i1 %z + } + + !llvm.dbg.cu = !{!0} + !llvm.module.flags = !{!5} + !0 = distinct !DICompileUnit(language: DW_LANG_Julia, file: !1, producer: "julia", isOptimized: true, runtimeVersion: 0, emissionKind: LineTablesOnly) + !1 = !DIFile(filename: "julia", directory: ".") + !2 = !DISubroutineType(types: !{}) + !3 = distinct !DISubprogram(name: "caller", linkageName: "caller", scope: null, file: !1, line: 1, type: !2, scopeLine: 1, spFlags: DISPFlagDefinition, unit: !0) + !4 = !DILocation(line: 2, scope: !3) + !5 = !{i32 2, !"Debug Info Version", i32 3} + """) + legalize_module!(mod) + @test LLVM.verify(mod) === nothing + caller = LLVM.functions(mod)["caller"] + calls = [inst for block in LLVM.blocks(caller) for inst in LLVM.instructions(block) + if inst isa LLVM.CallInst] + @test !isempty(calls) + locations = [call.debug_location for call in calls] + @test all(loc -> loc !== nothing && LLVM.scope(loc) == caller.subprogram, locations) + # the call replacing the `fcmp` gets a line-0 location in the caller + @test any(loc -> loc !== nothing && LLVM.line(loc) == 0, locations) + dispose(mod) + end +end + +@testset "external binary64 ABI is not silently changed" begin + LLVM.Context(; opaque_pointers=true) do ctx + mod = parse(LLVM.Module, """ + declare double @external(double) + define double @caller(double %x) { + %y = call double @external(double %x) + ret double %y + } + """) + @test_throws "unsupported binary64 call" legalize_module!(mod) + dispose(mod) + end + for attribute in ("byval", "sret"), indirect in (false, true) + LLVM.Context(; opaque_pointers=true) do ctx + callee = indirect ? "%fn" : "@external" + mod = parse(LLVM.Module, """ + %payload = type { double } + declare void @external(ptr $attribute(%payload)) + define void @caller(ptr %fn, ptr %x) { + call void $callee(ptr $attribute(%payload) %x) + ret void + } + """) + @test LLVM.verify(mod) === nothing + @test_throws "unsupported binary64 call" legalize_module!(mod) + dispose(mod) + end + end +end + +@testset "unsupported binary64 intrinsic diagnostic" begin + LLVM.Context(; opaque_pointers=true) do ctx + mod = parse(LLVM.Module, """ + declare double @llvm.ldexp.f64.i32(double, i32) + define double @caller(double %value, i32 %power) { + %result = call double @llvm.ldexp.f64.i32(double %value, i32 %power) + ret double %result + } + """) + try + @test_throws "unsupported binary64" legalize_module!(mod) + finally + dispose(mod) + end + end +end + +@testset "opaque-pointer storage still requires legalization" begin + LLVM.Context(; opaque_pointers=true) do ctx + mod = parse(LLVM.Module, """ + define ptr @storage() { + %slot = alloca double + %next = getelementptr double, ptr %slot, i64 1 + ret ptr %next + } + """) + @test LLVM.verify(mod) === nothing + @test count(contains("retains double"), validate_module(mod)) == 2 + legalize_module!(mod) + @test occursin("alloca i64", string(mod)) + @test occursin("getelementptr i64", string(mod)) + @test isempty(validate_module(mod)) + dispose(mod) + end +end + +@testset "unrelated equal-width mapping remains generic" begin + LLVM.Context(; opaque_pointers=true) do ctx + remapper = GPUCompiler.SoftFloat.Float64TypeRemapper() + @test GPUCompiler.SoftFloat.remap_type(remapper, LLVM.HalfType()) == LLVM.HalfType() + @test GPUCompiler.SoftFloat.remap_type(remapper, + LLVM.ArrayType(LLVM.DoubleType(), 3)) == LLVM.ArrayType(LLVM.Int64Type(), 3) + end +end + +@testset "native floating-point operations survive" begin + LLVM.Context(; opaque_pointers=true) do ctx + mod = parse(LLVM.Module, """ + define float @native_float(i1 %flag, float %x) { + %converted = uitofp i1 %flag to float + %sum = fadd float %converted, %x + ret float %sum + } + """) + legalize_module!(mod) + output = string(mod) + @test isempty(validate_module(mod)) + @test occursin("uitofp i1", output) + @test occursin("fadd float", output) + dispose(mod) + end +end + +if LLVM.version() >= v"17" + @testset "floating-point class attributes" begin + LLVM.Context(; opaque_pointers=true) do ctx + mod = parse(LLVM.Module, """ + define nofpclass(nan inf) double @callee(double nofpclass(nan) %x, + float nofpclass(nan) %y) { + ret double %x + } + define double @caller(double %x, float %y) { + %z = call nofpclass(nan inf) double @callee(double nofpclass(nan) %x, + float nofpclass(nan) %y) + ret double %z + } + """) + @test LLVM.verify(mod) === nothing + legalize_module!(mod) + @test LLVM.verify(mod) === nothing + output = string(mod) + @test !occursin(r"nofpclass\([^)]*\) i64|i64 nofpclass", output) + @test count("float nofpclass(nan)", output) == 2 + dispose(mod) + end + end +end + +@testset "widening preserves input bits" begin + LLVM.Context(; opaque_pointers=true) do ctx + mod = parse(LLVM.Module, """ + define double @widen(half %x) { + %y = fpext half %x to double + ret double %y + } + """) + legalize_module!(mod) + @test LLVM.verify(mod) === nothing + @test occursin("call i64 @gpu_softfloat_f16_to_f64(i16", string(mod)) + @test !occursin("fpext", string(mod)) + dispose(mod) + end +end + +if LLVM.version() < v"19" + # LLVM 19 removed `icmp` constant expressions; before, Base's `sin` contained one + @testset "unrelated constant expressions" begin + LLVM.Context(; opaque_pointers=true) do ctx + mod = parse(LLVM.Module, """ + @flag = global i8 0 + + define double @select_by_address(double %x) { + %y = fadd double %x, 1.0 + %z = select i1 icmp ult (i64 ptrtoint (ptr @flag to i64), i64 4096), double %x, double %y + ret double %z + } + """) + legalize_module!(mod) + @test LLVM.verify(mod) === nothing + @test isempty(validate_module(mod)) + @test occursin("icmp ult (i64 ptrtoint (ptr @flag to i64), i64 4096)", string(mod)) + dispose(mod) + end + end +end diff --git a/test/softfloat/native.jl b/test/softfloat/native.jl new file mode 100644 index 000000000..bf5b10330 --- /dev/null +++ b/test/softfloat/native.jl @@ -0,0 +1,118 @@ +# execution of the emulated binary64 routines through the native back-end's JIT + +using GPUCompiler.SoftFloat: softfloat_libraries +using LLVM, Test + +module SoftFloatFixtures +soft_widen16(a::UInt64, b::UInt64) = + reinterpret(UInt64, Float64(reinterpret(Float16, a % UInt16))) +soft_widen32(a::UInt64, b::UInt64) = + reinterpret(UInt64, Float64(reinterpret(Float32, a % UInt32))) +soft_add(a::UInt64, b::UInt64) = + reinterpret(UInt64, reinterpret(Float64, a) + reinterpret(Float64, b)) +# (using `sqrt_llvm` directly, as the native target has no runtime to throw a DomainError) +soft_math(a::UInt64, b::UInt64) = + reinterpret(UInt64, Base.sqrt_llvm(abs(reinterpret(Float64, a))) * reinterpret(Float64, b) - + fma(reinterpret(Float64, a), 2.0, 0.5)) +end + +# the provider is only selected for the fixtures, not for any other native job +GPUCompiler.device_libraries(job::Native.NativeCompilerJob) = + job.config.params.softfloat ? + (GPUCompiler.runtime_library(TestRuntime), softfloat_libraries()...) : () + +function GPUCompiler.lower_intrinsics!(job::Native.NativeCompilerJob, mod::LLVM.Module) + job.config.params.softfloat && GPUCompiler.SoftFloat.legalize!(job, mod) + return +end + +function compile_and_load(f) + job, _ = Native.create_job(f, (UInt64, UInt64); softfloat=true) + obj, entry, relocs = GPUCompiler.JuliaContext() do ctx + obj, meta = GPUCompiler.compile(:obj, job) + obj, LLVM.name(meta.entry), meta.relocations + end + Native.load(Vector{UInt8}(codeunits(obj)), entry, relocs) +end + +@testset "emulated binary64 on the native back-end" begin + ir = sprint(io -> Native.code_llvm(io, SoftFloatFixtures.soft_add, Tuple{UInt64, UInt64}; + dump_module=true, softfloat=true)) + @test !occursin(r"\bdouble\b", ir) + @test !occursin("declare i64 @gpu_softfloat_add64", ir) + + for (f, reference) in ((SoftFloatFixtures.soft_add, (a, b) -> a + b), + (SoftFloatFixtures.soft_math, (a, b) -> sqrt(abs(a)) * b - fma(a, 2.0, 0.5))) + ptr, jit, _ = compile_and_load(f) + try + for (a, b) in ((1.0, 2.0), (floatmax(Float64), -floatmax(Float64)), + (nextfloat(0.0), nextfloat(0.0)), (-0.0, -0.0), (1e-300, 1e300)) + result = ccall(ptr, UInt64, (UInt64, UInt64), reinterpret(UInt64, a), reinterpret(UInt64, b)) + @test result == reinterpret(UInt64, reference(a, b)) + end + finally + LLVM.dispose(jit) + end + end + + plain, _ = Native.create_job(identity, (UInt64,)) + @test isempty(GPUCompiler.device_libraries(plain)) +end + +@testset "library routines compile to integer code" begin + # every routine must lower to plain integer arithmetic: no double-precision values, + # no 128-bit integers, and no calls into the Julia runtime (e.g. to throw) + for method in GPUCompiler.SoftFloat.METHODS + ir = sprint(io -> Native.code_llvm(io, method.def, method.types; dump_module=true)) + @test !occursin(r"\bdouble\b", ir) + @test !occursin(r"(? !startswith(String(lib.name), "softfloat64"), GPUCompiler.device_libraries(native_job)) +end