Re: Missing math functions
[
Date Prev][
Date Next][
Thread Prev][
Thread Next]
[
Date Index]
[
Thread Index]
- Subject: Re: Missing math functions
- From: Luis Carvalho <lexcarvalho@...>
- Date: 2013年1月19日 12:21:08 -0500
> 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+%.)","")))'