safe_pow_polynomial Function

public elemental function safe_pow_polynomial(base, exponent) result(y)

Public only for the unit-test suite (no production module imports it); ignore when developing production code in other modules. base^exponent via the identity base^exponent = exp(exponent * log(base)) Routes through the polynomial log and exp paths so the result is deterministic across builds.

Domain: * base > 0 — proper case; returns the polynomial value. * base == 0 — returns 0 if exponent > 0, huge() if exponent < 0, 1 if exponent == 0 (the conventional IEEE-754 pow(0, 0) = 1 rule). * base < 0 — non-integer exponent is undefined for reals; returns -huge() as a sentinel. Callers that need integer base**n for negative base should use repeated multiplication, not safe_pow.

Accuracy: composes two ~ULP-level polynomial passes plus one multiplication, so expect ~5-10 ULP relative error.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: base
real(kind=wp), intent(in) :: exponent

Return Value real(kind=wp)


Calls

proc~~safe_pow_polynomial~~CallsGraph proc~safe_pow_polynomial safe_pow_polynomial proc~safe_exp_polynomial safe_exp_polynomial proc~safe_pow_polynomial->proc~safe_exp_polynomial proc~safe_log_polynomial safe_log_polynomial proc~safe_pow_polynomial->proc~safe_log_polynomial

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: log_base

Source Code

   elemental function safe_pow_polynomial(base, exponent) result(y)
      !! Public only for the unit-test suite (no production module imports it);
      !! ignore when developing production code in other modules.
      !! `base^exponent` via the identity
      !!   base^exponent = exp(exponent * log(base))
      !! Routes through the polynomial log and exp paths so the
      !! result is deterministic across builds.
      !!
      !! Domain:
      !!   * `base > 0` — proper case; returns the polynomial value.
      !!   * `base == 0` — returns 0 if `exponent > 0`, `huge()` if
      !!     `exponent < 0`, `1` if `exponent == 0` (the conventional
      !!     IEEE-754 `pow(0, 0) = 1` rule).
      !!   * `base < 0` — non-integer exponent is undefined for
      !!     reals; returns `-huge()` as a sentinel.  Callers that
      !!     need integer `base**n` for negative `base` should use
      !!     repeated multiplication, not `safe_pow`.
      !!
      !! Accuracy: composes two ~ULP-level polynomial passes plus
      !! one multiplication, so expect ~5-10 ULP relative error.
      real(wp), intent(in) :: base, exponent
      real(wp) :: y
      real(wp) :: log_base

      if (base > 0.0_wp) then
         log_base = safe_log_polynomial(base)
         y = safe_exp_polynomial(exponent*log_base)
      else if (base == 0.0_wp) then
         if (exponent > 0.0_wp) then
            y = 0.0_wp
         else if (exponent < 0.0_wp) then
            y = huge(y)
         else
            y = 1.0_wp
         end if
      else
         y = -huge(y)
      end if
   end function safe_pow_polynomial