Re: Missing math functions

[Luis Carvalho <lexcarvalho@...>]

> That's not good enough.
> 
> The asinh and atanh formulas are known to be inaccurate near the
> origin because of smearing.  When x*x+1 is formed, digits get thrown
> away. When x is added to the sqrt of that, more digits get thrown away.
> When log of something close to 1 is calculated, 1 is subtracted in the
> process.  This cancellation does not by itself introduce a new error,
> but exposes the fact that the digits are already irretrievably gone.
<snip>

This is slightly off-topic, but here's a pure Lua implementation of inverse
hyperbolic trig functions. It also provides log1p since it's needed to keep
some precision. Unfortunately, log1p does not give "full" precision, but it
should be good enough (of course, it's possible to implement a more precise
log1p, but that would be more convoluted.)

Cheers,
Luis


-- invhyp.lua: inverse hyperbolic trig functions
-- Adapted from glibc
local abs, log, sqrt = math.abs, math.log, math.sqrt
local log2 = log(2)

-- good for IEEE754, double precision
local function islarge (x) return x > 2 ^ 28 end
local function issmall (x) return x < 2 ^ (-28) end

local INF = math.huge
local function isinfornan (x)
  return x ~= x or x == INF or x == -INF
end

local function log1p (x) -- not very precise, but works well
  local u = 1 + x
  if u == 1 then return x end -- x < eps?
  return log(u) * x / (u - 1)
end

local function acosh (x)
  if x < 1 then return (x - x) / (x - x) end -- nan
  if islarge(x) then
    if isinfornan(x) then return x + x end
    return log2 + log(x)
  end
  if x + 1 == 1 then return 0 end -- acosh(1) == 0
  if x > 2 then
    local x2 = x * x
    return log(2 * x - 1 / (x + sqrt(x2 - 1)))
  end
  -- 1 < x < 2:
  local t = x - 1
  return log1p(t + sqrt(2 * t + t * t))
end

local function asinh (x)
  local y = abs(x)
  if issmall(y) then return x end
  local a
  if islarge(y) then -- very large?
    if isinfornan(x) then return x + x end
    a = log2 + log(y)
  else
    if y > 2 then
      a = log(2 * y + 1 / (y + sqrt(1 + y * y)))
    else
      local y2 = y * y
      a = log1p(y + y2 / (1 + sqrt(1 + y2)))
    end
  end
  return x < 0 and -a or a -- transfer sign
end

local function atanh (x)
  local y = abs(x)
  local a
  if y < .5 then
    if issmall(y) then return x end
    a = 2 * y
    a = .5 * log1p(a + a * y / (1 - y))
  else
    if y < 1 then
      a = .5 * log1p(2 * y / (1 - y))
    elseif y > 1 then
      return (x - x) / (x - x) -- nan
    else -- y == 1
      return x / 0 -- inf with sign
    end
  end
  return x < 0 and -a or a -- transfer sign
end

return {log1p = log1p, asinh = asinh, acosh = acosh, atanh = atanh}


-- 
Computers are useless. They can only give you answers.
                -- Pablo Picasso

-- 
Luis Carvalho (Kozure)
lua -e 'print((("lexcarvalho@NO.gmail.SPAM.com"):gsub("(%u+%.)","")))'