# logApprox.py
# Andrew Davison, ad@coe.psu.ac.th, June 2026

"""
  Evaluate ln(1+x) using argument reduction and the 
  Taylor series for atanh().

  The reduction is based on:
      1+x = m * 2**k
  so that m is scaled into the interval
      sqrt(1/2) <= m <= sqrt(2)
  Set r = m-1, so that it lies in the interval
      [-0.292893..., 0.414214...]

  Rather than evaluating a logarithm series, use
  the identity
      ln(1+r) = 2*atanh(r/(2+r))
  to use a series for atanh(). Let z = r/(2+r),
  then
      atanh(z) = z + z^3/3 + z^5/5 + z^7/7 + ...
  so
      ln(1+r) = 2*(z + z^3/3 + z^5/5 + ...)

  Because the reduction guarantees |z| <= 0.171573...
  the series converges much faster than the standard series
  for ln(1+r).

  The result for the original input is constructed using
      ln(1+x) = k*ln(2) + ln(m)

  The test-rig compares the approximation and math.log1p(x)
  (i.e. ln(1+x)) over a wide range of values. 
"""

import math


LN2 = 0.6931471805599453094
SQRT2 = 1.4142135623730950488
SQRT_HALF = SQRT2/2        # sqrt(2/4) == sqrt(2)/2


def logApprox(x, maxIters=16, eps=1e-15):
  if x <= -1.0:
    raise ValueError("logApprox() requires x > -1")
  if x == 0.0:
    return 0.0

  m = 1.0 + x
  k = 0
  # m is scaled into the interval
  # sqrt(1/2) <= m <= sqrt(2)
  while m > SQRT2:
    m *= 0.5
    k += 1
  while m < SQRT_HALF:
    m *= 2
    k -= 1

  r = m - 1
  z = r/(2 + r)
  z2 = z*z

  term = z
  tot = z    # for atanh series
  denom = 3
  for _ in range(maxIters):
    term *= z2
    delta = term/denom
    tot += delta
    if abs(delta) <= eps * max(1.0, abs(tot)):
      break
    denom += 2
  return k*LN2 + 2.0*tot


# -------------------------------------
if __name__ == "__main__":
  print("       x        Approx   math.log1p(x)   Rel Error")
  print("-"*52)

  testVals = [ -0.999999, -0.999, -0.99, -0.9, -0.5, 
               -0.1, -0.01, -0.001, 0.001, 0.01, 0.1,  
                0.5, 1, 5, 10, 100, 1.0e6  ]
  for x in testVals:
    approx = logApprox(x)
    exact = math.log1p(x)
    relErr = 0.0
    if exact != 0.0:
      relErr = abs(approx - exact)/abs(exact)
    print(f"{x:10.6g} " f"{approx:13.5e} "       
          f"{exact:13.5e} " f"{relErr:11.3e}")