Back to Blog

Apple's GPU Has No 64-bit Floats. I Made It Sort 50 Million Doubles Anyway

Prateek SinghSeptember 20, 202615 min read
Apple's GPU Has No 64-bit Floats. I Made It Sort 50 Million Doubles Anyway

The Metal Shading Language has no double type, and float64 is the default number in Python, pandas and Apache Arrow. Building ArrowMetal meant getting past three walls: a GPU with no 64-bit floats, a missing 64-bit atomic add, and a Swift compiler bug that reports errors nobody threw. Here is how each one was solved, what it cost, and why a GPU that cannot add two doubles sorts 50,000,000 of them in 32 ms.

The Metal Shading Language has no double. Apple GPUs have no 64-bit floating-point hardware at all. That would be a footnote if float64 were exotic. It is the opposite: it is what Python calls a float, what pandas gives you for any column with a decimal point, and the default numeric type in Apache Arrow. I was building ArrowMetal, an Arrow compute library for the Apple silicon GPU, and the first real column anyone would hand it was a type the GPU cannot add.

Hardware doubles on an Apple GPU
0
Metal has no double type. float64 is the default number in Python, pandas and Arrow.
Error of the software square root
0 ulp
bit-identical to the CPU over a million random inputs, subnormals included
Sort 50,000,000 doubles, M4 Max
32 ms
without one floating-point instruction. Polars: 133 ms on 9.7 cores.

This is the story of three walls I hit building it and how each one was solved: a GPU with no doubles, a 64-bit atomic add that the documentation describes and the compiler refuses, and a Swift compiler bug that made a function report an error nobody threw. It is not a benchmark post. The numbers that appear are consequences of the decisions, every one is from an Apple M4 Max with 16 CPU cores and 64 GB at 50,000,000 rows unless stated, and the places where the decision costs speed are in here next to the places where it does not.

Wall one: the GPU cannot hold the most common number in data

There were three ways through, and I tried the first before settling on the third.

optionwhat you getwhy I did or did not take it
Narrow to float32, use Metal's own math, widen backHardware speed. About seven correct significant digits out of sixteen.This is what the transcendental kernels did at first. A float64 column that silently carries float32 answers is wrong in a way nobody notices until it matters.
Double-float: carry each value as a pair of float32More digits than float32, still not IEEE-754 binary64. Results differ from the CPU in the last places.float64 is the default numeric type in Python. An approximate result there would be a support burden forever.
Software IEEE-754 binary64 on 64-bit integersThe same bits the CPU produces. Tens of integer instructions per operation.Taken. And then avoided wherever an operation can be done without any float arithmetic at all.
The three options for float64 on a GPU with no binary64 hardware. The decision is recorded in the project's docs/DECISIONS.md, 6 September 2026.

Software floating point on a GPU sounds like a way to make a fast chip slow. The part that made it workable is the second half of that last row. Before emulating anything, I asked of every operation: does this need float arithmetic, or only float ordering?

Most of what you do to a column never adds two numbers

Sort, argsort, compare, min, max, filter, take, slice, top-k, median, rank. None of them computes a new floating-point value. They only need to know which of two doubles is larger, or they just move bytes. And a double's ordering can be read straight off its bits, as an integer, with one trick.

A double is 64 bits. Flip them the right way and integer order is float order. sign 11 exponent bits 52 significand bits key = (bits & SIGN) ? ~bits : (bits | SIGN) the double its unsigned 64-bit key key, drawn to scale -inf 0x000fffffffffffff -1.5 0x4007ffffffffffff -5e-324 0x7ffffffffffffffe 0.0 and -0.0 0x8000000000000000 5e-324 0x8000000000000001 1.5 0xbff8000000000000 +inf 0xfff0000000000000 NaN (every one) 0xfff0000000000001 -0.0 is folded onto 0.0 and every NaN onto one value after +inf first, which is the order Apache Arrow defines.
The order-preserving key ArrowMetal builds for every float64 before sorting. Keys shown are computed by the Python snippet below; the bar is the key's top 16 bits to scale.

An IEEE-754 double is a sign bit, eleven exponent bits and fifty-two significand bits, in that order. For positive numbers that layout is already sorted: a bigger exponent is a bigger number, and within one exponent a bigger significand is a bigger number, so comparing the raw bits as an unsigned integer gives the right answer. Negative numbers are the same thing mirrored. So: set the top bit on positives to lift them above every negative, and invert every bit of the negatives so that the larger magnitude sorts first. Two special cases are folded in first, because Arrow defines them: -0.0 equals 0.0, and every NaN is one value that sorts after +inf. This is the kernel, from the repository:

// Sources/ArrowMetal/Kernels/SortSource.swift, the float64 key, as Metal Shading Language
ulong b = (ulong)a[i];                                   // the double's 64 bits, read as an integer
if ((b & 0x7FFFFFFFFFFFFFFFul) == 0ul) b = 0ul;          // -0.0 == +0.0
if (nan) b = 0x7FF0000000000001ul;                       // every NaN is one value, just after +inf
ulong k = (b & 0x8000000000000000ul) ? ~b                // negative: invert, so bigger magnitude sorts first
                                     : (b | 0x8000000000000000ul);   // positive: lift above every negative

You do not have to take that on trust. Here is the same mapping in twelve lines of Python, checked against Python's own float sort on two hundred thousand random doubles plus the awkward ones:

import struct, math, random

SIGN, M = 1 << 63, (1 << 64) - 1
def key(x):
    b = struct.unpack("<Q", struct.pack("<d", x))[0]     # the double's 64 bits as an unsigned integer
    if b & (M >> 1) == 0: b = 0                          # -0.0 == +0.0
    if (b & (M >> 1)) > 0x7FF0000000000000: b = 0x7FF0000000000001
    return (~b & M) if b & SIGN else (b | SIGN)

xs = [random.uniform(-1e300, 1e300) for _ in range(200_000)]
xs += [0.0, -0.0, math.inf, -math.inf, 5e-324, -5e-324, 1.0, -1.0]
assert sorted(xs) == sorted(xs, key=key)                 # passes: integer order of the keys is float order

Once every double is an unsigned 64-bit key, sorting is a radix sort over integers, which is the kind of work a GPU is built for: eight passes of counting and scattering, thousands of lanes wide. The result is that a GPU with no float64 sorts 50,000,000 float64 values in 32.12 ms on an M4 Max. Polars, the fastest CPU library on that row, takes 133 ms on 9.7 cores. The median of the same column takes 4.61 ms against 178 ms for Polars lazy, because a median does not need a sort either: a radix select histograms the top byte of the key, keeps the one bin the answer lives in, and repeats. The sort runs without one floating-point instruction.

The rest: IEEE-754, rebuilt out of integers

Addition cannot be dodged that way. For arithmetic, DoubleMath.swift carries a software binary64 over 64-bit integers: unpack sign, exponent and significand, align, operate, normalise, round to nearest with ties to even, repack. Subnormals, signed zeros, infinities and NaN propagation included. A few things I learned writing it that I have not seen written down elsewhere:

  • Apple's GPU arithmetic units are 32 bits wide. A 64-bit add or compare is two instructions. So the exponent travels as a 32-bit int, and the 53 by 53 bit significand product in a multiply is done as four 32 by 32 products with the partial products shared rather than computed twice.
  • A data-dependent loop costs every lane, not one. Normalising after a subtraction that cancels used to be a shift loop that could run fifty times. On a GPU the whole SIMD group waits for the slowest lane, so that loop became one count-leading-zeros and one shift.
  • Division gets one hardware instruction as a hint. d_div seeds a Newton reciprocal with a single float32 hardware division, then the exact 128-bit remainder settles the last bit. Square root is extracted digit by digit in integers. Both are correctly rounded, not close.

"Correctly rounded" is a claim that can be checked, so the tests check it: the arithmetic is held to Swift's own Double bit for bit, and each transcendental is compared with the CPU over a million random inputs drawn across its whole domain.

functioninputs sampled against the CPUlargest error
add, subtract, multiply, divideheld to Swift's own Double, bit for bit0 (bit-identical)
sqrt10⁶ random bit patterns over the whole exponent range, subnormals, 1e-320 and 1e308 included, plus perfect squares and their neighbours0 (bit-identical)
exp-745.2 to 709.78, plus the subnormal-result and near-overflow edges1 ulp
ln, log2, log105e-324 to 1.8e308, plus near 1 and the subnormals1 ulp
power10⁶ random pairs, 5·10⁵ of them negative bases with integer exponents1 ulp
trigonometric familyseparate kernel and budget; the test asserts 62 to 5 ulp
Measured accuracy of ArrowMetal 0.1.0's software binary64 against the CPU, 6 September 2026, from docs/DESIGN.md. One ulp is one unit in the last place of a double. The test suite asserts a 2-ulp bound on the transcendentals and bit equality on sqrt, so a regression of one ulp changes a number here before it fails a test.

The hardest one was pow. A result good to one ulp needs the product y·log2(x), which can reach 1024 in magnitude, to be accurate to 2⁻⁶¹: more than a double can hold. So log2(x) is carried as an unevaluated high and low pair, and y is split the same way so that the leading product is exact. That is fdlibm's layout, and its coefficients and constants are reused verbatim; the three logarithms fall out of the same reduction.

Here it is running on the GPU, from the PyPI wheel on an M4 Max. The comments are the actual outputs:

import math, pyarrow as pa, arrowmetal as am
am.device_name()                                    # 'Apple M4 Max'

x = am.array(pa.array([0.1, 1e308, 2.0, -0.0], pa.float64()))
y = am.array(pa.array([0.2, 1e308, 3.0,  5.0], pa.float64()))

(x + y).to_arrow().to_pylist()      # [0.30000000000000004, inf, 5.0, 5.0]   the CPU's bits, overflow included
x.sqrt().to_arrow()[2].as_py() == math.sqrt(2.0)    # True: bit-identical, not "close"
x.ln().to_arrow()[2].as_py(), math.log(2.0)         # (0.6931471805599453, 0.6931471805599453)

z = am.array(pa.array([3.5, float("nan"), -0.0, None, float("-inf"), 1e-320], pa.float64()))
z.sort().to_arrow().to_pylist()     # [-inf, -0.0, 1e-320, 3.5, nan, None]   NaN after +inf, nulls last

What it cost

For add and multiply, almost nothing: those rows run at the roughly 390 GB/s that any single pass over the data reaches on this machine, so they are limited by memory, not by the forty integer instructions. Divide is 3.83 ms against 4.38 ms for Polars lazy on 13.2 cores, and sqrt is 3.92 ms against 3.42 ms.

The transcendentals are where the bill arrives. Forty-odd software operations per element, three to four GPU instructions each, against one vectorised hardware instruction on twelve to fifteen CPU cores: ln takes 81.31 ms where Polars lazy takes 8.46 ms, and sin takes 150 ms against 23.28 ms. The earlier float32-detour ln ran at 250 GB/s; the binary64 one runs at 10 GB/s. That is a 25x drop in throughput for nine more correct digits. I think it is the right trade for a library whose whole claim is that a float64 column means float64, and those rows are on the project's to-improve list with this cause written next to them. If your workload is logarithms over doubles, use Polars. The float32 kernels are untouched and still take the hardware path.

float64 on a GPU with no float64: 50,000,000 rows, Apple M4 Max, against the fastest CPU idiom 0.1x1x10x100x median (quantile 0.5) vs Polars lazy, 6.8 cores 38.7x argsort vs DuckDB, 12.6 cores 8.60x cumulative sum vs pyarrow, 1 core 7.63x sort vs Polars, 9.7 cores 4.14x divide vs Polars lazy, 13.2 cores 1.14x sqrt vs Polars lazy, 12.7 cores 0.87x sin vs Polars lazy, 15.1 cores 0.16x ln vs Polars lazy, 14.2 cores 0.10x
float64 operations in ArrowMetal 0.1.0 against the fastest CPU idiom on each row, log scale. Apple M4 Max, 16 CPU cores, 50,000,000 rows, best of up to five calls after a warm-up. Polars 1.44.1, pyarrow 25.0.1, matrix generated 7 September 2026; DuckDB 1.5.5, run 12 September 2026. Grey bars are rows to improve.

Wall two: the 64-bit atomic add that is in the documentation and not in the compiler

A GPU group-by is thousands of threads adding into the same table of sums at once, so every add has to be atomic. Arrow's sums are 64-bit. Apple's Metal feature set tables describe the Apple9 GPU family as having the full set of 64-bit atomic operations. The shading-language headers admit a 64-bit type only for atomic min and max: no 64-bit atomic add, compare-exchange, load or store compiles, whatever language version is asked for. That is measured on an M4 Max with macOS 26.6.2, I filed it with Apple through Feedback Assistant on 20 September 2026 as FB24858110, with the compile probe, its output and the three kernels that carry the workarounds attached, and the row is in the project's upstream tracker. I do not know which of the two is the intended behaviour. Either way, the group-by had to work without it.

  • Sums: a 64-bit sum is a pair of 32-bit atomic adds with an explicit carry. For up to 1,024 groups each threadgroup keeps a private table in fast threadgroup memory and merges it into the device table once.
  • Min and max: done with 32-bit atomics in two passes. Each value maps to the same order-preserving 64-bit key as the sort. Pass one takes the extreme of the key's high 32 bits with 32-bit atomics. Pass two takes the extreme of the low 32 bits among only the rows whose high word already equals the winner. Some row attains the winning high word, and among those the smallest low word is the overall minimum. No sort of the key column, no 64-bit atomic.

With those two in place, a sum of 50,000,000 rows into 100,000 groups takes 7.41 ms on the M4 Max, against 46.39 ms for pyarrow's Acero on 11.7 cores.

Wall three: the error nobody threw

The strangest one was a crash in the array initialiser, in release builds only. A function was returning an error. The function contained no reachable throw. In a debug build it passed.

It turned out to be Swift 6.3.3. In a generic throws function compiled at -O in its own module, the closure passed to withUnsafeBytes uses the register that carries Swift's error result as scratch space around an Objective-C message send, and then exits by tail-calling memmove without restoring it. The caller looks at the error register, finds a leftover value, believes an error was thrown, and crashes trying to retain it. The issue had already been reported upstream by someone else; I reduced it to a two-module reproducer with no Metal in it and added that to swiftlang/swift#90477 on 8 September 2026, where it is open. ArrowMetal's fix is unglamorous: MetalArray.init uses a plain element loop instead, and the findings notes now say to run every test in release as well as debug.

What stays with me is how that one was found. ArrowMetal's test suite is differential: 39,069 cases comparing its answers with pyarrow's over 45 column types, plus 769 Swift tests against a CPU oracle, all run in release. A suite built to catch my mistakes also catches everyone else's, and it has: the upstream tracker lists 18 findings in other projects, each with a test that will fail on the day upstream fixes it so the workaround can come out. That is a post of its own.

What the three walls add up to

A GPU that cannot represent a double sorts them 4.1x sooner than the fastest CPU library, finds their median 38.7x sooner, adds them with the CPU's exact bits, and takes ten times longer than Polars to compute their logarithm. All four of those are the same decision seen from different sides: never approximate, and never do float arithmetic when ordering will do.

ArrowMetal is an independent Apache-2.0 project that implements Apache Arrow, listed on Arrow's Powered By page; it is version 0.1.0 and every number here is from one machine. To try the float64 path yourself:

pip install arrowmetal          # macOS 14+ on Apple silicon, Python 3.10+; pyarrow comes with it
pip install polars duckdb       # optional: the two bridges used below

The launch post covers what the library is, and the benchmark post has the full sort and group-by tables against Polars and DuckDB.

FAQ: double precision, float64 and 64-bit atomics in Metal

Does Metal support double precision (float64)?

No. The Metal Shading Language has no double type and Apple GPUs have no binary64 hardware. A library that needs float64 on an Apple GPU has three choices: narrow to float32 and accept about seven correct digits, emulate with pairs of float32, or implement IEEE-754 binary64 in software on 64-bit integers. ArrowMetal does the third, correctly rounded and bit-exact against the CPU for add, subtract, multiply, divide and sqrt.

How do you sort float64 values on a GPU that has no float64?

By never doing float arithmetic. A double's 64 bits are read as an unsigned integer and mapped to an order-preserving key: negative values are bit-inverted and positive values get the top bit set, after folding -0.0 onto 0.0 and every NaN onto one value after +inf. A radix sort of those integer keys is a sort of the doubles. On an Apple M4 Max, 50,000,000 float64 values sort in 32 ms this way, against 133 ms for Polars on 9.7 cores.

How accurate is software float64 on the GPU?

In ArrowMetal 0.1.0, add, subtract, multiply, divide and sqrt are correctly rounded and bit-identical to the CPU. exp, ln, log2, log10 and power measure 1 ulp against a 2-ulp asserted bound over a million random inputs each, and the trigonometric family measures 2 to 5 ulp against an asserted 6.

Is float64 math slow on an Apple GPU?

It depends on the operation. Add and multiply are memory bound, so the software arithmetic is nearly invisible. Divide is 1.14x and sqrt 0.87x against Polars lazy at 50,000,000 rows on an M4 Max. Transcendentals cost forty-odd software operations per element: ln is 0.10x and sin 0.16x against Polars lazy on 14 to 15 cores. Operations that need only ordering, such as sort, argsort and median, are 4.1x to 38.7x ahead.

Does Metal have 64-bit atomic operations?

Only atomic min and max, as measured on an M4 Max with macOS 26.6.2: no 64-bit atomic add, compare-exchange, load or store compiles, although Apple's feature set tables describe the Apple9 family as having the full set. ArrowMetal builds 64-bit grouped sums from a pair of 32-bit atomic adds with an explicit carry, and grouped min and max from two 32-bit passes over the high and low words.

References & Citations

  • ArrowMetal (2026). docs/DECISIONS.md: "Float64 on the GPU: bit-pattern ordering plus software IEEE-754" and "Reductions without atomics", 6 September 2026.
  • ArrowMetal (2026). docs/DESIGN.md: "Float64 without hardware doubles", the measured ulp table, the group-by "Atomic accumulation" and "Sort-free extremes" sections, and the radix select.
  • ArrowMetal (2026). Sources/ArrowMetal/Kernels/SortSource.swift (the float64 key) and DoubleMath.swift (software binary64).
  • ArrowMetal (2026). docs/UPSTREAM.md: the Apple Metal 64-bit atomics row (measured on an M4 Max, macOS 26.6.2; filed with Apple through Feedback Assistant on 20 September 2026 as FB24858110, open) and the Swift 6.3.3 row. docs/FINDINGS.md, "Toolchain".
  • Swift. swiftlang/swift#90477, open; the two-module Metal-free reproducer was added on 8 September 2026.
  • ArrowMetal (2026). docs/BENCHMARKS_MATRIX.md, generated 7 September 2026 on an Apple M4 Max, 16 cores, 64 GB; Polars 1.44.1, pyarrow 25.0.1. docs/TO_IMPROVE.md §3. DuckDB argsort cell from docs/DUCKDB.md, DuckDB 1.5.5, 12 September 2026.
  • The two Python snippets were run on 20 September 2026 on an Apple M4 Max, Python 3.13.9, arrowmetal 0.1.0 and pyarrow 25.0.1 from PyPI; the comments are their actual outputs.

Subscribe to new posts from theaivibe.org

No spam — just new posts. One-click unsubscribe.
Share this article

Related Posts