Skip to content

Add Float64-free rem2pi for Float32 and Float16 - #40

Merged
maleadt merged 2 commits into
JuliaGPU:mainfrom
Beforerr:rem2pi-float32
Nov 5, 2026
Merged

maleadt merged 2 commits into
JuliaGPU:mainfrom
Beforerr:rem2pi-float32

Conversation

@Beforerr

@Beforerr Beforerr commented Nov 5, 2026 •

Copy link
Copy Markdown
Contributor

mod2pi and rem2pi on Float32 and Float16 don't compile on devices without double precision:

using Metal
mod2pi.(MtlArray(Float32[-1, 1, 7]))  # InvalidIRError: unsupported use of double value

Base reduces these types in Float64, on every supported Julia version including current master, and no back-end overrides them:

rem2pi(x::Float32, r::RoundingMode) = Float32(rem2pi(Float64(x), r))
rem2pi(x::Float16, r::RoundingMode) = Float16(rem2pi(Float64(x), r))

This adds overrides to float64_overrides for RoundDown, RoundUp, RoundToZero and RoundNearest. The Float32 path computes the fractional part of x / 2π with integer Payne-Hanek: the 24-bit mantissa times a 96-bit window of the bits of 1/2π, using 32×32→64-bit products, without Float64, UInt128 or a lookup table. The high 64 bits of that product are then scaled by 2π. The window needs 96 bits because x / 2π comes as close as 2^-29.9 to an integer (at x = 3.0916716f29); a 64-bit window gives wrong results. Float16 goes through the Float32 path. rem2pi(::Int32, r) is not covered, because Base returns Float64 for it.

Verification, against Base over all 2^32 Float32 inputs, in each of the four modes:

  • CPU: within 0.51 ulp. RoundToZero and RoundNearest differ from Base on 2 and 8 inputs, each by exactly 0.5 ulp. No signed zeros differ, and every Float16 result matches Base.
  • Metal (Apple M1, Julia 1.13, Metal.jl main): the same results as on the CPU, also for subnormal inputs. The example above matches the host.

Beforerr and others added 2 commits November 4, 2026 18:23
Base computes rem2pi(::Float32, r) and rem2pi(::Float16, r) in Float64, so
mod2pi fails to compile on devices without double precision. Reduce in
fixed point with integer Payne-Hanek instead, which needs no Float64 and
stays accurate up to floatmax(Float32).

Compared with Base over all 2^32 Float32 values on the CPU, results are
within 0.51 ulp in all four rounding modes, and all Float16 results match.
The low 32 bits of the 96-bit product never affect the result: over all
finite Float32 inputs, results are bit-identical without them, for both the
fraction and its negation. Dropping them removes the separate negation and
the extra word during normalization.

Also correct the stated worst case (x / 2π comes within 2^-29.9 of an
integer, not 2^-30) and test the next two closest inputs.
@maleadt

maleadt commented Nov 5, 2026

Copy link
Copy Markdown
Member

Thanks!

@maleadt
maleadt merged commit 5c4f197 into JuliaGPU:main Nov 5, 2026
9 of 10 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants