Is there a way to accurately evaluate f(x) as given below for moderate negative x (around -1e6 to -1e2), while sticking with float64 numbers?
def f(x):
return 0.5 + arctan(x/sqrt(3)) / pi + sqrt(3)*x/(pi*(x*x+3))
For example:
g(-10_000) == 1.1026577511479050e-12 # actual closest float64 (what I want)
f(-10_000) == 1.1026461116240595e-12 # current implementation
# ^^^^^^^^^^^^ <- wrong
I ended up with this:
def g(x):
# approximation to f(x)
# method: For x<0 we have the identity: atan(x) === -pi/2 - atan(1/x)
# substitute this into f(x) then expand about -infinity, simplify & truncate
assert x <= -100
k = 3**0.5 / (pi * x)
y = -3.0/(x*x)
return k*y*(2/3 + y*(4/5 + y*(6/7 + y*(8/9 + y*10/11))))
This keeps the relative error below ~3e-16 for the range I'm interested in and has about the same performance, or even a little faster:
| x | rel_err(f) | rel_err(g) | abs_err(f) | abs_err(g) |
|---|---|---|---|---|
| -100 | -2.012e-11 | -2.774e-16 | -2.217e-17 | -3.058e-22 |
| -1,000 | -5.660e-08 | -4.603e-18 | -6.241e-17 | -5.076e-27 |
| -10,000 | -1.056e-05 | -1.571e-16 | -1.164e-17 | -1.733e-28 |
| -100,000 | -6.418e-02 | -2.529e-17 | -7.077e-17 | -2.789e-32 |
| -1,000,000 | -4.585e+01 | -9.824e-17 | -5.056e-17 | -1.083e-34 |