Ozaki scheme, maybe support FP128? - DGEMM without FP64 Arithmetic – Using FP64 Emulation and FP8 Tensor Cores with Ozaki Scheme

Internals & Design
1

This scheme was new to me:
https://arxiv.org/pdf/2508.00441

It’s about GPUs. Got me thinking maybe this would also work for FP128, on CPUs (or GPUs). [EDIT: has been done, for matrices, see my follow-up post.]

It doesn’t supports scalars (or well), pays off with a large enough matrix. Probably Julia can’t do Float64 faster on CPUs, using Float32, I suppose all CPUs have as many units for, or half as fast bandwidth wise.

@photor

2

For CPU it will be a lot slower than Quadmath.jl / libquadmath. CPUs have fast 64*64->128 multiply which is a much faster way to do 128 bit math.

3

I’m was thinking of for BLAS like in the paper (I wasn’t clear on that), I thought this scheme would be possible for Float128 too, and I see I’m right since it has already been done for CPUs:

https://dl.acm.org/doi/fullHtml/10.1145/3472456.3472493

Compared to this previous study, this paper provides the following contributions.

  • We apply the Ozaki scheme to high-precision (binary128) computations in environments where high-precision operations are not supported on hardware.
  • We present specific optimization methods to compute binary128 using binary64.
  • We present an implementation and evaluation on x86 hardware, while the previous study targeted NVIDIA GPUs.

..
Oz-b128 outperformed MP-b128 in all cases, and MP-dd was comparable to Oz-b128 when R = 64. Note that the DD arithmetic in the QD library used in MP-dd is not SIMD-optimized,

In short from the other paper:

The Ozaki scheme [28], proposed by Ozaki et al. in 2012, is a method that enables high-precision matrix multiplication to be computed using low-precision matrix multiplications.

i.e. OpenBLAS and MKL do not support Float128, but maybe Julia should, and then it can be implemented with such libraries (on Float64, or lesser on GPUs).

Julia doesn’t strictly need to support this for matrices (or scalars), since a library can. I just got exited about this relatively new research I didn’t know of; and the scheme, yes 12 years old, I guess you knew of it.

I believe the scheme can split Float128 into 3 Float64, or 5 Float32, or 10 Float16. And you use k^2 BLAS operations, so at least 9 times slower done with Float64 as opposed to “2x to 5x slower” for scalars with libquadmath, but since the scalars don’t work with BLAS you will end up faster, given large enough matrix.

I think you have in mind only the scalar (integer) case::

64-bit to 128-bit multiplication actually takes exactly 1 instruction, but on ARM64 processors, it requires exactly 2 instructions.

But that excludes unpacking floats and repacking in software, killing performance, and can’t use BLAS libraries like OpenBLAS.