Closed SepandKashani closed 9 months ago
After some more digging, I believe you are right: the abs/rel errors obtained numerically are close to those listed in the tables.
from mpmath import mp
import drjit as dr
import drjit.llvm as drb
N = 1_000
a32 = dr.linspace(drb.Float, -20, 30, N)
a64 = dr.linspace(drb.Float64, -20, 30, N)
b32 = dr.exp(a32)
b64 = dr.exp(a64)
with mp.workdps(8):
gt32 = list(map(mp.exp, a32))
with mp.workdps(16):
gt64 = list(map(mp.exp, a64))
abs32, rel32 = [], []
abs64, rel64 = [], []
for _b, _gt in zip(b32, gt32):
abs32.append(abs(_b - _gt))
rel32.append(abs(_b - _gt) / _gt)
for _b, _gt in zip(b64, gt64):
abs64.append(abs(_b - _gt))
rel64.append(abs(_b - _gt) / _gt)
max(abs32), max(rel32), max(abs64), max(rel64)
# (mpf('770048.0'),
# mpf('1.078238811961534729e-7'),
# mpf('0.000732421875'),
# mpf('2.52219106387520338e-16'))
Sorry for the spam; this PR can be closed.
I think that's actually correct. Even a 1 ULP error can translate into a huge absolute error when working with the exponential function due to its large output.