safe_exp_polynomial Function

   elemental function safe_exp_polynomial(x) result(y)
      !! Public only for the unit-test suite (no production module imports it);
      !! ignore when developing production code in other modules.
      !! Cody-Waite range reduction `x = k*ln(2) + r`, |r| ≤ ln(2)/2,
      !! followed by a 14-term Horner Taylor on `exp(r)` and an
      !! IEEE-exact `2^k` multiply via `scale()`.  Max relative
      !! error ≈ 2-3 ULP across [-709, 709].
      !!
      !! Edge guards: overflow → `huge(y)`, underflow → `0`,
      !! `exp(0) == 1.0` exact.
      real(wp), intent(in) :: x
      real(wp) :: y
      real(wp) :: k_real, r, p
      integer :: k_int

      if (x >= EXP_OVERFLOW) then
         y = huge(y)
         return
      end if
      if (x <= EXP_UNDERFLOW) then
         y = 0.0_wp
         return
      end if

      k_real = anint(x*INV_LN2)
      k_int = nint(k_real)
      r = (x - k_real*LN2_HI) - k_real*LN2_LO

      p = EXP_C13
      p = p*r + EXP_C12
      p = p*r + EXP_C11
      p = p*r + EXP_C10
      p = p*r + EXP_C9
      p = p*r + EXP_C8
      p = p*r + EXP_C7
      p = p*r + EXP_C6
      p = p*r + EXP_C5
      p = p*r + EXP_C4
      p = p*r + EXP_C3
      p = p*r + EXP_C2
      p = p*r + EXP_C1
      p = p*r + EXP_C0

#ifdef LFORTRAN_PASSING
      ! LFortran 0.64 runtime SCALE(x,i) returns 0 for negative i (integer
      ! 2**i); use the exact real power-of-two multiply instead.  Bit-identical
      ! to scale() for powers of two on conforming compilers, which keep scale().
      y = p*2.0_wp**k_int
#else
      y = scale(p, k_int)
#endif
   end function safe_exp_polynomial