Software Float64 emulation, and job-scoped device libraries

[maleadt]

This PR makes ordinary Float64 code compile for GPU targets that have no double-precision
support, such as Metal and some oneAPI devices, by emulating the arithmetic in software.
It also generalizes the runtime-library machinery that this builds on, so that back-ends
can link additional Julia-implemented libraries into individual compilation jobs.

Heavily LLM assisted (Fable 5.1, Astra 6), so this will need a bunch of testing and review before it's viable.

How it works

The emulation operates on the final LLVM module, after all Julia code has been compiled and
linked, and is therefore transparent to Julia code: any Float64 code that works on the CPU
works on the device, including Base's pure-Julia math functions (exp, log, sin, ^,
cbrt, ...), ComplexF64, tuples and structs containing Float64, and so on.

It happens in two steps, right before optimization (src/softfloat/legalize.jl):

1. Outlining. Every floating-point operation on double values (fadd, fmul, fcmp,
   fptosi, llvm.sqrt.f64, ...) is replaced by 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).
2. Rebuilding. The module is cloned with every double type replaced by i64, in
   function signatures, aggregates, vectors, globals, constants and typed attributes
   (byval etc.). Constants keep their exact bit patterns. This turns the placeholders into
   declarations of the actual emulation routines, which are then linked in.

The routines themselves (src/softfloat/binary64/) are plain Julia functions operating on
UInt64 bit patterns, ported from metal-softfloat
(which derives from Berkeley SoftFloat). They cover addition, multiplication, division,
square root, FMA, comparisons, min/max, integral rounding, and conversions to and from
integers, Float32 and Float16. Subtraction is addition of a negated operand, and
negation, abs and copysign are single bit operations that inline into user code. The
expensive routines are @noinline, so every kernel contains at most one copy of each.

For example, Metal.code_llvm(+, (Float64, Float64)) gives:

define i64 @julia_+(i64 %"x::Float64", i64 %"y::Float64") {
top:
  %0 = call fastcc i64 @gpu_softfloat_add64(i64 %"x::Float64", i64 %"y::Float64")
  ret i64 %0
}

while -(x::Float64) compiles to xor i64 %x, -9223372036854775808, and x * 2.0 + 1.0 to
two calls with the constants' bit patterns as immediate operands.

Semantics: round-to-nearest-even, gradual underflow, signed zeros, infinities, quiet
comparisons, and a canonical quiet NaN wherever an operation produces NaN. Not supported:
floating-point exception flags, other rounding modes, Float64 atomics and frem (Julia's
rem does not use it).

One piece of Base does not compile as-is: the Payne-Hanek argument reduction used by
sin/cos/tan for arguments above 2^20·π/2 uses UInt128 arithmetic, which most GPU
back-ends cannot lower (this is unrelated to Float64 emulation). SoftFloat.paynehanek is
an equivalent implementation using pairs of UInt64, bit-for-bit identical to Base's, for
back-ends to override Base.Math.paynehanek with.

Device-library providers

The runtime library is compiled from Julia methods into per-function relocatable bitcode,
cached with their CodeInstances, and linked into every kernel. This PR generalizes that
machinery into providers: a back-end returns providers from
GPUCompiler.device_library_providers(job), each listing its methods as
DeviceLibraryMethods. A provider's library is linked after finish_linked_module!, right
after its prepare_device_library! hook has had a chance to rewrite the module (which is
where the legalization above runs). Libraries are loaded only when a kernel references one
of their exports, so enabling a provider for every job is cheap, and they are cached per
provider and runtime configuration, sharing the runtime's bitcode cache and invalidation.

Nothing is registered globally: loading a package that defines a provider has no effect on
jobs that do not select it. A back-end enables emulation with one method:

GPUCompiler.device_library_providers(::MetalCompilerJob) =
    (GPUCompiler.SoftFloat.SoftFloat64Provider(),)

plus a method-table override for Base.Math.paynehanek. Since the selection is a function
of the compiler configuration, cached compilation results are keyed correctly.

Validation

• test/softfloat.jl: over a million bit-exact comparisons of every routine against native
  CPU arithmetic, including edge cases (subnormals, signed zeros, infinities, NaNs, every
  Float16 and both sides of every Float16 rounding midpoint), and paynehanek against
  Base's.
• test/softfloat/legalize.jl: LLVM-level tests of the legalization (aggregates, vectors,
  globals, nested constant arrays, typed attributes, metadata, opaque-pointer storage), and
  that unknown external ABIs are rejected rather than silently changed.
• test/softfloat/native.jl, test/device_library.jl: execution through the native JIT,
  provider selection, and library invalidation on method redefinition.
• test/spirv/softfloat.jl: a kernel compiled with supports_fp64=false through LLVM's
  SPIR-V back-end validates and has neither OpTypeFloat 64 nor the Float64 capability.
  (The Khronos translator rejects the odd integer widths the optimizer introduces, e.g.
  i11; actual oneAPI execution has not been tested.)
• On Metal, the complete GPUArrays test suite passes with Float64 and ComplexF64,
  on Julia 1.10, 1.11 and 1.12 (LLVM 15 through 18, typed and opaque pointers).

Performance

See JuliaGPU/Metal.jl#955 for measurements. In short, emulated operations cost roughly 2-10x a
native Float32 broadcast on an M1, and dependent chains of additions run at about a
tenth of the native Float32 rate. Code size stays bounded thanks to the out-of-line
routines: a kernel using exp, sqrt and sin grows by a few hundred lines of AIR.

Provenance

src/softfloat/binary64/ is a Julia port of metal-softfloat (MIT), whose tables and
algorithms derive from Berkeley SoftFloat 3e (BSD-3-Clause); paynehanek.jl is adapted from
Julia's base/special/rem_pio2.jl (MIT). The notices are in LICENSES/.

[codecov]

Codecov Report

❌ Patch coverage is 94.09142% with 53 lines in your changes missing coverage. Please review.
✅ Project coverage is 87.25%. Comparing base (d5581b1) to head (3452534).

Files with missing lines	Patch %	Lines
src/softfloat/legalize.jl	90.24%	43 Missing ⚠️
src/softfloat/binary64/types.jl	87.80%	5 Missing ⚠️
src/rtlib.jl	92.59%	2 Missing ⚠️
src/softfloat/uint128.jl	93.93%	2 Missing ⚠️
src/driver.jl	94.11%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #926      +/-   ##
==========================================
+ Coverage   85.74%   87.25%   +1.51%     
==========================================
  Files          29       39      +10     
  Lines        5619     6506     +887     
==========================================
+ Hits         4818     5677     +859     
- Misses        801      829      +28     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[maleadt]

Depends on llvm/llvm-project#208026, which I'll add in JuliaPackaging/Yggdrasil#14727