efp_decompose_impl Subroutine

   pure subroutine efp_decompose_impl(r, e1, e2, e3, e4, e5, e6, epoison)
      !! In-module `!$acc routine seq` duplicate of `rdb_efp::efp_decompose`
      !! -- six SCALAR bin outputs (not an `int64(6)` array) so the call
      !! sites below can accumulate directly into
      !! `reduction(+:e1..e6,epoison)` clauses (OpenACC has no portable
      !! array reduction).  Duplicated rather than called from `rdb_efp`
      !! because NVHPC's device codegen does not inline a
      !! `pure !$acc routine seq` helper across a module boundary (CLAUDE.md
      !! Gotchas); `test_efp_impl_matches_canonical` (`RDB_ENABLE_TESTING`-
      !! gated) pins this copy bin-for-bit against the canonical procedure
      !! over the same magnitude table -- `compute_ice_totals`'s docstring
      !! documents the identical pattern for `ice_cell_concentration_impl`.
      !!
      !! `epoison` is the device-reduction twin of `rdb_efp::efp_t%poison`:
      !! 1 when `r` is NaN, +-Inf, or exceeds bin 1's representable
      !! ceiling, 0 otherwise -- summed by the caller's OWN
      !! `reduction(+:epoison)` across the k-slab, and from there folded
      !! into the running `efp_t%poison` counter (`compute_total_h_efp`
      !! etc.) so a non-finite console summand (a NaN'd `h_layer` or
      !! `u_face_x_layer`, the diagnosed vcm_rx0_040_lagrangian failure
      !! mode) makes the FINAL reduced total read as NaN via
      !! `efp_to_real`, instead of the zeroed/saturated bins below
      !! silently reconstructing a plausible finite number -- the
      !! console's own NaN panic (`rdb_console_stats.F90`'s
      !! `panic_on_nan`) sees the poisoned value it is supposed to.
      !! `ieee_is_finite` (not a `>=`/`<` comparison chain) is the guard,
      !! per CLAUDE.md's NaN-blind-if/else-clamp gotcha: it is the
      !! dedicated bit-pattern test and is not subject to `-fast`
      !! reassociation.
      !$acc routine seq
      real(real64), intent(in) :: r
      integer(int64), intent(out) :: e1, e2, e3, e4, e5, e6
      integer(int64), intent(out) :: epoison
      real(real64) :: rs, s
      real(real64), parameter :: PR1 = 2.0_real64**(2*EFP_PREC_WIDTH)
      real(real64), parameter :: PR2 = 2.0_real64**(1*EFP_PREC_WIDTH)
      real(real64), parameter :: PR3 = 1.0_real64
      real(real64), parameter :: PR4 = 2.0_real64**(-1*EFP_PREC_WIDTH)
      real(real64), parameter :: PR5 = 2.0_real64**(-2*EFP_PREC_WIDTH)
      real(real64), parameter :: PR6 = 2.0_real64**(-3*EFP_PREC_WIDTH)
      real(real64), parameter :: IPR1 = 1.0_real64/PR1
      real(real64), parameter :: IPR2 = 1.0_real64/PR2
      real(real64), parameter :: IPR3 = 1.0_real64/PR3
      real(real64), parameter :: IPR4 = 1.0_real64/PR4
      real(real64), parameter :: IPR5 = 1.0_real64/PR5
      real(real64), parameter :: IPR6 = 1.0_real64/PR6
      real(real64), parameter :: MAX_E1 = real(huge(0_int64), real64)
         !! Same bin-1 ceiling as `rdb_efp::efp_decompose`'s `MAX_E1`.

      e1 = 0_int64
      e2 = 0_int64
      e3 = 0_int64
      e4 = 0_int64
      e5 = 0_int64
      e6 = 0_int64
      epoison = 0_int64

      if (.not. ieee_is_finite(r)) then
         ! Covers NaN and +-Inf alike: zeroed bins (matching the canonical
         ! `efp_decompose`'s NaN branch) plus the poison flag -- never an
         ! `int(NaN, int64)` conversion, which is compiler-undefined and is
         ! exactly how a NaN summand used to launder into "0.000" on the
         ! console (see rdb_efp's module docstring / CLAUDE.md Gotchas).
         epoison = 1_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*IPR1 >= MAX_E1) then
         epoison = 1_int64
         return
      end if

      e1 = int(s*aint(rs*IPR1), int64)
      rs = rs - real(abs(e1), real64)*PR1
      e2 = int(s*aint(rs*IPR2), int64)
      rs = rs - real(abs(e2), real64)*PR2
      e3 = int(s*aint(rs*IPR3), int64)
      rs = rs - real(abs(e3), real64)*PR3
      e4 = int(s*aint(rs*IPR4), int64)
      rs = rs - real(abs(e4), real64)*PR4
      e5 = int(s*aint(rs*IPR5), int64)
      rs = rs - real(abs(e5), real64)*PR5
      e6 = int(s*aint(rs*IPR6), int64)
   end subroutine efp_decompose_impl