efp_decompose Subroutine

   pure subroutine efp_decompose(r, e, is_nan, is_ovf)
      !! Greedy sign-magnitude fixed-point decomposition of `r` into six
      !! `int64` bins of weight `pr(n)`, `n = 1..6`.  Unrolled (no loop, no
      !! array indexing over `pr`/`I_pr`) so this is a template a
      !! `!$acc routine seq` in-module duplicate can mirror exactly (see
      !! `rdb_ocean_console_stats::efp_decompose_impl`, which pins against
      !! this procedure in `test_efp_impl_matches_canonical`).
      !!
      !! `is_nan` is set (and `e = 0`) for a NaN input, using
      !! `ieee_is_nan`-equivalent structural comparison so the module stays
      !! independent of `ieee_arithmetic`'s import list elsewhere.
      !! `is_ovf` is set when `|r| >= pr(1) * huge(1_int64)` -- the largest
      !! magnitude bin 1 can represent -- and `e` is truncated to that
      !! bound rather than silently wrapping. `+-Inf` is caught by the SAME
      !! `is_ovf` flag via an explicit `ieee_is_finite` check (not by
      !! falling through the magnitude comparison below): a relaxed-FP
      !! build (`-fast`, no `-Kieee`) is not guaranteed to keep a plain
      !! `>=` comparison against Infinity well-behaved under aggressive
      !! reassociation, the same hazard CLAUDE.md's NaN-blind if/else-clamp
      !! gotcha documents for `if/else` chains -- `ieee_is_finite` is the
      !! dedicated bit-pattern test and is not subject to it.
      real(real64), intent(in) :: r
      integer(int64), intent(out) :: e(EFP_DIGITS)
      logical, intent(out) :: is_nan
      logical, intent(out) :: is_ovf
      real(real64) :: rs, s
      real(real64), parameter :: MAX_E1 = real(huge(0_int64), real64)
         !! `huge(1_int64)` widened to real64 (rounds to the nearest
         !! representable double, ~9.223372036854776e18) -- the largest
         !! magnitude bin 1 can safely hold; used as the overflow ceiling
         !! below.

      is_nan = ieee_is_nan(r)
      is_ovf = .false.
      e = 0_int64
      if (is_nan) return

      if (.not. ieee_is_finite(r)) then
         ! +-Inf: bin 1's ceiling is exceeded by construction. `sign(MAX_E1,
         ! r)` is well-defined for an infinity (its sign bit is real).
         ! Callers (`efp_from_real`) fold `is_ovf` into the poison counter,
         ! so the saturated bin below is never actually read back out --
         ! it exists only so `efp_carry`/`efp_regularize` see a normal
         ! int64, not an attempted NaN-to-integer conversion.
         is_ovf = .true.
         e(1) = int(sign(MAX_E1, r), int64)
         return
      end if

      s = 1.0_real64
      rs = r
      if (rs < 0.0_real64) then
         s = -1.0_real64
         rs = -rs
      end if

      if (rs*EFP_IPR1 >= MAX_E1) then
         is_ovf = .true.
         e(1) = int(sign(MAX_E1, s), int64)
         return
      end if

      e(1) = int(s*aint(rs*EFP_IPR1), int64)
      rs = rs - real(abs(e(1)), real64)*EFP_PR1

      e(2) = int(s*aint(rs*EFP_IPR2), int64)
      rs = rs - real(abs(e(2)), real64)*EFP_PR2

      e(3) = int(s*aint(rs*EFP_IPR3), int64)
      rs = rs - real(abs(e(3)), real64)*EFP_PR3

      e(4) = int(s*aint(rs*EFP_IPR4), int64)
      rs = rs - real(abs(e(4)), real64)*EFP_PR4

      e(5) = int(s*aint(rs*EFP_IPR5), int64)
      rs = rs - real(abs(e(5)), real64)*EFP_PR5

      e(6) = int(s*aint(rs*EFP_IPR6), int64)
      ! Remainder `rs - |e(6)|*pr(6)` is the deterministic dropped quantum,
      ! `0 <= remainder < pr(6) = 2**-3P` -- not retained (matches MOM6:
      ! only six bins are kept).
   end subroutine efp_decompose