safe_log_polynomial Function

public elemental function safe_log_polynomial(x) result(y)

Public only for the unit-test suite (no production module imports it); ignore when developing production code in other modules. log(x) via IEEE exponent extraction + atanh-form series.

  1. Write x = m * 2^k using fraction() / exponent(). m initially lives in [0.5, 1).
  2. √2-symmetric split: if m < √(1/2), halve k and double m, putting m in [√(1/2), √2) — keeps |m - 1| ≤ √2 - 1 ≈ 0.414.
  3. atanh substitution f = (m-1)/(m+1) (|f| ≤ 0.172), log(m) = 2*f*(1 + f²/3 + f⁴/5 + ... + f¹⁸/19).
  4. log(x) = k*ln(2) + log(m), with k*ln(2) evaluated via Cody-Waite LN2_HI / LN2_LO.

Accuracy: < 1e-15 relative error across positive doubles. log(1) == 0 exact (the polynomial collapses to 0 at f=0 and k = 0 after the √2 split).

Domain: x <= 0 returns -huge() (sentinel; callers that care about IEEE NaN semantics should guard upstream).

Arguments

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

Return Value real(kind=wp)


Called by

proc~~safe_log_polynomial~~CalledByGraph proc~safe_log_polynomial safe_log_polynomial proc~safe_pow_polynomial safe_pow_polynomial proc~safe_pow_polynomial->proc~safe_log_polynomial

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: f
integer, private :: k_int
real(kind=wp), private :: k_real
real(kind=wp), private :: m
real(kind=wp), private :: p
real(kind=wp), private :: u

Source Code

   elemental function safe_log_polynomial(x) result(y)
      !! Public only for the unit-test suite (no production module imports it);
      !! ignore when developing production code in other modules.
      !! `log(x)` via IEEE exponent extraction + atanh-form series.
      !!
      !!   1. Write `x = m * 2^k` using `fraction()` / `exponent()`.
      !!      `m` initially lives in [0.5, 1).
      !!   2. √2-symmetric split: if `m < √(1/2)`, halve k and
      !!      double m, putting `m` in [√(1/2), √2) — keeps
      !!      |m - 1| ≤ √2 - 1 ≈ 0.414.
      !!   3. atanh substitution `f = (m-1)/(m+1)` (|f| ≤ 0.172),
      !!      `log(m) = 2*f*(1 + f²/3 + f⁴/5 + ... + f¹⁸/19)`.
      !!   4. `log(x) = k*ln(2) + log(m)`, with `k*ln(2)` evaluated
      !!      via Cody-Waite `LN2_HI / LN2_LO`.
      !!
      !! Accuracy: < 1e-15 relative error across positive doubles.
      !! `log(1) == 0` exact (the polynomial collapses to 0 at f=0
      !! and `k = 0` after the √2 split).
      !!
      !! Domain: `x <= 0` returns `-huge()` (sentinel; callers that
      !! care about IEEE NaN semantics should guard upstream).
      real(wp), intent(in) :: x
      real(wp) :: y
      real(wp) :: m, f, u, p, k_real
      integer :: k_int

      if (x <= 0.0_wp) then
         y = -huge(y)
         return
      end if

      k_int = exponent(x)
      m = fraction(x)

      ! √2-symmetric range: shift [0.5, √(1/2)) up by factor 2.
      if (m < SQRT_HALF) then
         m = 2.0_wp*m
         k_int = k_int - 1
      end if
      ! m now in [√(1/2), √2);  k_int is `exponent(x)` (or one less).

      f = (m - 1.0_wp)/(m + 1.0_wp)
      u = f*f

      ! Horner over u = f²: 1 + u/3 + u²/5 + ... + u⁹/19.
      p = LOG_C9
      p = p*u + LOG_C8
      p = p*u + LOG_C7
      p = p*u + LOG_C6
      p = p*u + LOG_C5
      p = p*u + LOG_C4
      p = p*u + LOG_C3
      p = p*u + LOG_C2
      p = p*u + LOG_C1
      p = p*u + LOG_C0

      k_real = real(k_int, wp)
      y = (k_real*LN2_HI + 2.0_wp*f*p) + k_real*LN2_LO
   end function safe_log_polynomial