From 7cff0f8b9db51a69d1a7c1d0ecdbb284f5f07b13 Mon Sep 17 00:00:00 2001
From: Tim Besard <tim.besard@gmail.com>
Date: Wed, 23 Sep 2026 00:08:11 +0200
Subject: [PATCH] Use GPUToolbox's Float64-free math overrides on devices
 without fp64

Base computes some single-precision math in double precision, e.g.,
`div` on Float32 (on Julia 1.12 and 1.13), `sind` and `cosd`, or
inversion of ComplexF32. On devices that don't support Float64, stack
GPUToolbox.Overlays.float64_overrides underneath SPIRVIntrinsics'
method table to keep these computations in single precision.

This replaces the `div(::Float32, ::Float32)` quirk, which truncated
the rounded quotient and thus returned incorrect results (e.g.,
`div(1f0, 0.1f0) == 10`), for all devices.
---
 Project.toml                |  2 +-
 src/compiler/compilation.jl | 12 ++++++++++--
 src/device/quirks.jl        |  4 ----
 src/oneAPI.jl               |  1 +
 4 files changed, 12 insertions(+), 7 deletions(-)

diff --git a/Project.toml b/Project.toml
index 0d9b1798..55486501 100644
--- a/Project.toml
+++ b/Project.toml
@@ -40,7 +40,7 @@ CEnum = "0.4, 0.5"
 ExprTools = "0.1"
 GPUArrays = "11.5.14"
 GPUCompiler = "2.8"
-GPUToolbox = "0.1, 0.2, 0.3, 1, 3"
+GPUToolbox = "3.1"
 KernelAbstractions = "0.9.39"
 LLVM = "6, 7, 8, 9"
 NEO_jll = "=26.18.38308"
diff --git a/src/compiler/compilation.jl b/src/compiler/compilation.jl
index 92557313..f1472403 100644
--- a/src/compiler/compilation.jl
+++ b/src/compiler/compilation.jl
@@ -32,8 +32,16 @@ GPUCompiler.runtime_module(::oneAPICompilerJob) = oneAPI
 
 GPUCompiler.kernel_state_type(::oneAPICompilerJob) = KernelState
 
-GPUCompiler.method_table_view(job::oneAPICompilerJob) =
-    GPUCompiler.StackedMethodTable(job.world, method_table, SPIRVIntrinsics.method_table)
+function GPUCompiler.method_table_view(job::oneAPICompilerJob)
+    if job.config.target.supports_fp64
+        parent = SPIRVIntrinsics.method_table
+    else
+        # keep single-precision math that Base computes in Float64 out of double precision
+        parent = GPUCompiler.StackedMethodTable(job.world, SPIRVIntrinsics.method_table,
+                                                GPUToolbox.Overlays.float64_overrides)
+    end
+    GPUCompiler.StackedMethodTable(job.world, method_table, parent)
+end
 
 # filter out OpenCL built-ins
 # TODO: eagerly lower these using the translator API
diff --git a/src/device/quirks.jl b/src/device/quirks.jl
index 1d151f1b..c62bc57e 100644
--- a/src/device/quirks.jl
+++ b/src/device/quirks.jl
@@ -65,8 +65,4 @@ end
         @inline
         return checkindex(Bool, eachindex(IndexLinear(), v), i)
     end
-
-    # Less accurate division for Float32 than Base Julia which relies on Float64
-    # https://github.com/JuliaLang/julia/pull/49637
-    @device_override Base.div(x::Float32, y::Float32) = trunc(x / y)
 end
diff --git a/src/oneAPI.jl b/src/oneAPI.jl
index 7baebfa4..6b6db42c 100644
--- a/src/oneAPI.jl
+++ b/src/oneAPI.jl
@@ -4,6 +4,7 @@ using GPUArrays
 using Adapt
 
 using GPUCompiler
+import GPUToolbox
 
 import ExprTools